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