← Back to branch summary

~ubuntu-branches/ubuntu/maverick/commons-math/maverick

« back to all changes in this revision

Viewing changes to src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java

  • Committer: Bazaar Package Importer
  • Author(s): Damien Raude-Morvan
  • Date: 2009-08-22 01:13:25 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20090822011325-hi4peq1ua5weguwn
Tags: 2.0-1
* New upstream release.
* Set Maintainer field to Debian Java Team
* Add myself as Uploaders
* Switch to Quilt patch system:
  - Refresh all patchs
  - Remove B-D on dpatch, Add B-D on quilt
  - Include patchsys-quilt.mk in debian/rules
* Bump Standards-Version to 3.8.3:
  - Add a README.source to describe patch system
* Maven POMs:
  - Add a Build-Depends-Indep dependency on maven-repo-helper
  - Use mh_installpom and mh_installjar to install the POM and the jar to the
    Maven repository
* Use default-jdk/jre:
  - Depends on java5-runtime-headless
  - Build-Depends on default-jdk
  - Use /usr/lib/jvm/default-java as JAVA_HOME
* Move api documentation to /usr/share/doc/libcommons-math-java/api
* Build-Depends on junit4 instead of junit

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java
 
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
}