1 //$$evalue.cpp eigen-value decomposition
2
3 // Copyright (C) 1991,2,3,4: R B Davies
4
5 #define WANT_MATH
6
7 #include "include.h"
8 #include "newmatap.h"
9 #include "newmatrm.h"
10 #include "precisio.h"
11
12 #ifdef use_namespace
13 namespace NEWMAT {
14 #endif
15
16 #ifdef DO_REPORT
17 #define REPORT { static ExeCounter ExeCount(__LINE__,17); ++ExeCount; }
18 #else
19 #define REPORT {}
20 #endif
21
22
23
tred2(const SymmetricMatrix & A,DiagonalMatrix & D,DiagonalMatrix & E,Matrix & Z)24 static void tred2(const SymmetricMatrix& A, DiagonalMatrix& D,
25 DiagonalMatrix& E, Matrix& Z)
26 {
27 Tracer et("Evalue(tred2)");
28 REPORT
29 Real tol =
30 FloatingPointPrecision::Minimum()/FloatingPointPrecision::Epsilon();
31 int n = A.Nrows(); Z.ReSize(n,n); Z.Inject(A);
32 D.ReSize(n); E.ReSize(n);
33 Real* z = Z.Store(); int i;
34
35 for (i=n-1; i > 0; i--) // i=0 is excluded
36 {
37 Real f = Z.element(i,i-1); Real g = 0.0;
38 int k = i-1; Real* zik = z + i*n;
39 while (k--) g += square(*zik++);
40 Real h = g + square(f);
41 if (g <= tol) { REPORT E.element(i) = f; h = 0.0; }
42 else
43 {
44 REPORT
45 g = sign(-sqrt(h), f); E.element(i) = g; h -= f*g;
46 Z.element(i,i-1) = f-g; f = 0.0;
47 Real* zji = z + i; Real* zij = z + i*n; Real* ej = E.Store();
48 int j;
49 for (j=0; j<i; j++)
50 {
51 *zji = (*zij++)/h; g = 0.0;
52 Real* zjk = z + j*n; zik = z + i*n;
53 k = j; while (k--) g += *zjk++ * (*zik++);
54 k = i-j;
55 if (k) for(;;)
56 { g += *zjk * (*zik++); if (!(--k)) break; zjk += n; }
57 *ej++ = g/h; f += g * (*zji); zji += n;
58 }
59 Real hh = f / (h + h); zij = z + i*n; ej = E.Store();
60 for (j=0; j<i; j++)
61 {
62 f = *zij++; g = *ej - hh * f; *ej++ = g;
63 Real* zjk = z + j*n; Real* zik = z + i*n;
64 Real* ek = E.Store(); k = j+1;
65 while (k--) *zjk++ -= ( f*(*ek++) + g*(*zik++) );
66 }
67 }
68 D.element(i) = h;
69 }
70
71 D.element(0) = 0.0; E.element(0) = 0.0;
72 for (i=0; i<n; i++)
73 {
74 if (D.element(i) != 0.0)
75 {
76 REPORT
77 for (int j=0; j<i; j++)
78 {
79 Real g = 0.0;
80 Real* zik = z + i*n; Real* zkj = z + j;
81 int k = i;
82 if (k) for (;;)
83 { g += *zik++ * (*zkj); if (!(--k)) break; zkj += n; }
84 Real* zki = z + i; zkj = z + j;
85 k = i;
86 if (k) for (;;)
87 { *zkj -= g * (*zki); if (!(--k)) break; zkj += n; zki += n; }
88 }
89 }
90 Real* zij = z + i*n; Real* zji = z + i;
91 int j = i;
92 if (j) for (;;)
93 { *zij++ = 0.0; *zji = 0.0; if (!(--j)) break; zji += n; }
94 D.element(i) = *zij; *zij = 1.0;
95 }
96 }
97
tql2(DiagonalMatrix & D,DiagonalMatrix & E,Matrix & Z)98 static void tql2(DiagonalMatrix& D, DiagonalMatrix& E, Matrix& Z)
99 {
100 Tracer et("Evalue(tql2)");
101 REPORT
102 Real eps = FloatingPointPrecision::Epsilon();
103 int n = D.Nrows(); Real* z = Z.Store(); int l;
104 for (l=1; l<n; l++) E.element(l-1) = E.element(l);
105 Real b = 0.0; Real f = 0.0; E.element(n-1) = 0.0;
106 for (l=0; l<n; l++)
107 {
108 int i,j;
109 Real& dl = D.element(l); Real& el = E.element(l);
110 Real h = eps * ( fabs(dl) + fabs(el) );
111 if (b < h) { REPORT b = h; }
112 int m;
113 for (m=l; m<n; m++) if (fabs(E.element(m)) <= b) break;
114 bool test = false;
115 for (j=0; j<30; j++)
116 {
117 if (m==l) { REPORT test = true; break; }
118 Real& dl1 = D.element(l+1);
119 Real g = dl; Real p = (dl1-g) / (2.0*el); Real r = sqrt(p*p + 1.0);
120 dl = el / (p < 0.0 ? p-r : p+r); Real h = g - dl; f += h;
121 Real* dlx = &dl1; i = n-l-1; while (i--) *dlx++ -= h;
122
123 p = D.element(m); Real c = 1.0; Real s = 0.0;
124 for (i=m-1; i>=l; i--)
125 {
126 Real ei = E.element(i); Real di = D.element(i);
127 Real& ei1 = E.element(i+1);
128 g = c * ei; h = c * p;
129 if ( fabs(p) >= fabs(ei))
130 {
131 REPORT
132 c = ei / p; r = sqrt(c*c + 1.0);
133 ei1 = s*p*r; s = c/r; c = 1.0/r;
134 }
135 else
136 {
137 REPORT
138 c = p / ei; r = sqrt(c*c + 1.0);
139 ei1 = s * ei * r; s = 1.0/r; c /= r;
140 }
141 p = c * di - s*g; D.element(i+1) = h + s * (c*g + s*di);
142
143 Real* zki = z + i; Real* zki1 = zki + 1; int k = n;
144 if (k) for (;;)
145 {
146 REPORT
147 h = *zki1; *zki1 = s*(*zki) + c*h; *zki = c*(*zki) - s*h;
148 if (!(--k)) break;
149 zki += n; zki1 += n;
150 }
151 }
152 el = s*p; dl = c*p;
153 if (fabs(el) <= b) { REPORT; test = true; break; }
154 }
155 if (!test) Throw ( ConvergenceException(D) );
156 dl += f;
157 }
158 /*
159 for (int i=0; i<n; i++)
160 {
161 int k = i; Real p = D.element(i);
162 for (int j=i+1; j<n; j++)
163 { if (D.element(j) < p) { k = j; p = D.element(j); } }
164 if (k != i)
165 {
166 D.element(k) = D.element(i); D.element(i) = p; int j = n;
167 Real* zji = z + i; Real* zjk = z + k;
168 if (j) for(;;)
169 {
170 p = *zji; *zji = *zjk; *zjk = p;
171 if (!(--j)) break;
172 zji += n; zjk += n;
173 }
174 }
175 }
176 */
177 }
178
tred3(const SymmetricMatrix & X,DiagonalMatrix & D,DiagonalMatrix & E,SymmetricMatrix & A)179 static void tred3(const SymmetricMatrix& X, DiagonalMatrix& D,
180 DiagonalMatrix& E, SymmetricMatrix& A)
181 {
182 Tracer et("Evalue(tred3)");
183 REPORT
184 Real tol =
185 FloatingPointPrecision::Minimum()/FloatingPointPrecision::Epsilon();
186 int n = X.Nrows(); A = X; D.ReSize(n); E.ReSize(n);
187 Real* ei = E.Store() + n;
188 for (int i = n-1; i >= 0; i--)
189 {
190 Real h = 0.0; Real f = - FloatingPointPrecision::Maximum();
191 Real* d = D.Store(); Real* a = A.Store() + (i*(i+1))/2; int k = i;
192 while (k--) { f = *a++; *d++ = f; h += square(f); }
193 if (h <= tol) { REPORT *(--ei) = 0.0; h = 0.0; }
194 else
195 {
196 REPORT
197 Real g = sign(-sqrt(h), f); *(--ei) = g; h -= f*g;
198 f -= g; *(d-1) = f; *(a-1) = f; f = 0.0;
199 Real* dj = D.Store(); Real* ej = E.Store(); int j;
200 for (j = 0; j < i; j++)
201 {
202 Real* dk = D.Store(); Real* ak = A.Store()+(j*(j+1))/2;
203 Real g = 0.0; k = j;
204 while (k--) g += *ak++ * *dk++;
205 k = i-j; int l = j;
206 if (k) for (;;) { g += *ak * *dk++; if (!(--k)) break; ak += ++l; }
207 g /= h; *ej++ = g; f += g * *dj++;
208 }
209 Real hh = f / (2 * h); Real* ak = A.Store();
210 dj = D.Store(); ej = E.Store();
211 for (j = 0; j < i; j++)
212 {
213 f = *dj++; g = *ej - hh * f; *ej++ = g;
214 Real* dk = D.Store(); Real* ek = E.Store(); k = j+1;
215 while (k--) { *ak++ -= (f * *ek++ + g * *dk++); }
216 }
217 }
218 *d = *a; *a = h;
219 }
220 }
221
tql1(DiagonalMatrix & D,DiagonalMatrix & E)222 static void tql1(DiagonalMatrix& D, DiagonalMatrix& E)
223 {
224 Tracer et("Evalue(tql1)");
225 REPORT
226 Real eps = FloatingPointPrecision::Epsilon();
227 int n = D.Nrows(); int l;
228 for (l=1; l<n; l++) E.element(l-1) = E.element(l);
229 Real b = 0.0; Real f = 0.0; E.element(n-1) = 0.0;
230 for (l=0; l<n; l++)
231 {
232 int i,j;
233 Real& dl = D.element(l); Real& el = E.element(l);
234 Real h = eps * ( fabs(dl) + fabs(el) );
235 if (b < h) b = h;
236 int m;
237 for (m=l; m<n; m++) if (fabs(E.element(m)) <= b) break;
238 bool test = false;
239 for (j=0; j<30; j++)
240 {
241 if (m==l) { REPORT test = true; break; }
242 Real& dl1 = D.element(l+1);
243 Real g = dl; Real p = (dl1-g) / (2.0*el); Real r = sqrt(p*p + 1.0);
244 dl = el / (p < 0.0 ? p-r : p+r); Real h = g - dl; f += h;
245 Real* dlx = &dl1; i = n-l-1; while (i--) *dlx++ -= h;
246
247 p = D.element(m); Real c = 1.0; Real s = 0.0;
248 for (i=m-1; i>=l; i--)
249 {
250 Real ei = E.element(i); Real di = D.element(i);
251 Real& ei1 = E.element(i+1);
252 g = c * ei; h = c * p;
253 if ( fabs(p) >= fabs(ei))
254 {
255 REPORT
256 c = ei / p; r = sqrt(c*c + 1.0);
257 ei1 = s*p*r; s = c/r; c = 1.0/r;
258 }
259 else
260 {
261 REPORT
262 c = p / ei; r = sqrt(c*c + 1.0);
263 ei1 = s * ei * r; s = 1.0/r; c /= r;
264 }
265 p = c * di - s*g; D.element(i+1) = h + s * (c*g + s*di);
266 }
267 el = s*p; dl = c*p;
268 if (fabs(el) <= b) { REPORT test = true; break; }
269 }
270 if (!test) Throw ( ConvergenceException(D) );
271 Real p = dl + f;
272 test = false;
273 for (i=l; i>0; i--)
274 {
275 if (p < D.element(i-1)) { REPORT D.element(i) = D.element(i-1); }
276 else { REPORT test = true; break; }
277 }
278 if (!test) i=0;
279 D.element(i) = p;
280 }
281 }
282
EigenValues(const SymmetricMatrix & A,DiagonalMatrix & D,Matrix & Z)283 void EigenValues(const SymmetricMatrix& A, DiagonalMatrix& D, Matrix& Z)
284 { REPORT DiagonalMatrix E; tred2(A, D, E, Z); tql2(D, E, Z); SortSV(D,Z,true); }
285
EigenValues(const SymmetricMatrix & X,DiagonalMatrix & D)286 void EigenValues(const SymmetricMatrix& X, DiagonalMatrix& D)
287 { REPORT DiagonalMatrix E; SymmetricMatrix A; tred3(X,D,E,A); tql1(D,E); }
288
EigenValues(const SymmetricMatrix & X,DiagonalMatrix & D,SymmetricMatrix & A)289 void EigenValues(const SymmetricMatrix& X, DiagonalMatrix& D,
290 SymmetricMatrix& A)
291 { REPORT DiagonalMatrix E; tred3(X,D,E,A); tql1(D,E); }
292
293
294 #ifdef use_namespace
295 }
296 #endif
297
298