View Javadoc

1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  package org.apache.commons.math.linear;
19  
20  import java.util.ArrayList;
21  import java.util.Arrays;
22  import java.util.List;
23  
24  import org.apache.commons.math.ConvergenceException;
25  import org.apache.commons.math.MathRuntimeException;
26  import org.apache.commons.math.MaxIterationsExceededException;
27  import org.apache.commons.math.util.MathUtils;
28  
29  /**
30   * Calculates the eigen decomposition of a <strong>symmetric</strong> matrix.
31   * <p>The eigen decomposition of matrix A is a set of two matrices:
32   * V and D such that A = V D V<sup>T</sup>. A, V and D are all m &times; m
33   * matrices.</p>
34   * <p>As of 2.0, this class supports only <strong>symmetric</strong> matrices,
35   * and hence computes only real realEigenvalues. This implies the D matrix returned by
36   * {@link #getD()} is always diagonal and the imaginary values returned {@link
37   * #getImagEigenvalue(int)} and {@link #getImagEigenvalues()} are always null.</p>
38   * <p>When called with a {@link RealMatrix} argument, this implementation only uses
39   * the upper part of the matrix, the part below the diagonal is not accessed at all.</p>
40   * <p>Eigenvalues are computed as soon as the matrix is decomposed, but eigenvectors
41   * are computed only when required, i.e. only when one of the {@link #getEigenvector(int)},
42   * {@link #getV()}, {@link #getVT()}, {@link #getSolver()} methods is called.</p>
43   * <p>This implementation is based on Inderjit Singh Dhillon thesis
44   * <a href="http://www.cs.utexas.edu/users/inderjit/public_papers/thesis.pdf">A
45   * New O(n<sup>2</sup>) Algorithm for the Symmetric Tridiagonal Eigenvalue/Eigenvector
46   * Problem</a>, on Beresford N. Parlett and Osni A. Marques paper <a
47   * href="http://www.netlib.org/lapack/lawnspdf/lawn155.pdf">An Implementation of the
48   * dqds Algorithm (Positive Case)</a> and on the corresponding LAPACK routines (DLARRE,
49   * DLASQ2, DLAZQ3, DLAZQ4, DLASQ5 and DLASQ6).</p>
50   * @author Beresford Parlett, University of California, Berkeley, USA (fortran version)
51   * @author Jim Demmel, University of California, Berkeley, USA (fortran version)
52   * @author Inderjit Dhillon, University of Texas, Austin, USA(fortran version)
53   * @author Osni Marques, LBNL/NERSC, USA (fortran version)
54   * @author Christof Voemel, University of California, Berkeley, USA(fortran version)
55   * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
56   * @since 2.0
57   */
58  public class EigenDecompositionImpl implements EigenDecomposition {
59  
60      /** Tolerance. */
61      private static final double TOLERANCE = 100 * MathUtils.EPSILON;
62  
63      /** Squared tolerance. */
64      private static final double TOLERANCE_2 = TOLERANCE * TOLERANCE;
65  
66      /** Split tolerance. */
67      private double splitTolerance;
68  
69      /** Main diagonal of the tridiagonal matrix. */
70      private double[] main;
71  
72      /** Secondary diagonal of the tridiagonal matrix. */
73      private double[] secondary;
74  
75      /** Squared secondary diagonal of the tridiagonal matrix. */
76      private double[] squaredSecondary;
77  
78      /** Transformer to tridiagonal (may be null if matrix is already tridiagonal). */
79      private TriDiagonalTransformer transformer;
80  
81      /** Lower bound of spectra. */
82      private double lowerSpectra;
83  
84      /** Upper bound of spectra. */
85      private double upperSpectra;
86  
87      /** Minimum pivot in the Sturm sequence. */
88      private double minPivot;
89  
90      /** Current shift. */
91      private double sigma;
92  
93      /** Low part of the current shift. */
94      private double sigmaLow;
95  
96      /** Shift increment to apply. */
97      private double tau;
98  
99      /** 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 }