← Back to branch summary

~ubuntu-branches/ubuntu/quantal/commons-math/quantal

« 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: 2010-04-05 23:33:02 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20100405233302-gpqlceked76nw28a
Tags: 2.1-1
* New upstream release.
* Bump Standards-Version to 3.8.4: no changes needed
* Bump debhelper to >= 7
* Switch to 3.0 (quilt) source format:
  - Remove B-D on quilt
  - Add d/source/format
  - Remove d/README.source

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java
17
17
 
18
18
package org.apache.commons.math.linear;
19
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
20
import org.apache.commons.math.MathRuntimeException;
26
21
import org.apache.commons.math.MaxIterationsExceededException;
27
22
import org.apache.commons.math.util.MathUtils;
28
23
 
29
24
/**
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) $
 
25
 * Calculates the eigen decomposition of a real <strong>symmetric</strong>
 
26
 * matrix.
 
27
 * <p>
 
28
 * The eigen decomposition of matrix A is a set of two matrices: V and D such
 
29
 * that A = V D V<sup>T</sup>. A, V and D are all m &times; m matrices.
 
30
 * </p>
 
31
 * <p>
 
32
 * As of 2.0, this class supports only <strong>symmetric</strong> matrices, and
 
33
 * hence computes only real realEigenvalues. This implies the D matrix returned
 
34
 * by {@link #getD()} is always diagonal and the imaginary values returned
 
35
 * {@link #getImagEigenvalue(int)} and {@link #getImagEigenvalues()} are always
 
36
 * null.
 
37
 * </p>
 
38
 * <p>
 
39
 * When called with a {@link RealMatrix} argument, this implementation only uses
 
40
 * the upper part of the matrix, the part below the diagonal is not accessed at
 
41
 * all.
 
42
 * </p>
 
43
 * <p>
 
44
 * This implementation is based on the paper by A. Drubrulle, R.S. Martin and
 
45
 * J.H. Wilkinson 'The Implicit QL Algorithm' in Wilksinson and Reinsch (1971)
 
46
 * Handbook for automatic computation, vol. 2, Linear algebra, Springer-Verlag,
 
47
 * New-York
 
48
 * </p>
 
49
 * @version $Revision: 912413 $ $Date: 2010-02-21 16:46:12 -0500 (Sun, 21 Feb 2010) $
56
50
 * @since 2.0
57
51
 */
58
52
public class EigenDecompositionImpl implements EigenDecomposition {
59
53
 
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;
 
54
    /** Maximum number of iterations accepted in the implicit QL transformation */
 
55
    private byte maxIter = 30;
68
56
 
69
57
    /** Main diagonal of the tridiagonal matrix. */
70
58
    private double[] main;
72
60
    /** Secondary diagonal of the tridiagonal matrix. */
73
61
    private double[] secondary;
74
62
 
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). */
 
63
    /**
 
64
     * Transformer to tridiagonal (may be null if matrix is already
 
65
     * tridiagonal).
 
66
     */
79
67
    private TriDiagonalTransformer transformer;
80
68
 
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
69
    /** Real part of the realEigenvalues. */
136
70
    private double[] realEigenvalues;
137
71
 
151
85
    private RealMatrix cachedVt;
152
86
 
153
87
    /**
154
 
     * Calculates the eigen decomposition of the given symmetric matrix. 
 
88
     * Calculates the eigen decomposition of the given symmetric matrix.
155
89
     * @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
 
90
     * @param splitTolerance dummy parameter, present for backward compatibility only.
 
91
     * @exception InvalidMatrixException (wrapping a
 
92
     * {@link org.apache.commons.math.ConvergenceException} if algorithm
 
93
     * fails to converge
161
94
     */
162
 
    public EigenDecompositionImpl(final RealMatrix matrix,
163
 
                                  final double splitTolerance)
