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