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 * @version $Id: EigenDecomposition.java 1537616 2013-10-31 20:07:38Z tn $
075 * @since 2.0 (changed to concrete class in 3.0)
076 */
077public 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 positive
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            double largestEigenvalueNorm = 0.0;
517            // Looping over all values (in case they are not sorted in decreasing
518            // order of their norm).
519            for (int i = 0; i < realEigenvalues.length; ++i) {
520                largestEigenvalueNorm = FastMath.max(largestEigenvalueNorm, eigenvalueNorm(i));
521            }
522            // Corner case: zero matrix, all exactly 0 eigenvalues
523            if (largestEigenvalueNorm == 0.0) {
524                return false;
525            }
526            for (int i = 0; i < realEigenvalues.length; ++i) {
527                // Looking for eigenvalues that are 0, where we consider anything much much smaller
528                // than the largest eigenvalue to be effectively 0.
529                if (Precision.equals(eigenvalueNorm(i) / largestEigenvalueNorm, 0, EPSILON)) {
530                    return false;
531                }
532            }
533            return true;
534        }
535
536        /**
537         * @param i which eigenvalue to find the norm of
538         * @return the norm of ith (complex) eigenvalue.
539         */
540        private double eigenvalueNorm(int i) {
541            final double re = realEigenvalues[i];
542            final double im = imagEigenvalues[i];
543            return FastMath.sqrt(re * re + im * im);
544        }
545
546        /**
547         * Get the inverse of the decomposed matrix.
548         *
549         * @return the inverse matrix.
550         * @throws SingularMatrixException if the decomposed matrix is singular.
551         */
552        public RealMatrix getInverse() {
553            if (!isNonSingular()) {
554                throw new SingularMatrixException();
555            }
556
557            final int m = realEigenvalues.length;
558            final double[][] invData = new double[m][m];
559
560            for (int i = 0; i < m; ++i) {
561                final double[] invI = invData[i];
562                for (int j = 0; j < m; ++j) {
563                    double invIJ = 0;
564                    for (int k = 0; k < m; ++k) {
565                        final double[] vK = eigenvectors[k].getDataRef();
566                        invIJ += vK[i] * vK[j] / realEigenvalues[k];
567                    }
568                    invI[j] = invIJ;
569                }
570            }
571            return MatrixUtils.createRealMatrix(invData);
572        }
573    }
574
575    /**
576     * Transforms the matrix to tridiagonal form.
577     *
578     * @param matrix Matrix to transform.
579     */
580    private void transformToTridiagonal(final RealMatrix matrix) {
581        // transform the matrix to tridiagonal
582        transformer = new TriDiagonalTransformer(matrix);
583        main = transformer.getMainDiagonalRef();
584        secondary = transformer.getSecondaryDiagonalRef();
585    }
586
587    /**
588     * Find eigenvalues and eigenvectors (Dubrulle et al., 1971)
589     *
590     * @param householderMatrix Householder matrix of the transformation
591     * to tridiagonal form.
592     */
593    private void findEigenVectors(final double[][] householderMatrix) {
594        final double[][]z = householderMatrix.clone();
595        final int n = main.length;
596        realEigenvalues = new double[n];
597        imagEigenvalues = new double[n];
598        final double[] e = new double[n];
599        for (int i = 0; i < n - 1; i++) {
600            realEigenvalues[i] = main[i];
601            e[i] = secondary[i];
602        }
603        realEigenvalues[n - 1] = main[n - 1];
604        e[n - 1] = 0;
605
606        // Determine the largest main and secondary value in absolute term.
607        double maxAbsoluteValue = 0;
608        for (int i = 0; i < n; i++) {
609            if (FastMath.abs(realEigenvalues[i]) > maxAbsoluteValue) {
610                maxAbsoluteValue = FastMath.abs(realEigenvalues[i]);
611            }
612            if (FastMath.abs(e[i]) > maxAbsoluteValue) {
613                maxAbsoluteValue = FastMath.abs(e[i]);
614            }
615        }
616        // Make null any main and secondary value too small to be significant
617        if (maxAbsoluteValue != 0) {
618            for (int i=0; i < n; i++) {
619                if (FastMath.abs(realEigenvalues[i]) <= Precision.EPSILON * maxAbsoluteValue) {
620                    realEigenvalues[i] = 0;
621                }
622                if (FastMath.abs(e[i]) <= Precision.EPSILON * maxAbsoluteValue) {
623                    e[i]=0;
624                }
625            }
626        }
627
628        for (int j = 0; j < n; j++) {
629            int its = 0;
630            int m;
631            do {
632                for (m = j; m < n - 1; m++) {
633                    double delta = FastMath.abs(realEigenvalues[m]) +
634                        FastMath.abs(realEigenvalues[m + 1]);
635                    if (FastMath.abs(e[m]) + delta == delta) {
636                        break;
637                    }
638                }
639                if (m != j) {
640                    if (its == maxIter) {
641                        throw new MaxCountExceededException(LocalizedFormats.CONVERGENCE_FAILED,
642                                                            maxIter);
643                    }
644                    its++;
645                    double q = (realEigenvalues[j + 1] - realEigenvalues[j]) / (2 * e[j]);
646                    double t = FastMath.sqrt(1 + q * q);
647                    if (q < 0.0) {
648                        q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q - t);
649                    } else {
650                        q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q + t);
651                    }
652                    double u = 0.0;
653                    double s = 1.0;
654                    double c = 1.0;
655                    int i;
656                    for (i = m - 1; i >= j; i--) {
657                        double p = s * e[i];
658                        double h = c * e[i];
659                        if (FastMath.abs(p) >= FastMath.abs(q)) {
660                            c = q / p;
661                            t = FastMath.sqrt(c * c + 1.0);
662                            e[i + 1] = p * t;
663                            s = 1.0 / t;
664                            c = c * s;
665                        } else {
666                            s = p / q;
667                            t = FastMath.sqrt(s * s + 1.0);
668                            e[i + 1] = q * t;
669                            c = 1.0 / t;
670                            s = s * c;
671                        }
672                        if (e[i + 1] == 0.0) {
673                            realEigenvalues[i + 1] -= u;
674                            e[m] = 0.0;
675                            break;
676                        }
677                        q = realEigenvalues[i + 1] - u;
678                        t = (realEigenvalues[i] - q) * s + 2.0 * c * h;
679                        u = s * t;
680                        realEigenvalues[i + 1] = q + u;
681                        q = c * t - h;
682                        for (int ia = 0; ia < n; ia++) {
683                            p = z[ia][i + 1];
684                            z[ia][i + 1] = s * z[ia][i] + c * p;
685                            z[ia][i] = c * z[ia][i] - s * p;
686                        }
687                    }
688                    if (t == 0.0 && i >= j) {
689                        continue;
690                    }
691                    realEigenvalues[j] -= u;
692                    e[j] = q;
693                    e[m] = 0.0;
694                }
695            } while (m != j);
696        }
697
698        //Sort the eigen values (and vectors) in increase order
699        for (int i = 0; i < n; i++) {
700            int k = i;
701            double p = realEigenvalues[i];
702            for (int j = i + 1; j < n; j++) {
703                if (realEigenvalues[j] > p) {
704                    k = j;
705                    p = realEigenvalues[j];
706                }
707            }
708            if (k != i) {
709                realEigenvalues[k] = realEigenvalues[i];
710                realEigenvalues[i] = p;
711                for (int j = 0; j < n; j++) {
712                    p = z[j][i];
713                    z[j][i] = z[j][k];
714                    z[j][k] = p;
715                }
716            }
717        }
718
719        // Determine the largest eigen value in absolute term.
720        maxAbsoluteValue = 0;
721        for (int i = 0; i < n; i++) {
722            if (FastMath.abs(realEigenvalues[i]) > maxAbsoluteValue) {
723                maxAbsoluteValue=FastMath.abs(realEigenvalues[i]);
724            }
725        }
726        // Make null any eigen value too small to be significant
727        if (maxAbsoluteValue != 0.0) {
728            for (int i=0; i < n; i++) {
729                if (FastMath.abs(realEigenvalues[i]) < Precision.EPSILON * maxAbsoluteValue) {
730                    realEigenvalues[i] = 0;
731                }
732            }
733        }
734        eigenvectors = new ArrayRealVector[n];
735        final double[] tmp = new double[n];
736        for (int i = 0; i < n; i++) {
737            for (int j = 0; j < n; j++) {
738                tmp[j] = z[j][i];
739            }
740            eigenvectors[i] = new ArrayRealVector(tmp);
741        }
742    }
743
744    /**
745     * Transforms the matrix to Schur form and calculates the eigenvalues.
746     *
747     * @param matrix Matrix to transform.
748     * @return the {@link SchurTransformer Shur transform} for this matrix
749     */
750    private SchurTransformer transformToSchur(final RealMatrix matrix) {
751        final SchurTransformer schurTransform = new SchurTransformer(matrix);
752        final double[][] matT = schurTransform.getT().getData();
753
754        realEigenvalues = new double[matT.length];
755        imagEigenvalues = new double[matT.length];
756
757        for (int i = 0; i < realEigenvalues.length; i++) {
758            if (i == (realEigenvalues.length - 1) ||
759                Precision.equals(matT[i + 1][i], 0.0, EPSILON)) {
760                realEigenvalues[i] = matT[i][i];
761            } else {
762                final double x = matT[i + 1][i + 1];
763                final double p = 0.5 * (matT[i][i] - x);
764                final double z = FastMath.sqrt(FastMath.abs(p * p + matT[i + 1][i] * matT[i][i + 1]));
765                realEigenvalues[i] = x + p;
766                imagEigenvalues[i] = z;
767                realEigenvalues[i + 1] = x + p;
768                imagEigenvalues[i + 1] = -z;
769                i++;
770            }
771        }
772        return schurTransform;
773    }
774
775    /**
776     * Performs a division of two complex numbers.
777     *
778     * @param xr real part of the first number
779     * @param xi imaginary part of the first number
780     * @param yr real part of the second number
781     * @param yi imaginary part of the second number
782     * @return result of the complex division
783     */
784    private Complex cdiv(final double xr, final double xi,
785                         final double yr, final double yi) {
786        return new Complex(xr, xi).divide(new Complex(yr, yi));
787    }
788
789    /**
790     * Find eigenvectors from a matrix transformed to Schur form.
791     *
792     * @param schur the schur transformation of the matrix
793     * @throws MathArithmeticException if the Schur form has a norm of zero
794     */
795    private void findEigenVectorsFromSchur(final SchurTransformer schur)
796        throws MathArithmeticException {
797        final double[][] matrixT = schur.getT().getData();
798        final double[][] matrixP = schur.getP().getData();
799
800        final int n = matrixT.length;
801
802        // compute matrix norm
803        double norm = 0.0;
804        for (int i = 0; i < n; i++) {
805           for (int j = FastMath.max(i - 1, 0); j < n; j++) {
806              norm = norm + FastMath.abs(matrixT[i][j]);
807           }
808        }
809
810        // we can not handle a matrix with zero norm
811        if (Precision.equals(norm, 0.0, EPSILON)) {
812           throw new MathArithmeticException(LocalizedFormats.ZERO_NORM);
813        }
814
815        // Backsubstitute to find vectors of upper triangular form
816
817        double r = 0.0;
818        double s = 0.0;
819        double z = 0.0;
820
821        for (int idx = n - 1; idx >= 0; idx--) {
822            double p = realEigenvalues[idx];
823            double q = imagEigenvalues[idx];
824
825            if (Precision.equals(q, 0.0)) {
826                // Real vector
827                int l = idx;
828                matrixT[idx][idx] = 1.0;
829                for (int i = idx - 1; i >= 0; i--) {
830                    double w = matrixT[i][i] - p;
831                    r = 0.0;
832                    for (int j = l; j <= idx; j++) {
833                        r = r + matrixT[i][j] * matrixT[j][idx];
834                    }
835                    if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) < 0) {
836                        z = w;
837                        s = r;
838                    } else {
839                        l = i;
840                        if (Precision.equals(imagEigenvalues[i], 0.0)) {
841                            if (w != 0.0) {
842                                matrixT[i][idx] = -r / w;
843                            } else {
844                                matrixT[i][idx] = -r / (Precision.EPSILON * norm);
845                            }
846                        } else {
847                            // Solve real equations
848                            double x = matrixT[i][i + 1];
849                            double y = matrixT[i + 1][i];
850                            q = (realEigenvalues[i] - p) * (realEigenvalues[i] - p) +
851                                imagEigenvalues[i] * imagEigenvalues[i];
852                            double t = (x * s - z * r) / q;
853                            matrixT[i][idx] = t;
854                            if (FastMath.abs(x) > FastMath.abs(z)) {
855                                matrixT[i + 1][idx] = (-r - w * t) / x;
856                            } else {
857                                matrixT[i + 1][idx] = (-s - y * t) / z;
858                            }
859                        }
860
861                        // Overflow control
862                        double t = FastMath.abs(matrixT[i][idx]);
863                        if ((Precision.EPSILON * t) * t > 1) {
864                            for (int j = i; j <= idx; j++) {
865                                matrixT[j][idx] = matrixT[j][idx] / t;
866                            }
867                        }
868                    }
869                }
870            } else if (q < 0.0) {
871                // Complex vector
872                int l = idx - 1;
873
874                // Last vector component imaginary so matrix is triangular
875                if (FastMath.abs(matrixT[idx][idx - 1]) > FastMath.abs(matrixT[idx - 1][idx])) {
876                    matrixT[idx - 1][idx - 1] = q / matrixT[idx][idx - 1];
877                    matrixT[idx - 1][idx]     = -(matrixT[idx][idx] - p) / matrixT[idx][idx - 1];
878                } else {
879                    final Complex result = cdiv(0.0, -matrixT[idx - 1][idx],
880                                                matrixT[idx - 1][idx - 1] - p, q);
881                    matrixT[idx - 1][idx - 1] = result.getReal();
882                    matrixT[idx - 1][idx]     = result.getImaginary();
883                }
884
885                matrixT[idx][idx - 1] = 0.0;
886                matrixT[idx][idx]     = 1.0;
887
888                for (int i = idx - 2; i >= 0; i--) {
889                    double ra = 0.0;
890                    double sa = 0.0;
891                    for (int j = l; j <= idx; j++) {
892                        ra = ra + matrixT[i][j] * matrixT[j][idx - 1];
893                        sa = sa + matrixT[i][j] * matrixT[j][idx];
894                    }
895                    double w = matrixT[i][i] - p;
896
897                    if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) < 0) {
898                        z = w;
899                        r = ra;
900                        s = sa;
901                    } else {
902                        l = i;
903                        if (Precision.equals(imagEigenvalues[i], 0.0)) {
904                            final Complex c = cdiv(-ra, -sa, w, q);
905                            matrixT[i][idx - 1] = c.getReal();
906                            matrixT[i][idx] = c.getImaginary();
907                        } else {
908                            // Solve complex equations
909                            double x = matrixT[i][i + 1];
910                            double y = matrixT[i + 1][i];
911                            double vr = (realEigenvalues[i] - p) * (realEigenvalues[i] - p) +
912                                        imagEigenvalues[i] * imagEigenvalues[i] - q * q;
913                            final double vi = (realEigenvalues[i] - p) * 2.0 * q;
914                            if (Precision.equals(vr, 0.0) && Precision.equals(vi, 0.0)) {
915                                vr = Precision.EPSILON * norm *
916                                     (FastMath.abs(w) + FastMath.abs(q) + FastMath.abs(x) +
917                                      FastMath.abs(y) + FastMath.abs(z));
918                            }
919                            final Complex c     = cdiv(x * r - z * ra + q * sa,
920                                                       x * s - z * sa - q * ra, vr, vi);
921                            matrixT[i][idx - 1] = c.getReal();
922                            matrixT[i][idx]     = c.getImaginary();
923
924                            if (FastMath.abs(x) > (FastMath.abs(z) + FastMath.abs(q))) {
925                                matrixT[i + 1][idx - 1] = (-ra - w * matrixT[i][idx - 1] +
926                                                           q * matrixT[i][idx]) / x;
927                                matrixT[i + 1][idx]     = (-sa - w * matrixT[i][idx] -
928                                                           q * matrixT[i][idx - 1]) / x;
929                            } else {
930                                final Complex c2        = cdiv(-r - y * matrixT[i][idx - 1],
931                                                               -s - y * matrixT[i][idx], z, q);
932                                matrixT[i + 1][idx - 1] = c2.getReal();
933                                matrixT[i + 1][idx]     = c2.getImaginary();
934                            }
935                        }
936
937                        // Overflow control
938                        double t = FastMath.max(FastMath.abs(matrixT[i][idx - 1]),
939                                                FastMath.abs(matrixT[i][idx]));
940                        if ((Precision.EPSILON * t) * t > 1) {
941                            for (int j = i; j <= idx; j++) {
942                                matrixT[j][idx - 1] = matrixT[j][idx - 1] / t;
943                                matrixT[j][idx]     = matrixT[j][idx] / t;
944                            }
945                        }
946                    }
947                }
948            }
949        }
950
951        // Back transformation to get eigenvectors of original matrix
952        for (int j = n - 1; j >= 0; j--) {
953            for (int i = 0; i <= n - 1; i++) {
954                z = 0.0;
955                for (int k = 0; k <= FastMath.min(j, n - 1); k++) {
956                    z = z + matrixP[i][k] * matrixT[k][j];
957                }
958                matrixP[i][j] = z;
959            }
960        }
961
962        eigenvectors = new ArrayRealVector[n];
963        final double[] tmp = new double[n];
964        for (int i = 0; i < n; i++) {
965            for (int j = 0; j < n; j++) {
966                tmp[j] = matrixP[j][i];
967            }
968            eigenvectors[i] = new ArrayRealVector(tmp);
969        }
970    }
971}