home *** CD-ROM | disk | FTP | other *** search
/ QBasic & Borland Pascal & C / Delphi5.iso / C / Samples / C-SSP.ARJ / MATEIG.C < prev    next >
Encoding:
Text File  |  1984-06-18  |  3.2 KB  |  102 lines
  1.     mateig(lra,n,a,s)
  2.  
  3.       /*to calculate eigen values and corresponding eigen*/
  4.       /*vectors of a real symmetric matrix by the Jacobi*/
  5.       /*method.*/
  6.  
  7.     int lra,n;
  8.     float a[],s[];
  9.  
  10.   {
  11.     double fabs(),sqrt();
  12.     int i,ia,ii,iia,iib,ij,ik,il,in,ip,iq,is,j,ja,ka,la,m;
  13.     float anorm,app,aqq,b,c,constf,csn,csn2,e,fnorm,snn,snn2,stc,u,v,w;
  14.  
  15.       e = 0.000001;
  16.       constf = n;
  17.       in = 0;
  18.     for(j = 0; j <= n-1; j++)
  19.      for(i = 0; i <= n-1; i++)
  20.       {
  21.       is = j * lra +i;             /* move columnwise to new element */
  22.       if ((i-j) == 0) s[is] = 1.;  /* look for diagonal */
  23.         else s[is] = 0.0;          /* set non-diagonal to Z */
  24.       }            /* end of i loop */
  25.                    /* end of j loop */
  26.     anorm = 0.0;
  27.     for(j = 0; j <= n-1; j++)
  28.      for(i = 0; i <= n-1; i++)
  29.       {
  30.       ia = i * lra +j;            /* column */
  31.       if ((j-i) != 0)
  32.         anorm = anorm + a[ia] * a[ia];
  33.       } /* end of i loop */
  34.         /* end of j loop */
  35.  
  36.       anorm = sqrt(anorm);
  37.       fnorm = anorm * e;
  38. l90:;
  39.       anorm = anorm/constf;
  40. l80:;
  41.       for(iq = 1; iq <= n-1; iq++)
  42.        {
  43.        m = iq - 1;
  44.        for(ip = 0; ip <= m; ip++)
  45.         {
  46.  
  47.          ia = lra * iq + ip;          /* row iq , col ip */
  48.          ja = lra * ip + iq;          /* row ip , col iq */
  49.          ka = lra * ip + ip;          /* row ip , col ip */
  50.          la = lra * iq + iq;          /* row iq , col iq */
  51.          if ((fabs(a[ja]) - anorm) <  0.0) goto l21;
  52.          in = 1;
  53.          u = -a[ja];
  54.          v = (a[ka] - a[la]) / 2.0;
  55.          w = u/(sqrt(u*u + v*v));
  56.          if (v <  0.0) w = -w;
  57.          snn = w/(sqrt(2.0 * (1.0 + sqrt(1.0 - w*w))));
  58.          snn2 = snn * snn;
  59.          csn = sqrt(1.0 - snn2);
  60.  
  61.          for(i = 0; i <= n-1; i++)
  62.          {
  63.           iia = lra * i + ip;               /* row i, col ip */
  64.           iib = lra * i + iq;               /* row i, col iq */
  65.           if ((i - ip) == 0) goto l301;
  66.            if ((i - iq)  == 0) goto l301;
  67.             b = a[iia] * csn - a[iib] * snn;
  68.             a[iib] = a[iia] * snn + a[iib] * csn;
  69.             a[iia] = b;
  70.   l301:
  71.          c = s[iia] * csn - s[iib] * snn;
  72.          s[iib] = s[iia] * snn + s[iib] * csn;
  73.          s[iia] = c;
  74.  
  75.          }   /* end of i loop */
  76.          csn2 = csn * csn;
  77.          stc = snn * csn;
  78.          app = a[ka]*csn2 + a[la]*snn2 - 2.0*a[ja]*stc;
  79.          aqq = a[ka]*snn2 + a[la] * csn2 + 2.0* a[ja] *stc;
  80.          a[ja] = (a[ka] - a[la])*stc + a[ja]*(csn2-snn2);
  81.          a[ia] = a[ja];
  82.          a[ka] = app;
  83.          a[la] = aqq;
  84.  
  85.          for(i = 0; i <= n-1; i++)
  86.          {
  87.           ii = lra * ip + i;             /* row ip, col i */
  88.           ij = lra * i + ip;             /* row i, col ip */
  89.           ik = lra * iq + i;             /* row iq, col i */
  90.           il = lra * i + iq;             /* row i, col iq */
  91.           a[ii] = a[ij];
  92.           a[ik] = a[il];
  93.          }              /* end of i loop */
  94. l21: ;
  95.         }               /* end of ip loop */
  96.        }                /* end of iq loop */
  97.       if ((in-1) != 0) goto l600;
  98.       in = 0;
  99.       goto l80;
  100. l600: if((anorm-fnorm) > 0.0) goto l90;
  101.   }
  102.