Search moodle.org's
Developer Documentation
Developer Documentation
See Release Notes
Long Term Support Release
- Bug fixes for general core bugs in 4.1.x will end 13 November 2023 (12 months).
- Bug fixes for security issues in 4.1.x will end 10 November 2025 (36 months).
- PHP version: minimum PHP 7.4.0 Note: minimum PHP version has increased since Moodle 4.0. PHP 8.0.x is supported too.
Differences Between: [Versions 310 and 401] [Versions 311 and 401] [Versions 39 and 401]
1 <?php 2 3 declare(strict_types=1); 4 5 /** 6 * Class to obtain eigenvalues and eigenvectors of a real matrix. 7 * 8 * If A is symmetric, then A = V*D*V' where the eigenvalue matrix D 9 * is diagonal and the eigenvector matrix V is orthogonal (i.e. 10 * A = V.times(D.times(V.transpose())) and V.times(V.transpose()) 11 * equals the identity matrix). 12 * 13 * If A is not symmetric, then the eigenvalue matrix D is block diagonal 14 * with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues, 15 * lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda]. The 16 * columns of V represent the eigenvectors in the sense that A*V = V*D, 17 * i.e. A.times(V) equals V.times(D). The matrix V may be badly 18 * conditioned, or even singular, so the validity of the equation 19 * A = V*D*inverse(V) depends upon V.cond(). 20 * 21 * @author Paul Meagher 22 * @license PHP v3.0 23 * 24 * @version 1.1 25 * 26 * Slightly changed to adapt the original code to PHP-ML library 27 * @date 2017/04/11 28 * 29 * @author Mustafa Karabulut 30 */ 31 32 namespace Phpml\Math\LinearAlgebra; 33 34 use Phpml\Math\Matrix; 35 36 class EigenvalueDecomposition 37 { 38 /** 39 * Row and column dimension (square matrix). 40 * 41 * @var int 42 */ 43 private $n; 44 45 /** 46 * Arrays for internal storage of eigenvalues. 47 * 48 * @var array 49 */ 50 private $d = []; 51 52 /** 53 * @var array 54 */ 55 private $e = []; 56 57 /** 58 * Array for internal storage of eigenvectors. 59 * 60 * @var array 61 */ 62 private $V = []; 63 64 /** 65 * Array for internal storage of nonsymmetric Hessenberg form. 66 * 67 * @var array 68 */ 69 private $H = []; 70 71 /** 72 * Working storage for nonsymmetric algorithm. 73 * 74 * @var array 75 */ 76 private $ort = []; 77 78 /** 79 * Used for complex scalar division. 80 * 81 * @var float 82 */ 83 private $cdivr; 84 85 /** 86 * @var float 87 */ 88 private $cdivi; 89 90 /** 91 * Constructor: Check for symmetry, then construct the eigenvalue decomposition 92 */ 93 public function __construct(array $arg) 94 { 95 $this->n = count($arg[0]); 96 $symmetric = true; 97 98 for ($j = 0; ($j < $this->n) & $symmetric; ++$j) { 99 for ($i = 0; ($i < $this->n) & $symmetric; ++$i) { 100 $symmetric = $arg[$i][$j] == $arg[$j][$i]; 101 } 102 } 103 104 if ($symmetric) { 105 $this->V = $arg; 106 // Tridiagonalize. 107 $this->tred2(); 108 // Diagonalize. 109 $this->tql2(); 110 } else { 111 $this->H = $arg; 112 $this->ort = []; 113 // Reduce to Hessenberg form. 114 $this->orthes(); 115 // Reduce Hessenberg to real Schur form. 116 $this->hqr2(); 117 } 118 } 119 120 /** 121 * Return the eigenvector matrix 122 */ 123 public function getEigenvectors(): array 124 { 125 $vectors = $this->V; 126 127 // Always return the eigenvectors of length 1.0 128 $vectors = new Matrix($vectors); 129 $vectors = array_map(function ($vect) { 130 $sum = 0; 131 $count = count($vect); 132 for ($i = 0; $i < $count; ++$i) { 133 $sum += $vect[$i] ** 2; 134 } 135 136 $sum **= .5; 137 for ($i = 0; $i < $count; ++$i) { 138 $vect[$i] /= $sum; 139 } 140 141 return $vect; 142 }, $vectors->transpose()->toArray()); 143 144 return $vectors; 145 } 146 147 /** 148 * Return the real parts of the eigenvalues<br> 149 * d = real(diag(D)); 150 */ 151 public function getRealEigenvalues(): array 152 { 153 return $this->d; 154 } 155 156 /** 157 * Return the imaginary parts of the eigenvalues <br> 158 * d = imag(diag(D)) 159 */ 160 public function getImagEigenvalues(): array 161 { 162 return $this->e; 163 } 164 165 /** 166 * Return the block diagonal eigenvalue matrix 167 */ 168 public function getDiagonalEigenvalues(): array 169 { 170 $D = []; 171 172 for ($i = 0; $i < $this->n; ++$i) { 173 $D[$i] = array_fill(0, $this->n, 0.0); 174 $D[$i][$i] = $this->d[$i]; 175 if ($this->e[$i] == 0) { 176 continue; 177 } 178 179 $o = $this->e[$i] > 0 ? $i + 1 : $i - 1; 180 $D[$i][$o] = $this->e[$i]; 181 } 182 183 return $D; 184 } 185 186 /** 187 * Symmetric Householder reduction to tridiagonal form. 188 */ 189 private function tred2(): void 190 { 191 // This is derived from the Algol procedures tred2 by 192 // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for 193 // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding 194 // Fortran subroutine in EISPACK. 195 $this->d = $this->V[$this->n - 1]; 196 // Householder reduction to tridiagonal form. 197 for ($i = $this->n - 1; $i > 0; --$i) { 198 $i_ = $i - 1; 199 // Scale to avoid under/overflow. 200 $h = $scale = 0.0; 201 $scale += array_sum(array_map('abs', $this->d)); 202 if ($scale == 0.0) { 203 $this->e[$i] = $this->d[$i_]; 204 $this->d = array_slice($this->V[$i_], 0, $this->n - 1); 205 for ($j = 0; $j < $i; ++$j) { 206 $this->V[$j][$i] = $this->V[$i][$j] = 0.0; 207 } 208 } else { 209 // Generate Householder vector. 210 for ($k = 0; $k < $i; ++$k) { 211 $this->d[$k] /= $scale; 212 $h += $this->d[$k] ** 2; 213 } 214 215 $f = $this->d[$i_]; 216 $g = $h ** .5; 217 if ($f > 0) { 218 $g = -$g; 219 } 220 221 $this->e[$i] = $scale * $g; 222 $h -= $f * $g; 223 $this->d[$i_] = $f - $g; 224 225 for ($j = 0; $j < $i; ++$j) { 226 $this->e[$j] = 0.0; 227 } 228 229 // Apply similarity transformation to remaining columns. 230 for ($j = 0; $j < $i; ++$j) { 231 $f = $this->d[$j]; 232 $this->V[$j][$i] = $f; 233 $g = $this->e[$j] + $this->V[$j][$j] * $f; 234 235 for ($k = $j + 1; $k <= $i_; ++$k) { 236 $g += $this->V[$k][$j] * $this->d[$k]; 237 $this->e[$k] += $this->V[$k][$j] * $f; 238 } 239 240 $this->e[$j] = $g; 241 } 242 243 $f = 0.0; 244 245 if ($h == 0.0) { 246 $h = 1e-32; 247 } 248 249 for ($j = 0; $j < $i; ++$j) { 250 $this->e[$j] /= $h; 251 $f += $this->e[$j] * $this->d[$j]; 252 } 253 254 $hh = $f / (2 * $h); 255 for ($j = 0; $j < $i; ++$j) { 256 $this->e[$j] -= $hh * $this->d[$j]; 257 } 258 259 for ($j = 0; $j < $i; ++$j) { 260 $f = $this->d[$j]; 261 $g = $this->e[$j]; 262 for ($k = $j; $k <= $i_; ++$k) { 263 $this->V[$k][$j] -= ($f * $this->e[$k] + $g * $this->d[$k]); 264 } 265 266 $this->d[$j] = $this->V[$i - 1][$j]; 267 $this->V[$i][$j] = 0.0; 268 } 269 } 270 271 $this->d[$i] = $h; 272 } 273 274 // Accumulate transformations. 275 for ($i = 0; $i < $this->n - 1; ++$i) { 276 $this->V[$this->n - 1][$i] = $this->V[$i][$i]; 277 $this->V[$i][$i] = 1.0; 278 $h = $this->d[$i + 1]; 279 if ($h != 0.0) { 280 for ($k = 0; $k <= $i; ++$k) { 281 $this->d[$k] = $this->V[$k][$i + 1] / $h; 282 } 283 284 for ($j = 0; $j <= $i; ++$j) { 285 $g = 0.0; 286 for ($k = 0; $k <= $i; ++$k) { 287 $g += $this->V[$k][$i + 1] * $this->V[$k][$j]; 288 } 289 290 for ($k = 0; $k <= $i; ++$k) { 291 $this->V[$k][$j] -= $g * $this->d[$k]; 292 } 293 } 294 } 295 296 for ($k = 0; $k <= $i; ++$k) { 297 $this->V[$k][$i + 1] = 0.0; 298 } 299 } 300 301 $this->d = $this->V[$this->n - 1]; 302 $this->V[$this->n - 1] = array_fill(0, $this->n, 0.0); 303 $this->V[$this->n - 1][$this->n - 1] = 1.0; 304 $this->e[0] = 0.0; 305 } 306 307 /** 308 * Symmetric tridiagonal QL algorithm. 309 * 310 * This is derived from the Algol procedures tql2, by 311 * Bowdler, Martin, Reinsch, and Wilkinson, Handbook for 312 * Auto. Comp., Vol.ii-Linear Algebra, and the corresponding 313 * Fortran subroutine in EISPACK. 314 */ 315 private function tql2(): void 316 { 317 for ($i = 1; $i < $this->n; ++$i) { 318 $this->e[$i - 1] = $this->e[$i]; 319 } 320 321 $this->e[$this->n - 1] = 0.0; 322 $f = 0.0; 323 $tst1 = 0.0; 324 $eps = 2.0 ** -52.0; 325 326 for ($l = 0; $l < $this->n; ++$l) { 327 // Find small subdiagonal element 328 $tst1 = max($tst1, abs($this->d[$l]) + abs($this->e[$l])); 329 $m = $l; 330 while ($m < $this->n) { 331 if (abs($this->e[$m]) <= $eps * $tst1) { 332 break; 333 } 334 335 ++$m; 336 } 337 338 // If m == l, $this->d[l] is an eigenvalue, 339 // otherwise, iterate. 340 if ($m > $l) { 341 do { 342 // Compute implicit shift 343 $g = $this->d[$l]; 344 $p = ($this->d[$l + 1] - $g) / (2.0 * $this->e[$l]); 345 $r = hypot($p, 1.0); 346 if ($p < 0) { 347 $r *= -1; 348 } 349 350 $this->d[$l] = $this->e[$l] / ($p + $r); 351 $this->d[$l + 1] = $this->e[$l] * ($p + $r); 352 $dl1 = $this->d[$l + 1]; 353 $h = $g - $this->d[$l]; 354 for ($i = $l + 2; $i < $this->n; ++$i) { 355 $this->d[$i] -= $h; 356 } 357 358 $f += $h; 359 // Implicit QL transformation. 360 $p = $this->d[$m]; 361 $c = 1.0; 362 $c2 = $c3 = $c; 363 $el1 = $this->e[$l + 1]; 364 $s = $s2 = 0.0; 365 for ($i = $m - 1; $i >= $l; --$i) { 366 $c3 = $c2; 367 $c2 = $c; 368 $s2 = $s; 369 $g = $c * $this->e[$i]; 370 $h = $c * $p; 371 $r = hypot($p, $this->e[$i]); 372 $this->e[$i + 1] = $s * $r; 373 $s = $this->e[$i] / $r; 374 $c = $p / $r; 375 $p = $c * $this->d[$i] - $s * $g; 376 $this->d[$i + 1] = $h + $s * ($c * $g + $s * $this->d[$i]); 377 // Accumulate transformation. 378 for ($k = 0; $k < $this->n; ++$k) { 379 $h = $this->V[$k][$i + 1]; 380 $this->V[$k][$i + 1] = $s * $this->V[$k][$i] + $c * $h; 381 $this->V[$k][$i] = $c * $this->V[$k][$i] - $s * $h; 382 } 383 } 384 385 $p = -$s * $s2 * $c3 * $el1 * $this->e[$l] / $dl1; 386 $this->e[$l] = $s * $p; 387 $this->d[$l] = $c * $p; 388 // Check for convergence. 389 } while (abs($this->e[$l]) > $eps * $tst1); 390 } 391 392 $this->d[$l] += $f; 393 $this->e[$l] = 0.0; 394 } 395 396 // Sort eigenvalues and corresponding vectors. 397 for ($i = 0; $i < $this->n - 1; ++$i) { 398 $k = $i; 399 $p = $this->d[$i]; 400 for ($j = $i + 1; $j < $this->n; ++$j) { 401 if ($this->d[$j] < $p) { 402 $k = $j; 403 $p = $this->d[$j]; 404 } 405 } 406 407 if ($k != $i) { 408 $this->d[$k] = $this->d[$i]; 409 $this->d[$i] = $p; 410 for ($j = 0; $j < $this->n; ++$j) { 411 $p = $this->V[$j][$i]; 412 $this->V[$j][$i] = $this->V[$j][$k]; 413 $this->V[$j][$k] = $p; 414 } 415 } 416 } 417 } 418 419 /** 420 * Nonsymmetric reduction to Hessenberg form. 421 * 422 * This is derived from the Algol procedures orthes and ortran, 423 * by Martin and Wilkinson, Handbook for Auto. Comp., 424 * Vol.ii-Linear Algebra, and the corresponding 425 * Fortran subroutines in EISPACK. 426 */ 427 private function orthes(): void 428 { 429 $low = 0; 430 $high = $this->n - 1; 431 432 for ($m = $low + 1; $m <= $high - 1; ++$m) { 433 // Scale column. 434 $scale = 0.0; 435 for ($i = $m; $i <= $high; ++$i) { 436 $scale += abs($this->H[$i][$m - 1]); 437 } 438 439 if ($scale != 0.0) { 440 // Compute Householder transformation. 441 $h = 0.0; 442 for ($i = $high; $i >= $m; --$i) { 443 $this->ort[$i] = $this->H[$i][$m - 1] / $scale; 444 $h += $this->ort[$i] * $this->ort[$i]; 445 } 446 447 $g = $h ** .5; 448 if ($this->ort[$m] > 0) { 449 $g *= -1; 450 } 451 452 $h -= $this->ort[$m] * $g; 453 $this->ort[$m] -= $g; 454 // Apply Householder similarity transformation 455 // H = (I -u * u' / h) * H * (I -u * u') / h) 456 for ($j = $m; $j < $this->n; ++$j) { 457 $f = 0.0; 458 for ($i = $high; $i >= $m; --$i) { 459 $f += $this->ort[$i] * $this->H[$i][$j]; 460 } 461 462 $f /= $h; 463 for ($i = $m; $i <= $high; ++$i) { 464 $this->H[$i][$j] -= $f * $this->ort[$i]; 465 } 466 } 467 468 for ($i = 0; $i <= $high; ++$i) { 469 $f = 0.0; 470 for ($j = $high; $j >= $m; --$j) { 471 $f += $this->ort[$j] * $this->H[$i][$j]; 472 } 473 474 $f /= $h; 475 for ($j = $m; $j <= $high; ++$j) { 476 $this->H[$i][$j] -= $f * $this->ort[$j]; 477 } 478 } 479 480 $this->ort[$m] = $scale * $this->ort[$m]; 481 $this->H[$m][$m - 1] = $scale * $g; 482 } 483 } 484 485 // Accumulate transformations (Algol's ortran). 486 for ($i = 0; $i < $this->n; ++$i) { 487 for ($j = 0; $j < $this->n; ++$j) { 488 $this->V[$i][$j] = ($i == $j ? 1.0 : 0.0); 489 } 490 } 491 492 for ($m = $high - 1; $m >= $low + 1; --$m) { 493 if ($this->H[$m][$m - 1] != 0.0) { 494 for ($i = $m + 1; $i <= $high; ++$i) { 495 $this->ort[$i] = $this->H[$i][$m - 1]; 496 } 497 498 for ($j = $m; $j <= $high; ++$j) { 499 $g = 0.0; 500 for ($i = $m; $i <= $high; ++$i) { 501 $g += $this->ort[$i] * $this->V[$i][$j]; 502 } 503 504 // Double division avoids possible underflow 505 $g /= $this->ort[$m]; 506 $g /= $this->H[$m][$m - 1]; 507 for ($i = $m; $i <= $high; ++$i) { 508 $this->V[$i][$j] += $g * $this->ort[$i]; 509 } 510 } 511 } 512 } 513 } 514 515 /** 516 * Performs complex division. 517 * 518 * @param int|float $xr 519 * @param int|float $xi 520 * @param int|float $yr 521 * @param int|float $yi 522 */ 523 private function cdiv($xr, $xi, $yr, $yi): void 524 { 525 if (abs($yr) > abs($yi)) { 526 $r = $yi / $yr; 527 $d = $yr + $r * $yi; 528 $this->cdivr = ($xr + $r * $xi) / $d; 529 $this->cdivi = ($xi - $r * $xr) / $d; 530 } else { 531 $r = $yr / $yi; 532 $d = $yi + $r * $yr; 533 $this->cdivr = ($r * $xr + $xi) / $d; 534 $this->cdivi = ($r * $xi - $xr) / $d; 535 } 536 } 537 538 /** 539 * Nonsymmetric reduction from Hessenberg to real Schur form. 540 * 541 * Code is derived from the Algol procedure hqr2, 542 * by Martin and Wilkinson, Handbook for Auto. Comp., 543 * Vol.ii-Linear Algebra, and the corresponding 544 * Fortran subroutine in EISPACK. 545 */ 546 private function hqr2(): void 547 { 548 // Initialize 549 $nn = $this->n; 550 $n = $nn - 1; 551 $low = 0; 552 $high = $nn - 1; 553 $eps = 2.0 ** -52.0; 554 $exshift = 0.0; 555 $p = $q = $r = $s = $z = 0; 556 // Store roots isolated by balanc and compute matrix norm 557 $norm = 0.0; 558 559 for ($i = 0; $i < $nn; ++$i) { 560 if (($i < $low) or ($i > $high)) { 561 $this->d[$i] = $this->H[$i][$i]; 562 $this->e[$i] = 0.0; 563 } 564 565 for ($j = max($i - 1, 0); $j < $nn; ++$j) { 566 $norm += abs($this->H[$i][$j]); 567 } 568 } 569 570 // Outer loop over eigenvalue index 571 $iter = 0; 572 while ($n >= $low) { 573 // Look for single small sub-diagonal element 574 $l = $n; 575 while ($l > $low) { 576 $s = abs($this->H[$l - 1][$l - 1]) + abs($this->H[$l][$l]); 577 if ($s == 0.0) { 578 $s = $norm; 579 } 580 581 if (abs($this->H[$l][$l - 1]) < $eps * $s) { 582 break; 583 } 584 585 --$l; 586 } 587 588 // Check for convergence 589 // One root found 590 if ($l == $n) { 591 $this->H[$n][$n] += $exshift; 592 $this->d[$n] = $this->H[$n][$n]; 593 $this->e[$n] = 0.0; 594 --$n; 595 $iter = 0; 596 // Two roots found 597 } elseif ($l == $n - 1) { 598 $w = $this->H[$n][$n - 1] * $this->H[$n - 1][$n]; 599 $p = ($this->H[$n - 1][$n - 1] - $this->H[$n][$n]) / 2.0; 600 $q = $p * $p + $w; 601 $z = abs($q) ** .5; 602 $this->H[$n][$n] += $exshift; 603 $this->H[$n - 1][$n - 1] += $exshift; 604 $x = $this->H[$n][$n]; 605 // Real pair 606 if ($q >= 0) { 607 if ($p >= 0) { 608 $z = $p + $z; 609 } else { 610 $z = $p - $z; 611 } 612 613 $this->d[$n - 1] = $x + $z; 614 $this->d[$n] = $this->d[$n - 1]; 615 if ($z != 0.0) { 616 $this->d[$n] = $x - $w / $z; 617 } 618 619 $this->e[$n - 1] = 0.0; 620 $this->e[$n] = 0.0; 621 $x = $this->H[$n][$n - 1]; 622 $s = abs($x) + abs($z); 623 $p = $x / $s; 624 $q = $z / $s; 625 $r = ($p * $p + $q * $q) ** .5; 626 $p /= $r; 627 $q /= $r; 628 // Row modification 629 for ($j = $n - 1; $j < $nn; ++$j) { 630 $z = $this->H[$n - 1][$j]; 631 $this->H[$n - 1][$j] = $q * $z + $p * $this->H[$n][$j]; 632 $this->H[$n][$j] = $q * $this->H[$n][$j] - $p * $z; 633 } 634 635 // Column modification 636 for ($i = 0; $i <= $n; ++$i) { 637 $z = $this->H[$i][$n - 1]; 638 $this->H[$i][$n - 1] = $q * $z + $p * $this->H[$i][$n]; 639 $this->H[$i][$n] = $q * $this->H[$i][$n] - $p * $z; 640 } 641 642 // Accumulate transformations 643 for ($i = $low; $i <= $high; ++$i) { 644 $z = $this->V[$i][$n - 1]; 645 $this->V[$i][$n - 1] = $q * $z + $p * $this->V[$i][$n]; 646 $this->V[$i][$n] = $q * $this->V[$i][$n] - $p * $z; 647 } 648 649 // Complex pair 650 } else { 651 $this->d[$n - 1] = $x + $p; 652 $this->d[$n] = $x + $p; 653 $this->e[$n - 1] = $z; 654 $this->e[$n] = -$z; 655 } 656 657 $n -= 2; 658 $iter = 0; 659 // No convergence yet 660 } else { 661 // Form shift 662 $x = $this->H[$n][$n]; 663 $y = 0.0; 664 $w = 0.0; 665 if ($l < $n) { 666 $y = $this->H[$n - 1][$n - 1]; 667 $w = $this->H[$n][$n - 1] * $this->H[$n - 1][$n]; 668 } 669 670 // Wilkinson's original ad hoc shift 671 if ($iter == 10) { 672 $exshift += $x; 673 for ($i = $low; $i <= $n; ++$i) { 674 $this->H[$i][$i] -= $x; 675 } 676 677 $s = abs($this->H[$n][$n - 1]) + abs($this->H[$n - 1][$n - 2]); 678 $x = $y = 0.75 * $s; 679 $w = -0.4375 * $s * $s; 680 } 681 682 // MATLAB's new ad hoc shift 683 if ($iter == 30) { 684 $s = ($y - $x) / 2.0; 685 $s *= $s + $w; 686 if ($s > 0) { 687 $s **= .5; 688 if ($y < $x) { 689 $s = -$s; 690 } 691 692 $s = $x - $w / (($y - $x) / 2.0 + $s); 693 for ($i = $low; $i <= $n; ++$i) { 694 $this->H[$i][$i] -= $s; 695 } 696 697 $exshift += $s; 698 $x = $y = $w = 0.964; 699 } 700 } 701 702 // Could check iteration count here. 703 ++$iter; 704 // Look for two consecutive small sub-diagonal elements 705 $m = $n - 2; 706 while ($m >= $l) { 707 $z = $this->H[$m][$m]; 708 $r = $x - $z; 709 $s = $y - $z; 710 $p = ($r * $s - $w) / $this->H[$m + 1][$m] + $this->H[$m][$m + 1]; 711 $q = $this->H[$m + 1][$m + 1] - $z - $r - $s; 712 $r = $this->H[$m + 2][$m + 1]; 713 $s = abs($p) + abs($q) + abs($r); 714 $p /= $s; 715 $q /= $s; 716 $r /= $s; 717 if ($m == $l) { 718 break; 719 } 720 721 if (abs($this->H[$m][$m - 1]) * (abs($q) + abs($r)) < 722 $eps * (abs($p) * (abs($this->H[$m - 1][$m - 1]) + abs($z) + abs($this->H[$m + 1][$m + 1])))) { 723 break; 724 } 725 726 --$m; 727 } 728 729 for ($i = $m + 2; $i <= $n; ++$i) { 730 $this->H[$i][$i - 2] = 0.0; 731 if ($i > $m + 2) { 732 $this->H[$i][$i - 3] = 0.0; 733 } 734 } 735 736 // Double QR step involving rows l:n and columns m:n 737 for ($k = $m; $k <= $n - 1; ++$k) { 738 $notlast = $k != $n - 1; 739 if ($k != $m) { 740 $p = $this->H[$k][$k - 1]; 741 $q = $this->H[$k + 1][$k - 1]; 742 $r = ($notlast ? $this->H[$k + 2][$k - 1] : 0.0); 743 $x = abs($p) + abs($q) + abs($r); 744 if ($x != 0.0) { 745 $p /= $x; 746 $q /= $x; 747 $r /= $x; 748 } 749 } 750 751 if ($x == 0.0) { 752 break; 753 } 754 755 $s = ($p * $p + $q * $q + $r * $r) ** .5; 756 if ($p < 0) { 757 $s = -$s; 758 } 759 760 if ($s != 0) { 761 if ($k != $m) { 762 $this->H[$k][$k - 1] = -$s * $x; 763 } elseif ($l != $m) { 764 $this->H[$k][$k - 1] = -$this->H[$k][$k - 1]; 765 } 766 767 $p += $s; 768 $x = $p / $s; 769 $y = $q / $s; 770 $z = $r / $s; 771 $q /= $p; 772 $r /= $p; 773 // Row modification 774 for ($j = $k; $j < $nn; ++$j) { 775 $p = $this->H[$k][$j] + $q * $this->H[$k + 1][$j]; 776 if ($notlast) { 777 $p += $r * $this->H[$k + 2][$j]; 778 $this->H[$k + 2][$j] -= $p * $z; 779 } 780 781 $this->H[$k][$j] -= $p * $x; 782 $this->H[$k + 1][$j] -= $p * $y; 783 } 784 785 // Column modification 786 for ($i = 0; $i <= min($n, $k + 3); ++$i) { 787 $p = $x * $this->H[$i][$k] + $y * $this->H[$i][$k + 1]; 788 if ($notlast) { 789 $p += $z * $this->H[$i][$k + 2]; 790 $this->H[$i][$k + 2] -= $p * $r; 791 } 792 793 $this->H[$i][$k] -= $p; 794 $this->H[$i][$k + 1] -= $p * $q; 795 } 796 797 // Accumulate transformations 798 for ($i = $low; $i <= $high; ++$i) { 799 $p = $x * $this->V[$i][$k] + $y * $this->V[$i][$k + 1]; 800 if ($notlast) { 801 $p += $z * $this->V[$i][$k + 2]; 802 $this->V[$i][$k + 2] -= $p * $r; 803 } 804 805 $this->V[$i][$k] -= $p; 806 $this->V[$i][$k + 1] -= $p * $q; 807 } 808 } // ($s != 0) 809 } // k loop 810 } // check convergence 811 } // while ($n >= $low) 812 813 // Backsubstitute to find vectors of upper triangular form 814 if ($norm == 0.0) { 815 return; 816 } 817 818 for ($n = $nn - 1; $n >= 0; --$n) { 819 $p = $this->d[$n]; 820 $q = $this->e[$n]; 821 // Real vector 822 if ($q == 0) { 823 $l = $n; 824 $this->H[$n][$n] = 1.0; 825 for ($i = $n - 1; $i >= 0; --$i) { 826 $w = $this->H[$i][$i] - $p; 827 $r = 0.0; 828 for ($j = $l; $j <= $n; ++$j) { 829 $r += $this->H[$i][$j] * $this->H[$j][$n]; 830 } 831 832 if ($this->e[$i] < 0.0) { 833 $z = $w; 834 $s = $r; 835 } else { 836 $l = $i; 837 if ($this->e[$i] == 0.0) { 838 if ($w != 0.0) { 839 $this->H[$i][$n] = -$r / $w; 840 } else { 841 $this->H[$i][$n] = -$r / ($eps * $norm); 842 } 843 844 // Solve real equations 845 } else { 846 $x = $this->H[$i][$i + 1]; 847 $y = $this->H[$i + 1][$i]; 848 $q = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i]; 849 $t = ($x * $s - $z * $r) / $q; 850 $this->H[$i][$n] = $t; 851 if (abs($x) > abs($z)) { 852 $this->H[$i + 1][$n] = (-$r - $w * $t) / $x; 853 } else { 854 $this->H[$i + 1][$n] = (-$s - $y * $t) / $z; 855 } 856 } 857 858 // Overflow control 859 $t = abs($this->H[$i][$n]); 860 if (($eps * $t) * $t > 1) { 861 for ($j = $i; $j <= $n; ++$j) { 862 $this->H[$j][$n] /= $t; 863 } 864 } 865 } 866 } 867 868 // Complex vector 869 } elseif ($q < 0) { 870 $l = $n - 1; 871 // Last vector component imaginary so matrix is triangular 872 if (abs($this->H[$n][$n - 1]) > abs($this->H[$n - 1][$n])) { 873 $this->H[$n - 1][$n - 1] = $q / $this->H[$n][$n - 1]; 874 $this->H[$n - 1][$n] = -($this->H[$n][$n] - $p) / $this->H[$n][$n - 1]; 875 } else { 876 $this->cdiv(0.0, -$this->H[$n - 1][$n], $this->H[$n - 1][$n - 1] - $p, $q); 877 $this->H[$n - 1][$n - 1] = $this->cdivr; 878 $this->H[$n - 1][$n] = $this->cdivi; 879 } 880 881 $this->H[$n][$n - 1] = 0.0; 882 $this->H[$n][$n] = 1.0; 883 for ($i = $n - 2; $i >= 0; --$i) { 884 // double ra,sa,vr,vi; 885 $ra = 0.0; 886 $sa = 0.0; 887 for ($j = $l; $j <= $n; ++$j) { 888 $ra += $this->H[$i][$j] * $this->H[$j][$n - 1]; 889 $sa += $this->H[$i][$j] * $this->H[$j][$n]; 890 } 891 892 $w = $this->H[$i][$i] - $p; 893 if ($this->e[$i] < 0.0) { 894 $z = $w; 895 $r = $ra; 896 $s = $sa; 897 } else { 898 $l = $i; 899 if ($this->e[$i] == 0) { 900 $this->cdiv(-$ra, -$sa, $w, $q); 901 $this->H[$i][$n - 1] = $this->cdivr; 902 $this->H[$i][$n] = $this->cdivi; 903 } else { 904 // Solve complex equations 905 $x = $this->H[$i][$i + 1]; 906 $y = $this->H[$i + 1][$i]; 907 $vr = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i] - $q * $q; 908 $vi = ($this->d[$i] - $p) * 2.0 * $q; 909 if ($vr == 0.0 && $vi == 0.0) { 910 $vr = $eps * $norm * (abs($w) + abs($q) + abs($x) + abs($y) + abs($z)); 911 } 912 913 $this->cdiv($x * $r - $z * $ra + $q * $sa, $x * $s - $z * $sa - $q * $ra, $vr, $vi); 914 $this->H[$i][$n - 1] = $this->cdivr; 915 $this->H[$i][$n] = $this->cdivi; 916 if (abs($x) > (abs($z) + abs($q))) { 917 $this->H[$i + 1][$n - 1] = (-$ra - $w * $this->H[$i][$n - 1] + $q * $this->H[$i][$n]) / $x; 918 $this->H[$i + 1][$n] = (-$sa - $w * $this->H[$i][$n] - $q * $this->H[$i][$n - 1]) / $x; 919 } else { 920 $this->cdiv(-$r - $y * $this->H[$i][$n - 1], -$s - $y * $this->H[$i][$n], $z, $q); 921 $this->H[$i + 1][$n - 1] = $this->cdivr; 922 $this->H[$i + 1][$n] = $this->cdivi; 923 } 924 } 925 926 // Overflow control 927 $t = max(abs($this->H[$i][$n - 1]), abs($this->H[$i][$n])); 928 if (($eps * $t) * $t > 1) { 929 for ($j = $i; $j <= $n; ++$j) { 930 $this->H[$j][$n - 1] /= $t; 931 $this->H[$j][$n] /= $t; 932 } 933 } 934 } // end else 935 } // end for 936 } // end else for complex case 937 } // end for 938 939 // Vectors of isolated roots 940 for ($i = 0; $i < $nn; ++$i) { 941 if ($i < $low || $i > $high) { 942 for ($j = $i; $j < $nn; ++$j) { 943 $this->V[$i][$j] = $this->H[$i][$j]; 944 } 945 } 946 } 947 948 // Back transformation to get eigenvectors of original matrix 949 for ($j = $nn - 1; $j >= $low; --$j) { 950 for ($i = $low; $i <= $high; ++$i) { 951 $z = 0.0; 952 for ($k = $low; $k <= min($j, $high); ++$k) { 953 $z += $this->V[$i][$k] * $this->H[$k][$j]; 954 } 955 956 $this->V[$i][$j] = $z; 957 } 958 } 959 } 960 }
title
Description
Body
title
Description
Body
title
Description
Body
title
Body
