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