|
/Users/lyon/j4p/src/math/linearAlgebra/EigenvalueDecomposition.java
|
1 package math.linearAlgebra;
2
3 import math.linearAlgebra.util.Maths;
4
5 /** Eigenvalues and eigenvectors of a real matrix.
6 <P>
7 If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is
8 diagonal and the eigenvector matrix V is orthogonal.
9 I.e. A = V.times(D.times(V.transpose())) and
10 V.times(V.transpose()) equals the identity matrix.
11 <P>
12 If A is not symmetric, then the eigenvalue matrix D is block diagonal
13 with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues,
14 lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda]. The
15 columns of V represent the eigenvectors in the sense that A*V = V*D,
16 i.e. A.times(V) equals V.times(D). The matrix V may be badly
17 conditioned, or even singular, so the validity of the equation
18 A = V*D*inverse(V) depends upon V.cond().
19 **/
20
21 public class EigenvalueDecomposition implements java.io.Serializable {
22
23 /* ------------------------
24 Class variables
25 * ------------------------ */
26
27 /** Row and column dimension (square matrix).
28 @serial matrix dimension.
29 */
30 private int n;
31
32 /** Symmetry flag.
33 @serial internal symmetry flag.
34 */
35 private boolean issymmetric;
36
37 /** Arrays for internal storage of eigenvalues.
38 @serial internal storage of eigenvalues.
39 */
40 private double[] d, e;
41
42 /** Array for internal storage of eigenvectors.
43 @serial internal storage of eigenvectors.
44 */
45 private double[][] V;
46
47 /** Array for internal storage of nonsymmetric Hessenberg form.
48 @serial internal storage of nonsymmetric Hessenberg form.
49 */
50 private double[][] H;
51
52 /** Working storage for nonsymmetric algorithm.
53 @serial working storage for nonsymmetric algorithm.
54 */
55 private double[] ort;
56
57 /* ------------------------
58 Private Methods
59 * ------------------------ */
60
61 // Symmetric Householder reduction to tridiagonal form.
62
63 private void tred2() {
64
65 // This is derived from the Algol procedures tred2 by
66 // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
67 // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
68 // Fortran subroutine in EISPACK.
69
70 for (int j = 0; j < n; j++) {
71 d[j] = V[n - 1][j];
72 }
73
74 // Householder reduction to tridiagonal form.
75
76 for (int i = n - 1; i > 0; i--) {
77
78 // Scale to avoid under/overflow.
79
80 double scale = 0.0;
81 double h = 0.0;
82 for (int k = 0; k < i; k++) {
83 scale = scale + Math.abs(d[k]);
84 }
85 if (scale == 0.0) {
86 e[i] = d[i - 1];
87 for (int j = 0; j < i; j++) {
88 d[j] = V[i - 1][j];
89 V[i][j] = 0.0;
90 V[j][i] = 0.0;
91 }
92 } else {
93
94 // Generate Householder vector.
95
96 for (int k = 0; k < i; k++) {
97 d[k] /= scale;
98 h += d[k] * d[k];
99 }
100 double f = d[i - 1];
101 double g = Math.sqrt(h);
102 if (f > 0) {
103 g = -g;
104 }
105 e[i] = scale * g;
106 h = h - f * g;
107 d[i - 1] = f - g;
108 for (int j = 0; j < i; j++) {
109 e[j] = 0.0;
110 }
111
112 // Apply similarity transformation to remaining columns.
113
114 for (int j = 0; j < i; j++) {
115 f = d[j];
116 V[j][i] = f;
117 g = e[j] + V[j][j] * f;
118 for (int k = j + 1; k <= i - 1; k++) {
119 g += V[k][j] * d[k];
120 e[k] += V[k][j] * f;
121 }
122 e[j] = g;
123 }
124 f = 0.0;
125 for (int j = 0; j < i; j++) {
126 e[j] /= h;
127 f += e[j] * d[j];
128 }
129 double hh = f / (h + h);
130 for (int j = 0; j < i; j++) {
131 e[j] -= hh * d[j];
132 }
133 for (int j = 0; j < i; j++) {
134 f = d[j];
135 g = e[j];
136 for (int k = j; k <= i - 1; k++) {
137 V[k][j] -= (f * e[k] + g * d[k]);
138 }
139 d[j] = V[i - 1][j];
140 V[i][j] = 0.0;
141 }
142 }
143 d[i] = h;
144 }
145
146 // Accumulate transformations.
147
148 for (int i = 0; i < n - 1; i++) {
149 V[n - 1][i] = V[i][i];
150 V[i][i] = 1.0;
151 double h = d[i + 1];
152 if (h != 0.0) {
153 for (int k = 0; k <= i; k++) {
154 d[k] = V[k][i + 1] / h;
155 }
156 for (int j = 0; j <= i; j++) {
157 double g = 0.0;
158 for (int k = 0; k <= i; k++) {
159 g += V[k][i + 1] * V[k][j];
160 }
161 for (int k = 0; k <= i; k++) {
162 V[k][j] -= g * d[k];
163 }
164 }
165 }
166 for (int k = 0; k <= i; k++) {
167 V[k][i + 1] = 0.0;
168 }
169 }
170 for (int j = 0; j < n; j++) {
171 d[j] = V[n - 1][j];
172 V[n - 1][j] = 0.0;
173 }
174 V[n - 1][n - 1] = 1.0;
175 e[0] = 0.0;
176 }
177
178 // Symmetric tridiagonal QL algorithm.
179
180 private void tql2() {
181
182 // This is derived from the Algol procedures tql2, by
183 // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
184 // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
185 // Fortran subroutine in EISPACK.
186
187 for (int i = 1; i < n; i++) {
188 e[i - 1] = e[i];
189 }
190 e[n - 1] = 0.0;
191
192 double f = 0.0;
193 double tst1 = 0.0;
194 double eps = Math.pow(2.0, -52.0);
195 for (int l = 0; l < n; l++) {
196
197 // Find small subdiagonal element
198
199 tst1 = Math.max(tst1, Math.abs(d[l]) + Math.abs(e[l]));
200 int m = l;
201 while (m < n) {
202 if (Math.abs(e[m]) <= eps * tst1) {
203 break;
204 }
205 m++;
206 }
207
208 // If m == l, d[l] is an eigenvalue,
209 // otherwise, iterate.
210
211 if (m > l) {
212 int iter = 0;
213 do {
214 iter = iter + 1; // (Could check iteration count here.)
215
216 // Compute implicit shift
217
218 double g = d[l];
219 double p = (d[l + 1] - g) / (2.0 * e[l]);
220 double r = Maths.hypot(p, 1.0);
221 if (p < 0) {
222 r = -r;
223 }
224 d[l] = e[l] / (p + r);
225 d[l + 1] = e[l] * (p + r);
226 double dl1 = d[l + 1];
227 double h = g - d[l];
228 for (int i = l + 2; i < n; i++) {
229 d[i] -= h;
230 }
231 f = f + h;
232
233 // Implicit QL transformation.
234
235 p = d[m];
236 double c = 1.0;
237 double c2 = c;
238 double c3 = c;
239 double el1 = e[l + 1];
240 double s = 0.0;
241 double s2 = 0.0;
242 for (int i = m - 1; i >= l; i--) {
243 c3 = c2;
244 c2 = c;
245 s2 = s;
246 g = c * e[i];
247 h = c * p;
248 r = Maths.hypot(p, e[i]);
249 e[i + 1] = s * r;
250 s = e[i] / r;
251 c = p / r;
252 p = c * d[i] - s * g;
253 d[i + 1] = h + s * (c * g + s * d[i]);
254
255 // Accumulate transformation.
256
257 for (int k = 0; k < n; k++) {
258 h = V[k][i + 1];
259 V[k][i + 1] = s * V[k][i] + c * h;
260 V[k][i] = c * V[k][i] - s * h;
261 }
262 }
263 p = -s * s2 * c3 * el1 * e[l] / dl1;
264 e[l] = s * p;
265 d[l] = c * p;
266
267 // Check for convergence.
268
269 } while (Math.abs(e[l]) > eps * tst1);
270 }
271 d[l] = d[l] + f;
272 e[l] = 0.0;
273 }
274
275 // Sort eigenvalues and corresponding vectors.
276
277 for (int i = 0; i < n - 1; i++) {
278 int k = i;
279 double p = d[i];
280 for (int j = i + 1; j < n; j++) {
281 if (d[j] < p) {
282 k = j;
283 p = d[j];
284 }
285 }
286 if (k != i) {
287 d[k] = d[i];
288 d[i] = p;
289 for (int j = 0; j < n; j++) {
290 p = V[j][i];
291 V[j][i] = V[j][k];
292 V[j][k] = p;
293 }
294 }
295 }
296 }
297
298 // Nonsymmetric reduction to Hessenberg form.
299
300 private void orthes() {
301
302 // This is derived from the Algol procedures orthes and ortran,
303 // by Martin and Wilkinson, Handbook for Auto. Comp.,
304 // Vol.ii-Linear Algebra, and the corresponding
305 // Fortran subroutines in EISPACK.
306
307 int low = 0;
308 int high = n - 1;
309
310 for (int m = low + 1; m <= high - 1; m++) {
311
312 // Scale column.
313
314 double scale = 0.0;
315 for (int i = m; i <= high; i++) {
316 scale = scale + Math.abs(H[i][m - 1]);
317 }
318 if (scale != 0.0) {
319
320 // Compute Householder transformation.
321
322 double h = 0.0;
323 for (int i = high; i >= m; i--) {
324 ort[i] = H[i][m - 1] / scale;
325 h += ort[i] * ort[i];
326 }
327 double g = Math.sqrt(h);
328 if (ort[m] > 0) {
329 g = -g;
330 }
331 h = h - ort[m] * g;
332 ort[m] = ort[m] - g;
333
334 // Apply Householder similarity transformation
335 // H = (I-u*u'/h)*H*(I-u*u')/h)
336
337 for (int j = m; j < n; j++) {
338 double f = 0.0;
339 for (int i = high; i >= m; i--) {
340 f += ort[i] * H[i][j];
341 }
342 f = f / h;
343 for (int i = m; i <= high; i++) {
344 H[i][j] -= f * ort[i];
345 }
346 }
347
348 for (int i = 0; i <= high; i++) {
349 double f = 0.0;
350 for (int j = high; j >= m; j--) {
351 f += ort[j] * H[i][j];
352 }
353 f = f / h;
354 for (int j = m; j <= high; j++) {
355 H[i][j] -= f * ort[j];
356 }
357 }
358 ort[m] = scale * ort[m];
359 H[m][m - 1] = scale * g;
360 }
361 }
362
363 // Accumulate transformations (Algol's ortran).
364
365 for (int i = 0; i < n; i++) {
366 for (int j = 0; j < n; j++) {
367 V[i][j] = (i == j ? 1.0 : 0.0);
368 }
369 }
370
371 for (int m = high - 1; m >= low + 1; m--) {
372 if (H[m][m - 1] != 0.0) {
373 for (int i = m + 1; i <= high; i++) {
374 ort[i] = H[i][m - 1];
375 }
376 for (int j = m; j <= high; j++) {
377 double g = 0.0;
378 for (int i = m; i <= high; i++) {
379 g += ort[i] * V[i][j];
380 }
381 // Double division avoids possible underflow
382 g = (g / ort[m]) / H[m][m - 1];
383 for (int i = m; i <= high; i++) {
384 V[i][j] += g * ort[i];
385 }
386 }
387 }
388 }
389 }
390
391
392 // Complex scalar division.
393
394 private transient double cdivr, cdivi;
395
396 private void cdiv(double xr, double xi, double yr, double yi) {
397 double r,d;
398 if (Math.abs(yr) > Math.abs(yi)) {
399 r = yi / yr;
400 d = yr + r * yi;
401 cdivr = (xr + r * xi) / d;
402 cdivi = (xi - r * xr) / d;
403 } else {
404 r = yr / yi;
405 d = yi + r * yr;
406 cdivr = (r * xr + xi) / d;
407 cdivi = (r * xi - xr) / d;
408 }
409 }
410
411
412 // Nonsymmetric reduction from Hessenberg to real Schur form.
413
414 private void hqr2() {
415
416 // This is derived from the Algol procedure hqr2,
417 // by Martin and Wilkinson, Handbook for Auto. Comp.,
418 // Vol.ii-Linear Algebra, and the corresponding
419 // Fortran subroutine in EISPACK.
420
421 // Initialize
422
423 int nn = this.n;
424 int n = nn - 1;
425 int low = 0;
426 int high = nn - 1;
427 double eps = Math.pow(2.0, -52.0);
428 double exshift = 0.0;
429 double p = 0,q = 0,r = 0,s = 0,z = 0,t,w,x,y;
430
431 // Store roots isolated by balanc and compute matrix norm
432
433 double norm = 0.0;
434 for (int i = 0; i < nn; i++) {
435 if (i < low | i > high) {
436 d[i] = H[i][i];
437 e[i] = 0.0;
438 }
439 for (int j = Math.max(i - 1, 0); j < nn; j++) {
440 norm = norm + Math.abs(H[i][j]);
441 }
442 }
443
444 // Outer loop over eigenvalue index
445
446 int iter = 0;
447 while (n >= low) {
448
449 // Look for single small sub-diagonal element
450
451 int l = n;
452 while (l > low) {
453 s = Math.abs(H[l - 1][l - 1]) + Math.abs(H[l][l]);
454 if (s == 0.0) {
455 s = norm;
456 }
457 if (Math.abs(H[l][l - 1]) < eps * s) {
458 break;
459 }
460 l--;
461 }
462
463 // Check for convergence
464 // One root found
465
466 if (l == n) {
467 H[n][n] = H[n][n] + exshift;
468 d[n] = H[n][n];
469 e[n] = 0.0;
470 n--;
471 iter = 0;
472
473 // Two roots found
474
475 } else if (l == n - 1) {
476 w = H[n][n - 1] * H[n - 1][n];
477 p = (H[n - 1][n - 1] - H[n][n]) / 2.0;
478 q = p * p + w;
479 z = Math.sqrt(Math.abs(q));
480 H[n][n] = H[n][n] + exshift;
481 H[n - 1][n - 1] = H[n - 1][n - 1] + exshift;
482 x = H[n][n];
483
484 // Real pair
485
486 if (q >= 0) {
487 if (p >= 0) {
488 z = p + z;
489 } else {
490 z = p - z;
491 }
492 d[n - 1] = x + z;
493 d[n] = d[n - 1];
494 if (z != 0.0) {
495 d[n] = x - w / z;
496 }
497 e[n - 1] = 0.0;
498 e[n] = 0.0;
499 x = H[n][n - 1];
500 s = Math.abs(x) + Math.abs(z);
501 p = x / s;
502 q = z / s;
503 r = Math.sqrt(p * p + q * q);
504 p = p / r;
505 q = q / r;
506
507 // Row modification
508
509 for (int j = n - 1; j < nn; j++) {
510 z = H[n - 1][j];
511 H[n - 1][j] = q * z + p * H[n][j];
512 H[n][j] = q * H[n][j] - p * z;
513 }
514
515 // Column modification
516
517 for (int i = 0; i <= n; i++) {
518 z = H[i][n - 1];
519 H[i][n - 1] = q * z + p * H[i][n];
520 H[i][n] = q * H[i][n] - p * z;
521 }
522
523 // Accumulate transformations
524
525 for (int i = low; i <= high; i++) {
526 z = V[i][n - 1];
527 V[i][n - 1] = q * z + p * V[i][n];
528 V[i][n] = q * V[i][n] - p * z;
529 }
530
531 // Complex pair
532
533 } else {
534 d[n - 1] = x + p;
535 d[n] = x + p;
536 e[n - 1] = z;
537 e[n] = -z;
538 }
539 n = n - 2;
540 iter = 0;
541
542 // No convergence yet
543
544 } else {
545
546 // Form shift
547
548 x = H[n][n];
549 y = 0.0;
550 w = 0.0;
551 if (l < n) {
552 y = H[n - 1][n - 1];
553 w = H[n][n - 1] * H[n - 1][n];
554 }
555
556 // Wilkinson's original ad hoc shift
557
558 if (iter == 10) {
559 exshift += x;
560 for (int i = low; i <= n; i++) {
561 H[i][i] -= x;
562 }
563 s = Math.abs(H[n][n - 1]) + Math.abs(H[n - 1][n - 2]);
564 x = y = 0.75 * s;
565 w = -0.4375 * s * s;
566 }
567
568 // MATLAB's new ad hoc shift
569
570 if (iter == 30) {
571 s = (y - x) / 2.0;
572 s = s * s + w;
573 if (s > 0) {
574 s = Math.sqrt(s);
575 if (y < x) {
576 s = -s;
577 }
578 s = x - w / ((y - x) / 2.0 + s);
579 for (int i = low; i <= n; i++) {
580 H[i][i] -= s;
581 }
582 exshift += s;
583 x = y = w = 0.964;
584 }
585 }
586
587 iter = iter + 1; // (Could check iteration count here.)
588
589 // Look for two consecutive small sub-diagonal elements
590
591 int m = n - 2;
592 while (m >= l) {
593 z = H[m][m];
594 r = x - z;
595 s = y - z;
596 p = (r * s - w) / H[m + 1][m] + H[m][m + 1];
597 q = H[m + 1][m + 1] - z - r - s;
598 r = H[m + 2][m + 1];
599 s = Math.abs(p) + Math.abs(q) + Math.abs(r);
600 p = p / s;
601 q = q / s;
602 r = r / s;
603 if (m == l) {
604 break;
605 }
606 if (Math.abs(H[m][m - 1]) * (Math.abs(q) + Math.abs(r)) <
607 eps * (Math.abs(p) * (Math.abs(H[m - 1][m - 1]) + Math.abs(z) +
608 Math.abs(H[m + 1][m + 1])))) {
609 break;
610 }
611 m--;
612 }
613
614 for (int i = m + 2; i <= n; i++) {
615 H[i][i - 2] = 0.0;
616 if (i > m + 2) {
617 H[i][i - 3] = 0.0;
618 }
619 }
620
621 // Double QR step involving rows l:n and columns m:n
622
623 for (int k = m; k <= n - 1; k++) {
624 boolean notlast = (k != n - 1);
625 if (k != m) {
626 p = H[k][k - 1];
627 q = H[k + 1][k - 1];
628 r = (notlast ? H[k + 2][k - 1] : 0.0);
629 x = Math.abs(p) + Math.abs(q) + Math.abs(r);
630 if (x != 0.0) {
631 p = p / x;
632 q = q / x;
633 r = r / x;
634 }
635 }
636 if (x == 0.0) {
637 break;
638 }
639 s = Math.sqrt(p * p + q * q + r * r);
640 if (p < 0) {
641 s = -s;
642 }
643 if (s != 0) {
644 if (k != m) {
645 H[k][k - 1] = -s * x;
646 } else if (l != m) {
647 H[k][k - 1] = -H[k][k - 1];
648 }
649 p = p + s;
650 x = p / s;
651 y = q / s;
652 z = r / s;
653 q = q / p;
654 r = r / p;
655
656 // Row modification
657
658 for (int j = k; j < nn; j++) {
659 p = H[k][j] + q * H[k + 1][j];
660 if (notlast) {
661 p = p + r * H[k + 2][j];
662 H[k + 2][j] = H[k + 2][j] - p * z;
663 }
664 H[k][j] = H[k][j] - p * x;
665 H[k + 1][j] = H[k + 1][j] - p * y;
666 }
667
668 // Column modification
669
670 for (int i = 0; i <= Math.min(n, k + 3); i++) {
671 p = x * H[i][k] + y * H[i][k + 1];
672 if (notlast) {
673 p = p + z * H[i][k + 2];
674 H[i][k + 2] = H[i][k + 2] - p * r;
675 }
676 H[i][k] = H[i][k] - p;
677 H[i][k + 1] = H[i][k + 1] - p * q;
678 }
679
680 // Accumulate transformations
681
682 for (int i = low; i <= high; i++) {
683 p = x * V[i][k] + y * V[i][k + 1];
684 if (notlast) {
685 p = p + z * V[i][k + 2];
686 V[i][k + 2] = V[i][k + 2] - p * r;
687 }
688 V[i][k] = V[i][k] - p;
689 V[i][k + 1] = V[i][k + 1] - p * q;
690 }
691 } // (s != 0)
692 } // k loop
693 } // check convergence
694 } // while (n >= low)
695
696 // Backsubstitute to find vectors of upper triangular form
697
698 if (norm == 0.0) {
699 return;
700 }
701
702 for (n = nn - 1; n >= 0; n--) {
703 p = d[n];
704 q = e[n];
705
706 // Real vector
707
708 if (q == 0) {
709 int l = n;
710 H[n][n] = 1.0;
711 for (int i = n - 1; i >= 0; i--) {
712 w = H[i][i] - p;
713 r = 0.0;
714 for (int j = l; j <= n; j++) {
715 r = r + H[i][j] * H[j][n];
716 }
717 if (e[i] < 0.0) {
718 z = w;
719 s = r;
720 } else {
721 l = i;
722 if (e[i] == 0.0) {
723 if (w != 0.0) {
724 H[i][n] = -r / w;
725 } else {
726 H[i][n] = -r / (eps * norm);
727 }
728
729 // Solve real equations
730
731 } else {
732 x = H[i][i + 1];
733 y = H[i + 1][i];
734 q = (d[i] - p) * (d[i] - p) + e[i] * e[i];
735 t = (x * s - z * r) / q;
736 H[i][n] = t;
737 if (Math.abs(x) > Math.abs(z)) {
738 H[i + 1][n] = (-r - w * t) / x;
739 } else {
740 H[i + 1][n] = (-s - y * t) / z;
741 }
742 }
743
744 // Overflow control
745
746 t = Math.abs(H[i][n]);
747 if ((eps * t) * t > 1) {
748 for (int j = i; j <= n; j++) {
749 H[j][n] = H[j][n] / t;
750 }
751 }
752 }
753 }
754
755 // Complex vector
756
757 } else if (q < 0) {
758 int l = n - 1;
759
760 // Last vector component imaginary so matrix is triangular
761
762 if (Math.abs(H[n][n - 1]) > Math.abs(H[n - 1][n])) {
763 H[n - 1][n - 1] = q / H[n][n - 1];
764 H[n - 1][n] = -(H[n][n] - p) / H[n][n - 1];
765 } else {
766 cdiv(0.0, -H[n - 1][n], H[n - 1][n - 1] - p, q);
767 H[n - 1][n - 1] = cdivr;
768 H[n - 1][n] = cdivi;
769 }
770 H[n][n - 1] = 0.0;
771 H[n][n] = 1.0;
772 for (int i = n - 2; i >= 0; i--) {
773 double ra,sa,vr,vi;
774 ra = 0.0;
775 sa = 0.0;
776 for (int j = l; j <= n; j++) {
777 ra = ra + H[i][j] * H[j][n - 1];
778 sa = sa + H[i][j] * H[j][n];
779 }
780 w = H[i][i] - p;
781
782 if (e[i] < 0.0) {
783 z = w;
784 r = ra;
785 s = sa;
786 } else {
787 l = i;
788 if (e[i] == 0) {
789 cdiv(-ra, -sa, w, q);
790 H[i][n - 1] = cdivr;
791 H[i][n] = cdivi;
792 } else {
793
794 // Solve complex equations
795
796 x = H[i][i + 1];
797 y = H[i + 1][i];
798 vr = (d[i] - p) * (d[i] - p) + e[i] * e[i] - q * q;
799 vi = (d[i] - p) * 2.0 * q;
800 if (vr == 0.0 & vi == 0.0) {
801 vr = eps * norm * (Math.abs(w) + Math.abs(q) +
802 Math.abs(x) + Math.abs(y) + Math.abs(z));
803 }
804 cdiv(x * r - z * ra + q * sa, x * s - z * sa - q * ra, vr, vi);
805 H[i][n - 1] = cdivr;
806 H[i][n] = cdivi;
807 if (Math.abs(x) > (Math.abs(z) + Math.abs(q))) {
808 H[i + 1][n - 1] = (-ra - w * H[i][n - 1] + q * H[i][n]) / x;
809 H[i + 1][n] = (-sa - w * H[i][n] - q * H[i][n - 1]) / x;
810 } else {
811 cdiv(-r - y * H[i][n - 1], -s - y * H[i][n], z, q);
812 H[i + 1][n - 1] = cdivr;
813 H[i + 1][n] = cdivi;
814 }
815 }
816
817 // Overflow control
818
819 t = Math.max(Math.abs(H[i][n - 1]), Math.abs(H[i][n]));
820 if ((eps * t) * t > 1) {
821 for (int j = i; j <= n; j++) {
822 H[j][n - 1] = H[j][n - 1] / t;
823 H[j][n] = H[j][n] / t;
824 }
825 }
826 }
827 }
828 }
829 }
830
831 // Vectors of isolated roots
832
833 for (int i = 0; i < nn; i++) {
834 if (i < low | i > high) {
835 for (int j = i; j < nn; j++) {
836 V[i][j] = H[i][j];
837 }
838 }
839 }
840
841 // Back transformation to get eigenvectors of original matrix
842
843 for (int j = nn - 1; j >= low; j--) {
844 for (int i = low; i <= high; i++) {
845 z = 0.0;
846 for (int k = low; k <= Math.min(j, high); k++) {
847 z = z + V[i][k] * H[k][j];
848 }
849 V[i][j] = z;
850 }
851 }
852 }
853
854
855 /* ------------------------
856 Constructor
857 * ------------------------ */
858
859 /** Check for symmetry, then construct the eigenvalue decomposition
860 @param A Square matrix
861 @return Structure to access D and V.
862 */
863
864 public EigenvalueDecomposition(Matrix Arg) {
865 double[][] A = Arg.getArray();
866 n = Arg.getColumnDimension();
867 V = new double[n][n];
868 d = new double[n];
869 e = new double[n];
870
871 issymmetric = true;
872 for (int j = 0; (j < n) & issymmetric; j++) {
873 for (int i = 0; (i < n) & issymmetric; i++) {
874 issymmetric = (A[i][j] == A[j][i]);
875 }
876 }
877
878 if (issymmetric) {
879 for (int i = 0; i < n; i++) {
880 for (int j = 0; j < n; j++) {
881 V[i][j] = A[i][j];
882 }
883 }
884
885 // Tridiagonalize.
886 tred2();
887
888 // Diagonalize.
889 tql2();
890
891 } else {
892 H = new double[n][n];
893 ort = new double[n];
894
895 for (int j = 0; j < n; j++) {
896 for (int i = 0; i < n; i++) {
897 H[i][j] = A[i][j];
898 }
899 }
900
901 // Reduce to Hessenberg form.
902 orthes();
903
904 // Reduce Hessenberg to real Schur form.
905 hqr2();
906 }
907 }
908
909 /* ------------------------
910 Public Methods
911 * ------------------------ */
912
913 /** Return the eigenvector matrix
914 @return V
915 */
916
917 public Matrix getV() {
918 return new Matrix(V, n, n);
919 }
920
921 /** Return the real parts of the eigenvalues
922 @return real(diag(D))
923 */
924
925 public double[] getRealEigenvalues() {
926 return d;
927 }
928
929 /** Return the imaginary parts of the eigenvalues
930 @return imag(diag(D))
931 */
932
933 public double[] getImagEigenvalues() {
934 return e;
935 }
936
937 /** Return the block diagonal eigenvalue matrix
938 @return D
939 */
940
941 public Matrix getD() {
942 Matrix X = new Matrix(n, n);
943 double[][] D = X.getArray();
944 for (int i = 0; i < n; i++) {
945 for (int j = 0; j < n; j++) {
946 D[i][j] = 0.0;
947 }
948 D[i][i] = d[i];
949 if (e[i] > 0) {
950 D[i][i + 1] = e[i];
951 } else if (e[i] < 0) {
952 D[i][i - 1] = e[i];
953 }
954 }
955 return X;
956 }
957 }
958