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