$search
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