home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / newmat03 / jacobi.cxx < prev    next >
Encoding:
C/C++ Source or Header  |  1991-11-20  |  2.8 KB  |  92 lines
  1. //$$jacobi.cxx                           jacobi eigenvalue analysis
  2.  
  3. // Copyright (C) 1991: R B Davies and DSIR
  4.  
  5. #define WANT_MATH
  6.  
  7. #include "include.hxx"
  8. #include "newmat.hxx"
  9. #include "precisio.hxx"
  10. #include "newmatrm.hxx"
  11.  
  12.  
  13.  
  14. void Jacobi(const SymmetricMatrix& X, DiagonalMatrix& D, SymmetricMatrix& A,
  15.    Matrix& V, BOOL eivec)
  16. {
  17.    int n = X.Nrows(); DiagonalMatrix B(n), Z(n); D.ReDimension(n); A = X.c();
  18.    if (eivec) { V.ReDimension(n,n); D = 1.0; V = D; }
  19.    B << A; D = B; Z = 0.0; A.Inject(Z);
  20.    for (int i=1; i<=50; i++)
  21.    {
  22.       real sm=0.0; real* a = A.Store(); int p = A.Storage();
  23.       while (p--) sm += fabs(*a++);            // have previously zeroed diags
  24.       if (sm==0.0) return;
  25.       real tresh = (i<4) ? 0.2 * sm / square(n) : 0.0; a = A.Store();
  26.       for (p = 0; p < n; p++)
  27.       {
  28.      real* ap1 = a + (p*(p+1))/2;
  29.          real& zp = Z.element(p); real& dp = D.element(p);
  30.      for (int q = p+1; q < n; q++)
  31.      {
  32.         real* ap = ap1; real* aq = a + (q*(q+1))/2;
  33.             real& zq = Z.element(q); real& dq = D.element(q);
  34.             real& apq = A.element(q,p);
  35.             real g = 100 * fabs(apq); real adp = fabs(dp); real adq = fabs(dq);
  36.         if (i>4 && adp+g==adp && adq+g==adq) apq = 0.0;
  37.         else if (fabs(apq) > tresh)
  38.         {
  39.            real t; real h = dq - dp; real ah = fabs(h);
  40.            if (ah+g==ah) t = apq / h;
  41.            else
  42.            {
  43.           real theta = 0.5 * h / apq;
  44.           t = 1.0 / ( fabs(theta) + sqrt(1.0 + square(theta)) );
  45.           if (theta<0.0) t = -t;
  46.            }
  47.            real c = 1.0 / sqrt(1.0 + square(t)); real s = t * c;
  48.            real tau = s / (1.0 + c); h = t * apq;
  49.                zp -= h; zq += h; dp -= h; dq += h; apq = 0.0;
  50.            int j = p;
  51.            while (j--)
  52.            {
  53.           g = *ap; h = *aq;
  54.           *ap++ = g-s*(h+g*tau); *aq++ = h+s*(g-h*tau);
  55.            }
  56.            int ip = p+1; j = q-ip; ap += ip++; aq++;
  57.            while (j--)
  58.            {
  59.           g = *ap; h = *aq;
  60.                   *ap = g-s*(h+g*tau); *aq++ = h+s*(g-h*tau);
  61.           ap += ip++;
  62.            }
  63.            int iq = q+1; j = n-iq; ap += ip++; aq += iq++;
  64.            while (j--)
  65.            {
  66.           g = *ap; h = *aq;
  67.           *ap = g-s*(h+g*tau); *aq = h+s*(g-h*tau);
  68.           ap += ip++; aq += iq++;
  69.            }
  70.            if (eivec)
  71.            {
  72.                   RectMatrixCol VP(V,p); RectMatrixCol VQ(V,q);
  73.                   Rotate(VP, VQ, tau, s);
  74.            }
  75.         }
  76.      }
  77.       }
  78.       B = B + Z; D = B; Z = 0.0;
  79.    }
  80. }
  81.  
  82. void Jacobi(const SymmetricMatrix& X, DiagonalMatrix& D)
  83. { SymmetricMatrix A; Matrix V; Jacobi(X,D,A,V,FALSE); }
  84.  
  85. void Jacobi(const SymmetricMatrix& X, DiagonalMatrix& D, SymmetricMatrix& A)
  86. { Matrix V; Jacobi(X,D,A,V,FALSE); }
  87.  
  88. void Jacobi(const SymmetricMatrix& X, DiagonalMatrix& D, Matrix& V)
  89. { SymmetricMatrix A; Jacobi(X,D,A,V,TRUE); }
  90.  
  91.  
  92.