Tekkotsu Homepage
Demos
Overview
Downloads
Dev. Resources
Reference
Credits

evalue.cpp

Go to the documentation of this file.
00001 //$$evalue.cpp                           eigen-value decomposition
00002 
00003 // Copyright (C) 1991,2,3,4: R B Davies
00004 
00005 #ifndef WANT_MATH
00006 #define WANT_MATH
00007 #endif
00008 
00009 #include "include.h"
00010 #include "newmatap.h"
00011 #include "newmatrm.h"
00012 #include "precisio.h"
00013 
00014 #ifdef use_namespace
00015 namespace NEWMAT {
00016 #endif
00017 
00018 #ifdef DO_REPORT
00019 #define REPORT { static ExeCounter ExeCount(__LINE__,17); ++ExeCount; }
00020 #else
00021 #define REPORT {}
00022 #endif
00023 
00024 
00025 
00026 static void tred2(const SymmetricMatrix& A, DiagonalMatrix& D,
00027    DiagonalMatrix& E, Matrix& Z)
00028 {
00029    Tracer et("Evalue(tred2)");
00030    REPORT
00031    Real tol =
00032       FloatingPointPrecision::Minimum()/FloatingPointPrecision::Epsilon();
00033    int n = A.Nrows(); Z.ReSize(n,n); Z.Inject(A);
00034    D.ReSize(n); E.ReSize(n);
00035    Real* z = Z.Store(); int i;
00036 
00037    for (i=n-1; i > 0; i--)                   // i=0 is excluded
00038    {
00039       Real f = Z.element(i,i-1); Real g = 0.0;
00040       int k = i-1; Real* zik = z + i*n;
00041       while (k--) g += square(*zik++);
00042       Real h = g + square(f);
00043       if (g <= tol) { REPORT E.element(i) = f; h = 0.0; }
00044       else
00045       {
00046          REPORT
00047          g = sign(-sqrt(h), f); E.element(i) = g; h -= f*g;
00048          Z.element(i,i-1) = f-g; f = 0.0;
00049          Real* zji = z + i; Real* zij = z + i*n; Real* ej = E.Store();
00050          int j;
00051          for (j=0; j<i; j++)
00052          {
00053             *zji = (*zij++)/h; g = 0.0;
00054             Real* zjk = z + j*n; zik = z + i*n;
00055             k = j; while (k--) g += *zjk++ * (*zik++);
00056             k = i-j;
00057             if (k) for(;;)
00058                { g += *zjk * (*zik++); if (!(--k)) break; zjk += n; }
00059             *ej++ = g/h; f += g * (*zji); zji += n;
00060          }
00061          Real hh = f / (h + h); zij = z + i*n; ej = E.Store();
00062          for (j=0; j<i; j++)
00063          {
00064             f = *zij++; g = *ej - hh * f; *ej++ = g;
00065             Real* zjk = z + j*n; Real* zik = z + i*n;
00066             Real* ek = E.Store(); k = j+1;
00067             while (k--)  *zjk++ -= ( f*(*ek++) + g*(*zik++) ); 
00068          }
00069       }
00070       D.element(i) = h;
00071    }
00072 
00073    D.element(0) = 0.0; E.element(0) = 0.0;
00074    for (i=0; i<n; i++)
00075    {
00076       if (D.element(i) != 0.0)
00077       {
00078          REPORT
00079          for (int j=0; j<i; j++)
00080          {
00081             Real g = 0.0;
00082             Real* zik = z + i*n; Real* zkj = z + j;
00083             int k = i;
00084             if (k) for (;;)
00085                { g += *zik++ * (*zkj); if (!(--k)) break; zkj += n; }
00086             Real* zki = z + i; zkj = z + j;
00087             k = i;
00088             if (k) for (;;)
00089                { *zkj -= g * (*zki); if (!(--k)) break; zkj += n; zki += n; }
00090          }
00091       }
00092       Real* zij = z + i*n; Real* zji = z + i;
00093       int j = i;
00094       if (j) for (;;)
00095          { *zij++ = 0.0; *zji = 0.0; if (!(--j)) break; zji += n; }
00096       D.element(i) = *zij; *zij = 1.0;
00097    }
00098 }
00099 
00100 static void tql2(DiagonalMatrix& D, DiagonalMatrix& E, Matrix& Z)
00101 {
00102    Tracer et("Evalue(tql2)");
00103    REPORT
00104    Real eps = FloatingPointPrecision::Epsilon();
00105    int n = D.Nrows(); Real* z = Z.Store(); int l;
00106    for (l=1; l<n; l++) E.element(l-1) = E.element(l);
00107    Real b = 0.0; Real f = 0.0; E.element(n-1) = 0.0;
00108    for (l=0; l<n; l++)
00109    {
00110       int i,j;
00111       Real& dl = D.element(l); Real& el = E.element(l);
00112       Real h = eps * ( fabs(dl) + fabs(el) );
00113       if (b < h) { REPORT b = h; }
00114       int m;
00115       for (m=l; m<n; m++) if (fabs(E.element(m)) <= b) break;
00116       bool test = false;
00117       for (j=0; j<30; j++)
00118       {
00119          if (m==l) { REPORT test = true; break; }
00120          Real& dl1 = D.element(l+1);
00121          Real g = dl; Real p = (dl1-g) / (2.0*el); Real r = sqrt(p*p + 1.0);
00122          dl = el / (p < 0.0 ? p-r : p+r); Real h = g - dl; f += h;
00123          Real* dlx = &dl1; i = n-l-1; while (i--) *dlx++ -= h;
00124 
00125          p = D.element(m); Real c = 1.0; Real s = 0.0;
00126          for (i=m-1; i>=l; i--)
00127          {
00128             Real ei = E.element(i); Real di = D.element(i);
00129             Real& ei1 = E.element(i+1);
00130             g = c * ei; h = c * p;
00131             if ( fabs(p) >= fabs(ei))
00132             {
00133                REPORT
00134                c = ei / p; r = sqrt(c*c + 1.0);
00135                ei1 = s*p*r; s = c/r; c = 1.0/r;
00136             }
00137             else
00138             {
00139                REPORT
00140                c = p / ei; r = sqrt(c*c + 1.0);
00141                ei1 = s * ei * r; s = 1.0/r; c /= r;
00142             }
00143             p = c * di - s*g; D.element(i+1) = h + s * (c*g + s*di);
00144 
00145             Real* zki = z + i; Real* zki1 = zki + 1; int k = n;
00146             if (k) for (;;)
00147             {
00148                REPORT
00149                h = *zki1; *zki1 = s*(*zki) + c*h; *zki = c*(*zki) - s*h;
00150                if (!(--k)) break;
00151                zki += n; zki1 += n;
00152             }
00153          }
00154          el = s*p; dl = c*p;
00155          if (fabs(el) <= b) { REPORT; test = true; break; }
00156       }
00157       if (!test) Throw ( ConvergenceException(D) );
00158       dl += f;
00159    }
00160 /*
00161    for (int i=0; i<n; i++)
00162    {
00163       int k = i; Real p = D.element(i);
00164       for (int j=i+1; j<n; j++)
00165          { if (D.element(j) < p) { k = j; p = D.element(j); } }
00166       if (k != i)
00167       {
00168          D.element(k) = D.element(i); D.element(i) = p; int j = n;
00169          Real* zji = z + i; Real* zjk = z + k;
00170          if (j) for(;;)
00171          {
00172             p = *zji; *zji = *zjk; *zjk = p;
00173             if (!(--j)) break;
00174             zji += n; zjk += n;
00175          }
00176       }
00177    }
00178 */
00179 }
00180 
00181 static void tred3(const SymmetricMatrix& X, DiagonalMatrix& D,
00182    DiagonalMatrix& E, SymmetricMatrix& A)
00183 {
00184    Tracer et("Evalue(tred3)");
00185    REPORT
00186    Real tol =
00187       FloatingPointPrecision::Minimum()/FloatingPointPrecision::Epsilon();
00188    int n = X.Nrows(); A = X; D.ReSize(n); E.ReSize(n);
00189    Real* ei = E.Store() + n;
00190    for (int i = n-1; i >= 0; i--)
00191    {
00192       Real h = 0.0; Real f = - FloatingPointPrecision::Maximum();
00193       Real* d = D.Store(); Real* a = A.Store() + (i*(i+1))/2; int k = i;
00194       while (k--) { f = *a++; *d++ = f; h += square(f); }
00195       if (h <= tol) { REPORT *(--ei) = 0.0; h = 0.0; }
00196       else
00197       {
00198          REPORT
00199          Real g = sign(-sqrt(h), f); *(--ei) = g; h -= f*g;
00200          f -= g; *(d-1) = f; *(a-1) = f; f = 0.0;
00201          Real* dj = D.Store(); Real* ej = E.Store(); int j;
00202          for (j = 0; j < i; j++)
00203          {
00204             Real* dk = D.Store(); Real* ak = A.Store()+(j*(j+1))/2;
00205             Real g = 0.0; k = j;
00206             while (k--)  g += *ak++ * *dk++;
00207             k = i-j; int l = j; 
00208             if (k) for (;;) { g += *ak * *dk++; if (!(--k)) break; ak += ++l; }
00209             g /= h; *ej++ = g; f += g * *dj++;
00210          }  
00211          Real hh = f / (2 * h); Real* ak = A.Store();
00212          dj = D.Store(); ej = E.Store();
00213          for (j = 0; j < i; j++)
00214          {
00215             f = *dj++; g = *ej - hh * f; *ej++ = g;
00216             Real* dk = D.Store(); Real* ek = E.Store(); k = j+1;
00217             while (k--) { *ak++ -= (f * *ek++ + g * *dk++); }
00218          }
00219       }
00220       *d = *a; *a = h;
00221    }
00222 }
00223 
00224 static void tql1(DiagonalMatrix& D, DiagonalMatrix& E)
00225 {
00226    Tracer et("Evalue(tql1)");
00227    REPORT
00228    Real eps = FloatingPointPrecision::Epsilon();
00229    int n = D.Nrows(); int l;
00230    for (l=1; l<n; l++) E.element(l-1) = E.element(l);
00231    Real b = 0.0; Real f = 0.0; E.element(n-1) = 0.0;
00232    for (l=0; l<n; l++)
00233    {
00234       int i,j;
00235       Real& dl = D.element(l); Real& el = E.element(l);
00236       Real h = eps * ( fabs(dl) + fabs(el) );
00237       if (b < h) b = h;
00238       int m;
00239       for (m=l; m<n; m++) if (fabs(E.element(m)) <= b) break;
00240       bool test = false;
00241       for (j=0; j<30; j++)
00242       {
00243          if (m==l) { REPORT test = true; break; }
00244          Real& dl1 = D.element(l+1);
00245          Real g = dl; Real p = (dl1-g) / (2.0*el); Real r = sqrt(p*p + 1.0);
00246          dl = el / (p < 0.0 ? p-r : p+r); Real h = g - dl; f += h;
00247          Real* dlx = &dl1; i = n-l-1; while (i--) *dlx++ -= h;
00248 
00249          p = D.element(m); Real c = 1.0; Real s = 0.0;
00250          for (i=m-1; i>=l; i--)
00251          {
00252             Real ei = E.element(i); Real di = D.element(i);
00253             Real& ei1 = E.element(i+1);
00254             g = c * ei; h = c * p;
00255             if ( fabs(p) >= fabs(ei))
00256             {
00257                REPORT
00258                c = ei / p; r = sqrt(c*c + 1.0);
00259                ei1 = s*p*r; s = c/r; c = 1.0/r;
00260             }
00261             else
00262             {
00263                REPORT
00264                c = p / ei; r = sqrt(c*c + 1.0);
00265                ei1 = s * ei * r; s = 1.0/r; c /= r;
00266             }
00267             p = c * di - s*g; D.element(i+1) = h + s * (c*g + s*di);
00268          }
00269          el = s*p; dl = c*p;
00270          if (fabs(el) <= b) { REPORT test = true; break; }
00271       }
00272       if (!test) Throw ( ConvergenceException(D) );
00273       Real p = dl + f;
00274       test = false;
00275       for (i=l; i>0; i--)
00276       {
00277          if (p < D.element(i-1)) { REPORT D.element(i) = D.element(i-1); }
00278          else { REPORT test = true; break; }
00279       }
00280       if (!test) i=0;
00281       D.element(i) = p;
00282    }
00283 }
00284 
00285 void EigenValues(const SymmetricMatrix& A, DiagonalMatrix& D, Matrix& Z)
00286 { REPORT DiagonalMatrix E; tred2(A, D, E, Z); tql2(D, E, Z); SortSV(D,Z,true); }
00287 
00288 void EigenValues(const SymmetricMatrix& X, DiagonalMatrix& D)
00289 { REPORT DiagonalMatrix E; SymmetricMatrix A; tred3(X,D,E,A); tql1(D,E); }
00290 
00291 void EigenValues(const SymmetricMatrix& X, DiagonalMatrix& D,
00292    SymmetricMatrix& A)
00293 { REPORT DiagonalMatrix E; tred3(X,D,E,A); tql1(D,E); }
00294 
00295 
00296 #ifdef use_namespace
00297 }
00298 #endif
00299 

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