3 * $Date: 2005-11-10 09:52:44f -0600 (Thu, 10 Nov 2005) $
\r
6 * Copyright (C) 2003-2005 Miguel, Jmol Development, www.jmol.org
\r
8 * Contact: jmol-developers@lists.sf.net
\r
10 * This library is free software; you can redistribute it and/or
\r
11 * modify it under the terms of the GNU Lesser General Public
\r
12 * License as published by the Free Software Foundation; either
\r
13 * version 2.1 of the License, or (at your option) any later version.
\r
15 * This library is distributed in the hope that it will be useful,
\r
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
\r
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
\r
18 * Lesser General Public License for more details.
\r
20 * You should have received a copy of the GNU Lesser General Public
\r
21 * License along with this library; if not, write to the Free Software
\r
22 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
\r
25 package javajs.util;
\r
27 import javajs.api.EigenInterface;
\r
31 * Eigenvalues and eigenvectors of a real matrix.
\r
32 * See javajs.api.EigenInterface() as well.
\r
34 * adapted by Bob Hanson from http://math.nist.gov/javanumerics/jama/ (public
\r
35 * domain); adding quaternion superimposition capability; removing
\r
36 * nonsymmetric reduction to Hessenberg form, which we do not need in Jmol.
\r
38 * Output is as a set of double[n] columns, but for the EigenInterface
\r
39 * we return them as V3[3] and float[3] (or double[3]) values.
\r
41 * Eigenvalues and eigenvectors are sorted from smallest to largest eigenvalue.
\r
44 * If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is diagonal
\r
45 * and the eigenvector matrix V is orthogonal. I.e. A =
\r
46 * V.times(D.times(V.transpose())) and V.times(V.transpose()) equals the
\r
49 * If A is not symmetric, then the eigenvalue matrix D is block diagonal with
\r
50 * the real eigenvalues in 1-by-1 blocks and any complex eigenvalues, lambda +
\r
51 * i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda]. The columns of V represent
\r
52 * the eigenvectors in the sense that A*V = V*D, i.e. A.times(V) equals
\r
53 * V.times(D). The matrix V may be badly conditioned, or even singular, so the
\r
54 * validity of the equation A = V*D*inverse(V) depends upon V.cond().
\r
57 public class Eigen implements EigenInterface {
\r
59 /* ------------------------
\r
61 * ------------------------ */
\r
65 public Eigen set(int n) {
\r
67 V = new double[n][n];
\r
74 public Eigen setM(double[][] m) {
\r
81 * return values sorted from smallest to largest value.
\r
84 public double[] getEigenvalues() {
\r
89 * Specifically for 3x3 systems, returns eigenVectors as V3[3]
\r
90 * and values as float[3]; sorted from smallest to largest value.
\r
92 * @param eigenVectors returned vectors
\r
93 * @param eigenValues returned values
\r
97 public void fillFloatArrays(V3[] eigenVectors, float[] eigenValues) {
\r
98 for (int i = 0; i < 3; i++) {
\r
99 if (eigenVectors != null) {
\r
100 if (eigenVectors[i] == null)
\r
101 eigenVectors[i] = new V3();
\r
102 eigenVectors[i].set((float) V[0][i], (float) V[1][i], (float) V[2][i]);
\r
104 if (eigenValues != null)
\r
105 eigenValues[i] = (float) d[i];
\r
110 * Transpose V and turn into floats; sorted from smallest to largest value.
\r
112 * @return ROWS of eigenvectors f[0], f[1], f[2], etc.
\r
115 public float[][] getEigenvectorsFloatTransposed() {
\r
116 float[][] f = new float[n][n];
\r
117 for (int i = n; --i >= 0;)
\r
118 for (int j = n; --j >= 0;)
\r
119 f[j][i] = (float) V[i][j];
\r
125 * Check for symmetry, then construct the eigenvalue decomposition
\r
131 public void calc(double[][] A) {
\r
133 /* Jmol only has need of symmetric solutions
\r
135 issymmetric = true;
\r
137 for (int j = 0; (j < n) & issymmetric; j++) {
\r
138 for (int i = 0; (i < n) & issymmetric; i++) {
\r
139 issymmetric = (A[i][j] == A[j][i]);
\r
145 for (int i = 0; i < n; i++) {
\r
146 for (int j = 0; j < n; j++) {
\r
158 H = new double[n][n];
\r
159 ort = new double[n];
\r
161 for (int j = 0; j < n; j++) {
\r
162 for (int i = 0; i < n; i++) {
\r
167 // Reduce to Hessenberg form.
\r
170 // Reduce Hessenberg to real Schur form.
\r
178 * Return the real parts of the eigenvalues
\r
180 * @return real(diag(D))
\r
183 public double[] getRealEigenvalues() {
\r
188 * Return the imaginary parts of the eigenvalues
\r
190 * @return imag(diag(D))
\r
193 public double[] getImagEigenvalues() {
\r
197 /* ------------------------
\r
199 * ------------------------ */
\r
202 * Row and column dimension (square matrix).
\r
204 * @serial matrix dimension.
\r
211 * @serial internal symmetry flag.
\r
213 //private boolean issymmetric = true;
\r
216 * Arrays for internal storage of eigenvalues.
\r
218 * @serial internal storage of eigenvalues.
\r
220 private double[] d, e;
\r
223 * Array for internal storage of eigenvectors.
\r
225 * @serial internal storage of eigenvectors.
\r
227 private double[][] V;
\r
230 * Array for internal storage of nonsymmetric Hessenberg form.
\r
232 * @serial internal storage of nonsymmetric Hessenberg form.
\r
234 //private double[][] H;
\r
237 * Working storage for nonsymmetric algorithm.
\r
239 * @serial working storage for nonsymmetric algorithm.
\r
241 //private double[] ort;
\r
243 /* ------------------------
\r
245 * ------------------------ */
\r
247 // Symmetric Householder reduction to tridiagonal form.
\r
249 private void tred2() {
\r
251 // This is derived from the Algol procedures tred2 by
\r
252 // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
\r
253 // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
\r
254 // Fortran subroutine in EISPACK.
\r
256 for (int j = 0; j < n; j++) {
\r
257 d[j] = V[n - 1][j];
\r
260 // Householder reduction to tridiagonal form.
\r
262 for (int i = n - 1; i > 0; i--) {
\r
264 // Scale to avoid under/overflow.
\r
266 double scale = 0.0;
\r
268 for (int k = 0; k < i; k++) {
\r
269 scale = scale + Math.abs(d[k]);
\r
271 if (scale == 0.0) {
\r
273 for (int j = 0; j < i; j++) {
\r
274 d[j] = V[i - 1][j];
\r
280 // Generate Householder vector.
\r
282 for (int k = 0; k < i; k++) {
\r
286 double f = d[i - 1];
\r
287 double g = Math.sqrt(h);
\r
294 for (int j = 0; j < i; j++) {
\r
298 // Apply similarity transformation to remaining columns.
\r
300 for (int j = 0; j < i; j++) {
\r
303 g = e[j] + V[j][j] * f;
\r
304 for (int k = j + 1; k <= i - 1; k++) {
\r
305 g += V[k][j] * d[k];
\r
306 e[k] += V[k][j] * f;
\r
311 for (int j = 0; j < i; j++) {
\r
315 double hh = f / (h + h);
\r
316 for (int j = 0; j < i; j++) {
\r
319 for (int j = 0; j < i; j++) {
\r
322 for (int k = j; k <= i - 1; k++) {
\r
323 V[k][j] -= (f * e[k] + g * d[k]);
\r
325 d[j] = V[i - 1][j];
\r
332 // Accumulate transformations.
\r
334 for (int i = 0; i < n - 1; i++) {
\r
335 V[n - 1][i] = V[i][i];
\r
337 double h = d[i + 1];
\r
339 for (int k = 0; k <= i; k++) {
\r
340 d[k] = V[k][i + 1] / h;
\r
342 for (int j = 0; j <= i; j++) {
\r
344 for (int k = 0; k <= i; k++) {
\r
345 g += V[k][i + 1] * V[k][j];
\r
347 for (int k = 0; k <= i; k++) {
\r
348 V[k][j] -= g * d[k];
\r
352 for (int k = 0; k <= i; k++) {
\r
356 for (int j = 0; j < n; j++) {
\r
357 d[j] = V[n - 1][j];
\r
360 V[n - 1][n - 1] = 1.0;
\r
364 // Symmetric tridiagonal QL algorithm.
\r
366 private void tql2() {
\r
368 // This is derived from the Algol procedures tql2, by
\r
369 // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
\r
370 // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
\r
371 // Fortran subroutine in EISPACK.
\r
373 for (int i = 1; i < n; i++) {
\r
380 double eps = Math.pow(2.0, -52.0);
\r
381 for (int l = 0; l < n; l++) {
\r
383 // Find small subdiagonal element
\r
385 tst1 = Math.max(tst1, Math.abs(d[l]) + Math.abs(e[l]));
\r
388 if (Math.abs(e[m]) <= eps * tst1) {
\r
394 // If m == l, d[l] is an eigenvalue,
\r
395 // otherwise, iterate.
\r
400 iter = iter + 1; // (Could check iteration count here.)
\r
402 // Compute implicit shift
\r
405 double p = (d[l + 1] - g) / (2.0 * e[l]);
\r
406 double r = hypot(p, 1.0);
\r
410 d[l] = e[l] / (p + r);
\r
411 d[l + 1] = e[l] * (p + r);
\r
412 double dl1 = d[l + 1];
\r
413 double h = g - d[l];
\r
414 for (int i = l + 2; i < n; i++) {
\r
419 // Implicit QL transformation.
\r
425 double el1 = e[l + 1];
\r
428 for (int i = m - 1; i >= l; i--) {
\r
434 r = hypot(p, e[i]);
\r
438 p = c * d[i] - s * g;
\r
439 d[i + 1] = h + s * (c * g + s * d[i]);
\r
441 // Accumulate transformation.
\r
443 for (int k = 0; k < n; k++) {
\r
445 V[k][i + 1] = s * V[k][i] + c * h;
\r
446 V[k][i] = c * V[k][i] - s * h;
\r
449 p = -s * s2 * c3 * el1 * e[l] / dl1;
\r
453 // Check for convergence.
\r
455 } while (Math.abs(e[l]) > eps * tst1);
\r
461 // Sort eigenvalues and corresponding vectors.
\r
463 for (int i = 0; i < n - 1; i++) {
\r
466 for (int j = i + 1; j < n; j++) {
\r
475 for (int j = 0; j < n; j++) {
\r
484 private static double hypot(double a, double b) {
\r
486 // sqrt(a^2 + b^2) without under/overflow.
\r
489 if (Math.abs(a) > Math.abs(b)) {
\r
491 r = Math.abs(a) * Math.sqrt(1 + r * r);
\r
492 } else if (b != 0) {
\r
494 r = Math.abs(b) * Math.sqrt(1 + r * r);
\r
501 // Nonsymmetric reduction to Hessenberg form.
\r
504 private void orthes() {
\r
506 // This is derived from the Algol procedures orthes and ortran,
\r
507 // by Martin and Wilkinson, Handbook for Auto. Comp.,
\r
508 // Vol.ii-Linear Algebra, and the corresponding
\r
509 // Fortran subroutines in EISPACK.
\r
514 for (int m = low + 1; m <= high - 1; m++) {
\r
518 double scale = 0.0;
\r
519 for (int i = m; i <= high; i++) {
\r
520 scale = scale + Math.abs(H[i][m - 1]);
\r
522 if (scale != 0.0) {
\r
524 // Compute Householder transformation.
\r
527 for (int i = high; i >= m; i--) {
\r
528 ort[i] = H[i][m - 1] / scale;
\r
529 h += ort[i] * ort[i];
\r
531 double g = Math.sqrt(h);
\r
535 h = h - ort[m] * g;
\r
536 ort[m] = ort[m] - g;
\r
538 // Apply Householder similarity transformation
\r
539 // H = (I-u*u'/h)*H*(I-u*u')/h)
\r
541 for (int j = m; j < n; j++) {
\r
543 for (int i = high; i >= m; i--) {
\r
544 f += ort[i] * H[i][j];
\r
547 for (int i = m; i <= high; i++) {
\r
548 H[i][j] -= f * ort[i];
\r
552 for (int i = 0; i <= high; i++) {
\r
554 for (int j = high; j >= m; j--) {
\r
555 f += ort[j] * H[i][j];
\r
558 for (int j = m; j <= high; j++) {
\r
559 H[i][j] -= f * ort[j];
\r
562 ort[m] = scale * ort[m];
\r
563 H[m][m - 1] = scale * g;
\r
567 // Accumulate transformations (Algol's ortran).
\r
569 for (int i = 0; i < n; i++) {
\r
570 for (int j = 0; j < n; j++) {
\r
571 V[i][j] = (i == j ? 1.0 : 0.0);
\r
575 for (int m = high - 1; m >= low + 1; m--) {
\r
576 if (H[m][m - 1] != 0.0) {
\r
577 for (int i = m + 1; i <= high; i++) {
\r
578 ort[i] = H[i][m - 1];
\r
580 for (int j = m; j <= high; j++) {
\r
582 for (int i = m; i <= high; i++) {
\r
583 g += ort[i] * V[i][j];
\r
585 // Double division avoids possible underflow
\r
586 g = (g / ort[m]) / H[m][m - 1];
\r
587 for (int i = m; i <= high; i++) {
\r
588 V[i][j] += g * ort[i];
\r
595 // Complex scalar division.
\r
597 private transient double cdivr, cdivi;
\r
599 private void cdiv(double xr, double xi, double yr, double yi) {
\r
601 if (Math.abs(yr) > Math.abs(yi)) {
\r
604 cdivr = (xr + r * xi) / d;
\r
605 cdivi = (xi - r * xr) / d;
\r
609 cdivr = (r * xr + xi) / d;
\r
610 cdivi = (r * xi - xr) / d;
\r
614 // Nonsymmetric reduction from Hessenberg to real Schur form.
\r
616 private void hqr2() {
\r
618 // This is derived from the Algol procedure hqr2,
\r
619 // by Martin and Wilkinson, Handbook for Auto. Comp.,
\r
620 // Vol.ii-Linear Algebra, and the corresponding
\r
621 // Fortran subroutine in EISPACK.
\r
629 double eps = Math.pow(2.0, -52.0);
\r
630 double exshift = 0.0;
\r
631 double p = 0, q = 0, r = 0, s = 0, z = 0, t, w, x, y;
\r
633 // Store roots isolated by balanc and compute matrix norm
\r
636 for (int i = 0; i < nn; i++) {
\r
637 if (i < low || i > high) {
\r
641 for (int j = Math.max(i - 1, 0); j < nn; j++) {
\r
642 norm = norm + Math.abs(H[i][j]);
\r
646 // Outer loop over eigenvalue index
\r
651 // Look for single small sub-diagonal element
\r
655 s = Math.abs(H[l - 1][l - 1]) + Math.abs(H[l][l]);
\r
659 if (Math.abs(H[l][l - 1]) < eps * s) {
\r
665 // Check for convergence
\r
669 H[n][n] = H[n][n] + exshift;
\r
677 } else if (l == n - 1) {
\r
678 w = H[n][n - 1] * H[n - 1][n];
\r
679 p = (H[n - 1][n - 1] - H[n][n]) / 2.0;
\r
681 z = Math.sqrt(Math.abs(q));
\r
682 H[n][n] = H[n][n] + exshift;
\r
683 H[n - 1][n - 1] = H[n - 1][n - 1] + exshift;
\r
702 s = Math.abs(x) + Math.abs(z);
\r
705 r = Math.sqrt(p * p + q * q);
\r
709 // Row modification
\r
711 for (int j = n - 1; j < nn; j++) {
\r
713 H[n - 1][j] = q * z + p * H[n][j];
\r
714 H[n][j] = q * H[n][j] - p * z;
\r
717 // Column modification
\r
719 for (int i = 0; i <= n; i++) {
\r
721 H[i][n - 1] = q * z + p * H[i][n];
\r
722 H[i][n] = q * H[i][n] - p * z;
\r
725 // Accumulate transformations
\r
727 for (int i = low; i <= high; i++) {
\r
729 V[i][n - 1] = q * z + p * V[i][n];
\r
730 V[i][n] = q * V[i][n] - p * z;
\r
744 // No convergence yet
\r
754 y = H[n - 1][n - 1];
\r
755 w = H[n][n - 1] * H[n - 1][n];
\r
758 // Wilkinson's original ad hoc shift
\r
762 for (int i = low; i <= n; i++) {
\r
765 s = Math.abs(H[n][n - 1]) + Math.abs(H[n - 1][n - 2]);
\r
767 w = -0.4375 * s * s;
\r
770 // MATLAB's new ad hoc shift
\r
780 s = x - w / ((y - x) / 2.0 + s);
\r
781 for (int i = low; i <= n; i++) {
\r
789 iter = iter + 1; // (Could check iteration count here.)
\r
791 // Look for two consecutive small sub-diagonal elements
\r
798 p = (r * s - w) / H[m + 1][m] + H[m][m + 1];
\r
799 q = H[m + 1][m + 1] - z - r - s;
\r
800 r = H[m + 2][m + 1];
\r
801 s = Math.abs(p) + Math.abs(q) + Math.abs(r);
\r
808 if (Math.abs(H[m][m - 1]) * (Math.abs(q) + Math.abs(r)) < eps
\r
809 * (Math.abs(p) * (Math.abs(H[m - 1][m - 1]) + Math.abs(z) + Math
\r
810 .abs(H[m + 1][m + 1])))) {
\r
816 for (int i = m + 2; i <= n; i++) {
\r
823 // Double QR step involving rows l:n and columns m:n
\r
825 for (int k = m; k <= n - 1; k++) {
\r
826 boolean notlast = (k != n - 1);
\r
829 q = H[k + 1][k - 1];
\r
830 r = (notlast ? H[k + 2][k - 1] : 0.0);
\r
831 x = Math.abs(p) + Math.abs(q) + Math.abs(r);
\r
841 s = Math.sqrt(p * p + q * q + r * r);
\r
847 H[k][k - 1] = -s * x;
\r
848 } else if (l != m) {
\r
849 H[k][k - 1] = -H[k][k - 1];
\r
858 // Row modification
\r
860 for (int j = k; j < nn; j++) {
\r
861 p = H[k][j] + q * H[k + 1][j];
\r
863 p = p + r * H[k + 2][j];
\r
864 H[k + 2][j] = H[k + 2][j] - p * z;
\r
866 H[k][j] = H[k][j] - p * x;
\r
867 H[k + 1][j] = H[k + 1][j] - p * y;
\r
870 // Column modification
\r
872 for (int i = 0; i <= Math.min(n, k + 3); i++) {
\r
873 p = x * H[i][k] + y * H[i][k + 1];
\r
875 p = p + z * H[i][k + 2];
\r
876 H[i][k + 2] = H[i][k + 2] - p * r;
\r
878 H[i][k] = H[i][k] - p;
\r
879 H[i][k + 1] = H[i][k + 1] - p * q;
\r
882 // Accumulate transformations
\r
884 for (int i = low; i <= high; i++) {
\r
885 p = x * V[i][k] + y * V[i][k + 1];
\r
887 p = p + z * V[i][k + 2];
\r
888 V[i][k + 2] = V[i][k + 2] - p * r;
\r
890 V[i][k] = V[i][k] - p;
\r
891 V[i][k + 1] = V[i][k + 1] - p * q;
\r
895 } // check convergence
\r
896 } // while (n >= low)
\r
898 // Backsubstitute to find vectors of upper triangular form
\r
904 for (n = nn - 1; n >= 0; n--) {
\r
913 for (int i = n - 1; i >= 0; i--) {
\r
916 for (int j = l; j <= n; j++) {
\r
917 r = r + H[i][j] * H[j][n];
\r
928 H[i][n] = -r / (eps * norm);
\r
931 // Solve real equations
\r
936 q = (d[i] - p) * (d[i] - p) + e[i] * e[i];
\r
937 t = (x * s - z * r) / q;
\r
939 if (Math.abs(x) > Math.abs(z)) {
\r
940 H[i + 1][n] = (-r - w * t) / x;
\r
942 H[i + 1][n] = (-s - y * t) / z;
\r
946 // Overflow control
\r
948 t = Math.abs(H[i][n]);
\r
949 if ((eps * t) * t > 1) {
\r
950 for (int j = i; j <= n; j++) {
\r
951 H[j][n] = H[j][n] / t;
\r
959 } else if (q < 0) {
\r
962 // Last vector component imaginary so matrix is triangular
\r
964 if (Math.abs(H[n][n - 1]) > Math.abs(H[n - 1][n])) {
\r
965 H[n - 1][n - 1] = q / H[n][n - 1];
\r
966 H[n - 1][n] = -(H[n][n] - p) / H[n][n - 1];
\r
968 cdiv(0.0, -H[n - 1][n], H[n - 1][n - 1] - p, q);
\r
969 H[n - 1][n - 1] = cdivr;
\r
970 H[n - 1][n] = cdivi;
\r
974 for (int i = n - 2; i >= 0; i--) {
\r
975 double ra, sa, vr, vi;
\r
978 for (int j = l; j <= n; j++) {
\r
979 ra = ra + H[i][j] * H[j][n - 1];
\r
980 sa = sa + H[i][j] * H[j][n];
\r
991 cdiv(-ra, -sa, w, q);
\r
992 H[i][n - 1] = cdivr;
\r
996 // Solve complex equations
\r
1000 vr = (d[i] - p) * (d[i] - p) + e[i] * e[i] - q * q;
\r
1001 vi = (d[i] - p) * 2.0 * q;
\r
1002 if (vr == 0.0 & vi == 0.0) {
\r
1005 * (Math.abs(w) + Math.abs(q) + Math.abs(x) + Math.abs(y) + Math
\r
1008 cdiv(x * r - z * ra + q * sa, x * s - z * sa - q * ra, vr, vi);
\r
1009 H[i][n - 1] = cdivr;
\r
1011 if (Math.abs(x) > (Math.abs(z) + Math.abs(q))) {
\r
1012 H[i + 1][n - 1] = (-ra - w * H[i][n - 1] + q * H[i][n]) / x;
\r
1013 H[i + 1][n] = (-sa - w * H[i][n] - q * H[i][n - 1]) / x;
\r
1015 cdiv(-r - y * H[i][n - 1], -s - y * H[i][n], z, q);
\r
1016 H[i + 1][n - 1] = cdivr;
\r
1017 H[i + 1][n] = cdivi;
\r
1021 // Overflow control
\r
1023 t = Math.max(Math.abs(H[i][n - 1]), Math.abs(H[i][n]));
\r
1024 if ((eps * t) * t > 1) {
\r
1025 for (int j = i; j <= n; j++) {
\r
1026 H[j][n - 1] = H[j][n - 1] / t;
\r
1027 H[j][n] = H[j][n] / t;
\r
1035 // Vectors of isolated roots
\r
1037 for (int i = 0; i < nn; i++) {
\r
1038 if (i < low || i > high) {
\r
1039 for (int j = i; j < nn; j++) {
\r
1040 V[i][j] = H[i][j];
\r
1045 // Back transformation to get eigenvectors of original matrix
\r
1047 for (int j = nn - 1; j >= low; j--) {
\r
1048 for (int i = low; i <= high; i++) {
\r
1050 for (int k = low; k <= Math.min(j, high); k++) {
\r
1051 z = z + V[i][k] * H[k][j];
\r