164
 
        throws InvalidMatrixException {
 
95
    public EigenDecompositionImpl(final RealMatrix matrix,final double splitTolerance)
 
96
            throws InvalidMatrixException {
165
97
        if (isSymmetric(matrix)) {
166
 
            this.splitTolerance = splitTolerance;
167
98
            transformToTridiagonal(matrix);
168
 
            decompose();
 
99
            findEigenVectors(transformer.getQ().getData());
169
100
        } else {
170
 
            // as of 2.0, non-symmetric matrices (i.e. complex eigenvalues) are NOT supported
 
101
            // as of 2.0, non-symmetric matrices (i.e. complex eigenvalues) are
 
102
            // NOT supported
171
103
            // see issue https://issues.apache.org/jira/browse/MATH-235
172
 
            throw new InvalidMatrixException("eigen decomposition of assymetric matrices not supported yet");
 
104
            throw new InvalidMatrixException(
 
105
                    "eigen decomposition of assymetric matrices not supported yet");
173
106
        }
174
107
    }
175
108
 
176
109
    /**
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
 
110
     * Calculates the eigen decomposition of the symmetric tridiagonal
 
111
     * matrix.  The Householder matrix is assumed to be the identity matrix.
 
112
     * @param main Main diagonal of the symmetric triadiagonal form
 
113
     * @param secondary Secondary of the tridiagonal form
 
114
     * @param splitTolerance dummy parameter, present for backward compatibility only.
 
115
     * @exception InvalidMatrixException (wrapping a
 
116
     * {@link org.apache.commons.math.ConvergenceException} if algorithm
 
117
     * fails to converge
185
118
     */
186
 
    public EigenDecompositionImpl(final double[] main, double[] secondary,
 
119
    public EigenDecompositionImpl(final double[] main,final double[] secondary,
187
120
            final double splitTolerance)
188
 
        throws InvalidMatrixException {
189
 
 
 
121
            throws InvalidMatrixException {
190
122
        this.main      = main.clone();
191
123
        this.secondary = secondary.clone();
192
124
        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;
 
125
        final int size=main.length;
 
126
        double[][] z = new double[size][size];
 
127
        for (int i=0;i<size;i++) {
 
128
            z[i][i]=1.0;
199
129
        }
200
 
 
201
 
        this.splitTolerance = splitTolerance;
202
 
        decompose();
203
 
 
 
130
        findEigenVectors(z);
204
131
    }
205
132
 
206
133
    /**
207
134
     * Check if a matrix is symmetric.
208
 
     * @param matrix matrix to check
 
135
     * @param matrix
 
136
     *            matrix to check
209
137
     * @return true if matrix is symmetric
210
138
     */
211
139
    private boolean isSymmetric(final RealMatrix matrix) {
212
 
        final int rows    = matrix.getRowDimension();
 
140
        final int rows = matrix.getRowDimension();
213
141
        final int columns = matrix.getColumnDimension();
214
 
        final double eps  = 10 * rows * columns * MathUtils.EPSILON;
 
142
        final double eps = 10 * rows * columns * MathUtils.EPSILON;
215
143
        for (int i = 0; i < rows; ++i) {
216
144
            for (int j = i + 1; j < columns; ++j) {
217
145
                final double mij = matrix.getEntry(i, j);
218
146
                final double mji = matrix.getEntry(j, i);
219
 
                if (Math.abs(mij - mji) > (Math.max(Math.abs(mij), Math.abs(mji)) * eps)) {
 
147
                if (Math.abs(mij - mji) > (Math.max(Math.abs(mij), Math
 
148
                        .abs(mji)) * eps)) {
220
149
                    return false;
221
150
                }
222
151
            }
224
153
        return true;
225
154
    }
226
155
 
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
156
    /** {@inheritDoc} */
251
 
    public RealMatrix getV()
252
 
        throws InvalidMatrixException {
 
157
    public RealMatrix getV() throws InvalidMatrixException {
253
158
 
254
159
        if (cachedV == null) {
255
 
 
256
 
            if (eigenvectors == null) {
257
 
                findEigenVectors();
258
 
            }
259
 
 
260
160
            final int m = eigenvectors.length;
261
161
            cachedV = MatrixUtils.createRealMatrix(m, m);
262
162
            for (int k = 0; k < m; ++k) {
263
163
                cachedV.setColumnVector(k, eigenvectors[k]);
264
164
            }
265
 
 
266
165
        }
267
 
 
268
166
        // return the cached matrix
269
167
        return cachedV;
270
168
 
271
169
    }
272
170
 
273
171
    /** {@inheritDoc} */
274
 
    public RealMatrix getD()
275
 
        throws InvalidMatrixException {
 
172
    public RealMatrix getD() throws InvalidMatrixException {
276
173
        if (cachedD == null) {
277
174
            // cache the matrix for subsequent calls
278
175
            cachedD = MatrixUtils.createRealDiagonalMatrix(realEigenvalues);
281
178
    }
282
179
 
283
180
    /** {@inheritDoc} */
284
 
    public RealMatrix getVT()
285
 
        throws InvalidMatrixException {
 
181
    public RealMatrix getVT() throws InvalidMatrixException {
286
182
 
287
183
        if (cachedVt == null) {
288
 
 
289
 
            if (eigenvectors == null) {
290
 
                findEigenVectors();
291
 
            }
292
 
 
293
184
            final int m = eigenvectors.length;
294
185
            cachedVt = MatrixUtils.createRealMatrix(m, m);
295
186
            for (int k = 0; k < m; ++k) {
300
191
 
301
192
        // return the cached matrix
302
193
        return cachedVt;
303
 
 
304
194
    }
305
195
 
306
196
    /** {@inheritDoc} */
307
 
    public double[] getRealEigenvalues()
308
 
        throws InvalidMatrixException {
 
197
    public double[] getRealEigenvalues() throws InvalidMatrixException {
309
198
        return realEigenvalues.clone();
310
199
    }
311
200
 
312
201
    /** {@inheritDoc} */
313
 
    public double getRealEigenvalue(final int i)
314
 
        throws InvalidMatrixException, ArrayIndexOutOfBoundsException {
 
202
    public double getRealEigenvalue(final int i) throws InvalidMatrixException,
 
203
            ArrayIndexOutOfBoundsException {
315
204
        return realEigenvalues[i];
316
205
    }
317
206
 
318
207
    /** {@inheritDoc} */
319
 
    public double[] getImagEigenvalues()
320
 
        throws InvalidMatrixException {
 
208
    public double[] getImagEigenvalues() throws InvalidMatrixException {
321
209
        return imagEigenvalues.clone();
322
210
    }
323
211
 
324
212
    /** {@inheritDoc} */
325
 
    public double getImagEigenvalue(final int i)
326
 
        throws InvalidMatrixException, ArrayIndexOutOfBoundsException {
 
213
    public double getImagEigenvalue(final int i) throws InvalidMatrixException,
 
214
            ArrayIndexOutOfBoundsException {
327
215
        return imagEigenvalues[i];
328
216
    }
329
217
 
330
218
    /** {@inheritDoc} */
331
219
    public RealVector getEigenvector(final int i)
332
 
        throws InvalidMatrixException, ArrayIndexOutOfBoundsException {
333
 
        if (eigenvectors == null) {
334
 
            findEigenVectors();
335
 
        }
 
220
            throws InvalidMatrixException, ArrayIndexOutOfBoundsException {
336
221
        return eigenvectors[i].copy();
337
222
    }
338
223
 
350
235
 
351
236
    /** {@inheritDoc} */
352
237
    public DecompositionSolver getSolver() {
353
 
        if (eigenvectors == null) {
354
 
            findEigenVectors();
355
 
        }
356
238
        return new Solver(realEigenvalues, imagEigenvalues, eigenvectors);
357
239
    }
358
240
 
359
241
    /** Specialized solver. */
360
242
    private static class Solver implements DecompositionSolver {
361
 
    
 
243
 
362
244
        /** Real part of the realEigenvalues. */
363
245
        private double[] realEigenvalues;
364
246
 
370
252
 
371
253
        /**
372
254
         * 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
 
255
         * @param realEigenvalues
 
256
         *            real parts of the eigenvalues
 
257
         * @param imagEigenvalues
 
258
         *            imaginary parts of the eigenvalues
 
259
         * @param eigenvectors
 
260
         *            eigenvectors
376
261
         */
377
 
        private Solver(final double[] realEigenvalues, final double[] imagEigenvalues,
378
 
                       final ArrayRealVector[] eigenvectors) {
 
262
        private Solver(final double[] realEigenvalues,
 
263
                final double[] imagEigenvalues,
 
264
                final ArrayRealVector[] eigenvectors) {
379
265
            this.realEigenvalues = realEigenvalues;
380
266
            this.imagEigenvalues = imagEigenvalues;
381
 
            this.eigenvectors    = eigenvectors; 
 
267
            this.eigenvectors = eigenvectors;
382
268
        }
383
269
 
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
 
270
        /**
 
271
         * Solve the linear equation A &times; X = B for symmetric matrices A.
 
272
         * <p>
 
273
         * This method only find exact linear solutions, i.e. solutions for
 
274
         * which ||A &times; X - B|| is exactly 0.
 
275
         * </p>
 
276
         * @param b
 
277
         *            right-hand side of the equation A &times; X = B
388
278
         * @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
 
279
         * @exception IllegalArgumentException
 
280
         *                if matrices dimensions don't match
 
281
         * @exception InvalidMatrixException
 
282
         *                if decomposed matrix is singular
391
283
         */
392
284
        public double[] solve(final double[] b)
393
 
            throws IllegalArgumentException, InvalidMatrixException {
 
285
                throws IllegalArgumentException, InvalidMatrixException {
394
286
 
395
287
            if (!isNonSingular()) {
396
288
                throw new SingularMatrixException();
417
309
 
418
310
        }
419
311
 
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
 
312
        /**
 
313
         * Solve the linear equation A &times; X = B for symmetric matrices A.
 
314
         * <p>
 
315
         * This method only find exact linear solutions, i.e. solutions for
 
316
         * which ||A &times; X - B|| is exactly 0.
 
317
         * </p>
 
318
         * @param b
 
319
         *            right-hand side of the equation A &times; X = B
424
320
         * @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
 
321
         * @exception IllegalArgumentException
 
322
         *                if matrices dimensions don't match
 
323
         * @exception InvalidMatrixException
 
324
         *                if decomposed matrix is singular
427
325
         */
428
326
        public RealVector solve(final RealVector b)
429
 
            throws IllegalArgumentException, InvalidMatrixException {
 
327
                throws IllegalArgumentException, InvalidMatrixException {
430
328
 
431
329
            if (!isNonSingular()) {
432
330
                throw new SingularMatrixException();
435
333
            final int m = realEigenvalues.length;
436
334
            if (b.getDimension() != m) {
437
335
                throw MathRuntimeException.createIllegalArgumentException(
438
 
                        "vector length mismatch: got {0} but expected {1}",
439
 
                        b.getDimension(), m);
 
336
                        "vector length mismatch: got {0} but expected {1}", b
 
337
                                .getDimension(), m);
440
338
            }
441
339
 
442
340
            final double[] bp = new double[m];
453
351
 
454
352
        }
455
353
 
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
 
354
        /**
 
355
         * Solve the linear equation A &times; X = B for symmetric matrices A.
 
356
         * <p>
 
357
         * This method only find exact linear solutions, i.e. solutions for
 
358
         * which ||A &times; X - B|| is exactly 0.
 
359
         * </p>
 
360
         * @param b
 
361
         *            right-hand side of the equation A &times; X = B
460
362
         * @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
 
363
         * @exception IllegalArgumentException
 
364
         *                if matrices dimensions don't match
 
365
         * @exception InvalidMatrixException
 
366
         *                if decomposed matrix is singular
463
367
         */
464
368
        public RealMatrix solve(final RealMatrix b)
465
 
            throws IllegalArgumentException, InvalidMatrixException {
 
369
                throws IllegalArgumentException, InvalidMatrixException {
466
370
 
467
371
            if (!isNonSingular()) {
468
372
                throw new SingularMatrixException();
470
374
 
471
375
            final int m = realEigenvalues.length;
472
376
            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");
 
377
                throw MathRuntimeException
 
378
                        .createIllegalArgumentException(
 
379
                                "dimensions mismatch: got {0}x{1} but expected {2}x{3}",
 
380
                                b.getRowDimension(), b.getColumnDimension(), m,
 
381
                                "n");
476
382
            }
477
383
 
478
384
            final int nColB = b.getColumnDimension();
509
415
            return true;
510
416
        }
511
417
 
512
 
        /** Get the inverse of the decomposed matrix.
 
418
        /**
 
419
         * Get the inverse of the decomposed matrix.
513
420
         * @return inverse matrix
514
 
         * @throws InvalidMatrixException if decomposed matrix is singular
 
421
         * @throws InvalidMatrixException
 
422
         *             if decomposed matrix is singular
515
423
         */
516
 
        public RealMatrix getInverse()
517
 
            throws InvalidMatrixException {
 
424
        public RealMatrix getInverse() throws InvalidMatrixException {
518
425
 
519
426
            if (!isNonSingular()) {
520
427
                throw new SingularMatrixException();
542
449
 
543
450
    /**
544
451
     * Transform matrix to tridiagonal.
545
 
     * @param matrix matrix to transform
 
452
     * @param matrix
 
453
     *            matrix to transform
546
454
     */
547
455
    private void transformToTridiagonal(final RealMatrix matrix) {
548
456
 
549
457
        // transform the matrix to tridiagonal
550
458
        transformer = new TriDiagonalTransformer(matrix);
551
 
        main      = transformer.getMainDiagonalRef();
 
459
        main = transformer.getMainDiagonalRef();
552
460
        secondary = transformer.getSecondaryDiagonalRef();
553
461
 
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];
 
462
    }
 
463
 
 
464
    /**
 
465
     * Find eigenvalues and eigenvectors (Dubrulle et al., 1971)
 
466
     * @param householderMatrix Householder matrix of the transformation
 
467
     *  to tri-diagonal form.
 
468
     */
 
469
    private void findEigenVectors(double[][] householderMatrix) {
 
470
 
 
471
        double[][]z = householderMatrix.clone();
 
472
        final int n = main.length;
 
473
        realEigenvalues = new double[n];
 
474
        imagEigenvalues = new double[n];
 
475
        double[] e = new double[n];
 
476
        for (int i = 0; i < n - 1; i++) {
 
477
            realEigenvalues[i] = main[i];
 
478
            e[i] = secondary[i];
 
479
        }
 
480
        realEigenvalues[n - 1] = main[n - 1];
 
481
        e[n - 1] = 0.0;
 
482
 
 
483
        // Determine the largest main and secondary value in absolute term.
 
484
        double maxAbsoluteValue=0.0;
 
485
        for (int i = 0; i < n; i++) {
 
486
            if (Math.abs(realEigenvalues[i])>maxAbsoluteValue) {
 
487
                maxAbsoluteValue=Math.abs(realEigenvalues[i]);
 
488
            }
 
489
            if (Math.abs(e[i])>maxAbsoluteValue) {
 
490
                maxAbsoluteValue=Math.abs(e[i]);
 
491
            }
 
492
        }
 
493
        // Make null any main and secondary value too small to be significant
 
494
        if (maxAbsoluteValue!=0.0) {
 
495
            for (int i=0; i < n; i++) {
 
496
                if (Math.abs(realEigenvalues[i])<=MathUtils.EPSILON*maxAbsoluteValue) {
 
497
                    realEigenvalues[i]=0.0;
 
498
                }
 
499
                if (Math.abs(e[i])<=MathUtils.EPSILON*maxAbsoluteValue) {
 
500
                    e[i]=0.0;
 
501
                }
 
502
            }
 
503
        }
 
504
 
 
505
        for (int j = 0; j < n; j++) {
 
506
            int its = 0;
 
507
            int m;
 
508
            do {
 
509
                for (m = j; m < n - 1; m++) {
 
510
                    double delta = Math.abs(realEigenvalues[m]) + Math.abs(realEigenvalues[m + 1]);
 
511
                    if (Math.abs(e[m]) + delta == delta) {
 
512
                        break;
 
513
                    }
 
514
                }
 
515
                if (m != j) {
 
516
                    if (its == maxIter)
 
517
                        throw new InvalidMatrixException(
 
518
                                new MaxIterationsExceededException(maxIter));
 
519
                    its++;
 
520
                    double q = (realEigenvalues[j + 1] - realEigenvalues[j]) / (2 * e[j]);
 
521
                    double t = Math.sqrt(1 + q * q);
 
522
                    if (q < 0.0) {
 
523
                        q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q - t);
 
524
                    } else {
 
525
                        q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q + t);
 
526
                    }
 
527
                    double u = 0.0;
 
528
                    double s = 1.0;
 
529
                    double c = 1.0;
 
530
                    int i;
 
531
                    for (i = m - 1; i >= j; i--) {
 
532
                        double p = s * e[i];
 
533
                        double h = c * e[i];
 
534
                        if (Math.abs(p) >= Math.abs(q)) {
 
535
                            c = q / p;
 
536
                            t = Math.sqrt(c * c + 1.0);
 
537
                            e[i + 1] = p * t;
 
538
                            s = 1.0 / t;
 
539
                            c = c * s;
906
540
                        } 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
 
 
 
541
                            s = p / q;
 
542
                            t = Math.sqrt(s * s + 1.0);
 
543
                            e[i + 1] = q * t;
 
544
                            c = 1.0 / t;
 
545
                            s = s * c;
 
546
                        }
 
547
                        if (e[i + 1] == 0.0) {
 
548
                            realEigenvalues[i + 1] -= u;
 
549
                            e[m] = 0.0;
 
550
                            break;
 
551
                        }
 
552
                        q = realEigenvalues[i + 1] - u;
 
553
                        t = (realEigenvalues[i] - q) * s + 2.0 * c * h;
 
554
                        u = s * t;
 
555
                        realEigenvalues[i + 1] = q + u;
 
556
                        q = c * t - h;
 
557
                        for (int ia = 0; ia < n; ia++) {
 
558
                            p = z[ia][i + 1];
 
559
                            z[ia][i + 1] = s * z[ia][i] + c * p;
 
560
                            z[ia][i] = c * z[ia][i] - s * p;
 
561
                        }
 
562
                    }
 
563
                    if (e[i + 1] == 0.0 && i >= j)
 
564
                        continue;
 
565
                    realEigenvalues[j] -= u;
 
566
                    e[j] = q;
 
567
                    e[m] = 0.0;
 
568
                }
 
569
            } while (m != j);
 
570
        }
 
571
 
 
572
        //Sort the eigen values (and vectors) in increase order
 
573
        for (int i = 0; i < n; i++) {
 
574
            int k = i;
 
575
            double p = realEigenvalues[i];
 
576
            for (int j = i + 1; j < n; j++) {
 
577
                if (realEigenvalues[j] > p) {
 
578
                    k = j;
 
579
                    p = realEigenvalues[j];
 
580
                }
 
581
            }
 
582
            if (k != i) {
 
583
                realEigenvalues[k] = realEigenvalues[i];
 
584
                realEigenvalues[i] = p;
 
585
                for (int j = 0; j < n; j++) {
 
586
                    p = z[j][i];
 
587
                    z[j][i] = z[j][k];
 
588
                    z[j][k] = p;
 
589
                }
 
590
            }
 
591
        }
 
592
 
 
593
        // Determine the largest eigen value in absolute term.
 
594
        maxAbsoluteValue=0.0;
 
595
        for (int i = 0; i < n; i++) {
 
596
            if (Math.abs(realEigenvalues[i])>maxAbsoluteValue) {
 
597
                maxAbsoluteValue=Math.abs(realEigenvalues[i]);
 
598
            }
 
599
        }
 
600
        // Make null any eigen value too small to be significant
 
601
        if (maxAbsoluteValue!=0.0) {
 
602
            for (int i=0; i < n; i++) {
 
603
                if (Math.abs(realEigenvalues[i])<MathUtils.EPSILON*maxAbsoluteValue) {
 
604
                    realEigenvalues[i]=0.0;
 
605
                }
 
606
            }
 
607
        }
 
608
        eigenvectors = new ArrayRealVector[n];
 
609
        double[] tmp = new double[n];
 
610
        for (int i = 0; i < n; i++) {
 
611
            for (int j = 0; j < n; j++) {
 
612
                tmp[j] = z[j][i];
 
613
            }
 
614
            eigenvectors[i] = new ArrayRealVector(tmp);
 
615
        }
 
616
    }
1825
617
}