home *** CD-ROM | disk | FTP | other *** search
/ Aminet 18 / aminetcdnumber181997.iso / Aminet / dev / c / math_classes.lha / math_classes / matrix / jacobi.cc < prev    next >
Encoding:
C/C++ Source or Header  |  1996-09-19  |  1.2 KB  |  119 lines
  1.  
  2. #include <iostream.h>
  3. #include "Matrix.h"
  4. #include <math.h>
  5.  
  6.  
  7. double abs(double a);    // defined in matrix.cc
  8.  
  9.  
  10. Matrix jacobi(Matrix a0)
  11.     {
  12.     Matrix a(a0);
  13.     mtrxtype aii,aij,ajj;
  14.     mtrxtype t,s,c,tau;
  15.     mtrxtype r1,r2,temp;
  16.     mtrxtype p;
  17.     int i,j,k;
  18.     int f,g;
  19.     
  20.     mtrxtype aik,ajk;
  21.     
  22.     if(a.rows() != a.cols()) 
  23.         {
  24.         cout << "Matrix must be square" << endl;
  25.         exit(1);
  26.         }
  27.     
  28.     temp=0;
  29.     for(f=1;f<=a.rows();f++)
  30.         {
  31.         for(g=f+1;g<=a.cols();g++)
  32.             {
  33.             if(abs(a(f,g)) > temp)
  34.                 {
  35.                 temp=abs(a(f,g));
  36.                 i=f;
  37.                 j=g;
  38.                 }
  39.             }
  40.         }
  41.         
  42.         
  43.     aii=a(i,i);
  44.     ajj=a(j,j);
  45.     aij=a(i,j);
  46.     
  47.     p=-(ajj-aii)/aij/2;
  48.     
  49.     r1=p+sqrt(p*p+1);
  50.     r2=p-sqrt(p*p+1);
  51.     
  52.     if(abs(r1) < abs(r2))
  53.         {
  54.         t=r1;
  55.         }
  56.     else
  57.         {
  58.         t=r2;
  59.         }
  60.         
  61.     c=1/sqrt(1+t*t);
  62.     
  63.     s=c*t;
  64.     
  65.     tau=s/(1+c);
  66.     
  67.     
  68.     for(k=1;k<=a.rows();k++)
  69.         {
  70.         if(k!=i && k!=j)
  71.             {
  72.             aik=a(i,k)-s*(a(j,k)+tau*a(i,k));
  73.             ajk=a(j,k)+s*(a(i,k)-tau*a(j,k));
  74.             a(i,k)=aik;
  75.             a(j,k)=ajk;
  76.             a(k,j)=ajk;
  77.             a(k,i)=aik;
  78.             }
  79.         }
  80.         
  81.     a(i,i)=aii-t*aij;
  82.     a(j,j)=ajj+t*aij;
  83.     a(i,j)=0;
  84.     a(j,i)=0;
  85.     
  86.     
  87.     return a;
  88.     
  89.     }
  90.     
  91.     
  92.     
  93. int main()
  94.     {
  95.     mtrxtype b[]={1,2,3,4,2,1,0,3,3,0,1,2,4,3,2,1};
  96.     Matrix a(4,4,b);
  97.     int i;    
  98.     
  99.     print(a);
  100.     cout << endl;
  101.     
  102.     for(i=1;i<=8;i++)
  103.         {
  104.         a=jacobi(a);
  105.         print(a);
  106.         cout << endl;
  107.         }
  108.     return 0;
  109.     }
  110.     
  111.  
  112.             
  113.             
  114.             
  115.             
  116.     
  117.     
  118.     
  119.