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