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 × 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 × X = B for symmetric matrices A.
385 * <p>This method only find exact linear solutions, i.e. solutions for
386 * which ||A × X - B|| is exactly 0.</p>
387 * @param b right-hand side of the equation A × X = B
388 * @return a vector X that minimizes the two norm of A × 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 × X = B for symmetric matrices A.
421 * <p>This method only find exact linear solutions, i.e. solutions for
422 * which ||A × X - B|| is exactly 0.</p>
423 * @param b right-hand side of the equation A × X = B
424 * @return a vector X that minimizes the two norm of A × 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 × X = B for symmetric matrices A.
457 * <p>This method only find exact linear solutions, i.e. solutions for
458 * which ||A × X - B|| is exactly 0.</p>
459 * @param b right-hand side of the equation A × X = B
460 * @return a matrix X that minimizes the two norm of A × 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λ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-λ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> - λ 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> - λ 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 }