00001
00002
00003
00004
00005
00006
00007
00008 #define WANT_MATH
00009
00010 #include "include.h"
00011 #include "newmat.h"
00012 #include "precisio.h"
00013 #include "newmatrm.h"
00014
00015 #ifdef use_namespace
00016 namespace NEWMAT {
00017 #endif
00018
00019 void Jacobi(const SymmetricMatrix& X, DiagonalMatrix& D, SymmetricMatrix& A,
00020 Matrix& V, bool eivec) {
00021 Real epsilon = FloatingPointPrecision::Epsilon();
00022 Tracer et("Jacobi");
00023 int n = X.Nrows(); DiagonalMatrix B(n), Z(n); D.ReSize(n); A = X;
00024 if (eivec) { V.ReSize(n,n); D = 1.0; V = D; }
00025 B << A; D = B; Z = 0.0; A.Inject(Z);
00026 for (int i=1; i<=50; i++) {
00027 Real sm=0.0; Real* a = A.Store(); int p = A.Storage();
00028 while (p--) sm += fabs(*a++);
00029 if (sm==0.0) return;
00030 Real tresh = (i<4) ? 0.2 * sm / square(n) : 0.0; a = A.Store();
00031 for (p = 0; p < n; p++) {
00032 Real* ap1 = a + (p*(p+1))/2;
00033 Real& zp = Z.element(p); Real& dp = D.element(p);
00034 for (int q = p+1; q < n; q++) {
00035 Real* ap = ap1; Real* aq = a + (q*(q+1))/2;
00036 Real& zq = Z.element(q); Real& dq = D.element(q);
00037 Real& apq = A.element(q,p);
00038 Real g = 100 * fabs(apq); Real adp = fabs(dp); Real adq = fabs(dq);
00039
00040 if (i>4 && g < epsilon*adp && g < epsilon*adq) apq = 0.0;
00041 else if (fabs(apq) > tresh) {
00042 Real t; Real h = dq - dp; Real ah = fabs(h);
00043 if (g < epsilon*ah) t = apq / h;
00044 else {
00045 Real theta = 0.5 * h / apq;
00046 t = 1.0 / ( fabs(theta) + sqrt(1.0 + square(theta)) );
00047 if (theta<0.0) t = -t;
00048 }
00049 Real c = 1.0 / sqrt(1.0 + square(t)); Real s = t * c;
00050 Real tau = s / (1.0 + c); h = t * apq;
00051 zp -= h; zq += h; dp -= h; dq += h; apq = 0.0;
00052 int j = p;
00053 while (j--) {
00054 g = *ap; h = *aq;
00055 *ap++ = g-s*(h+g*tau); *aq++ = h+s*(g-h*tau);
00056 }
00057 int ip = p+1; j = q-ip; ap += ip++; aq++;
00058 while (j--) {
00059 g = *ap; h = *aq;
00060 *ap = g-s*(h+g*tau); *aq++ = h+s*(g-h*tau);
00061 ap += ip++;
00062 }
00063 int iq = q+1; j = n-iq; ap += ip++; aq += iq++;
00064 while (j--) {
00065 g = *ap; h = *aq;
00066 *ap = g-s*(h+g*tau); *aq = h+s*(g-h*tau);
00067 ap += ip++; aq += iq++;
00068 }
00069 if (eivec) {
00070 RectMatrixCol VP(V,p); RectMatrixCol VQ(V,q);
00071 Rotate(VP, VQ, tau, s);
00072 }
00073 }
00074 }
00075 }
00076 B = B + Z; D = B; Z = 0.0;
00077 }
00078 Throw(ConvergenceException(X));
00079 }
00080
00081 void Jacobi(const SymmetricMatrix& X, DiagonalMatrix& D)
00082 { SymmetricMatrix A; Matrix V; Jacobi(X,D,A,V,false); }
00083
00084 void Jacobi(const SymmetricMatrix& X, DiagonalMatrix& D, SymmetricMatrix& A)
00085 { Matrix V; Jacobi(X,D,A,V,false); }
00086
00087 void Jacobi(const SymmetricMatrix& X, DiagonalMatrix& D, Matrix& V)
00088 { SymmetricMatrix A; Jacobi(X,D,A,V,true); }
00089
00090 #ifdef use_namespace
00091 }
00092 #endif
