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