Tekkotsu Homepage
Demos
Overview
Downloads
Dev. Resources
Reference
Credits

jacobi.cpp

Go to the documentation of this file.
00001 //$$jacobi.cpp                           jacobi eigenvalue analysis
00002 
00003 // Copyright (C) 1991,2,3,4: R B Davies
00004 
00005 
00006 //#define WANT_STREAM
00007 
00008 #ifndef WANT_MATH
00009 #define WANT_MATH
00010 #endif
00011 
00012 #include "include.h"
00013 #include "newmatap.h"
00014 #include "precisio.h"
00015 #include "newmatrm.h"
00016 
00017 #ifdef use_namespace
00018 namespace NEWMAT {
00019 #endif
00020 
00021 #ifdef DO_REPORT
00022 #define REPORT { static ExeCounter ExeCount(__LINE__,18); ++ExeCount; }
00023 #else
00024 #define REPORT {}
00025 #endif
00026 
00027 
00028 void Jacobi(const SymmetricMatrix& X, DiagonalMatrix& D, SymmetricMatrix& A,
00029    Matrix& V, bool eivec)
00030 {
00031    Real epsilon = FloatingPointPrecision::Epsilon();
00032    Tracer et("Jacobi");
00033    REPORT
00034    int n = X.Nrows(); DiagonalMatrix B(n), Z(n); D.ReSize(n); A = X;
00035    if (eivec) { REPORT V.ReSize(n,n); D = 1.0; V = D; }
00036    B << A; D = B; Z = 0.0; A.Inject(Z);
00037    bool converged = false;
00038    for (int i=1; i<=50; i++)
00039    {
00040       Real sm=0.0; Real* a = A.Store(); int p = A.Storage();
00041       while (p--) sm += fabs(*a++);            // have previously zeroed diags
00042       if (sm==0.0) { REPORT converged = true; break; }
00043       Real tresh = (i<4) ? 0.2 * sm / square(n) : 0.0; a = A.Store();
00044       for (p = 0; p < n; p++)
00045       {
00046          Real* ap1 = a + (p*(p+1))/2;
00047          Real& zp = Z.element(p); Real& dp = D.element(p);
00048          for (int q = p+1; q < n; q++)
00049          {
00050             Real* ap = ap1; Real* aq = a + (q*(q+1))/2;
00051             Real& zq = Z.element(q); Real& dq = D.element(q);
00052             Real& apq = A.element(q,p);
00053             Real g = 100 * fabs(apq); Real adp = fabs(dp); Real adq = fabs(dq);
00054 
00055             if (i>4 && g < epsilon*adp && g < epsilon*adq) { REPORT apq = 0.0; }
00056             else if (fabs(apq) > tresh)
00057             {
00058                REPORT
00059                Real t; Real h = dq - dp; Real ah = fabs(h);
00060                if (g < epsilon*ah) { REPORT t = apq / h; }
00061                else
00062                {
00063                   REPORT
00064                   Real theta = 0.5 * h / apq;
00065                   t = 1.0 / ( fabs(theta) + sqrt(1.0 + square(theta)) );
00066                   if (theta<0.0) { REPORT t = -t; }
00067                }
00068                Real c = 1.0 / sqrt(1.0 + square(t)); Real s = t * c;
00069                Real tau = s / (1.0 + c); h = t * apq;
00070                zp -= h; zq += h; dp -= h; dq += h; apq = 0.0;
00071                int j = p;
00072                while (j--)
00073                {
00074                   g = *ap; h = *aq;
00075                   *ap++ = g-s*(h+g*tau); *aq++ = h+s*(g-h*tau);
00076                }
00077                int ip = p+1; j = q-ip; ap += ip++; aq++;
00078                while (j--)
00079                {
00080                   g = *ap; h = *aq;
00081                   *ap = g-s*(h+g*tau); *aq++ = h+s*(g-h*tau);
00082                   ap += ip++;
00083                }
00084                if (q < n-1)             // last loop is non-empty
00085                {
00086                   int iq = q+1; j = n-iq; ap += ip++; aq += iq++;
00087                   for (;;)
00088                   {
00089                      g = *ap; h = *aq;
00090                      *ap = g-s*(h+g*tau); *aq = h+s*(g-h*tau);
00091                      if (!(--j)) break;
00092                      ap += ip++; aq += iq++;
00093                   }
00094                }
00095                if (eivec)
00096                {
00097                   REPORT
00098                   RectMatrixCol VP(V,p); RectMatrixCol VQ(V,q);
00099                   Rotate(VP, VQ, tau, s);
00100                }
00101             }
00102          }
00103       }
00104       B = B + Z; D = B; Z = 0.0;
00105    }
00106    if (!converged) Throw(ConvergenceException(X));
00107    if (eivec) SortSV(D, V, true);
00108    else SortAscending(D);
00109 }
00110 
00111 void Jacobi(const SymmetricMatrix& X, DiagonalMatrix& D)
00112 { REPORT SymmetricMatrix A; Matrix V; Jacobi(X,D,A,V,false); }
00113 
00114 void Jacobi(const SymmetricMatrix& X, DiagonalMatrix& D, SymmetricMatrix& A)
00115 { REPORT Matrix V; Jacobi(X,D,A,V,false); }
00116 
00117 void Jacobi(const SymmetricMatrix& X, DiagonalMatrix& D, Matrix& V)
00118 { REPORT SymmetricMatrix A; Jacobi(X,D,A,V,true); }
00119 
00120 
00121 #ifdef use_namespace
00122 }
00123 #endif
00124 

newmat11b
Generated Mon May 9 04:54:17 2016 by Doxygen 1.6.3