001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    
018    package org.apache.commons.math3.linear;
019    
020    import org.apache.commons.math3.complex.Complex;
021    import org.apache.commons.math3.exception.MathArithmeticException;
022    import org.apache.commons.math3.exception.MathUnsupportedOperationException;
023    import org.apache.commons.math3.exception.MaxCountExceededException;
024    import org.apache.commons.math3.exception.DimensionMismatchException;
025    import org.apache.commons.math3.exception.util.LocalizedFormats;
026    import org.apache.commons.math3.util.Precision;
027    import org.apache.commons.math3.util.FastMath;
028    
029    /**
030     * Calculates the eigen decomposition of a real matrix.
031     * <p>The eigen decomposition of matrix A is a set of two matrices:
032     * V and D such that A = V &times; D &times; V<sup>T</sup>.
033     * A, V and D are all m &times; m matrices.</p>
034     * <p>This class is similar in spirit to the <code>EigenvalueDecomposition</code>
035     * class from the <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a>
036     * library, with the following changes:</p>
037     * <ul>
038     *   <li>a {@link #getVT() getVt} method has been added,</li>
039     *   <li>two {@link #getRealEigenvalue(int) getRealEigenvalue} and {@link #getImagEigenvalue(int)
040     *   getImagEigenvalue} methods to pick up a single eigenvalue have been added,</li>
041     *   <li>a {@link #getEigenvector(int) getEigenvector} method to pick up a single
042     *   eigenvector has been added,</li>
043     *   <li>a {@link #getDeterminant() getDeterminant} method has been added.</li>
044     *   <li>a {@link #getSolver() getSolver} method has been added.</li>
045     * </ul>
046     * <p>
047     * As of 3.1, this class supports general real matrices (both symmetric and non-symmetric):
048     * </p>
049     * <p>
050     * If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is diagonal and the eigenvector
051     * matrix V is orthogonal, i.e. A = V.multiply(D.multiply(V.transpose())) and
052     * V.multiply(V.transpose()) equals the identity matrix.
053     * </p>
054     * <p>
055     * If A is not symmetric, then the eigenvalue matrix D is block diagonal with the real eigenvalues
056     * in 1-by-1 blocks and any complex eigenvalues, lambda + i*mu, in 2-by-2 blocks:
057     * <pre>
058     *    [lambda, mu    ]
059     *    [   -mu, lambda]
060     * </pre>
061     * The columns of V represent the eigenvectors in the sense that A*V = V*D,
062     * i.e. A.multiply(V) equals V.multiply(D).
063     * The matrix V may be badly conditioned, or even singular, so the validity of the equation
064     * A = V*D*inverse(V) depends upon the condition of V.
065     * </p>
066     * <p>
067     * This implementation is based on the paper by A. Drubrulle, R.S. Martin and
068     * J.H. Wilkinson "The Implicit QL Algorithm" in Wilksinson and Reinsch (1971)
069     * Handbook for automatic computation, vol. 2, Linear algebra, Springer-Verlag,
070     * New-York
071     * </p>
072     * @see <a href="http://mathworld.wolfram.com/EigenDecomposition.html">MathWorld</a>
073     * @see <a href="http://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix">Wikipedia</a>
074     * @version $Id: EigenDecomposition.java 1422195 2012-12-15 06:45:18Z psteitz $
075     * @since 2.0 (changed to concrete class in 3.0)
076     */
077    public class EigenDecomposition {
078        /** Internally used epsilon criteria. */
079        private static final double EPSILON = 1e-12;
080        /** Maximum number of iterations accepted in the implicit QL transformation */
081        private byte maxIter = 30;
082        /** Main diagonal of the tridiagonal matrix. */
083        private double[] main;
084        /** Secondary diagonal of the tridiagonal matrix. */
085        private double[] secondary;
086        /**
087         * Transformer to tridiagonal (may be null if matrix is already
088         * tridiagonal).
089         */
090        private TriDiagonalTransformer transformer;
091        /** Real part of the realEigenvalues. */
092        private double[] realEigenvalues;
093        /** Imaginary part of the realEigenvalues. */
094        private double[] imagEigenvalues;
095        /** Eigenvectors. */
096        private ArrayRealVector[] eigenvectors;
097        /** Cached value of V. */
098        private RealMatrix cachedV;
099        /** Cached value of D. */
100        private RealMatrix cachedD;
101        /** Cached value of Vt. */
102        private RealMatrix cachedVt;
103        /** Whether the matrix is symmetric. */
104        private final boolean isSymmetric;
105    
106        /**
107         * Calculates the eigen decomposition of the given real matrix.
108         * <p>
109         * Supports decomposition of a general matrix since 3.1.
110         *
111         * @param matrix Matrix to decompose.
112         * @throws MaxCountExceededException if the algorithm fails to converge.
113         * @throws MathArithmeticException if the decomposition of a general matrix
114         * results in a matrix with zero norm
115         * @since 3.1
116         */
117        public EigenDecomposition(final RealMatrix matrix)
118            throws MathArithmeticException {
119            final double symTol = 10 * matrix.getRowDimension() * matrix.getColumnDimension() * Precision.EPSILON;
120            isSymmetric = MatrixUtils.isSymmetric(matrix, symTol);
121            if (isSymmetric) {
122                transformToTridiagonal(matrix);
123                findEigenVectors(transformer.getQ().getData());
124            } else {
125                final SchurTransformer t = transformToSchur(matrix);
126                findEigenVectorsFromSchur(t);
127            }
128        }
129    
130        /**
131         * Calculates the eigen decomposition of the given real matrix.
132         *
133         * @param matrix Matrix to decompose.
134         * @param splitTolerance Dummy parameter (present for backward
135         * compatibility only).
136         * @throws MathArithmeticException  if the decomposition of a general matrix
137         * results in a matrix with zero norm
138         * @throws MaxCountExceededException if the algorithm fails to converge.
139         * @deprecated in 3.1 (to be removed in 4.0) due to unused parameter
140         */
141        @Deprecated
142        public EigenDecomposition(final RealMatrix matrix,
143                                  final double splitTolerance)
144            throws MathArithmeticException {
145            this(matrix);
146        }
147    
148        /**
149         * Calculates the eigen decomposition of the symmetric tridiagonal
150         * matrix.  The Householder matrix is assumed to be the identity matrix.
151         *
152         * @param main Main diagonal of the symmetric tridiagonal form.
153         * @param secondary Secondary of the tridiagonal form.
154         * @throws MaxCountExceededException if the algorithm fails to converge.
155         * @since 3.1
156         */
157        public EigenDecomposition(final double[] main, final double[] secondary) {
158            isSymmetric = true;
159            this.main      = main.clone();
160            this.secondary = secondary.clone();
161            transformer    = null;
162            final int size = main.length;
163            final double[][] z = new double[size][size];
164            for (int i = 0; i < size; i++) {
165                z[i][i] = 1.0;
166            }
167            findEigenVectors(z);
168        }
169    
170        /**
171         * Calculates the eigen decomposition of the symmetric tridiagonal
172         * matrix.  The Householder matrix is assumed to be the identity matrix.
173         *
174         * @param main Main diagonal of the symmetric tridiagonal form.
175         * @param secondary Secondary of the tridiagonal form.
176         * @param splitTolerance Dummy parameter (present for backward
177         * compatibility only).
178         * @throws MaxCountExceededException if the algorithm fails to converge.
179         * @deprecated in 3.1 (to be removed in 4.0) due to unused parameter
180         */
181        @Deprecated
182        public EigenDecomposition(final double[] main, final double[] secondary,
183                                  final double splitTolerance) {
184            this(main, secondary);
185        }
186    
187        /**
188         * Gets the matrix V of the decomposition.
189         * V is an orthogonal matrix, i.e. its transpose is also its inverse.
190         * The columns of V are the eigenvectors of the original matrix.
191         * No assumption is made about the orientation of the system axes formed
192         * by the columns of V (e.g. in a 3-dimension space, V can form a left-
193         * or right-handed system).
194         *
195         * @return the V matrix.
196         */
197        public RealMatrix getV() {
198    
199            if (cachedV == null) {
200                final int m = eigenvectors.length;
201                cachedV = MatrixUtils.createRealMatrix(m, m);
202                for (int k = 0; k < m; ++k) {
203                    cachedV.setColumnVector(k, eigenvectors[k]);
204                }
205            }
206            // return the cached matrix
207            return cachedV;
208        }
209    
210        /**
211         * Gets the block diagonal matrix D of the decomposition.
212         * D is a block diagonal matrix.
213         * Real eigenvalues are on the diagonal while complex values are on
214         * 2x2 blocks { {real +imaginary}, {-imaginary, real} }.
215         *
216         * @return the D matrix.
217         *
218         * @see #getRealEigenvalues()
219         * @see #getImagEigenvalues()
220         */
221        public RealMatrix getD() {
222    
223            if (cachedD == null) {
224                // cache the matrix for subsequent calls
225                cachedD = MatrixUtils.createRealDiagonalMatrix(realEigenvalues);
226    
227                for (int i = 0; i < imagEigenvalues.length; i++) {
228                    if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) > 0) {
229                        cachedD.setEntry(i, i+1, imagEigenvalues[i]);
230                    } else if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) < 0) {
231                        cachedD.setEntry(i, i-1, imagEigenvalues[i]);
232                    }
233                }
234            }
235            return cachedD;
236        }
237    
238        /**
239         * Gets the transpose of the matrix V of the decomposition.
240         * V is an orthogonal matrix, i.e. its transpose is also its inverse.
241         * The columns of V are the eigenvectors of the original matrix.
242         * No assumption is made about the orientation of the system axes formed
243         * by the columns of V (e.g. in a 3-dimension space, V can form a left-
244         * or right-handed system).
245         *
246         * @return the transpose of the V matrix.
247         */
248        public RealMatrix getVT() {
249    
250            if (cachedVt == null) {
251                final int m = eigenvectors.length;
252                cachedVt = MatrixUtils.createRealMatrix(m, m);
253                for (int k = 0; k < m; ++k) {
254                    cachedVt.setRowVector(k, eigenvectors[k]);
255                }
256            }
257    
258            // return the cached matrix
259            return cachedVt;
260        }
261    
262        /**
263         * Returns whether the calculated eigen values are complex or real.
264         * <p>The method performs a zero check for each element of the
265         * {@link #getImagEigenvalues()} array and returns {@code true} if any
266         * element is not equal to zero.
267         *
268         * @return {@code true} if the eigen values are complex, {@code false} otherwise
269         * @since 3.1
270         */
271        public boolean hasComplexEigenvalues() {
272            for (int i = 0; i < imagEigenvalues.length; i++) {
273                if (!Precision.equals(imagEigenvalues[i], 0.0, EPSILON)) {
274                    return true;
275                }
276            }
277            return false;
278        }
279    
280        /**
281         * Gets a copy of the real parts of the eigenvalues of the original matrix.
282         *
283         * @return a copy of the real parts of the eigenvalues of the original matrix.
284         *
285         * @see #getD()
286         * @see #getRealEigenvalue(int)
287         * @see #getImagEigenvalues()
288         */
289        public double[] getRealEigenvalues() {
290            return realEigenvalues.clone();
291        }
292    
293        /**
294         * Returns the real part of the i<sup>th</sup> eigenvalue of the original
295         * matrix.
296         *
297         * @param i index of the eigenvalue (counting from 0)
298         * @return real part of the i<sup>th</sup> eigenvalue of the original
299         * matrix.
300         *
301         * @see #getD()
302         * @see #getRealEigenvalues()
303         * @see #getImagEigenvalue(int)
304         */
305        public double getRealEigenvalue(final int i) {
306            return realEigenvalues[i];
307        }
308    
309        /**
310         * Gets a copy of the imaginary parts of the eigenvalues of the original
311         * matrix.
312         *
313         * @return a copy of the imaginary parts of the eigenvalues of the original
314         * matrix.
315         *
316         * @see #getD()
317         * @see #getImagEigenvalue(int)
318         * @see #getRealEigenvalues()
319         */
320        public double[] getImagEigenvalues() {
321            return imagEigenvalues.clone();
322        }
323    
324        /**
325         * Gets the imaginary part of the i<sup>th</sup> eigenvalue of the original
326         * matrix.
327         *
328         * @param i Index of the eigenvalue (counting from 0).
329         * @return the imaginary part of the i<sup>th</sup> eigenvalue of the original
330         * matrix.
331         *
332         * @see #getD()
333         * @see #getImagEigenvalues()
334         * @see #getRealEigenvalue(int)
335         */
336        public double getImagEigenvalue(final int i) {
337            return imagEigenvalues[i];
338        }
339    
340        /**
341         * Gets a copy of the i<sup>th</sup> eigenvector of the original matrix.
342         *
343         * @param i Index of the eigenvector (counting from 0).
344         * @return a copy of the i<sup>th</sup> eigenvector of the original matrix.
345         * @see #getD()
346         */
347        public RealVector getEigenvector(final int i) {
348            return eigenvectors[i].copy();
349        }
350    
351        /**
352         * Computes the determinant of the matrix.
353         *
354         * @return the determinant of the matrix.
355         */
356        public double getDeterminant() {
357            double determinant = 1;
358            for (double lambda : realEigenvalues) {
359                determinant *= lambda;
360            }
361            return determinant;
362        }
363    
364        /**
365         * Computes the square-root of the matrix.
366         * This implementation assumes that the matrix is symmetric and postive
367         * definite.
368         *
369         * @return the square-root of the matrix.
370         * @throws MathUnsupportedOperationException if the matrix is not
371         * symmetric or not positive definite.
372         * @since 3.1
373         */
374        public RealMatrix getSquareRoot() {
375            if (!isSymmetric) {
376                throw new MathUnsupportedOperationException();
377            }
378    
379            final double[] sqrtEigenValues = new double[realEigenvalues.length];
380            for (int i = 0; i < realEigenvalues.length; i++) {
381                final double eigen = realEigenvalues[i];
382                if (eigen <= 0) {
383                    throw new MathUnsupportedOperationException();
384                }
385                sqrtEigenValues[i] = FastMath.sqrt(eigen);
386            }
387            final RealMatrix sqrtEigen = MatrixUtils.createRealDiagonalMatrix(sqrtEigenValues);
388            final RealMatrix v = getV();
389            final RealMatrix vT = getVT();
390    
391            return v.multiply(sqrtEigen).multiply(vT);
392        }
393    
394        /**
395         * Gets a solver for finding the A &times; X = B solution in exact
396         * linear sense.
397         * <p>
398         * Since 3.1, eigen decomposition of a general matrix is supported,
399         * but the {@link DecompositionSolver} only supports real eigenvalues.
400         *
401         * @return a solver
402         * @throws MathUnsupportedOperationException if the decomposition resulted in
403         * complex eigenvalues
404         */
405        public DecompositionSolver getSolver() {
406            if (hasComplexEigenvalues()) {
407                throw new MathUnsupportedOperationException();
408            }
409            return new Solver(realEigenvalues, imagEigenvalues, eigenvectors);
410        }
411    
412        /** Specialized solver. */
413        private static class Solver implements DecompositionSolver {
414            /** Real part of the realEigenvalues. */
415            private double[] realEigenvalues;
416            /** Imaginary part of the realEigenvalues. */
417            private double[] imagEigenvalues;
418            /** Eigenvectors. */
419            private final ArrayRealVector[] eigenvectors;
420    
421            /**
422             * Builds a solver from decomposed matrix.
423             *
424             * @param realEigenvalues Real parts of the eigenvalues.
425             * @param imagEigenvalues Imaginary parts of the eigenvalues.
426             * @param eigenvectors Eigenvectors.
427             */
428            private Solver(final double[] realEigenvalues,
429                    final double[] imagEigenvalues,
430                    final ArrayRealVector[] eigenvectors) {
431                this.realEigenvalues = realEigenvalues;
432                this.imagEigenvalues = imagEigenvalues;
433                this.eigenvectors = eigenvectors;
434            }
435    
436            /**
437             * Solves the linear equation A &times; X = B for symmetric matrices A.
438             * <p>
439             * This method only finds exact linear solutions, i.e. solutions for
440             * which ||A &times; X - B|| is exactly 0.
441             * </p>
442             *
443             * @param b Right-hand side of the equation A &times; X = B.
444             * @return a Vector X that minimizes the two norm of A &times; X - B.
445             *
446             * @throws DimensionMismatchException if the matrices dimensions do not match.
447             * @throws SingularMatrixException if the decomposed matrix is singular.
448             */
449            public RealVector solve(final RealVector b) {
450                if (!isNonSingular()) {
451                    throw new SingularMatrixException();
452                }
453    
454                final int m = realEigenvalues.length;
455                if (b.getDimension() != m) {
456                    throw new DimensionMismatchException(b.getDimension(), m);
457                }
458    
459                final double[] bp = new double[m];
460                for (int i = 0; i < m; ++i) {
461                    final ArrayRealVector v = eigenvectors[i];
462                    final double[] vData = v.getDataRef();
463                    final double s = v.dotProduct(b) / realEigenvalues[i];
464                    for (int j = 0; j < m; ++j) {
465                        bp[j] += s * vData[j];
466                    }
467                }
468    
469                return new ArrayRealVector(bp, false);
470            }
471    
472            /** {@inheritDoc} */
473            public RealMatrix solve(RealMatrix b) {
474    
475                if (!isNonSingular()) {
476                    throw new SingularMatrixException();
477                }
478    
479                final int m = realEigenvalues.length;
480                if (b.getRowDimension() != m) {
481                    throw new DimensionMismatchException(b.getRowDimension(), m);
482                }
483    
484                final int nColB = b.getColumnDimension();
485                final double[][] bp = new double[m][nColB];
486                final double[] tmpCol = new double[m];
487                for (int k = 0; k < nColB; ++k) {
488                    for (int i = 0; i < m; ++i) {
489                        tmpCol[i] = b.getEntry(i, k);
490                        bp[i][k]  = 0;
491                    }
492                    for (int i = 0; i < m; ++i) {
493                        final ArrayRealVector v = eigenvectors[i];
494                        final double[] vData = v.getDataRef();
495                        double s = 0;
496                        for (int j = 0; j < m; ++j) {
497                            s += v.getEntry(j) * tmpCol[j];
498                        }
499                        s /= realEigenvalues[i];
500                        for (int j = 0; j < m; ++j) {
501                            bp[j][k] += s * vData[j];
502                        }
503                    }
504                }
505    
506                return new Array2DRowRealMatrix(bp, false);
507    
508            }
509    
510            /**
511             * Checks whether the decomposed matrix is non-singular.
512             *
513             * @return true if the decomposed matrix is non-singular.
514             */
515            public boolean isNonSingular() {
516                for (int i = 0; i < realEigenvalues.length; ++i) {
517                    if (realEigenvalues[i] == 0 &&
518                        imagEigenvalues[i] == 0) {
519                        return false;
520                    }
521                }
522                return true;
523            }
524    
525            /**
526             * Get the inverse of the decomposed matrix.
527             *
528             * @return the inverse matrix.
529             * @throws SingularMatrixException if the decomposed matrix is singular.
530             */
531            public RealMatrix getInverse() {
532                if (!isNonSingular()) {
533                    throw new SingularMatrixException();
534                }
535    
536                final int m = realEigenvalues.length;
537                final double[][] invData = new double[m][m];
538    
539                for (int i = 0; i < m; ++i) {
540                    final double[] invI = invData[i];
541                    for (int j = 0; j < m; ++j) {
542                        double invIJ = 0;
543                        for (int k = 0; k < m; ++k) {
544                            final double[] vK = eigenvectors[k].getDataRef();
545                            invIJ += vK[i] * vK[j] / realEigenvalues[k];
546                        }
547                        invI[j] = invIJ;
548                    }
549                }
550                return MatrixUtils.createRealMatrix(invData);
551            }
552        }
553    
554        /**
555         * Transforms the matrix to tridiagonal form.
556         *
557         * @param matrix Matrix to transform.
558         */
559        private void transformToTridiagonal(final RealMatrix matrix) {
560            // transform the matrix to tridiagonal
561            transformer = new TriDiagonalTransformer(matrix);
562            main = transformer.getMainDiagonalRef();
563            secondary = transformer.getSecondaryDiagonalRef();
564        }
565    
566        /**
567         * Find eigenvalues and eigenvectors (Dubrulle et al., 1971)
568         *
569         * @param householderMatrix Householder matrix of the transformation
570         * to tridiagonal form.
571         */
572        private void findEigenVectors(final double[][] householderMatrix) {
573            final double[][]z = householderMatrix.clone();
574            final int n = main.length;
575            realEigenvalues = new double[n];
576            imagEigenvalues = new double[n];
577            final double[] e = new double[n];
578            for (int i = 0; i < n - 1; i++) {
579                realEigenvalues[i] = main[i];
580                e[i] = secondary[i];
581            }
582            realEigenvalues[n - 1] = main[n - 1];
583            e[n - 1] = 0;
584    
585            // Determine the largest main and secondary value in absolute term.
586            double maxAbsoluteValue = 0;
587            for (int i = 0; i < n; i++) {
588                if (FastMath.abs(realEigenvalues[i]) > maxAbsoluteValue) {
589                    maxAbsoluteValue = FastMath.abs(realEigenvalues[i]);
590                }
591                if (FastMath.abs(e[i]) > maxAbsoluteValue) {
592                    maxAbsoluteValue = FastMath.abs(e[i]);
593                }
594            }
595            // Make null any main and secondary value too small to be significant
596            if (maxAbsoluteValue != 0) {
597                for (int i=0; i < n; i++) {
598                    if (FastMath.abs(realEigenvalues[i]) <= Precision.EPSILON * maxAbsoluteValue) {
599                        realEigenvalues[i] = 0;
600                    }
601                    if (FastMath.abs(e[i]) <= Precision.EPSILON * maxAbsoluteValue) {
602                        e[i]=0;
603                    }
604                }
605            }
606    
607            for (int j = 0; j < n; j++) {
608                int its = 0;
609                int m;
610                do {
611                    for (m = j; m < n - 1; m++) {
612                        double delta = FastMath.abs(realEigenvalues[m]) +
613                            FastMath.abs(realEigenvalues[m + 1]);
614                        if (FastMath.abs(e[m]) + delta == delta) {
615                            break;
616                        }
617                    }
618                    if (m != j) {
619                        if (its == maxIter) {
620                            throw new MaxCountExceededException(LocalizedFormats.CONVERGENCE_FAILED,
621                                                                maxIter);
622                        }
623                        its++;
624                        double q = (realEigenvalues[j + 1] - realEigenvalues[j]) / (2 * e[j]);
625                        double t = FastMath.sqrt(1 + q * q);
626                        if (q < 0.0) {
627                            q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q - t);
628                        } else {
629                            q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q + t);
630                        }
631                        double u = 0.0;
632                        double s = 1.0;
633                        double c = 1.0;
634                        int i;
635                        for (i = m - 1; i >= j; i--) {
636                            double p = s * e[i];
637                            double h = c * e[i];
638                            if (FastMath.abs(p) >= FastMath.abs(q)) {
639                                c = q / p;
640                                t = FastMath.sqrt(c * c + 1.0);
641                                e[i + 1] = p * t;
642                                s = 1.0 / t;
643                                c = c * s;
644                            } else {
645                                s = p / q;
646                                t = FastMath.sqrt(s * s + 1.0);
647                                e[i + 1] = q * t;
648                                c = 1.0 / t;
649                                s = s * c;
650                            }
651                            if (e[i + 1] == 0.0) {
652                                realEigenvalues[i + 1] -= u;
653                                e[m] = 0.0;
654                                break;
655                            }
656                            q = realEigenvalues[i + 1] - u;
657                            t = (realEigenvalues[i] - q) * s + 2.0 * c * h;
658                            u = s * t;
659                            realEigenvalues[i + 1] = q + u;
660                            q = c * t - h;
661                            for (int ia = 0; ia < n; ia++) {
662                                p = z[ia][i + 1];
663                                z[ia][i + 1] = s * z[ia][i] + c * p;
664                                z[ia][i] = c * z[ia][i] - s * p;
665                            }
666                        }
667                        if (t == 0.0 && i >= j) {
668                            continue;
669                        }
670                        realEigenvalues[j] -= u;
671                        e[j] = q;
672                        e[m] = 0.0;
673                    }
674                } while (m != j);
675            }
676    
677            //Sort the eigen values (and vectors) in increase order
678            for (int i = 0; i < n; i++) {
679                int k = i;
680                double p = realEigenvalues[i];
681                for (int j = i + 1; j < n; j++) {
682                    if (realEigenvalues[j] > p) {
683                        k = j;
684                        p = realEigenvalues[j];
685                    }
686                }
687                if (k != i) {
688                    realEigenvalues[k] = realEigenvalues[i];
689                    realEigenvalues[i] = p;
690                    for (int j = 0; j < n; j++) {
691                        p = z[j][i];
692                        z[j][i] = z[j][k];
693                        z[j][k] = p;
694                    }
695                }
696            }
697    
698            // Determine the largest eigen value in absolute term.
699            maxAbsoluteValue = 0;
700            for (int i = 0; i < n; i++) {
701                if (FastMath.abs(realEigenvalues[i]) > maxAbsoluteValue) {
702                    maxAbsoluteValue=FastMath.abs(realEigenvalues[i]);
703                }
704            }
705            // Make null any eigen value too small to be significant
706            if (maxAbsoluteValue != 0.0) {
707                for (int i=0; i < n; i++) {
708                    if (FastMath.abs(realEigenvalues[i]) < Precision.EPSILON * maxAbsoluteValue) {
709                        realEigenvalues[i] = 0;
710                    }
711                }
712            }
713            eigenvectors = new ArrayRealVector[n];
714            final double[] tmp = new double[n];
715            for (int i = 0; i < n; i++) {
716                for (int j = 0; j < n; j++) {
717                    tmp[j] = z[j][i];
718                }
719                eigenvectors[i] = new ArrayRealVector(tmp);
720            }
721        }
722    
723        /**
724         * Transforms the matrix to Schur form and calculates the eigenvalues.
725         *
726         * @param matrix Matrix to transform.
727         * @return the {@link SchurTransformer Shur transform} for this matrix
728         */
729        private SchurTransformer transformToSchur(final RealMatrix matrix) {
730            final SchurTransformer schurTransform = new SchurTransformer(matrix);
731            final double[][] matT = schurTransform.getT().getData();
732    
733            realEigenvalues = new double[matT.length];
734            imagEigenvalues = new double[matT.length];
735    
736            for (int i = 0; i < realEigenvalues.length; i++) {
737                if (i == (realEigenvalues.length - 1) ||
738                    Precision.equals(matT[i + 1][i], 0.0, EPSILON)) {
739                    realEigenvalues[i] = matT[i][i];
740                } else {
741                    final double x = matT[i + 1][i + 1];
742                    final double p = 0.5 * (matT[i][i] - x);
743                    final double z = FastMath.sqrt(FastMath.abs(p * p + matT[i + 1][i] * matT[i][i + 1]));
744                    realEigenvalues[i] = x + p;
745                    imagEigenvalues[i] = z;
746                    realEigenvalues[i + 1] = x + p;
747                    imagEigenvalues[i + 1] = -z;
748                    i++;
749                }
750            }
751            return schurTransform;
752        }
753    
754        /**
755         * Performs a division of two complex numbers.
756         *
757         * @param xr real part of the first number
758         * @param xi imaginary part of the first number
759         * @param yr real part of the second number
760         * @param yi imaginary part of the second number
761         * @return result of the complex division
762         */
763        private Complex cdiv(final double xr, final double xi,
764                             final double yr, final double yi) {
765            return new Complex(xr, xi).divide(new Complex(yr, yi));
766        }
767    
768        /**
769         * Find eigenvectors from a matrix transformed to Schur form.
770         *
771         * @param schur the schur transformation of the matrix
772         * @throws MathArithmeticException if the Schur form has a norm of zero
773         */
774        private void findEigenVectorsFromSchur(final SchurTransformer schur)
775            throws MathArithmeticException {
776            final double[][] matrixT = schur.getT().getData();
777            final double[][] matrixP = schur.getP().getData();
778    
779            final int n = matrixT.length;
780    
781            // compute matrix norm
782            double norm = 0.0;
783            for (int i = 0; i < n; i++) {
784               for (int j = FastMath.max(i - 1, 0); j < n; j++) {
785                  norm = norm + FastMath.abs(matrixT[i][j]);
786               }
787            }
788    
789            // we can not handle a matrix with zero norm
790            if (Precision.equals(norm, 0.0, EPSILON)) {
791               throw new MathArithmeticException(LocalizedFormats.ZERO_NORM);
792            }
793    
794            // Backsubstitute to find vectors of upper triangular form
795    
796            double r = 0.0;
797            double s = 0.0;
798            double z = 0.0;
799    
800            for (int idx = n - 1; idx >= 0; idx--) {
801                double p = realEigenvalues[idx];
802                double q = imagEigenvalues[idx];
803    
804                if (Precision.equals(q, 0.0)) {
805                    // Real vector
806                    int l = idx;
807                    matrixT[idx][idx] = 1.0;
808                    for (int i = idx - 1; i >= 0; i--) {
809                        double w = matrixT[i][i] - p;
810                        r = 0.0;
811                        for (int j = l; j <= idx; j++) {
812                            r = r + matrixT[i][j] * matrixT[j][idx];
813                        }
814                        if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) < 0.0) {
815                            z = w;
816                            s = r;
817                        } else {
818                            l = i;
819                            if (Precision.equals(imagEigenvalues[i], 0.0)) {
820                                if (w != 0.0) {
821                                    matrixT[i][idx] = -r / w;
822                                } else {
823                                    matrixT[i][idx] = -r / (Precision.EPSILON * norm);
824                                }
825                            } else {
826                                // Solve real equations
827                                double x = matrixT[i][i + 1];
828                                double y = matrixT[i + 1][i];
829                                q = (realEigenvalues[i] - p) * (realEigenvalues[i] - p) +
830                                    imagEigenvalues[i] * imagEigenvalues[i];
831                                double t = (x * s - z * r) / q;
832                                matrixT[i][idx] = t;
833                                if (FastMath.abs(x) > FastMath.abs(z)) {
834                                    matrixT[i + 1][idx] = (-r - w * t) / x;
835                                } else {
836                                    matrixT[i + 1][idx] = (-s - y * t) / z;
837                                }
838                            }
839    
840                            // Overflow control
841                            double t = FastMath.abs(matrixT[i][idx]);
842                            if ((Precision.EPSILON * t) * t > 1) {
843                                for (int j = i; j <= idx; j++) {
844                                    matrixT[j][idx] = matrixT[j][idx] / t;
845                                }
846                            }
847                        }
848                    }
849                } else if (q < 0.0) {
850                    // Complex vector
851                    int l = idx - 1;
852    
853                    // Last vector component imaginary so matrix is triangular
854                    if (FastMath.abs(matrixT[idx][idx - 1]) > FastMath.abs(matrixT[idx - 1][idx])) {
855                        matrixT[idx - 1][idx - 1] = q / matrixT[idx][idx - 1];
856                        matrixT[idx - 1][idx]     = -(matrixT[idx][idx] - p) / matrixT[idx][idx - 1];
857                    } else {
858                        final Complex result = cdiv(0.0, -matrixT[idx - 1][idx],
859                                                    matrixT[idx - 1][idx - 1] - p, q);
860                        matrixT[idx - 1][idx - 1] = result.getReal();
861                        matrixT[idx - 1][idx]     = result.getImaginary();
862                    }
863    
864                    matrixT[idx][idx - 1] = 0.0;
865                    matrixT[idx][idx]     = 1.0;
866    
867                    for (int i = idx - 2; i >= 0; i--) {
868                        double ra = 0.0;
869                        double sa = 0.0;
870                        for (int j = l; j <= idx; j++) {
871                            ra = ra + matrixT[i][j] * matrixT[j][idx - 1];
872                            sa = sa + matrixT[i][j] * matrixT[j][idx];
873                        }
874                        double w = matrixT[i][i] - p;
875    
876                        if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) < 0.0) {
877                            z = w;
878                            r = ra;
879                            s = sa;
880                        } else {
881                            l = i;
882                            if (Precision.equals(imagEigenvalues[i], 0.0)) {
883                                final Complex c = cdiv(-ra, -sa, w, q);
884                                matrixT[i][idx - 1] = c.getReal();
885                                matrixT[i][idx] = c.getImaginary();
886                            } else {
887                                // Solve complex equations
888                                double x = matrixT[i][i + 1];
889                                double y = matrixT[i + 1][i];
890                                double vr = (realEigenvalues[i] - p) * (realEigenvalues[i] - p) +
891                                            imagEigenvalues[i] * imagEigenvalues[i] - q * q;
892                                final double vi = (realEigenvalues[i] - p) * 2.0 * q;
893                                if (Precision.equals(vr, 0.0) && Precision.equals(vi, 0.0)) {
894                                    vr = Precision.EPSILON * norm *
895                                         (FastMath.abs(w) + FastMath.abs(q) + FastMath.abs(x) +
896                                          FastMath.abs(y) + FastMath.abs(z));
897                                }
898                                final Complex c     = cdiv(x * r - z * ra + q * sa,
899                                                           x * s - z * sa - q * ra, vr, vi);
900                                matrixT[i][idx - 1] = c.getReal();
901                                matrixT[i][idx]     = c.getImaginary();
902    
903                                if (FastMath.abs(x) > (FastMath.abs(z) + FastMath.abs(q))) {
904                                    matrixT[i + 1][idx - 1] = (-ra - w * matrixT[i][idx - 1] +
905                                                               q * matrixT[i][idx]) / x;
906                                    matrixT[i + 1][idx]     = (-sa - w * matrixT[i][idx] -
907                                                               q * matrixT[i][idx - 1]) / x;
908                                } else {
909                                    final Complex c2        = cdiv(-r - y * matrixT[i][idx - 1],
910                                                                   -s - y * matrixT[i][idx], z, q);
911                                    matrixT[i + 1][idx - 1] = c2.getReal();
912                                    matrixT[i + 1][idx]     = c2.getImaginary();
913                                }
914                            }
915    
916                            // Overflow control
917                            double t = FastMath.max(FastMath.abs(matrixT[i][idx - 1]),
918                                                    FastMath.abs(matrixT[i][idx]));
919                            if ((Precision.EPSILON * t) * t > 1) {
920                                for (int j = i; j <= idx; j++) {
921                                    matrixT[j][idx - 1] = matrixT[j][idx - 1] / t;
922                                    matrixT[j][idx]     = matrixT[j][idx] / t;
923                                }
924                            }
925                        }
926                    }
927                }
928            }
929    
930            // Vectors of isolated roots
931            for (int i = 0; i < n; i++) {
932                if (i < 0 | i > n - 1) {
933                    for (int j = i; j < n; j++) {
934                        matrixP[i][j] = matrixT[i][j];
935                    }
936                }
937            }
938    
939            // Back transformation to get eigenvectors of original matrix
940            for (int j = n - 1; j >= 0; j--) {
941                for (int i = 0; i <= n - 1; i++) {
942                    z = 0.0;
943                    for (int k = 0; k <= FastMath.min(j, n - 1); k++) {
944                        z = z + matrixP[i][k] * matrixT[k][j];
945                    }
946                    matrixP[i][j] = z;
947                }
948            }
949    
950            eigenvectors = new ArrayRealVector[n];
951            final double[] tmp = new double[n];
952            for (int i = 0; i < n; i++) {
953                for (int j = 0; j < n; j++) {
954                    tmp[j] = matrixP[j][i];
955                }
956                eigenvectors[i] = new ArrayRealVector(tmp);
957            }
958        }
959    }