    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    }