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


kni
Author(s): Martin Günther
autogenerated on Thu Aug 27 2015 13:40:06