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.math.linear;
019    
020    import java.util.ArrayList;
021    import java.util.Arrays;
022    import java.util.List;
023    
024    import org.apache.commons.math.ConvergenceException;
025    import org.apache.commons.math.MathRuntimeException;
026    import org.apache.commons.math.MaxIterationsExceededException;
027    import org.apache.commons.math.util.MathUtils;
028    
029    /**
030     * Calculates the eigen decomposition of a <strong>symmetric</strong> matrix.
031     * <p>The eigen decomposition of matrix A is a set of two matrices:
032     * V and D such that A = V D V<sup>T</sup>. A, V and D are all m &times; m
033     * matrices.</p>
034     * <p>As of 2.0, this class supports only <strong>symmetric</strong> matrices,
035     * and hence computes only real realEigenvalues. This implies the D matrix returned by
036     * {@link #getD()} is always diagonal and the imaginary values returned {@link
037     * #getImagEigenvalue(int)} and {@link #getImagEigenvalues()} are always null.</p>
038     * <p>When called with a {@link RealMatrix} argument, this implementation only uses
039     * the upper part of the matrix, the part below the diagonal is not accessed at all.</p>
040     * <p>Eigenvalues are computed as soon as the matrix is decomposed, but eigenvectors
041     * are computed only when required, i.e. only when one of the {@link #getEigenvector(int)},
042     * {@link #getV()}, {@link #getVT()}, {@link #getSolver()} methods is called.</p>
043     * <p>This implementation is based on Inderjit Singh Dhillon thesis
044     * <a href="http://www.cs.utexas.edu/users/inderjit/public_papers/thesis.pdf">A
045     * New O(n<sup>2</sup>) Algorithm for the Symmetric Tridiagonal Eigenvalue/Eigenvector
046     * Problem</a>, on Beresford N. Parlett and Osni A. Marques paper <a
047     * href="http://www.netlib.org/lapack/lawnspdf/lawn155.pdf">An Implementation of the
048     * dqds Algorithm (Positive Case)</a> and on the corresponding LAPACK routines (DLARRE,
049     * DLASQ2, DLAZQ3, DLAZQ4, DLASQ5 and DLASQ6).</p>
050     * @author Beresford Parlett, University of California, Berkeley, USA (fortran version)
051     * @author Jim Demmel, University of California, Berkeley, USA (fortran version)
052     * @author Inderjit Dhillon, University of Texas, Austin, USA(fortran version)
053     * @author Osni Marques, LBNL/NERSC, USA (fortran version)
054     * @author Christof Voemel, University of California, Berkeley, USA(fortran version)
055     * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
056     * @since 2.0
057     */
058    public class EigenDecompositionImpl implements EigenDecomposition {
059    
060        /** Tolerance. */
061        private static final double TOLERANCE = 100 * MathUtils.EPSILON;
062    
063        /** Squared tolerance. */
064        private static final double TOLERANCE_2 = TOLERANCE * TOLERANCE;
065    
066        /** Split tolerance. */
067        private double splitTolerance;
068    
069        /** Main diagonal of the tridiagonal matrix. */
070        private double[] main;
071    
072        /** Secondary diagonal of the tridiagonal matrix. */
073        private double[] secondary;
074    
075        /** Squared secondary diagonal of the tridiagonal matrix. */
076        private double[] squaredSecondary;
077    
078        /** Transformer to tridiagonal (may be null if matrix is already tridiagonal). */
079        private TriDiagonalTransformer transformer;
080    
081        /** Lower bound of spectra. */
082        private double lowerSpectra;
083    
084        /** Upper bound of spectra. */
085        private double upperSpectra;
086    
087        /** Minimum pivot in the Sturm sequence. */
088        private double minPivot;
089    
090        /** Current shift. */
091        private double sigma;
092    
093        /** Low part of the current shift. */
094        private double sigmaLow;
095    
096        /** Shift increment to apply. */
097        private double tau;
098    
099        /** Work array for all decomposition algorithms. */
100        private double[] work;
101    
102        /** Shift within qd array for ping-pong implementation. */
103        private int pingPong;
104    
105        /** Max value of diagonal elements in current segment. */
106        private double qMax;
107    
108        /** Min value of off-diagonal elements in current segment. */
109        private double eMin;
110    
111        /** Type of the last dqds shift. */
112        private int    tType;
113    
114        /** Minimal value on current state of the diagonal. */
115        private double dMin;
116    
117        /** Minimal value on current state of the diagonal, excluding last element. */
118        private double dMin1;
119    
120        /** Minimal value on current state of the diagonal, excluding last two elements. */
121        private double dMin2;
122    
123        /** Last value on current state of the diagonal. */
124        private double dN;
125    
126        /** Last but one value on current state of the diagonal. */
127        private double dN1;
128    
129        /** Last but two on current state of the diagonal. */
130        private double dN2;
131    
132        /** Shift ratio with respect to dMin used when tType == 6. */
133        private double g;
134    
135        /** Real part of the realEigenvalues. */
136        private double[] realEigenvalues;
137    
138        /** Imaginary part of the realEigenvalues. */
139        private double[] imagEigenvalues;
140    
141        /** Eigenvectors. */
142        private ArrayRealVector[] eigenvectors;
143    
144        /** Cached value of V. */
145        private RealMatrix cachedV;
146    
147        /** Cached value of D. */
148        private RealMatrix cachedD;
149    
150        /** Cached value of Vt. */
151        private RealMatrix cachedVt;
152    
153        /**
154         * Calculates the eigen decomposition of the given symmetric matrix. 
155         * @param matrix The <strong>symmetric</strong> matrix to decompose.
156         * @param splitTolerance tolerance on the off-diagonal elements relative to the
157         * geometric mean to split the tridiagonal matrix (a suggested value is
158         * {@link MathUtils#SAFE_MIN})
159         * @exception InvalidMatrixException (wrapping a {@link ConvergenceException}
160         * if algorithm fails to converge
161         */
162        public EigenDecompositionImpl(final RealMatrix matrix,
163                                      final double splitTolerance)
164            throws InvalidMatrixException {
165            if (isSymmetric(matrix)) {
166                this.splitTolerance = splitTolerance;
167                transformToTridiagonal(matrix);
168                decompose();
169            } else {
170                // as of 2.0, non-symmetric matrices (i.e. complex eigenvalues) are NOT supported
171                // see issue https://issues.apache.org/jira/browse/MATH-235
172                throw new InvalidMatrixException("eigen decomposition of assymetric matrices not supported yet");
173            }
174        }
175    
176        /**
177         * Calculates the eigen decomposition of the given tridiagonal symmetric matrix. 
178         * @param main the main diagonal of the matrix (will be copied)
179         * @param secondary the secondary diagonal of the matrix (will be copied)
180         * @param splitTolerance tolerance on the off-diagonal elements relative to the
181         * geometric mean to split the tridiagonal matrix (a suggested value is
182         * {@link MathUtils#SAFE_MIN})
183         * @exception InvalidMatrixException (wrapping a {@link ConvergenceException}
184         * if algorithm fails to converge
185         */
186        public EigenDecompositionImpl(final double[] main, double[] secondary,
187                final double splitTolerance)
188            throws InvalidMatrixException {
189    
190            this.main      = main.clone();
191            this.secondary = secondary.clone();
192            transformer    = null;
193    
194            // pre-compute some elements
195            squaredSecondary = new double[secondary.length];
196            for (int i = 0; i < squaredSecondary.length; ++i) {
197                final double s = secondary[i];
198                squaredSecondary[i] = s * s;
199            }
200    
201            this.splitTolerance = splitTolerance;
202            decompose();
203    
204        }
205    
206        /**
207         * Check if a matrix is symmetric.
208         * @param matrix matrix to check
209         * @return true if matrix is symmetric
210         */
211        private boolean isSymmetric(final RealMatrix matrix) {
212            final int rows    = matrix.getRowDimension();
213            final int columns = matrix.getColumnDimension();
214            final double eps  = 10 * rows * columns * MathUtils.EPSILON;
215            for (int i = 0; i < rows; ++i) {
216                for (int j = i + 1; j < columns; ++j) {
217                    final double mij = matrix.getEntry(i, j);
218                    final double mji = matrix.getEntry(j, i);
219                    if (Math.abs(mij - mji) > (Math.max(Math.abs(mij), Math.abs(mji)) * eps)) {
220                        return false;
221                    }
222                }
223            }
224            return true;
225        }
226    
227        /**
228         * Decompose a tridiagonal symmetric matrix. 
229         * @exception InvalidMatrixException (wrapping a {@link ConvergenceException}
230         * if algorithm fails to converge
231         */
232        private void decompose() {
233    
234            cachedV  = null;
235            cachedD  = null;
236            cachedVt = null;
237            work     = new double[6 * main.length];
238    
239            // compute the Gershgorin circles
240            computeGershgorinCircles();
241    
242            // find all the realEigenvalues
243            findEigenvalues();
244    
245            // we will search for eigenvectors only if required
246            eigenvectors = null;
247    
248        }
249    
250        /** {@inheritDoc} */
251        public RealMatrix getV()
252            throws InvalidMatrixException {
253    
254            if (cachedV == null) {
255    
256                if (eigenvectors == null) {
257                    findEigenVectors();
258                }
259    
260                final int m = eigenvectors.length;
261                cachedV = MatrixUtils.createRealMatrix(m, m);
262                for (int k = 0; k < m; ++k) {
263                    cachedV.setColumnVector(k, eigenvectors[k]);
264                }
265    
266            }
267    
268            // return the cached matrix
269            return cachedV;
270    
271        }
272    
273        /** {@inheritDoc} */
274        public RealMatrix getD()
275            throws InvalidMatrixException {
276            if (cachedD == null) {
277                // cache the matrix for subsequent calls
278                cachedD = MatrixUtils.createRealDiagonalMatrix(realEigenvalues);
279            }
280            return cachedD;
281        }
282    
283        /** {@inheritDoc} */
284        public RealMatrix getVT()
285            throws InvalidMatrixException {
286    
287            if (cachedVt == null) {
288    
289                if (eigenvectors == null) {
290                    findEigenVectors();
291                }
292    
293                final int m = eigenvectors.length;
294                cachedVt = MatrixUtils.createRealMatrix(m, m);
295                for (int k = 0; k < m; ++k) {
296                    cachedVt.setRowVector(k, eigenvectors[k]);
297                }
298    
299            }
300    
301            // return the cached matrix
302            return cachedVt;
303    
304        }
305    
306        /** {@inheritDoc} */
307        public double[] getRealEigenvalues()
308            throws InvalidMatrixException {
309            return realEigenvalues.clone();
310        }
311    
312        /** {@inheritDoc} */
313        public double getRealEigenvalue(final int i)
314            throws InvalidMatrixException, ArrayIndexOutOfBoundsException {
315            return realEigenvalues[i];
316        }
317    
318        /** {@inheritDoc} */
319        public double[] getImagEigenvalues()
320            throws InvalidMatrixException {
321            return imagEigenvalues.clone();
322        }
323    
324        /** {@inheritDoc} */
325        public double getImagEigenvalue(final int i)
326            throws InvalidMatrixException, ArrayIndexOutOfBoundsException {
327            return imagEigenvalues[i];
328        }
329    
330        /** {@inheritDoc} */
331        public RealVector getEigenvector(final int i)
332            throws InvalidMatrixException, ArrayIndexOutOfBoundsException {
333            if (eigenvectors == null) {
334                findEigenVectors();
335            }
336            return eigenvectors[i].copy();
337        }
338    
339        /**
340         * Return the determinant of the matrix
341         * @return determinant of the matrix
342         */
343        public double getDeterminant() {
344            double determinant = 1;
345            for (double lambda : realEigenvalues) {
346                determinant *= lambda;
347            }
348            return determinant;
349        }
350    
351        /** {@inheritDoc} */
352        public DecompositionSolver getSolver() {
353            if (eigenvectors == null) {
354                findEigenVectors();
355            }
356            return new Solver(realEigenvalues, imagEigenvalues, eigenvectors);
357        }
358    
359        /** Specialized solver. */
360        private static class Solver implements DecompositionSolver {
361        
362            /** Real part of the realEigenvalues. */
363            private double[] realEigenvalues;
364    
365            /** Imaginary part of the realEigenvalues. */
366            private double[] imagEigenvalues;
367    
368            /** Eigenvectors. */
369            private final ArrayRealVector[] eigenvectors;
370    
371            /**
372             * Build a solver from decomposed matrix.
373             * @param realEigenvalues real parts of the eigenvalues
374             * @param imagEigenvalues imaginary parts of the eigenvalues
375             * @param eigenvectors eigenvectors
376             */
377            private Solver(final double[] realEigenvalues, final double[] imagEigenvalues,
378                           final ArrayRealVector[] eigenvectors) {
379                this.realEigenvalues = realEigenvalues;
380                this.imagEigenvalues = imagEigenvalues;
381                this.eigenvectors    = eigenvectors; 
382            }
383    
384            /** Solve the linear equation A &times; X = B for symmetric matrices A.
385             * <p>This method only find exact linear solutions, i.e. solutions for
386             * which ||A &times; X - B|| is exactly 0.</p>
387             * @param b right-hand side of the equation A &times; X = B
388             * @return a vector X that minimizes the two norm of A &times; X - B
389             * @exception IllegalArgumentException if matrices dimensions don't match
390             * @exception InvalidMatrixException if decomposed matrix is singular
391             */
392            public double[] solve(final double[] b)
393                throws IllegalArgumentException, InvalidMatrixException {
394    
395                if (!isNonSingular()) {
396                    throw new SingularMatrixException();
397                }
398    
399                final int m = realEigenvalues.length;
400                if (b.length != m) {
401                    throw MathRuntimeException.createIllegalArgumentException(
402                            "vector length mismatch: got {0} but expected {1}",
403                            b.length, m);
404                }
405    
406                final double[] bp = new double[m];
407                for (int i = 0; i < m; ++i) {
408                    final ArrayRealVector v = eigenvectors[i];
409                    final double[] vData = v.getDataRef();
410                    final double s = v.dotProduct(b) / realEigenvalues[i];
411                    for (int j = 0; j < m; ++j) {
412                        bp[j] += s * vData[j];
413                    }
414                }
415    
416                return bp;
417    
418            }
419    
420            /** Solve the linear equation A &times; X = B for symmetric matrices A.
421             * <p>This method only find exact linear solutions, i.e. solutions for
422             * which ||A &times; X - B|| is exactly 0.</p>
423             * @param b right-hand side of the equation A &times; X = B
424             * @return a vector X that minimizes the two norm of A &times; X - B
425             * @exception IllegalArgumentException if matrices dimensions don't match
426             * @exception InvalidMatrixException if decomposed matrix is singular
427             */
428            public RealVector solve(final RealVector b)
429                throws IllegalArgumentException, InvalidMatrixException {
430    
431                if (!isNonSingular()) {
432                    throw new SingularMatrixException();
433                }
434    
435                final int m = realEigenvalues.length;
436                if (b.getDimension() != m) {
437                    throw MathRuntimeException.createIllegalArgumentException(
438                            "vector length mismatch: got {0} but expected {1}",
439                            b.getDimension(), m);
440                }
441    
442                final double[] bp = new double[m];
443                for (int i = 0; i < m; ++i) {
444                    final ArrayRealVector v = eigenvectors[i];
445                    final double[] vData = v.getDataRef();
446                    final double s = v.dotProduct(b) / realEigenvalues[i];
447                    for (int j = 0; j < m; ++j) {
448                        bp[j] += s * vData[j];
449                    }
450                }
451    
452                return new ArrayRealVector(bp, false);
453    
454            }
455    
456            /** Solve the linear equation A &times; X = B for symmetric matrices A.
457             * <p>This method only find exact linear solutions, i.e. solutions for
458             * which ||A &times; X - B|| is exactly 0.</p>
459             * @param b right-hand side of the equation A &times; X = B
460             * @return a matrix X that minimizes the two norm of A &times; X - B
461             * @exception IllegalArgumentException if matrices dimensions don't match
462             * @exception InvalidMatrixException if decomposed matrix is singular
463             */
464            public RealMatrix solve(final RealMatrix b)
465                throws IllegalArgumentException, InvalidMatrixException {
466    
467                if (!isNonSingular()) {
468                    throw new SingularMatrixException();
469                }
470    
471                final int m = realEigenvalues.length;
472                if (b.getRowDimension() != m) {
473                    throw MathRuntimeException.createIllegalArgumentException(
474                            "dimensions mismatch: got {0}x{1} but expected {2}x{3}",
475                            b.getRowDimension(), b.getColumnDimension(), m, "n");
476                }
477    
478                final int nColB = b.getColumnDimension();
479                final double[][] bp = new double[m][nColB];
480                for (int k = 0; k < nColB; ++k) {
481                    for (int i = 0; i < m; ++i) {
482                        final ArrayRealVector v = eigenvectors[i];
483                        final double[] vData = v.getDataRef();
484                        double s = 0;
485                        for (int j = 0; j < m; ++j) {
486                            s += v.getEntry(j) * b.getEntry(j, k);
487                        }
488                        s /= realEigenvalues[i];
489                        for (int j = 0; j < m; ++j) {
490                            bp[j][k] += s * vData[j];
491                        }
492                    }
493                }
494    
495                return MatrixUtils.createRealMatrix(bp);
496    
497            }
498    
499            /**
500             * Check if the decomposed matrix is non-singular.
501             * @return true if the decomposed matrix is non-singular
502             */
503            public boolean isNonSingular() {
504                for (int i = 0; i < realEigenvalues.length; ++i) {
505                    if ((realEigenvalues[i] == 0) && (imagEigenvalues[i] == 0)) {
506                        return false;
507                    }
508                }
509                return true;
510            }
511    
512            /** Get the inverse of the decomposed matrix.
513             * @return inverse matrix
514             * @throws InvalidMatrixException if decomposed matrix is singular
515             */
516            public RealMatrix getInverse()
517                throws InvalidMatrixException {
518    
519                if (!isNonSingular()) {
520                    throw new SingularMatrixException();
521                }
522    
523                final int m = realEigenvalues.length;
524                final double[][] invData = new double[m][m];
525    
526                for (int i = 0; i < m; ++i) {
527                    final double[] invI = invData[i];
528                    for (int j = 0; j < m; ++j) {
529                        double invIJ = 0;
530                        for (int k = 0; k < m; ++k) {
531                            final double[] vK = eigenvectors[k].getDataRef();
532                            invIJ += vK[i] * vK[j] / realEigenvalues[k];
533                        }
534                        invI[j] = invIJ;
535                    }
536                }
537                return MatrixUtils.createRealMatrix(invData);
538    
539            }
540    
541        }
542    
543        /**
544         * Transform matrix to tridiagonal.
545         * @param matrix matrix to transform
546         */
547        private void transformToTridiagonal(final RealMatrix matrix) {
548    
549            // transform the matrix to tridiagonal
550            transformer = new TriDiagonalTransformer(matrix);
551            main      = transformer.getMainDiagonalRef();
552            secondary = transformer.getSecondaryDiagonalRef();
553    
554            // pre-compute some elements
555            squaredSecondary = new double[secondary.length];
556            for (int i = 0; i < squaredSecondary.length; ++i) {
557                final double s = secondary[i];
558                squaredSecondary[i] = s * s;
559            }
560    
561        }
562    
563        /**
564         * Compute the Gershgorin circles for all rows.
565         */
566        private void computeGershgorinCircles() {
567    
568            final int m     = main.length;
569            final int lowerStart = 4 * m;
570            final int upperStart = 5 * m;
571            lowerSpectra = Double.POSITIVE_INFINITY;
572            upperSpectra = Double.NEGATIVE_INFINITY;
573            double eMax = 0;
574    
575            double eCurrent = 0;
576            for (int i = 0; i < m - 1; ++i) {
577    
578                final double dCurrent = main[i];
579                final double ePrevious = eCurrent;
580                eCurrent = Math.abs(secondary[i]);
581                eMax = Math.max(eMax, eCurrent);
582                final double radius = ePrevious + eCurrent;
583    
584                final double lower = dCurrent - radius;
585                work[lowerStart + i] = lower;
586                lowerSpectra = Math.min(lowerSpectra, lower);
587    
588                final double upper = dCurrent + radius;
589                work[upperStart + i] = upper;
590                upperSpectra = Math.max(upperSpectra, upper);
591                
592            }
593    
594            final double dCurrent = main[m - 1];
595            work[lowerStart + m - 1] = dCurrent - eCurrent;
596            work[upperStart + m - 1] = dCurrent + eCurrent;
597            minPivot = MathUtils.SAFE_MIN * Math.max(1.0, eMax * eMax);
598    
599        }
600    
601        /**
602         * Find the realEigenvalues.
603         * @exception InvalidMatrixException if a block cannot be diagonalized
604         */
605        private void findEigenvalues()
606            throws InvalidMatrixException {
607    
608            // compute splitting points
609            List<Integer> splitIndices = computeSplits();
610    
611            // find realEigenvalues in each block
612            realEigenvalues = new double[main.length];
613            imagEigenvalues = new double[main.length];
614            int begin = 0;
615            for (final int end : splitIndices) {
616                final int n = end - begin;
617                switch (n) {
618    
619                case 1:
620                    // apply dedicated method for dimension 1
621                    process1RowBlock(begin);
622                    break;
623    
624                case 2:
625                    // apply dedicated method for dimension 2
626                    process2RowsBlock(begin);
627                    break;
628    
629                case 3:
630                    // apply dedicated method for dimension 3
631                    process3RowsBlock(begin);
632                    break;
633    
634                default:
635    
636                    // choose an initial shift for LDL<sup>T</sup> decomposition
637                    final double[] range       = eigenvaluesRange(begin, n);
638                    final double oneFourth     = 0.25 * (3 * range[0] + range[1]);
639                    final int oneFourthCount   = countEigenValues(oneFourth, begin, n);
640                    final double threeFourth   = 0.25 * (range[0] + 3 * range[1]);
641                    final int threeFourthCount = countEigenValues(threeFourth, begin, n);
642                    final boolean chooseLeft   = (oneFourthCount - 1) >= (n - threeFourthCount);
643                    final double lambda        = chooseLeft ? range[0] : range[1];
644    
645                    tau = (range[1] - range[0]) * MathUtils.EPSILON * n + 2 * minPivot;
646    
647                    // decompose T&lambda;I as LDL<sup>T</sup>
648                    ldlTDecomposition(lambda, begin, n);
649    
650                    // apply general dqd/dqds method
651                    processGeneralBlock(n);
652    
653                    // extract realEigenvalues
654                    if (chooseLeft) {
655                        for (int i = 0; i < n; ++i) {
656                            realEigenvalues[begin + i] = lambda + work[4 * i];
657                        }
658                    } else {
659                        for (int i = 0; i < n; ++i) {
660                            realEigenvalues[begin + i] = lambda - work[4 * i];
661                        }                    
662                    }
663    
664                }
665                begin = end;
666            }
667    
668            // sort the realEigenvalues in decreasing order
669            Arrays.sort(realEigenvalues);
670            for (int i = 0, j = realEigenvalues.length - 1; i < j; ++i, --j) {
671                final double tmp = realEigenvalues[i];
672                realEigenvalues[i] = realEigenvalues[j];
673                realEigenvalues[j] = tmp;
674            }
675    
676        }
677    
678        /**
679         * Compute splitting points.
680         * @return list of indices after matrix can be split
681         */
682        private List<Integer> computeSplits() {
683    
684            final List<Integer> list = new ArrayList<Integer>();
685    
686            // splitting preserving relative accuracy
687            double absDCurrent = Math.abs(main[0]);
688            for (int i = 0; i < secondary.length; ++i) {
689                final double absDPrevious = absDCurrent;
690                absDCurrent = Math.abs(main[i + 1]);
691                final double max = splitTolerance * Math.sqrt(absDPrevious * absDCurrent);
692                if (Math.abs(secondary[i]) <= max) {
693                    list.add(i + 1);
694                    secondary[i] = 0;
695                    squaredSecondary[i] = 0;
696                }
697            }
698    
699            list.add(secondary.length + 1);
700            return list;
701    
702        }
703    
704        /**
705         * Find eigenvalue in a block with 1 row.
706         * <p>In low dimensions, we simply solve the characteristic polynomial.</p>
707         * @param index index of the first row of the block
708         */
709        private void process1RowBlock(final int index) {
710            realEigenvalues[index] = main[index];
711        }
712    
713        /**
714         * Find realEigenvalues in a block with 2 rows.
715         * <p>In low dimensions, we simply solve the characteristic polynomial.</p>
716         * @param index index of the first row of the block
717         * @exception InvalidMatrixException if characteristic polynomial cannot be solved
718         */
719        private void process2RowsBlock(final int index)
720            throws InvalidMatrixException {
721    
722            // the characteristic polynomial is
723            // X^2 - (q0 + q1) X + q0 q1 - e1^2
724            final double q0   = main[index];
725            final double q1   = main[index + 1];
726            final double e12  = squaredSecondary[index];
727    
728            final double s     = q0 + q1;
729            final double p     = q0 * q1 - e12;
730            final double delta = s * s - 4 * p;
731            if (delta < 0) {
732                throw new InvalidMatrixException("cannot solve degree {0} equation", 2);
733            }
734    
735            final double largestRoot = 0.5 * (s + Math.sqrt(delta));
736            realEigenvalues[index]     = largestRoot;
737            realEigenvalues[index + 1] = p / largestRoot;
738    
739        }
740    
741        /**
742         * Find realEigenvalues in a block with 3 rows.
743         * <p>In low dimensions, we simply solve the characteristic polynomial.</p>
744         * @param index index of the first row of the block
745         * @exception InvalidMatrixException if diagonal elements are not positive
746         */
747        private void process3RowsBlock(final int index)
748            throws InvalidMatrixException {
749    
750            // the characteristic polynomial is
751            // X^3 - (q0 + q1 + q2) X^2 + (q0 q1 + q0 q2 + q1 q2 - e1^2 - e2^2) X + q0 e2^2 + q2 e1^2 - q0 q1 q2
752            final double q0       = main[index];
753            final double q1       = main[index + 1];
754            final double q2       = main[index + 2];
755            final double e12      = squaredSecondary[index];
756            final double q1q2Me22 = q1 * q2 - squaredSecondary[index + 1];
757    
758            // compute coefficients of the cubic equation as: x^3 + b x^2 + c x + d = 0
759            final double b        = -(q0 + q1 + q2);
760            final double c        = q0 * q1 + q0 * q2 + q1q2Me22 - e12;
761            final double d        = q2 * e12 - q0 * q1q2Me22;
762    
763            // solve cubic equation
764            final double b2       = b * b;
765            final double q        = (3 * c - b2) / 9;
766            final double r        = ((9 * c - 2 * b2) * b - 27 * d) / 54;
767            final double delta    = q * q * q + r * r;
768            if (delta >= 0) {
769                // in fact, there are solutions to the equation, but in the context
770                // of symmetric realEigenvalues problem, there should be three distinct
771                // real roots, so we throw an error if this condition is not met
772                throw new InvalidMatrixException("cannot solve degree {0} equation", 3);           
773            }
774            final double sqrtMq = Math.sqrt(-q);
775            final double theta  = Math.acos(r / (-q * sqrtMq));
776            final double alpha  = 2 * sqrtMq;
777            final double beta   = b / 3;
778    
779            double z0 = alpha * Math.cos(theta / 3) - beta;
780            double z1 = alpha * Math.cos((theta + 2 * Math.PI) / 3) - beta;
781            double z2 = alpha * Math.cos((theta + 4 * Math.PI) / 3) - beta;
782            if (z0 < z1) {
783                final double t = z0;
784                z0 = z1;
785                z1 = t;
786            }
787            if (z1 < z2) {
788                final double t = z1;
789                z1 = z2;
790                z2 = t;
791            }
792            if (z0 < z1) {
793                final double t = z0;
794                z0 = z1;
795                z1 = t;
796            }
797            realEigenvalues[index]     = z0;
798            realEigenvalues[index + 1] = z1;
799            realEigenvalues[index + 2] = z2;
800    
801        }
802    
803        /**
804         * Find realEigenvalues using dqd/dqds algorithms.
805         * <p>This implementation is based on Beresford N. Parlett
806         * and Osni A. Marques paper <a
807         * href="http://www.netlib.org/lapack/lawnspdf/lawn155.pdf">An
808         * Implementation of the dqds Algorithm (Positive Case)</a> and on the
809         * corresponding LAPACK routine DLASQ2.</p>
810         * @param n number of rows of the block
811         * @exception InvalidMatrixException if block cannot be diagonalized
812         * after 30 * n iterations
813         */
814        private void processGeneralBlock(final int n)
815            throws InvalidMatrixException {
816    
817            // check decomposed matrix data range
818            double sumOffDiag = 0;
819            for (int i = 0; i < n - 1; ++i) {
820                final int fourI = 4 * i;
821                final double ei = work[fourI + 2];
822                sumOffDiag += ei;
823            }
824    
825            if (sumOffDiag == 0) {
826                // matrix is already diagonal
827                return;
828            }
829    
830            // initial checks for splits (see Parlett & Marques section 3.3)
831            flipIfWarranted(n, 2);
832    
833            // two iterations with Li's test for initial splits
834            initialSplits(n);
835    
836            // initialize parameters used by goodStep
837            tType = 0;
838            dMin1 = 0;
839            dMin2 = 0;
840            dN    = 0;
841            dN1   = 0;
842            dN2   = 0;
843            tau   = 0;
844    
845            // process split segments
846            int i0 = 0;
847            int n0 = n;
848            while (n0 > 0) {
849    
850                // retrieve shift that was temporarily stored as a negative off-diagonal element
851                sigma    = (n0 == n) ? 0 : -work[4 * n0 - 2];
852                sigmaLow = 0;
853    
854                // find start of a new split segment to process
855                double eMin = (i0 == n0) ? 0 : work[4 * n0 - 6];
856                double eMax = 0;
857                double qMax = work[4 * n0 - 4];
858                double qMin = qMax;
859                i0 = 0;
860                for (int i = 4 * (n0 - 2); i >= 0; i -= 4) {
861                    if (work[i + 2] <= 0) {
862                        i0 = 1 + i / 4;
863                        break;
864                    }
865                    if (qMin >= 4 * eMax) {
866                        qMin = Math.min(qMin, work[i + 4]);
867                        eMax = Math.max(eMax, work[i + 2]);
868                    }
869                    qMax = Math.max(qMax, work[i] + work[i + 2]);
870                    eMin = Math.min(eMin, work[i + 2]);
871                }
872                work[4 * n0 - 2] = eMin;
873    
874                // lower bound of Gershgorin disk
875                dMin = -Math.max(0, qMin - 2 * Math.sqrt(qMin * eMax));
876    
877                pingPong = 0;
878                int maxIter = 30 * (n0 - i0);
879                for (int k = 0; i0 < n0; ++k) {
880                    if (k >= maxIter) {
881                        throw new InvalidMatrixException(new MaxIterationsExceededException(maxIter));
882                    }
883    
884                    // perform one step
885                    n0 = goodStep(i0, n0);
886                    pingPong = 1 - pingPong;
887    
888                    // check for new splits after "ping" steps
889                    // when the last elements of qd array are very small
890                    if ((pingPong == 0) && (n0 - i0 > 3) &&
891                        (work[4 * n0 - 1] <= TOLERANCE_2 * qMax) &&
892                        (work[4 * n0 - 2] <= TOLERANCE_2 * sigma)) {
893                        int split = i0 - 1;
894                        qMax = work[4 * i0];
895                        eMin = work[4 * i0 + 2];
896                        double previousEMin = work[4 * i0 + 3];
897                        for (int i = 4 * i0; i < 4 * n0 - 11; i += 4) {
898                            if ((work[i + 3] <= TOLERANCE_2 * work[i]) &&
899                                (work[i + 2] <= TOLERANCE_2 * sigma)) {
900                                // insert a split
901                                work[i + 2]  = -sigma;
902                                split        = i / 4;
903                                qMax         = 0;
904                                eMin         = work[i + 6];
905                                previousEMin = work[i + 7];
906                            } else {
907                                qMax         = Math.max(qMax, work[i + 4]);
908                                eMin         = Math.min(eMin, work[i + 2]);
909                                previousEMin = Math.min(previousEMin, work[i + 3]);
910                            }
911                        }
912                        work[4 * n0 - 2] = eMin;
913                        work[4 * n0 - 1] = previousEMin;
914                        i0 = split + 1;
915                    }
916                }
917    
918            }
919    
920        }
921    
922        /**
923         * Perform two iterations with Li's tests for initial splits.
924         * @param n number of rows of the matrix to process
925         */
926        private void initialSplits(final int n) {
927    
928            pingPong = 0;
929            for (int k = 0; k < 2; ++k) {
930    
931                // apply Li's reverse test
932                double d = work[4 * (n - 1) + pingPong];
933                for (int i = 4 * (n - 2) + pingPong; i >= 0; i -= 4) {
934                    if (work[i + 2] <= TOLERANCE_2 * d) {
935                        work[i + 2] = -0.0;
936                        d = work[i];
937                    } else {
938                        d *= work[i] / (d + work[i + 2]);
939                    }
940                }
941    
942                // apply dqd plus Li's forward test.
943                d = work[pingPong];
944                for (int i = 2 + pingPong; i < 4 * n - 2; i += 4) {
945                    final int j = i - 2 * pingPong - 1;
946                    work[j] = d + work[i];
947                    if (work[i] <= TOLERANCE_2 * d) {
948                        work[i]     = -0.0;
949                        work[j]     = d;
950                        work[j + 2] = 0.0;
951                        d = work[i + 2];
952                    } else if ((MathUtils.SAFE_MIN * work[i + 2] < work[j]) &&
953                               (MathUtils.SAFE_MIN * work[j] < work[i + 2])) {
954                        final double tmp = work[i + 2] / work[j];
955                        work[j + 2] = work[i] * tmp;
956                        d *= tmp;
957                    } else {
958                        work[j + 2] = work[i + 2] * (work[i] / work[j]);
959                        d *= work[i + 2] / work[j];
960                   }
961                }
962                work[4 * n - 3 - pingPong] = d;
963    
964                // from ping to pong
965                pingPong = 1 - pingPong;
966    
967            }
968    
969        }
970    
971        /**
972         * Perform one "good" dqd/dqds step.
973         * <p>This implementation is based on Beresford N. Parlett
974         * and Osni A. Marques paper <a
975         * href="http://www.netlib.org/lapack/lawnspdf/lawn155.pdf">An
976         * Implementation of the dqds Algorithm (Positive Case)</a> and on the
977         * corresponding LAPACK routine DLAZQ3.</p>
978         * @param start start index
979         * @param end end index
980         * @return new end (maybe deflated)
981         */
982        private int goodStep(final int start, final int end) {
983    
984            g = 0.0;
985    
986            // step 1: accepting realEigenvalues
987            int deflatedEnd = end;
988            for (boolean deflating = true; deflating;) {
989    
990                if (start >= deflatedEnd) {
991                    // the array has been completely deflated
992                    return deflatedEnd;
993                }
994    
995                final int k = 4 * deflatedEnd + pingPong - 1;
996    
997                if ((start == deflatedEnd - 1) ||
998                    ((start != deflatedEnd - 2) &&
999                     ((work[k - 5] <= TOLERANCE_2 * (sigma + work[k - 3])) ||
1000                      (work[k - 2 * pingPong - 4] <= TOLERANCE_2 * work[k - 7])))) {
1001    
1002                    // one eigenvalue found, deflate array
1003                    work[4 * deflatedEnd - 4] = sigma + work[4 * deflatedEnd - 4 + pingPong];
1004                    deflatedEnd -= 1;
1005    
1006                } else if ((start == deflatedEnd - 2) ||
1007                    (work[k - 9] <= TOLERANCE_2 * sigma) ||
1008                    (work[k - 2 * pingPong - 8] <= TOLERANCE_2 * work[k - 11])) {
1009    
1010                    // two realEigenvalues found, deflate array
1011                    if (work[k - 3] > work[k - 7]) {
1012                        final double tmp = work[k - 3];
1013                        work[k - 3] = work[k - 7];
1014                        work[k - 7] = tmp;
1015                    }
1016    
1017                    if (work[k - 5] > TOLERANCE_2 * work[k - 3]) {
1018                        double t = 0.5 * ((work[k - 7] - work[k - 3]) + work[k - 5]);
1019                        double s = work[k - 3] * (work[k - 5] / t);
1020                        if (s <= t) {
1021                            s = work[k - 3] * work[k - 5] / (t * (1 + Math.sqrt(1 + s / t)));
1022                        } else {
1023                            s = work[k - 3] * work[k - 5] / (t + Math.sqrt(t * (t + s)));                      
1024                        }
1025                        t = work[k - 7] + (s + work[k - 5]);
1026                        work[k - 3] *= work[k - 7] / t;
1027                        work[k - 7]  = t;
1028                    }
1029                    work[4 * deflatedEnd - 8] = sigma + work[k - 7];
1030                    work[4 * deflatedEnd - 4] = sigma + work[k - 3];
1031                    deflatedEnd -= 2;
1032                } else {
1033    
1034                    // no more realEigenvalues found, we need to iterate
1035                    deflating = false;
1036    
1037                }
1038    
1039            }
1040    
1041            final int l = 4 * deflatedEnd + pingPong - 1;
1042    
1043            // step 2: flip array if needed
1044            if ((dMin <= 0) || (deflatedEnd < end)) {
1045                if (flipIfWarranted(deflatedEnd, 1)) {
1046                    dMin2 = Math.min(dMin2, work[l - 1]);
1047                    work[l - 1] =
1048                        Math.min(work[l - 1],
1049                                 Math.min(work[3 + pingPong], work[7 + pingPong]));
1050                    work[l - 2 * pingPong] =
1051                        Math.min(work[l - 2 * pingPong],
1052                                 Math.min(work[6 + pingPong], work[6 + pingPong]));
1053                    qMax  = Math.max(qMax, Math.max(work[3 + pingPong], work[7 + pingPong]));
1054                    dMin  = -0.0;
1055                }
1056            }
1057    
1058            if ((dMin < 0) ||
1059                (MathUtils.SAFE_MIN * qMax < Math.min(work[l - 1],
1060                                                      Math.min(work[l - 9],
1061                                                               dMin2 + work[l - 2 * pingPong])))) {
1062                // step 3: choose a shift
1063                computeShiftIncrement(start, deflatedEnd, end - deflatedEnd);
1064    
1065                // step 4a: dqds
1066                for (boolean loop = true; loop;) {
1067    
1068                    // perform one dqds step with the chosen shift
1069                    dqds(start, deflatedEnd);
1070    
1071                    // check result of the dqds step
1072                    if ((dMin >= 0) && (dMin1 > 0)) {
1073                        // the shift was good
1074                        updateSigma(tau);
1075                        return deflatedEnd;
1076                    } else if ((dMin < 0.0) &&
1077                               (dMin1 > 0.0) &&
1078                               (work[4 * deflatedEnd - 5 - pingPong] < TOLERANCE * (sigma + dN1)) &&
1079                               (Math.abs(dN) < TOLERANCE * sigma)) {
1080                       // convergence hidden by negative DN.
1081                        work[4 * deflatedEnd - 3 - pingPong] = 0.0;
1082                        dMin = 0.0;
1083                        updateSigma(tau);
1084                        return deflatedEnd;
1085                    } else if (dMin < 0.0) {
1086                        // tau too big. Select new tau and try again.
1087                        if (tType < -22) {
1088                            // failed twice. Play it safe.
1089                            tau = 0.0;
1090                        } else if (dMin1 > 0.0) {
1091                            // late failure. Gives excellent shift.
1092                            tau = (tau + dMin) * (1.0 - 2.0 * MathUtils.EPSILON);
1093                            tType -= 11;
1094                        } else {
1095                            // early failure. Divide by 4.
1096                            tau *= 0.25;
1097                            tType -= 12;
1098                        }
1099                    } else if (Double.isNaN(dMin)) {
1100                        tau = 0.0;
1101                    } else {
1102                        // possible underflow. Play it safe.
1103                        loop = false;
1104                    }
1105                }
1106    
1107            }
1108    
1109            // perform a dqd step (i.e. no shift)
1110            dqd(start, deflatedEnd);
1111    
1112            return deflatedEnd;
1113    
1114        }
1115    
1116        /**
1117         * Flip qd array if warranted.
1118         * @param n number of rows in the block
1119         * @param step within the array (1 for flipping all elements, 2 for flipping
1120         * only every other element)
1121         * @return true if qd array was flipped
1122         */
1123        private boolean flipIfWarranted(final int n, final int step) {
1124            if (1.5 * work[pingPong] < work[4 * (n - 1) + pingPong]) {
1125                // flip array
1126                for (int i = 0, j = 4 * n - 1; i < j; i += 4, j -= 4) {
1127                    for (int k = 0; k < 4; k += step) {
1128                        final double tmp = work[i + k];
1129                        work[i + k] = work[j - k];
1130                        work[j - k] = tmp;
1131                    }
1132                }
1133                return true;
1134            }
1135            return false;
1136        }
1137    
1138        /**
1139         * Compute an interval containing all realEigenvalues of a block.
1140         * @param index index of the first row of the block
1141         * @param n number of rows of the block
1142         * @return an interval containing the realEigenvalues
1143         */
1144        private double[] eigenvaluesRange(final int index, final int n) {
1145    
1146            // find the bounds of the spectra of the local block
1147            final int lowerStart = 4 * main.length;
1148            final int upperStart = 5 * main.length;
1149            double lower = Double.POSITIVE_INFINITY;
1150            double upper = Double.NEGATIVE_INFINITY;
1151            for (int i = 0; i < n; ++i) {
1152                lower = Math.min(lower, work[lowerStart + index +i]);
1153                upper = Math.max(upper, work[upperStart + index +i]);
1154            }
1155    
1156            // set thresholds
1157            final double tNorm = Math.max(Math.abs(lower), Math.abs(upper));
1158            final double relativeTolerance = Math.sqrt(MathUtils.EPSILON);
1159            final double absoluteTolerance = 4 * minPivot;
1160            final int maxIter =
1161                2 + (int) ((Math.log(tNorm + minPivot) - Math.log(minPivot)) / Math.log(2.0));
1162            final double margin = 2 * (tNorm * MathUtils.EPSILON * n + 2 * minPivot);
1163    
1164            // search lower eigenvalue
1165            double left  = lower - margin;
1166            double right = upper + margin;
1167            for (int i = 0; i < maxIter; ++i) {
1168    
1169                final double range = right - left;
1170                if ((range < absoluteTolerance) ||
1171                    (range < relativeTolerance * Math.max(Math.abs(left), Math.abs(right)))) {
1172                    // search has converged
1173                    break;
1174                }
1175    
1176                final double middle = 0.5 * (left + right);
1177                if (countEigenValues(middle, index, n) >= 1) {
1178                    right = middle;
1179                } else {
1180                    left = middle;
1181                }
1182    
1183            }
1184            lower = Math.max(lower, left - 100 * MathUtils.EPSILON * Math.abs(left));
1185    
1186            // search upper eigenvalue
1187            left  = lower - margin;
1188            right = upper + margin;
1189            for (int i = 0; i < maxIter; ++i) {
1190    
1191                final double range = right - left;
1192                if ((range < absoluteTolerance) ||
1193                    (range < relativeTolerance * Math.max(Math.abs(left), Math.abs(right)))) {
1194                    // search has converged
1195                    break;
1196                }
1197    
1198                final double middle = 0.5 * (left + right);
1199                if (countEigenValues(middle, index, n) >= n) {
1200                    right = middle;
1201                } else {
1202                    left = middle;
1203                }
1204    
1205            }
1206            upper = Math.min(upper, right + 100 * MathUtils.EPSILON * Math.abs(right));
1207    
1208            return new double[] { lower, upper };
1209    
1210        }
1211    
1212        /**
1213         * Count the number of realEigenvalues below a point.
1214         * @param t value below which we must count the number of realEigenvalues
1215         * @param index index of the first row of the block
1216         * @param n number of rows of the block
1217         * @return number of realEigenvalues smaller than t
1218         */
1219        private int countEigenValues(final double t, final int index, final int n) {
1220            double ratio = main[index] - t;
1221            int count = (ratio > 0) ? 0 : 1;
1222            for (int i = 1; i < n; ++i) {
1223                ratio = main[index + i] - squaredSecondary[index + i - 1] / ratio - t;
1224                if (ratio <= 0) {
1225                    ++count;
1226                }
1227            }
1228            return count;
1229        }
1230    
1231        /**
1232         * Decompose the shifted tridiagonal matrix T-&lambda;I as LDL<sup>T</sup>.
1233         * <p>A shifted symmetric tridiagonal matrix T can be decomposed as
1234         * LDL<sup>T</sup> where L is a lower bidiagonal matrix with unit diagonal
1235         * and D is a diagonal matrix. This method is an implementation of
1236         * algorithm 4.4.7 from Dhillon's thesis.</p>
1237         * @param lambda shift to add to the matrix before decomposing it
1238         * to ensure it is positive definite
1239         * @param index index of the first row of the block
1240         * @param n number of rows of the block
1241         */
1242        private void ldlTDecomposition(final double lambda, final int index, final int n) {
1243            double di = main[index] - lambda;
1244            work[0] = Math.abs(di);
1245            for (int i = 1; i < n; ++i) {
1246                final int    fourI = 4 * i;
1247                final double eiM1  = secondary[index + i - 1];
1248                final double ratio = eiM1 / di;
1249                work[fourI - 2] = ratio * ratio * Math.abs(di);
1250                di = (main[index + i] - lambda) - eiM1 * ratio;
1251                work[fourI] = Math.abs(di);
1252            }
1253        }
1254    
1255        /**
1256         * Perform a dqds step, using current shift increment.
1257         * <p>This implementation is a translation of the LAPACK routine DLASQ5.</p>
1258         * @param start start index
1259         * @param end end index
1260         */
1261        private void dqds(final int start, final int end) {
1262    
1263            eMin = work[4 * start + pingPong + 4];
1264            double d = work[4 * start + pingPong] - tau;
1265            dMin = d;
1266            dMin1 = -work[4 * start + pingPong];
1267    
1268            if (pingPong == 0) {
1269                for (int j4 = 4 * start + 3; j4 <= 4 * (end - 3); j4 += 4) {
1270                    work[j4 - 2] = d + work[j4 - 1];
1271                    final double tmp = work[j4 + 1] / work[j4 - 2];
1272                    d = d * tmp - tau;
1273                    dMin = Math.min(dMin, d);
1274                    work[j4] = work[j4 - 1] * tmp;
1275                    eMin = Math.min(work[j4], eMin);
1276                }
1277            } else {
1278                for (int j4 = 4 * start + 3; j4 <= 4 * (end - 3); j4 += 4) {
1279                    work[j4 - 3] = d + work[j4];
1280                    final double tmp = work[j4 + 2] / work[j4 - 3];
1281                    d = d * tmp - tau;
1282                    dMin = Math.min(dMin, d);
1283                    work[j4 - 1] = work[j4] * tmp;
1284                    eMin = Math.min(work[j4 - 1], eMin);
1285                }
1286            }
1287    
1288            // unroll last two steps.
1289            dN2 = d;
1290            dMin2 = dMin;
1291            int j4 = 4 * (end - 2) - pingPong - 1;
1292            int j4p2 = j4 + 2 * pingPong - 1;
1293            work[j4 - 2] = dN2 + work[j4p2];
1294            work[j4] = work[j4p2 + 2] * (work[j4p2] / work[j4 - 2]);
1295            dN1 = work[j4p2 + 2] * (dN2 / work[j4 - 2]) - tau;
1296            dMin = Math.min(dMin, dN1);
1297    
1298            dMin1 = dMin;
1299            j4 = j4 + 4;
1300            j4p2 = j4 + 2 * pingPong - 1;
1301            work[j4 - 2] = dN1 + work[j4p2];
1302            work[j4] = work[j4p2 + 2] * (work[j4p2] / work[j4 - 2]);
1303            dN = work[j4p2 + 2] * (dN1 / work[j4 - 2]) - tau;
1304            dMin = Math.min(dMin, dN);
1305    
1306            work[j4 + 2] = dN;
1307            work[4 * end - pingPong - 1] = eMin;
1308    
1309        }
1310    
1311    
1312        /**
1313         * Perform a dqd step.
1314         * <p>This implementation is a translation of the LAPACK routine DLASQ6.</p>
1315         * @param start start index
1316         * @param end end index
1317         */
1318        private void dqd(final int start, final int end) {
1319    
1320            eMin = work[4 * start + pingPong + 4];
1321            double d = work[4 * start + pingPong];
1322            dMin = d;
1323    
1324            if (pingPong == 0) {
1325                for (int j4 = 4 * start + 3; j4 < 4 * (end - 3); j4 += 4) {
1326                    work[j4 - 2] = d + work[j4 - 1];
1327                    if (work[j4 - 2] == 0.0) {
1328                        work[j4] = 0.0;
1329                        d = work[j4 + 1];
1330                        dMin = d;
1331                        eMin = 0.0;
1332                    } else if ((MathUtils.SAFE_MIN * work[j4 + 1] < work[j4 - 2]) &&
1333                               (MathUtils.SAFE_MIN * work[j4 - 2] < work[j4 + 1])) {
1334                        final double tmp = work[j4 + 1] / work[j4 - 2];
1335                        work[j4] = work[j4 - 1] * tmp;
1336                        d *= tmp;
1337                    } else {
1338                        work[j4] = work[j4 + 1] * (work[j4 - 1] / work[j4 - 2]);
1339                        d *= work[j4 + 1] / work[j4 - 2];
1340                    }
1341                    dMin = Math.min(dMin, d);
1342                    eMin = Math.min(eMin, work[j4]);
1343                }
1344            } else {
1345                for (int j4 = 4 * start + 3; j4 < 4 * (end - 3); j4 += 4) {
1346                    work[j4 - 3] = d + work[j4];
1347                    if (work[j4 - 3] == 0.0) {
1348                        work[j4 - 1] = 0.0;
1349                        d = work[j4 + 2];
1350                        dMin = d;
1351                        eMin = 0.0;
1352                    } else if ((MathUtils.SAFE_MIN * work[j4 + 2] < work[j4 - 3]) &&
1353                               (MathUtils.SAFE_MIN * work[j4 - 3] < work[j4 + 2])) {
1354                        final double tmp = work[j4 + 2] / work[j4 - 3];
1355                        work[j4 - 1] = work[j4] * tmp;
1356                        d *= tmp;
1357                    } else {
1358                        work[j4 - 1] = work[j4 + 2] * (work[j4] / work[j4 - 3]);
1359                        d *= work[j4 + 2] / work[j4 - 3];
1360                    }
1361                    dMin = Math.min(dMin, d);
1362                    eMin = Math.min(eMin, work[j4 - 1]);
1363                }
1364            }
1365    
1366            // Unroll last two steps
1367            dN2   = d;
1368            dMin2 = dMin;
1369            int j4 = 4 * (end - 2) - pingPong - 1;
1370            int j4p2 = j4 + 2 * pingPong - 1;
1371            work[j4 - 2] = dN2 + work[j4p2];
1372            if (work[j4 - 2] == 0.0) {
1373                work[j4] = 0.0;
1374                dN1  = work[j4p2 + 2];
1375                dMin = dN1;
1376                eMin = 0.0;
1377            } else if ((MathUtils.SAFE_MIN * work[j4p2 + 2] < work[j4 - 2]) &&
1378                       (MathUtils.SAFE_MIN * work[j4 - 2] < work[j4p2 + 2])) {
1379                final double tmp = work[j4p2 + 2] / work[j4 - 2];
1380                work[j4] = work[j4p2] * tmp;
1381                dN1 = dN2 * tmp;
1382            } else {
1383                work[j4] = work[j4p2 + 2] * (work[j4p2] / work[j4 - 2]);
1384                dN1 = work[j4p2 + 2] * (dN2 / work[j4 - 2]);
1385            }
1386            dMin = Math.min(dMin, dN1);
1387    
1388            dMin1 = dMin;
1389            j4 = j4 + 4;
1390            j4p2 = j4 + 2 * pingPong - 1;
1391            work[j4 - 2] = dN1 + work[j4p2];
1392            if (work[j4 - 2] == 0.0) {
1393                work[j4] = 0.0;
1394                dN   = work[j4p2 + 2];
1395                dMin = dN;
1396                eMin = 0.0;
1397            } else if ((MathUtils.SAFE_MIN * work[j4p2 + 2] < work[j4 - 2]) &&
1398                       (MathUtils.SAFE_MIN * work[j4 - 2] < work[j4p2 + 2])) {
1399                final double tmp = work[j4p2 + 2] / work[j4 - 2];
1400                work[j4] = work[j4p2] * tmp;
1401                dN = dN1 * tmp;
1402            } else {
1403                work[j4] = work[j4p2 + 2] * (work[j4p2] / work[j4 - 2]);
1404                dN = work[j4p2 + 2] * (dN1 / work[j4 - 2]);
1405            }
1406            dMin = Math.min(dMin, dN);
1407    
1408            work[j4 + 2] = dN;
1409            work[4 * end - pingPong - 1] = eMin;
1410    
1411        }
1412    
1413        /**
1414         * Compute the shift increment as an estimate of the smallest eigenvalue.
1415         * <p>This implementation is a translation of the LAPACK routine DLAZQ4.</p>
1416         * @param start start index
1417         * @param end end index
1418         * @param deflated number of realEigenvalues just deflated
1419         */
1420        private void computeShiftIncrement(final int start, final int end, final int deflated) {
1421    
1422            final double cnst1 = 0.563;
1423            final double cnst2 = 1.010;
1424            final double cnst3 = 1.05;
1425    
1426            // a negative dMin forces the shift to take that absolute value
1427            // tType records the type of shift.
1428            if (dMin <= 0.0) {
1429                tau = -dMin;
1430                tType = -1;
1431                return;
1432            }
1433    
1434            int nn = 4 * end + pingPong - 1;
1435            switch (deflated) {
1436    
1437            case 0 : // no realEigenvalues deflated. 
1438                if (dMin == dN || dMin == dN1) {
1439    
1440                    double b1 = Math.sqrt(work[nn - 3]) * Math.sqrt(work[nn - 5]);
1441                    double b2 = Math.sqrt(work[nn - 7]) * Math.sqrt(work[nn - 9]);
1442                    double a2 = work[nn - 7] + work[nn - 5];
1443    
1444                    if (dMin == dN && dMin1 == dN1) {
1445                        // cases 2 and 3. 
1446                        final double gap2 = dMin2 - a2 - dMin2 * 0.25;
1447                        final double gap1 = a2 - dN - ((gap2 > 0.0 && gap2 > b2) ? (b2 / gap2) * b2 : (b1 + b2));
1448                        if (gap1 > 0.0 && gap1 > b1) {
1449                            tau   = Math.max(dN - (b1 / gap1) * b1, 0.5 * dMin);
1450                            tType = -2;
1451                        } else {
1452                            double s = 0.0;
1453                            if (dN > b1) {
1454                                s = dN - b1;
1455                            }
1456                            if (a2 > (b1 + b2)) {
1457                                s = Math.min(s, a2 - (b1 + b2));
1458                            }
1459                            tau   = Math.max(s, 0.333 * dMin);
1460                            tType = -3;
1461                        }
1462                    } else {
1463                        // case 4.
1464                        tType = -4;
1465                        double s = 0.25 * dMin;
1466                        double gam;
1467                        int np;
1468                        if (dMin == dN) {
1469                            gam = dN;
1470                            a2 = 0.0;
1471                            if (work[nn - 5]  >  work[nn - 7]) {
1472                                return;
1473                            }
1474                            b2 = work[nn - 5] / work[nn - 7];
1475                            np = nn - 9;
1476                        } else {
1477                            np = nn - 2 * pingPong;
1478                            b2 = work[np - 2];
1479                            gam = dN1;
1480                            if (work[np - 4]  >  work[np - 2]) {
1481                                return;
1482                            }
1483                            a2 = work[np - 4] / work[np - 2];
1484                            if (work[nn - 9]  >  work[nn - 11]) {
1485                                return;
1486                            }
1487                            b2 = work[nn - 9] / work[nn - 11];
1488                            np = nn - 13;
1489                        }
1490    
1491                        // approximate contribution to norm squared from i < nn-1.
1492                        a2 = a2 + b2;
1493                        for (int i4 = np; i4 >= 4 * start + 2 + pingPong; i4 -= 4) {
1494                            if(b2 == 0.0) {
1495                                break;
1496                            }
1497                            b1 = b2;
1498                            if (work[i4]  >  work[i4 - 2]) {
1499                                return;
1500                            }
1501                            b2 = b2 * (work[i4] / work[i4 - 2]);
1502                            a2 = a2 + b2;
1503                            if (100 * Math.max(b2, b1) < a2 || cnst1 < a2) {
1504                                break;
1505                            }
1506                        }
1507                        a2 = cnst3 * a2;
1508    
1509                        // rayleigh quotient residual bound.
1510                        if (a2 < cnst1) {
1511                            s = gam * (1 - Math.sqrt(a2)) / (1 + a2);
1512                        }
1513                        tau = s;
1514    
1515                    }
1516                } else if (dMin == dN2) {
1517    
1518                    // case 5.
1519                    tType = -5;
1520                    double s = 0.25 * dMin;
1521    
1522                    // compute contribution to norm squared from i > nn-2.
1523                    final int np = nn - 2 * pingPong;
1524                    double b1 = work[np - 2];
1525                    double b2 = work[np - 6];
1526                    final double gam = dN2;
1527                    if (work[np - 8] > b2 || work[np - 4] > b1) {
1528                        return;
1529                    }
1530                    double a2 = (work[np - 8] / b2) * (1 + work[np - 4] / b1);
1531    
1532                    // approximate contribution to norm squared from i < nn-2.
1533                    if (end - start > 2) {
1534                        b2 = work[nn - 13] / work[nn - 15];
1535                        a2 = a2 + b2;
1536                        for (int i4 = nn - 17; i4 >= 4 * start + 2 + pingPong; i4 -= 4) {
1537                            if (b2 == 0.0) {
1538                                break;
1539                            }
1540                            b1 = b2;
1541                            if (work[i4]  >  work[i4 - 2]) {
1542                                return;
1543                            }
1544                            b2 = b2 * (work[i4] / work[i4 - 2]);
1545                            a2 = a2 + b2;
1546                            if (100 * Math.max(b2, b1) < a2 || cnst1 < a2)  {
1547                                break;
1548                            }
1549                        }
1550                        a2 = cnst3 * a2;
1551                    }
1552    
1553                    if (a2 < cnst1) {
1554                        tau = gam * (1 - Math.sqrt(a2)) / (1 + a2);
1555                    } else {
1556                        tau = s;
1557                    }
1558    
1559                } else {
1560    
1561                    // case 6, no information to guide us.
1562                    if (tType == -6) {
1563                        g += 0.333 * (1 - g);
1564                    } else if (tType == -18) {
1565                        g = 0.25 * 0.333;
1566                    } else {
1567                        g = 0.25;
1568                    }
1569                    tau   = g * dMin;
1570                    tType = -6;
1571    
1572                }
1573                break;
1574    
1575            case 1 : // one eigenvalue just deflated. use dMin1, dN1 for dMin and dN.
1576                if (dMin1 == dN1 && dMin2 == dN2) { 
1577    
1578                    // cases 7 and 8.
1579                    tType = -7;
1580                    double s = 0.333 * dMin1;
1581                    if (work[nn - 5] > work[nn - 7]) {
1582                        return;
1583                    }
1584                    double b1 = work[nn - 5] / work[nn - 7];
1585                    double b2 = b1;
1586                    if (b2 != 0.0) {
1587                        for (int i4 = 4 * end - 10 + pingPong; i4 >= 4 * start + 2 + pingPong; i4 -= 4) {
1588                            final double oldB1 = b1;
1589                            if (work[i4] > work[i4 - 2]) {
1590                                return;
1591                            }
1592                            b1 = b1 * (work[i4] / work[i4 - 2]);
1593                            b2 = b2 + b1;
1594                            if (100 * Math.max(b1, oldB1) < b2) {
1595                                break;
1596                            }
1597                        }
1598                    }
1599                    b2 = Math.sqrt(cnst3 * b2);
1600                    final double a2 = dMin1 / (1 + b2 * b2);
1601                    final double gap2 = 0.5 * dMin2 - a2;
1602                    if (gap2 > 0.0 && gap2 > b2 * a2) {
1603                        tau = Math.max(s, a2 * (1 - cnst2 * a2 * (b2 / gap2) * b2));
1604                    } else {
1605                        tau = Math.max(s, a2 * (1 - cnst2 * b2));
1606                        tType = -8;
1607                    }
1608                } else {
1609    
1610                    // case 9.
1611                    tau = 0.25 * dMin1;
1612                    if (dMin1 == dN1) {
1613                        tau = 0.5 * dMin1;
1614                    }
1615                    tType = -9;
1616                }
1617                break;
1618    
1619            case 2 : // two realEigenvalues deflated. use dMin2, dN2 for dMin and dN.
1620    
1621                // cases 10 and 11.
1622                if (dMin2 == dN2 && 2 * work[nn - 5] < work[nn - 7]) { 
1623                    tType = -10;
1624                    final double s = 0.333 * dMin2;
1625                    if (work[nn - 5] > work[nn - 7]) {
1626                        return;
1627                    }
1628                    double b1 = work[nn - 5] / work[nn - 7];
1629                    double b2 = b1;
1630                    if (b2 != 0.0){
1631                        for (int i4 = 4 * end - 9 + pingPong; i4 >= 4 * start + 2 + pingPong; i4 -= 4) {
1632                            if (work[i4] > work[i4 - 2]) {
1633                                return;
1634                            }
1635                            b1 *= work[i4] / work[i4 - 2];
1636                            b2 += b1;
1637                            if (100 * b1 < b2) {
1638                                break;
1639                            }
1640                        }
1641                    }
1642                    b2 = Math.sqrt(cnst3 * b2);
1643                    final double a2 = dMin2 / (1 + b2 * b2);
1644                    final double gap2 = work[nn - 7] + work[nn - 9] -
1645                    Math.sqrt(work[nn - 11]) * Math.sqrt(work[nn - 9]) - a2;
1646                    if (gap2 > 0.0 && gap2 > b2 * a2) {
1647                        tau = Math.max(s, a2 * (1 - cnst2 * a2 * (b2 / gap2) * b2));
1648                    } else {
1649                        tau = Math.max(s, a2 * (1 - cnst2 * b2));
1650                    }
1651                } else {
1652                    tau   = 0.25 * dMin2;
1653                    tType = -11;
1654                }
1655                break;
1656    
1657            default : // case 12, more than two realEigenvalues deflated. no information.
1658                tau   = 0.0;
1659                tType = -12;
1660            }
1661    
1662        }
1663    
1664        /**
1665         * Update sigma.
1666         * @param tau shift to apply to sigma
1667         */
1668        private void updateSigma(final double tau) {
1669            // BEWARE: do NOT attempt to simplify the following statements
1670            // the expressions below take care to accumulate the part of sigma
1671            // that does not fit within a double variable into sigmaLow
1672            if (tau < sigma) {
1673                sigmaLow += tau;
1674                final double t = sigma + sigmaLow;
1675                sigmaLow -= t - sigma;
1676                sigma = t;
1677            } else {
1678                final double t = sigma + tau;
1679                sigmaLow += sigma - (t - tau);
1680                sigma = t;
1681            }
1682        }
1683    
1684        /**
1685         * Find eigenvectors.
1686         */
1687        private void findEigenVectors() {
1688    
1689            final int m = main.length;
1690            eigenvectors = new ArrayRealVector[m];
1691    
1692            // perform an initial non-shifted LDLt decomposition
1693            final double[] d = new double[m];
1694            final double[] l = new double[m - 1];
1695            double di = main[0];
1696            d[0] = di;
1697            for (int i = 1; i < m; ++i) {
1698                final double eiM1  = secondary[i - 1];
1699                final double ratio = eiM1 / di;
1700                di       = main[i] - eiM1 * ratio;
1701                l[i - 1] = ratio;
1702                d[i]     = di;
1703            }
1704    
1705            // compute eigenvectors
1706            for (int i = 0; i < m; ++i) {
1707                eigenvectors[i] = findEigenvector(realEigenvalues[i], d, l);
1708            }
1709    
1710        }
1711    
1712        /**
1713         * Find an eigenvector corresponding to an eigenvalue, using bidiagonals.
1714         * <p>This method corresponds to algorithm X from Dhillon's thesis.</p>
1715         * 
1716         * @param eigenvalue eigenvalue for which eigenvector is desired
1717         * @param d diagonal elements of the initial non-shifted D matrix
1718         * @param l off-diagonal elements of the initial non-shifted L matrix
1719         * @return an eigenvector
1720         */
1721        private ArrayRealVector findEigenvector(final double eigenvalue,
1722                                               final double[] d, final double[] l) {
1723    
1724            // compute the LDLt and UDUt decompositions of the
1725            // perfectly shifted tridiagonal matrix
1726            final int m = main.length;
1727            stationaryQuotientDifferenceWithShift(d, l, eigenvalue);
1728            progressiveQuotientDifferenceWithShift(d, l, eigenvalue);
1729    
1730            // select the twist index leading to
1731            // the least diagonal element in the twisted factorization
1732            int r = m - 1;
1733            double minG = Math.abs(work[6 * r] + work[6 * r + 3] + eigenvalue);
1734            for (int i = 0, sixI = 0; i < m - 1; ++i, sixI += 6) {
1735                final double g = work[sixI] + d[i] * work[sixI + 9] / work[sixI + 10];
1736                final double absG = Math.abs(g);
1737                if (absG < minG) {
1738                    r = i;
1739                    minG = absG;
1740                }
1741            }
1742    
1743            // solve the singular system by ignoring the equation
1744            // at twist index and propagating upwards and downwards
1745            double[] eigenvector = new double[m];
1746            double n2 = 1;
1747            eigenvector[r] = 1;
1748            double z = 1;
1749            for (int i = r - 1; i >= 0; --i) {
1750                z *= -work[6 * i + 2];
1751                eigenvector[i] = z;
1752                n2 += z * z;
1753            }
1754            z = 1;
1755            for (int i = r + 1; i < m; ++i) {
1756                z *= -work[6 * i - 1];
1757                eigenvector[i] = z;
1758                n2 += z * z;
1759            }
1760    
1761            // normalize vector
1762            final double inv = 1.0 / Math.sqrt(n2);
1763            for (int i = 0; i < m; ++i) {
1764                eigenvector[i] *= inv;
1765            }
1766    
1767            return (transformer == null) ?
1768                   new ArrayRealVector(eigenvector, false) :
1769                   new ArrayRealVector(transformer.getQ().operate(eigenvector), false);
1770    
1771        }
1772    
1773        /**
1774         * Decompose matrix LDL<sup>T</sup> - &lambda; I as
1775         * L<sub>+</sub>D<sub>+</sub>L<sub>+</sub><sup>T</sup>.
1776         * <p>This method corresponds to algorithm 4.4.3 (dstqds) from Dhillon's thesis.</p>
1777         * @param d diagonal elements of D,
1778         * @param l off-diagonal elements of L
1779         * @param lambda shift to apply
1780         */
1781        private void stationaryQuotientDifferenceWithShift(final double[] d, final double[] l,
1782                                                           final double lambda) {
1783            final int nM1 = d.length - 1;
1784            double si = -lambda;
1785            for (int i = 0, sixI = 0; i < nM1; ++i, sixI += 6) {
1786                final double di   = d[i];
1787                final double li   = l[i];
1788                final double diP1 = di + si;
1789                final double liP1 = li * di / diP1;
1790                work[sixI]        = si;
1791                work[sixI + 1]    = diP1;
1792                work[sixI + 2]    = liP1;
1793                si = li * liP1 * si - lambda;
1794            }
1795            work[6 * nM1 + 1] = d[nM1] + si;
1796            work[6 * nM1]     = si;
1797        }
1798    
1799        /**
1800         * Decompose matrix LDL<sup>T</sup> - &lambda; I as
1801         * U<sub>-</sub>D<sub>-</sub>U<sub>-</sub><sup>T</sup>.
1802         * <p>This method corresponds to algorithm 4.4.5 (dqds) from Dhillon's thesis.</p>
1803         * @param d diagonal elements of D
1804         * @param l off-diagonal elements of L
1805         * @param lambda shift to apply
1806         */
1807        private void progressiveQuotientDifferenceWithShift(final double[] d, final double[] l,
1808                                                            final double lambda) {
1809            final int nM1 = d.length - 1;
1810            double pi = d[nM1] - lambda;
1811            for (int i = nM1 - 1, sixI = 6 * i; i >= 0; --i, sixI -= 6) {
1812                final double di   = d[i];
1813                final double li   = l[i];
1814                final double diP1 = di * li * li + pi;
1815                final double t    = di / diP1;
1816                work[sixI +  9]   = pi;
1817                work[sixI + 10]   = diP1;
1818                work[sixI +  5]   = li * t;
1819                pi = pi * t - lambda;
1820            }
1821            work[3] = pi;
1822            work[4] = pi;
1823        }
1824    
1825    }