Source code: jmat/data/matrixDecompositions/EigenvalueDecomposition.java
1 package jmat.data.matrixDecompositions;
2
3 import jmat.data.Matrix;
4
5
6 /** Eigenvalues and eigenvectors of a real matrix.
7 <P>
8 If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is
9 diagonal and the eigenvector matrix V is orthogonal.
10 I.e. A = V.times(D.times(V.transpose())) and
11 V.times(V.transpose()) equals the identity matrix.
12 <P>
13 If A is not symmetric, then the eigenvalue matrix D is block diagonal
14 with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues,
15 lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda]. The
16 columns of V represent the eigenvectors in the sense that A*V = V*D,
17 i.e. A.times(V) equals V.times(D). The matrix V may be badly
18 conditioned, or even singular, so the validity of the equation
19 A = V*D*inverse(V) depends upon V.cond().
20 **/
21 public class EigenvalueDecomposition implements java.io.Serializable
22 {
23 //~ Instance fields ////////////////////////////////////////////////////////
24
25 /** Array for internal storage of nonsymmetric Hessenberg form.
26 @serial internal storage of nonsymmetric Hessenberg form.
27 */
28 private double[][] H;
29
30 /** Array for internal storage of eigenvectors.
31 @serial internal storage of eigenvectors.
32 */
33 private double[][] V;
34
35 /** Arrays for internal storage of eigenvalues.
36 @serial internal storage of eigenvalues.
37 */
38 private double[] d;
39
40 /** Arrays for internal storage of eigenvalues.
41 @serial internal storage of eigenvalues.
42 */
43 private double[] e;
44
45 /** Working storage for nonsymmetric algorithm.
46 @serial working storage for nonsymmetric algorithm.
47 */
48 private double[] ort;
49
50 /** Symmetry flag.
51 @serial internal symmetry flag.
52 */
53 private boolean issymmetric;
54
55 // Complex scalar division.
56 private transient double cdivi;
57
58 // Complex scalar division.
59 private transient double cdivr;
60
61 /* ------------------------
62 Class variables
63 * ------------------------ */
64
65 /** Row and column dimension (square matrix).
66 @serial matrix dimension.
67 */
68 private int n;
69
70 //~ Constructors ///////////////////////////////////////////////////////////
71
72 /* ------------------------
73 Constructor
74 * ------------------------ */
75
76 /** Check for symmetry, then construct the eigenvalue decomposition
77 @param Arg Square matrix
78 @return Structure to access D and V.
79 */
80 public EigenvalueDecomposition(Matrix Arg)
81 {
82 double[][] A = Arg.getArray();
83 n = Arg.getColumnDimension();
84 V = new double[n][n];
85 d = new double[n];
86 e = new double[n];
87
88 issymmetric = true;
89
90 for (int j = 0; (j < n) & issymmetric; j++)
91 {
92 for (int i = 0; (i < n) & issymmetric; i++)
93 {
94 issymmetric = (A[i][j] == A[j][i]);
95 }
96 }
97
98 if (issymmetric)
99 {
100 for (int i = 0; i < n; i++)
101 {
102 for (int j = 0; j < n; j++)
103 {
104 V[i][j] = A[i][j];
105 }
106 }
107
108 // Tridiagonalize.
109 tred2();
110
111 // Diagonalize.
112 tql2();
113 }
114 else
115 {
116 H = new double[n][n];
117 ort = new double[n];
118
119 for (int j = 0; j < n; j++)
120 {
121 for (int i = 0; i < n; i++)
122 {
123 H[i][j] = A[i][j];
124 }
125 }
126
127 // Reduce to Hessenberg form.
128 orthes();
129
130 // Reduce Hessenberg to real Schur form.
131 hqr2();
132 }
133 }
134
135 //~ Methods ////////////////////////////////////////////////////////////////
136
137 /** Return the block diagonal eigenvalue matrix
138 @return D
139 */
140 public Matrix getD()
141 {
142 Matrix X = new Matrix(n, n);
143 double[][] D = X.getArray();
144
145 for (int i = 0; i < n; i++)
146 {
147 for (int j = 0; j < n; j++)
148 {
149 D[i][j] = 0.0;
150 }
151
152 D[i][i] = d[i];
153
154 if (e[i] > 0)
155 {
156 D[i][i + 1] = e[i];
157 }
158 else if (e[i] < 0)
159 {
160 D[i][i - 1] = e[i];
161 }
162 }
163
164 return X;
165 }
166
167 /* ------------------------
168 Public Methods
169 * ------------------------ */
170
171 /** Return the eigenvector matrix
172 @return V
173 */
174 public boolean isDReal()
175 {
176 boolean b = false;
177
178 for (int i = 0; i < e.length; i++)
179 {
180 if (e[i] > 0)
181 {
182 b = true;
183 }
184 }
185
186 return b;
187 }
188
189 /** Return the imaginary diagonal eigenvalue matrix
190 @return Imag(D)
191 */
192 public Matrix getD_Imag()
193 {
194 Matrix X = new Matrix(n, n);
195 double[][] E = X.getArray();
196
197 for (int i = 0; i < n; i++)
198 {
199 for (int j = 0; j < n; j++)
200 {
201 E[i][j] = 0.0;
202 }
203
204 E[i][i] = e[i];
205 }
206
207 return X;
208 }
209
210 /** Return the real diagonal eigenvalue matrix
211 @return Real(D)
212 */
213 public Matrix getD_Real()
214 {
215 Matrix X = new Matrix(n, n);
216 double[][] D = X.getArray();
217
218 for (int i = 0; i < n; i++)
219 {
220 for (int j = 0; j < n; j++)
221 {
222 D[i][j] = 0.0;
223 }
224
225 D[i][i] = d[i];
226 }
227
228 return X;
229 }
230
231 public double[] getRealEigenvalues()
232 {
233 double[] D = new double[n];
234
235 for (int i = 0; i < n; i++)
236 {
237 D[i] = d[i];
238 }
239
240 return D;
241 }
242
243 /** Return the eigenvector matrix
244 @return V
245 */
246 public Matrix getV()
247 {
248 return new Matrix(V, n, n);
249 }
250
251 private void cdiv(double xr, double xi, double yr, double yi)
252 {
253 double r;
254 double d;
255
256 if (Math.abs(yr) > Math.abs(yi))
257 {
258 r = yi / yr;
259 d = yr + (r * yi);
260 cdivr = (xr + (r * xi)) / d;
261 cdivi = (xi - (r * xr)) / d;
262 }
263 else
264 {
265 r = yr / yi;
266 d = yi + (r * yr);
267 cdivr = ((r * xr) + xi) / d;
268 cdivi = ((r * xi) - xr) / d;
269 }
270 }
271
272 // Nonsymmetric reduction from Hessenberg to real Schur form.
273 private void hqr2()
274 {
275 // This is derived from the Algol procedure hqr2,
276 // by Martin and Wilkinson, Handbook for Auto. Comp.,
277 // Vol.ii-Linear Algebra, and the corresponding
278 // Fortran subroutine in EISPACK.
279 // Initialize
280 int nn = this.n;
281 int n = nn - 1;
282 int low = 0;
283 int high = nn - 1;
284 double eps = Math.pow(2.0, -52.0);
285 double exshift = 0.0;
286 double p = 0;
287 double q = 0;
288 double r = 0;
289 double s = 0;
290 double z = 0;
291 double t;
292 double w;
293 double x;
294 double y;
295
296 // Store roots isolated by balanc and compute matrix norm
297 double norm = 0.0;
298
299 for (int i = 0; i < nn; i++)
300 {
301 if ((i < low) | (i > high))
302 {
303 d[i] = H[i][i];
304 e[i] = 0.0;
305 }
306
307 for (int j = Math.max(i - 1, 0); j < nn; j++)
308 {
309 norm = norm + Math.abs(H[i][j]);
310 }
311 }
312
313 // Outer loop over eigenvalue index
314 int iter = 0;
315
316 while (n >= low)
317 {
318 // Look for single small sub-diagonal element
319 int l = n;
320
321 while (l > low)
322 {
323 s = Math.abs(H[l - 1][l - 1]) + Math.abs(H[l][l]);
324
325 if (s == 0.0)
326 {
327 s = norm;
328 }
329
330 if (Math.abs(H[l][l - 1]) < (eps * s))
331 {
332 break;
333 }
334
335 l--;
336 }
337
338 // Check for convergence
339 // One root found
340 if (l == n)
341 {
342 H[n][n] = H[n][n] + exshift;
343 d[n] = H[n][n];
344 e[n] = 0.0;
345 n--;
346 iter = 0;
347
348 // Two roots found
349 }
350 else if (l == (n - 1))
351 {
352 w = H[n][n - 1] * H[n - 1][n];
353 p = (H[n - 1][n - 1] - H[n][n]) / 2.0;
354 q = (p * p) + w;
355 z = Math.sqrt(Math.abs(q));
356 H[n][n] = H[n][n] + exshift;
357 H[n - 1][n - 1] = H[n - 1][n - 1] + exshift;
358 x = H[n][n];
359
360 // Real pair
361 if (q >= 0)
362 {
363 if (p >= 0)
364 {
365 z = p + z;
366 }
367 else
368 {
369 z = p - z;
370 }
371
372 d[n - 1] = x + z;
373 d[n] = d[n - 1];
374
375 if (z != 0.0)
376 {
377 d[n] = x - (w / z);
378 }
379
380 e[n - 1] = 0.0;
381 e[n] = 0.0;
382 x = H[n][n - 1];
383 s = Math.abs(x) + Math.abs(z);
384 p = x / s;
385 q = z / s;
386 r = Math.sqrt((p * p) + (q * q));
387 p = p / r;
388 q = q / r;
389
390 // Row modification
391 for (int j = n - 1; j < nn; j++)
392 {
393 z = H[n - 1][j];
394 H[n - 1][j] = (q * z) + (p * H[n][j]);
395 H[n][j] = (q * H[n][j]) - (p * z);
396 }
397
398 // Column modification
399 for (int i = 0; i <= n; i++)
400 {
401 z = H[i][n - 1];
402 H[i][n - 1] = (q * z) + (p * H[i][n]);
403 H[i][n] = (q * H[i][n]) - (p * z);
404 }
405
406 // Accumulate transformations
407 for (int i = low; i <= high; i++)
408 {
409 z = V[i][n - 1];
410 V[i][n - 1] = (q * z) + (p * V[i][n]);
411 V[i][n] = (q * V[i][n]) - (p * z);
412 }
413
414 // Complex pair
415 }
416 else
417 {
418 d[n - 1] = x + p;
419 d[n] = x + p;
420 e[n - 1] = z;
421 e[n] = -z;
422 }
423
424 n = n - 2;
425 iter = 0;
426
427 // No convergence yet
428 }
429 else
430 {
431 // Form shift
432 x = H[n][n];
433 y = 0.0;
434 w = 0.0;
435
436 if (l < n)
437 {
438 y = H[n - 1][n - 1];
439 w = H[n][n - 1] * H[n - 1][n];
440 }
441
442 // Wilkinson's original ad hoc shift
443 if (iter == 10)
444 {
445 exshift += x;
446
447 for (int i = low; i <= n; i++)
448 {
449 H[i][i] -= x;
450 }
451
452 s = Math.abs(H[n][n - 1]) + Math.abs(H[n - 1][n - 2]);
453 x = y = 0.75 * s;
454 w = -0.4375 * s * s;
455 }
456
457 // MATLAB's new ad hoc shift
458 if (iter == 30)
459 {
460 s = (y - x) / 2.0;
461 s = (s * s) + w;
462
463 if (s > 0)
464 {
465 s = Math.sqrt(s);
466
467 if (y < x)
468 {
469 s = -s;
470 }
471
472 s = x - (w / (((y - x) / 2.0) + s));
473
474 for (int i = low; i <= n; i++)
475 {
476 H[i][i] -= s;
477 }
478
479 exshift += s;
480 x = y = w = 0.964;
481 }
482 }
483
484 iter = iter + 1; // (Could check iteration count here.)
485
486 // Look for two consecutive small sub-diagonal elements
487 int m = n - 2;
488
489 while (m >= l)
490 {
491 z = H[m][m];
492 r = x - z;
493 s = y - z;
494 p = (((r * s) - w) / H[m + 1][m]) + H[m][m + 1];
495 q = H[m + 1][m + 1] - z - r - s;
496 r = H[m + 2][m + 1];
497 s = Math.abs(p) + Math.abs(q) + Math.abs(r);
498 p = p / s;
499 q = q / s;
500 r = r / s;
501
502 if (m == l)
503 {
504 break;
505 }
506
507 if ((Math.abs(H[m][m - 1]) * (Math.abs(q) + Math.abs(r))) < (eps * (Math.abs(
508 p) * (Math.abs(H[m - 1][m - 1]) + Math.abs(z) +
509 Math.abs(H[m + 1][m + 1])))))
510 {
511 break;
512 }
513
514 m--;
515 }
516
517 for (int i = m + 2; i <= n; i++)
518 {
519 H[i][i - 2] = 0.0;
520
521 if (i > (m + 2))
522 {
523 H[i][i - 3] = 0.0;
524 }
525 }
526
527 // Double QR step involving rows l:n and columns m:n
528 for (int k = m; k <= (n - 1); k++)
529 {
530 boolean notlast = (k != (n - 1));
531
532 if (k != m)
533 {
534 p = H[k][k - 1];
535 q = H[k + 1][k - 1];
536 r = (notlast ? H[k + 2][k - 1] : 0.0);
537 x = Math.abs(p) + Math.abs(q) + Math.abs(r);
538
539 if (x != 0.0)
540 {
541 p = p / x;
542 q = q / x;
543 r = r / x;
544 }
545 }
546
547 if (x == 0.0)
548 {
549 break;
550 }
551
552 s = Math.sqrt((p * p) + (q * q) + (r * r));
553
554 if (p < 0)
555 {
556 s = -s;
557 }
558
559 if (s != 0)
560 {
561 if (k != m)
562 {
563 H[k][k - 1] = -s * x;
564 }
565 else if (l != m)
566 {
567 H[k][k - 1] = -H[k][k - 1];
568 }
569
570 p = p + s;
571 x = p / s;
572 y = q / s;
573 z = r / s;
574 q = q / p;
575 r = r / p;
576
577 // Row modification
578 for (int j = k; j < nn; j++)
579 {
580 p = H[k][j] + (q * H[k + 1][j]);
581
582 if (notlast)
583 {
584 p = p + (r * H[k + 2][j]);
585 H[k + 2][j] = H[k + 2][j] - (p * z);
586 }
587
588 H[k][j] = H[k][j] - (p * x);
589 H[k + 1][j] = H[k + 1][j] - (p * y);
590 }
591
592 // Column modification
593 for (int i = 0; i <= Math.min(n, k + 3); i++)
594 {
595 p = (x * H[i][k]) + (y * H[i][k + 1]);
596
597 if (notlast)
598 {
599 p = p + (z * H[i][k + 2]);
600 H[i][k + 2] = H[i][k + 2] - (p * r);
601 }
602
603 H[i][k] = H[i][k] - p;
604 H[i][k + 1] = H[i][k + 1] - (p * q);
605 }
606
607 // Accumulate transformations
608 for (int i = low; i <= high; i++)
609 {
610 p = (x * V[i][k]) + (y * V[i][k + 1]);
611
612 if (notlast)
613 {
614 p = p + (z * V[i][k + 2]);
615 V[i][k + 2] = V[i][k + 2] - (p * r);
616 }
617
618 V[i][k] = V[i][k] - p;
619 V[i][k + 1] = V[i][k + 1] - (p * q);
620 }
621 }
622
623 // (s != 0)
624 }
625
626 // k loop
627 }
628
629 // check convergence
630 }
631
632 // while (n >= low)
633 // Backsubstitute to find vectors of upper triangular form
634 if (norm == 0.0)
635 {
636 return;
637 }
638
639 for (n = nn - 1; n >= 0; n--)
640 {
641 p = d[n];
642 q = e[n];
643
644 // Real vector
645 if (q == 0)
646 {
647 int l = n;
648 H[n][n] = 1.0;
649
650 for (int i = n - 1; i >= 0; i--)
651 {
652 w = H[i][i] - p;
653 r = 0.0;
654
655 for (int j = l; j <= n; j++)
656 {
657 r = r + (H[i][j] * H[j][n]);
658 }
659
660 if (e[i] < 0.0)
661 {
662 z = w;
663 s = r;
664 }
665 else
666 {
667 l = i;
668
669 if (e[i] == 0.0)
670 {
671 if (w != 0.0)
672 {
673 H[i][n] = -r / w;
674 }
675 else
676 {
677 H[i][n] = -r / (eps * norm);
678 }
679
680 // Solve real equations
681 }
682 else
683 {
684 x = H[i][i + 1];
685 y = H[i + 1][i];
686 q = ((d[i] - p) * (d[i] - p)) + (e[i] * e[i]);
687 t = ((x * s) - (z * r)) / q;
688 H[i][n] = t;
689
690 if (Math.abs(x) > Math.abs(z))
691 {
692 H[i + 1][n] = (-r - (w * t)) / x;
693 }
694 else
695 {
696 H[i + 1][n] = (-s - (y * t)) / z;
697 }
698 }
699
700 // Overflow control
701 t = Math.abs(H[i][n]);
702
703 if (((eps * t) * t) > 1)
704 {
705 for (int j = i; j <= n; j++)
706 {
707 H[j][n] = H[j][n] / t;
708 }
709 }
710 }
711 }
712
713 // Complex vector
714 }
715 else if (q < 0)
716 {
717 int l = n - 1;
718
719 // Last vector component imaginary so matrix is triangular
720 if (Math.abs(H[n][n - 1]) > Math.abs(H[n - 1][n]))
721 {
722 H[n - 1][n - 1] = q / H[n][n - 1];
723 H[n - 1][n] = -(H[n][n] - p) / H[n][n - 1];
724 }
725 else
726 {
727 cdiv(0.0, -H[n - 1][n], H[n - 1][n - 1] - p, q);
728 H[n - 1][n - 1] = cdivr;
729 H[n - 1][n] = cdivi;
730 }
731
732 H[n][n - 1] = 0.0;
733 H[n][n] = 1.0;
734
735 for (int i = n - 2; i >= 0; i--)
736 {
737 double ra;
738 double sa;
739 double vr;
740 double vi;
741 ra = 0.0;
742 sa = 0.0;
743
744 for (int j = l; j <= n; j++)
745 {
746 ra = ra + (H[i][j] * H[j][n - 1]);
747 sa = sa + (H[i][j] * H[j][n]);
748 }
749
750 w = H[i][i] - p;
751
752 if (e[i] < 0.0)
753 {
754 z = w;
755 r = ra;
756 s = sa;
757 }
758 else
759 {
760 l = i;
761
762 if (e[i] == 0)
763 {
764 cdiv(-ra, -sa, w, q);
765 H[i][n - 1] = cdivr;
766 H[i][n] = cdivi;
767 }
768 else
769 {
770 // Solve complex equations
771 x = H[i][i + 1];
772 y = H[i + 1][i];
773 vr = (((d[i] - p) * (d[i] - p)) + (e[i] * e[i])) -
774 (q * q);
775 vi = (d[i] - p) * 2.0 * q;
776
777 if ((vr == 0.0) & (vi == 0.0))
778 {
779 vr = eps * norm * (Math.abs(w) + Math.abs(q) +
780 Math.abs(x) + Math.abs(y) + Math.abs(z));
781 }
782
783 cdiv((x * r) - (z * ra) + (q * sa),
784 (x * s) - (z * sa) - (q * ra), vr, vi);
785 H[i][n - 1] = cdivr;
786 H[i][n] = cdivi;
787
788 if (Math.abs(x) > (Math.abs(z) + Math.abs(q)))
789 {
790 H[i + 1][n - 1] = (-ra - (w * H[i][n - 1]) +
791 (q * H[i][n])) / x;
792 H[i + 1][n] = (-sa - (w * H[i][n]) -
793 (q * H[i][n - 1])) / x;
794 }
795 else
796 {
797 cdiv(-r - (y * H[i][n - 1]),
798 -s - (y * H[i][n]), z, q);
799 H[i + 1][n - 1] = cdivr;
800 H[i + 1][n] = cdivi;
801 }
802 }
803
804 // Overflow control
805 t = Math.max(Math.abs(H[i][n - 1]), Math.abs(H[i][n]));
806
807 if (((eps * t) * t) > 1)
808 {
809 for (int j = i; j <= n; j++)
810 {
811 H[j][n - 1] = H[j][n - 1] / t;
812 H[j][n] = H[j][n] / t;
813 }
814 }
815 }
816 }
817 }
818 }
819
820 // Vectors of isolated roots
821 for (int i = 0; i < nn; i++)
822 {
823 if ((i < low) | (i > high))
824 {
825 for (int j = i; j < nn; j++)
826 {
827 V[i][j] = H[i][j];
828 }
829 }
830 }
831
832 // Back transformation to get eigenvectors of original matrix
833 for (int j = nn - 1; j >= low; j--)
834 {
835 for (int i = low; i <= high; i++)
836 {
837 z = 0.0;
838
839 for (int k = low; k <= Math.min(j, high); k++)
840 {
841 z = z + (V[i][k] * H[k][j]);
842 }
843
844 V[i][j] = z;
845 }
846 }
847 }
848
849 // Nonsymmetric reduction to Hessenberg form.
850 private void orthes()
851 {
852 // This is derived from the Algol procedures orthes and ortran,
853 // by Martin and Wilkinson, Handbook for Auto. Comp.,
854 // Vol.ii-Linear Algebra, and the corresponding
855 // Fortran subroutines in EISPACK.
856 int low = 0;
857 int high = n - 1;
858
859 for (int m = low + 1; m <= (high - 1); m++)
860 {
861 // Scale column.
862 double scale = 0.0;
863
864 for (int i = m; i <= high; i++)
865 {
866 scale = scale + Math.abs(H[i][m - 1]);
867 }
868
869 if (scale != 0.0)
870 {
871 // Compute Householder transformation.
872 double h = 0.0;
873
874 for (int i = high; i >= m; i--)
875 {
876 ort[i] = H[i][m - 1] / scale;
877 h += (ort[i] * ort[i]);
878 }
879
880 double g = Math.sqrt(h);
881
882 if (ort[m] > 0)
883 {
884 g = -g;
885 }
886
887 h = h - (ort[m] * g);
888 ort[m] = ort[m] - g;
889
890 // Apply Householder similarity transformation
891 // H = (I-u*u'/h)*H*(I-u*u')/h)
892 for (int j = m; j < n; j++)
893 {
894 double f = 0.0;
895
896 for (int i = high; i >= m; i--)
897 {
898 f += (ort[i] * H[i][j]);
899 }
900
901 f = f / h;
902
903 for (int i = m; i <= high; i++)
904 {
905 H[i][j] -= (f * ort[i]);
906 }
907 }
908
909 for (int i = 0; i <= high; i++)
910 {
911 double f = 0.0;
912
913 for (int j = high; j >= m; j--)
914 {
915 f += (ort[j] * H[i][j]);
916 }
917
918 f = f / h;
919
920 for (int j = m; j <= high; j++)
921 {
922 H[i][j] -= (f * ort[j]);
923 }
924 }
925
926 ort[m] = scale * ort[m];
927 H[m][m - 1] = scale * g;
928 }
929 }
930
931 // Accumulate transformations (Algol's ortran).
932 for (int i = 0; i < n; i++)
933 {
934 for (int j = 0; j < n; j++)
935 {
936 V[i][j] = ((i == j) ? 1.0 : 0.0);
937 }
938 }
939
940 for (int m = high - 1; m >= (low + 1); m--)
941 {
942 if (H[m][m - 1] != 0.0)
943 {
944 for (int i = m + 1; i <= high; i++)
945 {
946 ort[i] = H[i][m - 1];
947 }
948
949 for (int j = m; j <= high; j++)
950 {
951 double g = 0.0;
952
953 for (int i = m; i <= high; i++)
954 {
955 g += (ort[i] * V[i][j]);
956 }
957
958 // Double division avoids possible underflow
959 g = (g / ort[m]) / H[m][m - 1];
960
961 for (int i = m; i <= high; i++)
962 {
963 V[i][j] += (g * ort[i]);
964 }
965 }
966 }
967 }
968 }
969
970 // Symmetric tridiagonal QL algorithm.
971 private void tql2()
972 {
973 // This is derived from the Algol procedures tql2, by
974 // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
975 // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
976 // Fortran subroutine in EISPACK.
977 for (int i = 1; i < n; i++)
978 {
979 e[i - 1] = e[i];
980 }
981
982 e[n - 1] = 0.0;
983
984 double f = 0.0;
985 double tst1 = 0.0;
986 double eps = Math.pow(2.0, -52.0);
987
988 for (int l = 0; l < n; l++)
989 {
990 // Find small subdiagonal element
991 tst1 = Math.max(tst1, Math.abs(d[l]) + Math.abs(e[l]));
992
993 int m = l;
994
995 while (m < n)
996 {
997 if (Math.abs(e[m]) <= (eps * tst1))
998 {
999 break;
1000 }
1001
1002 m++;
1003 }
1004
1005 // If m == l, d[l] is an eigenvalue,
1006 // otherwise, iterate.
1007 if (m > l)
1008 {
1009 int iter = 0;
1010
1011 do
1012 {
1013 iter = iter + 1; // (Could check iteration count here.)
1014
1015 // Compute implicit shift
1016 double g = d[l];
1017 double p = (d[l + 1] - g) / (2.0 * e[l]);
1018 double r = Mathfun.hypot(p, 1.0);
1019
1020 if (p < 0)
1021 {
1022 r = -r;
1023 }
1024
1025 d[l] = e[l] / (p + r);
1026 d[l + 1] = e[l] * (p + r);
1027
1028 double dl1 = d[l + 1];
1029 double h = g - d[l];
1030
1031 for (int i = l + 2; i < n; i++)
1032 {
1033 d[i] -= h;
1034 }
1035
1036 f = f + h;
1037
1038 // Implicit QL transformation.
1039 p = d[m];
1040
1041 double c = 1.0;
1042 double c2 = c;
1043 double c3 = c;
1044 double el1 = e[l + 1];
1045 double s = 0.0;
1046 double s2 = 0.0;
1047
1048 for (int i = m - 1; i >= l; i--)
1049 {
1050 c3 = c2;
1051 c2 = c;
1052 s2 = s;
1053 g = c * e[i];
1054 h = c * p;
1055 r = Mathfun.hypot(p, e[i]);
1056 e[i + 1] = s * r;
1057 s = e[i] / r;
1058 c = p / r;
1059 p = (c * d[i]) - (s * g);
1060 d[i + 1] = h + (s * ((c * g) + (s * d[i])));
1061
1062 // Accumulate transformation.
1063 for (int k = 0; k < n; k++)
1064 {
1065 h = V[k][i + 1];
1066 V[k][i + 1] = (s * V[k][i]) + (c * h);
1067 V[k][i] = (c * V[k][i]) - (s * h);
1068 }
1069 }
1070
1071 p = (-s * s2 * c3 * el1 * e[l]) / dl1;
1072 e[l] = s * p;
1073 d[l] = c * p;
1074
1075 // Check for convergence.
1076 }
1077 while (Math.abs(e[l]) > (eps * tst1));
1078 }
1079
1080 d[l] = d[l] + f;
1081 e[l] = 0.0;
1082 }
1083
1084 // Sort eigenvalues and corresponding vectors.
1085 for (int i = 0; i < (n - 1); i++)
1086 {
1087 int k = i;
1088 double p = d[i];
1089
1090 for (int j = i + 1; j < n; j++)
1091 {
1092 if (d[j] < p)
1093 {
1094 k = j;
1095 p = d[j];
1096 }
1097 }
1098
1099 if (k != i)
1100 {
1101 d[k] = d[i];
1102 d[i] = p;
1103
1104 for (int j = 0; j < n; j++)
1105 {
1106 p = V[j][i];
1107 V[j][i] = V[j][k];
1108 V[j][k] = p;
1109 }
1110 }
1111 }
1112 }
1113
1114 /* ------------------------
1115 Private Methods
1116 * ------------------------ */
1117
1118 // Symmetric Householder reduction to tridiagonal form.
1119 private void tred2()
1120 {
1121 // This is derived from the Algol procedures tred2 by
1122 // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
1123 // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
1124 // Fortran subroutine in EISPACK.
1125 for (int j = 0; j < n; j++)
1126 {
1127 d[j] = V[n - 1][j];
1128 }
1129
1130 // Householder reduction to tridiagonal form.
1131 for (int i = n - 1; i > 0; i--)
1132 {
1133 // Scale to avoid under/overflow.
1134 double scale = 0.0;
1135 double h = 0.0;
1136
1137 for (int k = 0; k < i; k++)
1138 {
1139 scale = scale + Math.abs(d[k]);
1140 }
1141
1142 if (scale == 0.0)
1143 {
1144 e[i] = d[i - 1];
1145
1146 for (int j = 0; j < i; j++)
1147 {
1148 d[j] = V[i - 1][j];
1149 V[i][j] = 0.0;
1150 V[j][i] = 0.0;
1151 }
1152 }
1153 else
1154 {
1155 // Generate Householder vector.
1156 for (int k = 0; k < i; k++)
1157 {
1158 d[k] /= scale;
1159 h += (d[k] * d[k]);
1160 }
1161
1162 double f = d[i - 1];
1163 double g = Math.sqrt(h);
1164
1165 if (f > 0)
1166 {
1167 g = -g;
1168 }
1169
1170 e[i] = scale * g;
1171 h = h - (f * g);
1172 d[i - 1] = f - g;
1173
1174 for (int j = 0; j < i; j++)
1175 {
1176 e[j] = 0.0;
1177 }
1178
1179 // Apply similarity transformation to remaining columns.
1180 for (int j = 0; j < i; j++)
1181 {
1182 f = d[j];
1183 V[j][i] = f;
1184 g = e[j] + (V[j][j] * f);
1185
1186 for (int k = j + 1; k <= (i - 1); k++)
1187 {
1188 g += (V[k][j] * d[k]);
1189 e[k] += (V[k][j] * f);
1190 }
1191
1192 e[j] = g;
1193 }
1194
1195 f = 0.0;
1196
1197 for (int j = 0; j < i; j++)
1198 {
1199 e[j] /= h;
1200 f += (e[j] * d[j]);
1201 }
1202
1203 double hh = f / (h + h);
1204
1205 for (int j = 0; j < i; j++)
1206 {
1207 e[j] -= (hh * d[j]);
1208 }
1209
1210 for (int j = 0; j < i; j++)
1211 {
1212 f = d[j];
1213 g = e[j];
1214
1215 for (int k = j; k <= (i - 1); k++)
1216 {
1217 V[k][j] -= ((f * e[k]) + (g * d[k]));
1218 }
1219
1220 d[j] = V[i - 1][j];
1221 V[i][j] = 0.0;
1222 }
1223 }
1224
1225 d[i] = h;
1226 }
1227
1228 // Accumulate transformations.
1229 for (int i = 0; i < (n - 1); i++)
1230 {
1231 V[n - 1][i] = V[i][i];
1232 V[i][i] = 1.0;
1233
1234 double h = d[i + 1];
1235
1236 if (h != 0.0)
1237 {
1238 for (int k = 0; k <= i; k++)
1239 {
1240 d[k] = V[k][i + 1] / h;
1241 }
1242
1243 for (int j = 0; j <= i; j++)
1244 {
1245 double g = 0.0;
1246
1247 for (int k = 0; k <= i; k++)
1248 {
1249 g += (V[k][i + 1] * V[k][j]);
1250 }
1251
1252 for (int k = 0; k <= i; k++)
1253 {
1254 V[k][j] -= (g * d[k]);
1255 }
1256 }
1257 }
1258
1259 for (int k = 0; k <= i; k++)
1260 {
1261 V[k][i + 1] = 0.0;
1262 }
1263 }
1264
1265 for (int j = 0; j < n; j++)
1266 {
1267 d[j] = V[n - 1][j];
1268 V[n - 1][j] = 0.0;
1269 }
1270
1271 V[n - 1][n - 1] = 1.0;
1272 e[0] = 0.0;
1273 }
1274}
1275///////////////////////////////////////////////////////////////////////////////
1276// END OF FILE.
1277///////////////////////////////////////////////////////////////////////////////