001 /* 002 * Copyright (C) 2008-2010 by Holger Arndt 003 * 004 * This file is part of the Universal Java Matrix Package (UJMP). 005 * See the NOTICE file distributed with this work for additional 006 * information regarding copyright ownership and licensing. 007 * 008 * UJMP is free software; you can redistribute it and/or modify 009 * it under the terms of the GNU Lesser General Public License as 010 * published by the Free Software Foundation; either version 2 011 * of the License, or (at your option) any later version. 012 * 013 * UJMP is distributed in the hope that it will be useful, 014 * but WITHOUT ANY WARRANTY; without even the implied warranty of 015 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 016 * GNU Lesser General Public License for more details. 017 * 018 * You should have received a copy of the GNU Lesser General Public 019 * License along with UJMP; if not, write to the 020 * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, 021 * Boston, MA 02110-1301 USA 022 */ 023 024 package org.ujmp.core.doublematrix.calculation.general.decomposition; 025 026 import org.ujmp.core.Matrix; 027 import org.ujmp.core.MatrixFactory; 028 import org.ujmp.core.util.DecompositionOps; 029 import org.ujmp.core.util.MathUtil; 030 import org.ujmp.core.util.UJMPSettings; 031 032 /** 033 * Eigenvalues and eigenvectors of a real matrix. 034 * <P> 035 * If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is diagonal 036 * and the eigenvector matrix V is orthogonal. I.e. A = 037 * V.times(D.times(V.transpose())) and V.times(V.transpose()) equals the 038 * identity matrix. 039 * <P> 040 * If A is not symmetric, then the eigenvalue matrix D is block diagonal with 041 * the real eigenvalues in 1-by-1 blocks and any complex eigenvalues, lambda + 042 * i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda]. The columns of V represent 043 * the eigenvectors in the sense that A*V = V*D, i.e. A.times(V) equals 044 * V.times(D). The matrix V may be badly conditioned, or even singular, so the 045 * validity of the equation A = V*D*inverse(V) depends upon V.cond(). 046 **/ 047 048 public interface Eig<T> { 049 050 public static int THRESHOLD = 100; 051 052 public T[] calc(T source); 053 054 public static final Eig<Matrix> MATRIX = new Eig<Matrix>() { 055 056 public final Matrix[] calc(Matrix source) { 057 if (UJMPSettings.getNumberOfThreads() == 1) { 058 if (source.getRowCount() >= THRESHOLD && source.getColumnCount() >= THRESHOLD) { 059 return MATRIXLARGESINGLETHREADED.calc(source); 060 } else { 061 return MATRIXSMALLSINGLETHREADED.calc(source); 062 } 063 } else { 064 if (source.getRowCount() >= THRESHOLD && source.getColumnCount() >= THRESHOLD) { 065 return MATRIXLARGEMULTITHREADED.calc(source); 066 } else { 067 return MATRIXSMALLMULTITHREADED.calc(source); 068 } 069 } 070 } 071 }; 072 073 public static final Eig<Matrix> MATRIXLARGESINGLETHREADED = new Eig<Matrix>() { 074 public Matrix[] calc(Matrix source) { 075 Eig<Matrix> eig = null; 076 if (UJMPSettings.isUseJBlas()) { 077 eig = DecompositionOps.EIG_JBLAS; 078 } 079 if (eig == null && UJMPSettings.isUseEJML()) { 080 eig = DecompositionOps.EIG_EJML; 081 } 082 if (eig == null) { 083 eig = UJMP; 084 } 085 return eig.calc(source); 086 } 087 }; 088 089 public static final Eig<Matrix> MATRIXLARGEMULTITHREADED = new Eig<Matrix>() { 090 public Matrix[] calc(Matrix source) { 091 Eig<Matrix> eig = null; 092 if (UJMPSettings.isUseJBlas()) { 093 eig = DecompositionOps.EIG_JBLAS; 094 } 095 if (eig == null && UJMPSettings.isUseOjalgo()) { 096 eig = DecompositionOps.EIG_OJALGO; 097 } 098 if (eig == null && UJMPSettings.isUseEJML()) { 099 eig = DecompositionOps.EIG_EJML; 100 } 101 if (eig == null) { 102 eig = UJMP; 103 } 104 return eig.calc(source); 105 } 106 }; 107 108 public static final Eig<Matrix> INSTANCE = MATRIX; 109 110 public static final Eig<Matrix> UJMP = new Eig<Matrix>() { 111 112 public final Matrix[] calc(Matrix source) { 113 EigMatrix qr = new EigMatrix(source); 114 return new Matrix[] { qr.getV(), qr.getD() }; 115 } 116 }; 117 118 public static final Eig<Matrix> MATRIXSMALLMULTITHREADED = UJMP; 119 120 public static final Eig<Matrix> MATRIXSMALLSINGLETHREADED = UJMP; 121 122 final class EigMatrix { 123 private static final long serialVersionUID = -4312402808395971553L; 124 125 private static final double EPSILON = Math.pow(2.0, -52.0); 126 127 /** 128 * Row and column dimension (square matrix). 129 * 130 * @serial matrix dimension. 131 */ 132 private final int n; 133 134 /** 135 * Symmetry flag. 136 * 137 * @serial internal symmetry flag. 138 */ 139 private boolean issymmetric; 140 141 /** 142 * Arrays for internal storage of eigenvalues. 143 * 144 * @serial internal storage of eigenvalues. 145 */ 146 private final double[] d, e; 147 148 /** 149 * Array for internal storage of eigenvectors. 150 * 151 * @serial internal storage of eigenvectors. 152 */ 153 private final double[][] V; 154 155 /** 156 * Array for internal storage of nonsymmetric Hessenberg form. 157 * 158 * @serial internal storage of nonsymmetric Hessenberg form. 159 */ 160 private final double[][] H; 161 162 /** 163 * Working storage for nonsymmetric algorithm. 164 * 165 * @serial working storage for nonsymmetric algorithm. 166 */ 167 private final double[] ort; 168 169 /* 170 * ------------------------ Private Methods ------------------------ 171 */ 172 173 // Symmetric Householder reduction to tridiagonal form. 174 private final void tred2() { 175 176 // This is derived from the Algol procedures tred2 by 177 // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for 178 // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding 179 // Fortran subroutine in EISPACK. 180 181 final double[] Vn1 = V[n - 1]; 182 183 for (int j = 0; j < n; j++) { 184 d[j] = Vn1[j]; 185 } 186 187 // Householder reduction to tridiagonal form. 188 189 for (int i = n - 1; i > 0; i--) { 190 191 // Scale to avoid under/overflow. 192 193 double scale = 0.0; 194 double h = 0.0; 195 for (int k = 0; k < i; k++) { 196 scale = scale + Math.abs(d[k]); 197 } 198 if (scale == 0.0) { 199 e[i] = d[i - 1]; 200 for (int j = 0; j < i; j++) { 201 d[j] = V[i - 1][j]; 202 V[i][j] = 0.0; 203 V[j][i] = 0.0; 204 } 205 } else { 206 207 // Generate Householder vector. 208 209 for (int k = 0; k < i; k++) { 210 d[k] /= scale; 211 h += d[k] * d[k]; 212 } 213 double f = d[i - 1]; 214 double g = Math.sqrt(h); 215 if (f > 0) { 216 g = -g; 217 } 218 e[i] = scale * g; 219 h = h - f * g; 220 d[i - 1] = f - g; 221 for (int j = 0; j < i; j++) { 222 e[j] = 0.0; 223 } 224 225 // Apply similarity transformation to remaining columns. 226 227 for (int j = 0; j < i; j++) { 228 f = d[j]; 229 V[j][i] = f; 230 g = e[j] + V[j][j] * f; 231 for (int k = j + 1; k <= i - 1; k++) { 232 g += V[k][j] * d[k]; 233 e[k] += V[k][j] * f; 234 } 235 e[j] = g; 236 } 237 f = 0.0; 238 for (int j = 0; j < i; j++) { 239 e[j] /= h; 240 f += e[j] * d[j]; 241 } 242 final double hh = f / (h + h); 243 for (int j = 0; j < i; j++) { 244 e[j] -= hh * d[j]; 245 } 246 for (int j = 0; j < i; j++) { 247 f = d[j]; 248 g = e[j]; 249 for (int k = j; k <= i - 1; k++) { 250 V[k][j] -= (f * e[k] + g * d[k]); 251 } 252 d[j] = V[i - 1][j]; 253 V[i][j] = 0.0; 254 } 255 } 256 d[i] = h; 257 } 258 259 // Accumulate transformations. 260 261 for (int i = 0; i < n - 1; i++) { 262 Vn1[i] = V[i][i]; 263 V[i][i] = 1.0; 264 final double h = d[i + 1]; 265 if (h != 0.0) { 266 for (int k = 0; k <= i; k++) { 267 d[k] = V[k][i + 1] / h; 268 } 269 for (int j = 0; j <= i; j++) { 270 double g = 0.0; 271 for (int k = 0; k <= i; k++) { 272 g += V[k][i + 1] * V[k][j]; 273 } 274 for (int k = 0; k <= i; k++) { 275 V[k][j] -= g * d[k]; 276 } 277 } 278 } 279 for (int k = 0; k <= i; k++) { 280 V[k][i + 1] = 0.0; 281 } 282 } 283 for (int j = 0; j < n; j++) { 284 d[j] = Vn1[j]; 285 Vn1[j] = 0.0; 286 } 287 Vn1[n - 1] = 1.0; 288 e[0] = 0.0; 289 } 290 291 // Symmetric tridiagonal QL algorithm. 292 293 private final void tql2() { 294 295 // This is derived from the Algol procedures tql2, by 296 // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for 297 // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding 298 // Fortran subroutine in EISPACK. 299 300 for (int i = 1; i < n; i++) { 301 e[i - 1] = e[i]; 302 } 303 e[n - 1] = 0.0; 304 305 double f = 0.0; 306 double tst1 = 0.0; 307 for (int l = 0; l < n; l++) { 308 309 // Find small subdiagonal element 310 311 tst1 = Math.max(tst1, Math.abs(d[l]) + Math.abs(e[l])); 312 int m = l; 313 while (m < n) { 314 if (Math.abs(e[m]) <= EPSILON * tst1) { 315 break; 316 } 317 m++; 318 } 319 320 // If m == l, d[l] is an eigenvalue, 321 // otherwise, iterate. 322 323 if (m > l) { 324 int iter = 0; 325 do { 326 iter = iter + 1; // (Could check iteration count here.) 327 328 // Compute implicit shift 329 330 double g = d[l]; 331 double p = (d[l + 1] - g) / (2.0 * e[l]); 332 double r = MathUtil.hypot(p, 1.0); 333 if (p < 0) { 334 r = -r; 335 } 336 d[l] = e[l] / (p + r); 337 d[l + 1] = e[l] * (p + r); 338 final double dl1 = d[l + 1]; 339 double h = g - d[l]; 340 for (int i = l + 2; i < n; i++) { 341 d[i] -= h; 342 } 343 f = f + h; 344 345 // Implicit QL transformation. 346 347 p = d[m]; 348 double c = 1.0; 349 double c2 = c; 350 double c3 = c; 351 double el1 = e[l + 1]; 352 double s = 0.0; 353 double s2 = 0.0; 354 for (int i = m - 1; i >= l; i--) { 355 c3 = c2; 356 c2 = c; 357 s2 = s; 358 g = c * e[i]; 359 h = c * p; 360 r = MathUtil.hypot(p, e[i]); 361 e[i + 1] = s * r; 362 s = e[i] / r; 363 c = p / r; 364 p = c * d[i] - s * g; 365 d[i + 1] = h + s * (c * g + s * d[i]); 366 367 // Accumulate transformation. 368 369 for (int k = 0; k < n; k++) { 370 h = V[k][i + 1]; 371 V[k][i + 1] = s * V[k][i] + c * h; 372 V[k][i] = c * V[k][i] - s * h; 373 } 374 } 375 p = -s * s2 * c3 * el1 * e[l] / dl1; 376 e[l] = s * p; 377 d[l] = c * p; 378 379 // Check for convergence. 380 381 } while (Math.abs(e[l]) > EPSILON * tst1); 382 } 383 d[l] = d[l] + f; 384 e[l] = 0.0; 385 } 386 387 // Sort eigenvalues and corresponding vectors. 388 389 for (int i = 0; i < n - 1; i++) { 390 int k = i; 391 double p = d[i]; 392 for (int j = i + 1; j < n; j++) { 393 if (d[j] < p) { 394 k = j; 395 p = d[j]; 396 } 397 } 398 if (k != i) { 399 d[k] = d[i]; 400 d[i] = p; 401 for (int j = 0; j < n; j++) { 402 p = V[j][i]; 403 V[j][i] = V[j][k]; 404 V[j][k] = p; 405 } 406 } 407 } 408 } 409 410 // Nonsymmetric reduction to Hessenberg form. 411 private final void orthes() { 412 413 // This is derived from the Algol procedures orthes and ortran, 414 // by Martin and Wilkinson, Handbook for Auto. Comp., 415 // Vol.ii-Linear Algebra, and the corresponding 416 // Fortran subroutines in EISPACK. 417 418 final int high = n - 1; 419 420 for (int m = 1; m <= high - 1; m++) { 421 422 // Scale column. 423 424 double scale = 0.0; 425 for (int i = m; i <= high; i++) { 426 scale = scale + Math.abs(H[i][m - 1]); 427 } 428 if (scale != 0.0) { 429 430 // Compute Householder transformation. 431 432 double h = 0.0; 433 for (int i = high; i >= m; i--) { 434 ort[i] = H[i][m - 1] / scale; 435 h += ort[i] * ort[i]; 436 } 437 double g = Math.sqrt(h); 438 if (ort[m] > 0) { 439 g = -g; 440 } 441 h = h - ort[m] * g; 442 ort[m] = ort[m] - g; 443 444 // Apply Householder similarity transformation 445 // H = (I-u*u'/h)*H*(I-u*u')/h) 446 447 for (int j = m; j < n; j++) { 448 double f = 0.0; 449 for (int i = high; i >= m; i--) { 450 f += ort[i] * H[i][j]; 451 } 452 f = f / h; 453 for (int i = m; i <= high; i++) { 454 H[i][j] -= f * ort[i]; 455 } 456 } 457 458 for (int i = 0; i <= high; i++) { 459 double f = 0.0; 460 for (int j = high; j >= m; j--) { 461 f += ort[j] * H[i][j]; 462 } 463 f = f / h; 464 for (int j = m; j <= high; j++) { 465 H[i][j] -= f * ort[j]; 466 } 467 } 468 ort[m] = scale * ort[m]; 469 H[m][m - 1] = scale * g; 470 } 471 } 472 473 // Accumulate transformations (Algol's ortran). 474 475 for (int i = 0; i < n; i++) { 476 for (int j = 0; j < n; j++) { 477 V[i][j] = (i == j ? 1.0 : 0.0); 478 } 479 } 480 481 for (int m = high - 1; m >= 1; m--) { 482 if (H[m][m - 1] != 0.0) { 483 for (int i = m + 1; i <= high; i++) { 484 ort[i] = H[i][m - 1]; 485 } 486 for (int j = m; j <= high; j++) { 487 double g = 0.0; 488 for (int i = m; i <= high; i++) { 489 g += ort[i] * V[i][j]; 490 } 491 // Double division avoids possible underflow 492 g = (g / ort[m]) / H[m][m - 1]; 493 for (int i = m; i <= high; i++) { 494 V[i][j] += g * ort[i]; 495 } 496 } 497 } 498 } 499 } 500 501 // Complex scalar division. 502 503 private transient double cdivr, cdivi; 504 505 private final void cdiv(double xr, double xi, double yr, double yi) { 506 double r, d; 507 if (Math.abs(yr) > Math.abs(yi)) { 508 r = yi / yr; 509 d = yr + r * yi; 510 cdivr = (xr + r * xi) / d; 511 cdivi = (xi - r * xr) / d; 512 } else { 513 r = yr / yi; 514 d = yi + r * yr; 515 cdivr = (r * xr + xi) / d; 516 cdivi = (r * xi - xr) / d; 517 } 518 } 519 520 // Nonsymmetric reduction from Hessenberg to real Schur form. 521 522 private final void hqr2() { 523 524 // This is derived from the Algol procedure hqr2, 525 // by Martin and Wilkinson, Handbook for Auto. Comp., 526 // Vol.ii-Linear Algebra, and the corresponding 527 // Fortran subroutine in EISPACK. 528 529 // Initialize 530 531 final int nn = this.n; 532 int n = nn - 1; 533 int low = 0; 534 int high = nn - 1; 535 536 double exshift = 0.0; 537 double p = 0, q = 0, r = 0, s = 0, z = 0, t, w, x, y; 538 539 // Store roots isolated by balanc and compute matrix norm 540 541 double norm = 0.0; 542 for (int i = 0; i < nn; i++) { 543 if (i < low | i > high) { 544 d[i] = H[i][i]; 545 e[i] = 0.0; 546 } 547 for (int j = Math.max(i - 1, 0); j < nn; j++) { 548 norm = norm + Math.abs(H[i][j]); 549 } 550 } 551 552 // Outer loop over eigenvalue index 553 554 int iter = 0; 555 while (n >= low) { 556 557 // Look for single small sub-diagonal element 558 559 int l = n; 560 while (l > low) { 561 s = Math.abs(H[l - 1][l - 1]) + Math.abs(H[l][l]); 562 if (s == 0.0) { 563 s = norm; 564 } 565 if (Math.abs(H[l][l - 1]) < EPSILON * s) { 566 break; 567 } 568 l--; 569 } 570 571 // Check for convergence 572 // One root found 573 574 if (l == n) { 575 H[n][n] = H[n][n] + exshift; 576 d[n] = H[n][n]; 577 e[n] = 0.0; 578 n--; 579 iter = 0; 580 581 // Two roots found 582 583 } else if (l == n - 1) { 584 w = H[n][n - 1] * H[n - 1][n]; 585 p = (H[n - 1][n - 1] - H[n][n]) / 2.0; 586 q = p * p + w; 587 z = Math.sqrt(Math.abs(q)); 588 H[n][n] = H[n][n] + exshift; 589 H[n - 1][n - 1] = H[n - 1][n - 1] + exshift; 590 x = H[n][n]; 591 592 // Real pair 593 594 if (q >= 0) { 595 if (p >= 0) { 596 z = p + z; 597 } else { 598 z = p - z; 599 } 600 d[n - 1] = x + z; 601 d[n] = d[n - 1]; 602 if (z != 0.0) { 603 d[n] = x - w / z; 604 } 605 e[n - 1] = 0.0; 606 e[n] = 0.0; 607 x = H[n][n - 1]; 608 s = Math.abs(x) + Math.abs(z); 609 p = x / s; 610 q = z / s; 611 r = Math.sqrt(p * p + q * q); 612 p = p / r; 613 q = q / r; 614 615 // Row modification 616 617 for (int j = n - 1; j < nn; j++) { 618 z = H[n - 1][j]; 619 H[n - 1][j] = q * z + p * H[n][j]; 620 H[n][j] = q * H[n][j] - p * z; 621 } 622 623 // Column modification 624 625 for (int i = 0; i <= n; i++) { 626 z = H[i][n - 1]; 627 H[i][n - 1] = q * z + p * H[i][n]; 628 H[i][n] = q * H[i][n] - p * z; 629 } 630 631 // Accumulate transformations 632 633 for (int i = low; i <= high; i++) { 634 z = V[i][n - 1]; 635 V[i][n - 1] = q * z + p * V[i][n]; 636 V[i][n] = q * V[i][n] - p * z; 637 } 638 639 // Complex pair 640 641 } else { 642 d[n - 1] = x + p; 643 d[n] = x + p; 644 e[n - 1] = z; 645 e[n] = -z; 646 } 647 n = n - 2; 648 iter = 0; 649 650 // No convergence yet 651 652 } else { 653 654 // Form shift 655 656 x = H[n][n]; 657 y = 0.0; 658 w = 0.0; 659 if (l < n) { 660 y = H[n - 1][n - 1]; 661 w = H[n][n - 1] * H[n - 1][n]; 662 } 663 664 // Wilkinson's original ad hoc shift 665 666 if (iter == 10) { 667 exshift += x; 668 for (int i = low; i <= n; i++) { 669 H[i][i] -= x; 670 } 671 s = Math.abs(H[n][n - 1]) + Math.abs(H[n - 1][n - 2]); 672 x = y = 0.75 * s; 673 w = -0.4375 * s * s; 674 } 675 676 // MATLAB's new ad hoc shift 677 678 if (iter == 30) { 679 s = (y - x) / 2.0; 680 s = s * s + w; 681 if (s > 0) { 682 s = Math.sqrt(s); 683 if (y < x) { 684 s = -s; 685 } 686 s = x - w / ((y - x) / 2.0 + s); 687 for (int i = low; i <= n; i++) { 688 H[i][i] -= s; 689 } 690 exshift += s; 691 x = y = w = 0.964; 692 } 693 } 694 695 iter = iter + 1; // (Could check iteration count here.) 696 697 // Look for two consecutive small sub-diagonal elements 698 699 int m = n - 2; 700 while (m >= l) { 701 z = H[m][m]; 702 r = x - z; 703 s = y - z; 704 p = (r * s - w) / H[m + 1][m] + H[m][m + 1]; 705 q = H[m + 1][m + 1] - z - r - s; 706 r = H[m + 2][m + 1]; 707 s = Math.abs(p) + Math.abs(q) + Math.abs(r); 708 p = p / s; 709 q = q / s; 710 r = r / s; 711 if (m == l) { 712 break; 713 } 714 if (Math.abs(H[m][m - 1]) * (Math.abs(q) + Math.abs(r)) < EPSILON 715 * (Math.abs(p) * (Math.abs(H[m - 1][m - 1]) + Math.abs(z) + Math 716 .abs(H[m + 1][m + 1])))) { 717 break; 718 } 719 m--; 720 } 721 722 for (int i = m + 2; i <= n; i++) { 723 H[i][i - 2] = 0.0; 724 if (i > m + 2) { 725 H[i][i - 3] = 0.0; 726 } 727 } 728 729 // Double QR step involving rows l:n and columns m:n 730 731 for (int k = m; k <= n - 1; k++) { 732 boolean notlast = (k != n - 1); 733 if (k != m) { 734 p = H[k][k - 1]; 735 q = H[k + 1][k - 1]; 736 r = (notlast ? H[k + 2][k - 1] : 0.0); 737 x = Math.abs(p) + Math.abs(q) + Math.abs(r); 738 if (x != 0.0) { 739 p = p / x; 740 q = q / x; 741 r = r / x; 742 } 743 } 744 if (x == 0.0) { 745 break; 746 } 747 s = Math.sqrt(p * p + q * q + r * r); 748 if (p < 0) { 749 s = -s; 750 } 751 if (s != 0) { 752 if (k != m) { 753 H[k][k - 1] = -s * x; 754 } else if (l != m) { 755 H[k][k - 1] = -H[k][k - 1]; 756 } 757 p = p + s; 758 x = p / s; 759 y = q / s; 760 z = r / s; 761 q = q / p; 762 r = r / p; 763 764 // Row modification 765 766 for (int j = k; j < nn; j++) { 767 p = H[k][j] + q * H[k + 1][j]; 768 if (notlast) { 769 p = p + r * H[k + 2][j]; 770 H[k + 2][j] = H[k + 2][j] - p * z; 771 } 772 H[k][j] = H[k][j] - p * x; 773 H[k + 1][j] = H[k + 1][j] - p * y; 774 } 775 776 // Column modification 777 778 for (int i = 0; i <= Math.min(n, k + 3); i++) { 779 p = x * H[i][k] + y * H[i][k + 1]; 780 if (notlast) { 781 p = p + z * H[i][k + 2]; 782 H[i][k + 2] = H[i][k + 2] - p * r; 783 } 784 H[i][k] = H[i][k] - p; 785 H[i][k + 1] = H[i][k + 1] - p * q; 786 } 787 788 // Accumulate transformations 789 790 for (int i = low; i <= high; i++) { 791 p = x * V[i][k] + y * V[i][k + 1]; 792 if (notlast) { 793 p = p + z * V[i][k + 2]; 794 V[i][k + 2] = V[i][k + 2] - p * r; 795 } 796 V[i][k] = V[i][k] - p; 797 V[i][k + 1] = V[i][k + 1] - p * q; 798 } 799 } // (s != 0) 800 } // k loop 801 } // check convergence 802 } // while (n >= low) 803 804 // Backsubstitute to find vectors of upper triangular form 805 806 if (norm == 0.0) { 807 return; 808 } 809 810 for (n = nn - 1; n >= 0; n--) { 811 p = d[n]; 812 q = e[n]; 813 814 // Real vector 815 816 if (q == 0) { 817 int l = n; 818 H[n][n] = 1.0; 819 for (int i = n - 1; i >= 0; i--) { 820 w = H[i][i] - p; 821 r = 0.0; 822 for (int j = l; j <= n; j++) { 823 r = r + H[i][j] * H[j][n]; 824 } 825 if (e[i] < 0.0) { 826 z = w; 827 s = r; 828 } else { 829 l = i; 830 if (e[i] == 0.0) { 831 if (w != 0.0) { 832 H[i][n] = -r / w; 833 } else { 834 H[i][n] = -r / (EPSILON * norm); 835 } 836 837 // Solve real equations 838 839 } else { 840 x = H[i][i + 1]; 841 y = H[i + 1][i]; 842 q = (d[i] - p) * (d[i] - p) + e[i] * e[i]; 843 t = (x * s - z * r) / q; 844 H[i][n] = t; 845 if (Math.abs(x) > Math.abs(z)) { 846 H[i + 1][n] = (-r - w * t) / x; 847 } else { 848 H[i + 1][n] = (-s - y * t) / z; 849 } 850 } 851 852 // Overflow control 853 854 t = Math.abs(H[i][n]); 855 if ((EPSILON * t) * t > 1) { 856 for (int j = i; j <= n; j++) { 857 H[j][n] = H[j][n] / t; 858 } 859 } 860 } 861 } 862 863 // Complex vector 864 865 } else if (q < 0) { 866 int l = n - 1; 867 868 // Last vector component imaginary so matrix is triangular 869 870 if (Math.abs(H[n][n - 1]) > Math.abs(H[n - 1][n])) { 871 H[n - 1][n - 1] = q / H[n][n - 1]; 872 H[n - 1][n] = -(H[n][n] - p) / H[n][n - 1]; 873 } else { 874 cdiv(0.0, -H[n - 1][n], H[n - 1][n - 1] - p, q); 875 H[n - 1][n - 1] = cdivr; 876 H[n - 1][n] = cdivi; 877 } 878 H[n][n - 1] = 0.0; 879 H[n][n] = 1.0; 880 for (int i = n - 2; i >= 0; i--) { 881 double ra, sa, vr, vi; 882 ra = 0.0; 883 sa = 0.0; 884 for (int j = l; j <= n; j++) { 885 ra = ra + H[i][j] * H[j][n - 1]; 886 sa = sa + H[i][j] * H[j][n]; 887 } 888 w = H[i][i] - p; 889 890 if (e[i] < 0.0) { 891 z = w; 892 r = ra; 893 s = sa; 894 } else { 895 l = i; 896 if (e[i] == 0) { 897 cdiv(-ra, -sa, w, q); 898 H[i][n - 1] = cdivr; 899 H[i][n] = cdivi; 900 } else { 901 902 // Solve complex equations 903 904 x = H[i][i + 1]; 905 y = H[i + 1][i]; 906 vr = (d[i] - p) * (d[i] - p) + e[i] * e[i] - q * q; 907 vi = (d[i] - p) * 2.0 * q; 908 if (vr == 0.0 & vi == 0.0) { 909 vr = EPSILON 910 * norm 911 * (Math.abs(w) + Math.abs(q) + Math.abs(x) 912 + Math.abs(y) + Math.abs(z)); 913 } 914 cdiv(x * r - z * ra + q * sa, x * s - z * sa - q * ra, vr, vi); 915 H[i][n - 1] = cdivr; 916 H[i][n] = cdivi; 917 if (Math.abs(x) > (Math.abs(z) + Math.abs(q))) { 918 H[i + 1][n - 1] = (-ra - w * H[i][n - 1] + q * H[i][n]) / x; 919 H[i + 1][n] = (-sa - w * H[i][n] - q * H[i][n - 1]) / x; 920 } else { 921 cdiv(-r - y * H[i][n - 1], -s - y * H[i][n], z, q); 922 H[i + 1][n - 1] = cdivr; 923 H[i + 1][n] = cdivi; 924 } 925 } 926 927 // Overflow control 928 929 t = Math.max(Math.abs(H[i][n - 1]), Math.abs(H[i][n])); 930 if ((EPSILON * t) * t > 1) { 931 for (int j = i; j <= n; j++) { 932 H[j][n - 1] = H[j][n - 1] / t; 933 H[j][n] = H[j][n] / t; 934 } 935 } 936 } 937 } 938 } 939 } 940 941 // Vectors of isolated roots 942 943 for (int i = 0; i < nn; i++) { 944 if (i < low | i > high) { 945 for (int j = i; j < nn; j++) { 946 V[i][j] = H[i][j]; 947 } 948 } 949 } 950 951 // Back transformation to get eigenvectors of original matrix 952 953 for (int j = nn - 1; j >= low; j--) { 954 for (int i = low; i <= high; i++) { 955 z = 0.0; 956 for (int k = low; k <= Math.min(j, high); k++) { 957 z = z + V[i][k] * H[k][j]; 958 } 959 V[i][j] = z; 960 } 961 } 962 } 963 964 /* 965 * ------------------------ Constructor ------------------------ 966 */ 967 968 /** 969 * Check for symmetry, then construct the eigenvalue decomposition 970 * 971 * @param A 972 * Square matrix 973 * @return Structure to access D and V. 974 */ 975 976 public EigMatrix(Matrix Arg) { 977 final double[][] A = Arg.toDoubleArray(); 978 n = (int) Arg.getColumnCount(); 979 V = new double[n][n]; 980 d = new double[n]; 981 e = new double[n]; 982 H = new double[n][n]; 983 ort = new double[n]; 984 985 issymmetric = true; 986 for (int j = 0; (j < n) & issymmetric; j++) { 987 for (int i = 0; (i < n) & issymmetric; i++) { 988 issymmetric = (A[i][j] == A[j][i]); 989 } 990 } 991 992 if (issymmetric) { 993 for (int i = 0; i < n; i++) { 994 for (int j = 0; j < n; j++) { 995 V[i][j] = A[i][j]; 996 } 997 } 998 999 // Tridiagonalize. 1000 tred2(); 1001 1002 // Diagonalize. 1003 tql2(); 1004 1005 } else { 1006 1007 for (int j = 0; j < n; j++) { 1008 for (int i = 0; i < n; i++) { 1009 H[i][j] = A[i][j]; 1010 } 1011 } 1012 1013 // Reduce to Hessenberg form. 1014 orthes(); 1015 1016 // Reduce Hessenberg to real Schur form. 1017 hqr2(); 1018 } 1019 } 1020 1021 /* 1022 * ------------------------ Public Methods ------------------------ 1023 */ 1024 1025 /** 1026 * Return the eigenvector matrix 1027 * 1028 * @return V 1029 */ 1030 1031 public final Matrix getV() { 1032 return MatrixFactory.linkToArray(V); 1033 } 1034 1035 /** 1036 * Return the real parts of the eigenvalues 1037 * 1038 * @return real(diag(D)) 1039 */ 1040 1041 public final double[] getRealEigenvalues() { 1042 return d; 1043 } 1044 1045 /** 1046 * Return the imaginary parts of the eigenvalues 1047 * 1048 * @return imag(diag(D)) 1049 */ 1050 1051 public final double[] getImagEigenvalues() { 1052 return e; 1053 } 1054 1055 /** 1056 * Return the block diagonal eigenvalue matrix 1057 * 1058 * @return D 1059 */ 1060 1061 public final Matrix getD() { 1062 final double[][] D = new double[n][n]; 1063 for (int i = 0; i < n; i++) { 1064 for (int j = 0; j < n; j++) { 1065 D[i][j] = 0.0; 1066 } 1067 D[i][i] = d[i]; 1068 if (e[i] > 0) { 1069 D[i][i + 1] = e[i]; 1070 } else if (e[i] < 0) { 1071 D[i][i - 1] = e[i]; 1072 } 1073 } 1074 return MatrixFactory.linkToArray(D); 1075 } 1076 1077 } 1078 }