Search moodle.org's
Developer Documentation
Developer Documentation
- Bug fixes for general core bugs in 3.11.x will end 14 Nov 2022 (12 months plus 6 months extension).
- Bug fixes for security issues in 3.11.x will end 13 Nov 2023 (18 months plus 12 months extension).
- PHP version: minimum PHP 7.3.0 Note: minimum PHP version has increased since Moodle 3.10. PHP 7.4.x is supported too.
Differences Between: [Versions 311 and 400] [Versions 311 and 401] [Versions 311 and 402] [Versions 311 and 403]
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 = ($g / $this->ort[$m]) / $this->H[$m][$m - 1]; 506 for ($i = $m; $i <= $high; ++$i) { 507 $this->V[$i][$j] += $g * $this->ort[$i]; 508 } 509 } 510 } 511 } 512 } 513 514 /** 515 * Performs complex division. 516 * 517 * @param int|float $xr 518 * @param int|float $xi 519 * @param int|float $yr 520 * @param int|float $yi 521 */ 522 private function cdiv($xr, $xi, $yr, $yi): void 523 { 524 if (abs($yr) > abs($yi)) { 525 $r = $yi / $yr; 526 $d = $yr + $r * $yi; 527 $this->cdivr = ($xr + $r * $xi) / $d; 528 $this->cdivi = ($xi - $r * $xr) / $d; 529 } else { 530 $r = $yr / $yi; 531 $d = $yi + $r * $yr; 532 $this->cdivr = ($r * $xr + $xi) / $d; 533 $this->cdivi = ($r * $xi - $xr) / $d; 534 } 535 } 536 537 /** 538 * Nonsymmetric reduction from Hessenberg to real Schur form. 539 * 540 * Code is derived from the Algol procedure hqr2, 541 * by Martin and Wilkinson, Handbook for Auto. Comp., 542 * Vol.ii-Linear Algebra, and the corresponding 543 * Fortran subroutine in EISPACK. 544 */ 545 private function hqr2(): void 546 { 547 // Initialize 548 $nn = $this->n; 549 $n = $nn - 1; 550 $low = 0; 551 $high = $nn - 1; 552 $eps = 2.0 ** -52.0; 553 $exshift = 0.0; 554 $p = $q = $r = $s = $z = 0; 555 // Store roots isolated by balanc and compute matrix norm 556 $norm = 0.0; 557 558 for ($i = 0; $i < $nn; ++$i) { 559 if (($i < $low) or ($i > $high)) { 560 $this->d[$i] = $this->H[$i][$i]; 561 $this->e[$i] = 0.0; 562 } 563 564 for ($j = max($i - 1, 0); $j < $nn; ++$j) { 565 $norm += abs($this->H[$i][$j]); 566 } 567 } 568 569 // Outer loop over eigenvalue index 570 $iter = 0; 571 while ($n >= $low) { 572 // Look for single small sub-diagonal element 573 $l = $n; 574 while ($l > $low) { 575 $s = abs($this->H[$l - 1][$l - 1]) + abs($this->H[$l][$l]); 576 if ($s == 0.0) { 577 $s = $norm; 578 } 579 580 if (abs($this->H[$l][$l - 1]) < $eps * $s) { 581 break; 582 } 583 584 --$l; 585 } 586 587 // Check for convergence 588 // One root found 589 if ($l == $n) { 590 $this->H[$n][$n] += $exshift; 591 $this->d[$n] = $this->H[$n][$n]; 592 $this->e[$n] = 0.0; 593 --$n; 594 $iter = 0; 595 // Two roots found 596 } elseif ($l == $n - 1) { 597 $w = $this->H[$n][$n - 1] * $this->H[$n - 1][$n]; 598 $p = ($this->H[$n - 1][$n - 1] - $this->H[$n][$n]) / 2.0; 599 $q = $p * $p + $w; 600 $z = abs($q) ** .5; 601 $this->H[$n][$n] += $exshift; 602 $this->H[$n - 1][$n - 1] += $exshift; 603 $x = $this->H[$n][$n]; 604 // Real pair 605 if ($q >= 0) { 606 if ($p >= 0) { 607 $z = $p + $z; 608 } else { 609 $z = $p - $z; 610 } 611 612 $this->d[$n - 1] = $x + $z; 613 $this->d[$n] = $this->d[$n - 1]; 614 if ($z != 0.0) { 615 $this->d[$n] = $x - $w / $z; 616 } 617 618 $this->e[$n - 1] = 0.0; 619 $this->e[$n] = 0.0; 620 $x = $this->H[$n][$n - 1]; 621 $s = abs($x) + abs($z); 622 $p = $x / $s; 623 $q = $z / $s; 624 $r = ($p * $p + $q * $q) ** .5; 625 $p /= $r; 626 $q /= $r; 627 // Row modification 628 for ($j = $n - 1; $j < $nn; ++$j) { 629 $z = $this->H[$n - 1][$j]; 630 $this->H[$n - 1][$j] = $q * $z + $p * $this->H[$n][$j]; 631 $this->H[$n][$j] = $q * $this->H[$n][$j] - $p * $z; 632 } 633 634 // Column modification 635 for ($i = 0; $i <= $n; ++$i) { 636 $z = $this->H[$i][$n - 1]; 637 $this->H[$i][$n - 1] = $q * $z + $p * $this->H[$i][$n]; 638 $this->H[$i][$n] = $q * $this->H[$i][$n] - $p * $z; 639 } 640 641 // Accumulate transformations 642 for ($i = $low; $i <= $high; ++$i) { 643 $z = $this->V[$i][$n - 1]; 644 $this->V[$i][$n - 1] = $q * $z + $p * $this->V[$i][$n]; 645 $this->V[$i][$n] = $q * $this->V[$i][$n] - $p * $z; 646 } 647 648 // Complex pair 649 } else { 650 $this->d[$n - 1] = $x + $p; 651 $this->d[$n] = $x + $p; 652 $this->e[$n - 1] = $z; 653 $this->e[$n] = -$z; 654 } 655 656 $n -= 2; 657 $iter = 0; 658 // No convergence yet 659 } else { 660 // Form shift 661 $x = $this->H[$n][$n]; 662 $y = 0.0; 663 $w = 0.0; 664 if ($l < $n) { 665 $y = $this->H[$n - 1][$n - 1]; 666 $w = $this->H[$n][$n - 1] * $this->H[$n - 1][$n]; 667 } 668 669 // Wilkinson's original ad hoc shift 670 if ($iter == 10) { 671 $exshift += $x; 672 for ($i = $low; $i <= $n; ++$i) { 673 $this->H[$i][$i] -= $x; 674 } 675 676 $s = abs($this->H[$n][$n - 1]) + abs($this->H[$n - 1][$n - 2]); 677 $x = $y = 0.75 * $s; 678 $w = -0.4375 * $s * $s; 679 } 680 681 // MATLAB's new ad hoc shift 682 if ($iter == 30) { 683 $s = ($y - $x) / 2.0; 684 $s *= $s + $w; 685 if ($s > 0) { 686 $s **= .5; 687 if ($y < $x) { 688 $s = -$s; 689 } 690 691 $s = $x - $w / (($y - $x) / 2.0 + $s); 692 for ($i = $low; $i <= $n; ++$i) { 693 $this->H[$i][$i] -= $s; 694 } 695 696 $exshift += $s; 697 $x = $y = $w = 0.964; 698 } 699 } 700 701 // Could check iteration count here. 702 ++$iter; 703 // Look for two consecutive small sub-diagonal elements 704 $m = $n - 2; 705 while ($m >= $l) { 706 $z = $this->H[$m][$m]; 707 $r = $x - $z; 708 $s = $y - $z; 709 $p = ($r * $s - $w) / $this->H[$m + 1][$m] + $this->H[$m][$m + 1]; 710 $q = $this->H[$m + 1][$m + 1] - $z - $r - $s; 711 $r = $this->H[$m + 2][$m + 1]; 712 $s = abs($p) + abs($q) + abs($r); 713 $p /= $s; 714 $q /= $s; 715 $r /= $s; 716 if ($m == $l) { 717 break; 718 } 719 720 if (abs($this->H[$m][$m - 1]) * (abs($q) + abs($r)) < 721 $eps * (abs($p) * (abs($this->H[$m - 1][$m - 1]) + abs($z) + abs($this->H[$m + 1][$m + 1])))) { 722 break; 723 } 724 725 --$m; 726 } 727 728 for ($i = $m + 2; $i <= $n; ++$i) { 729 $this->H[$i][$i - 2] = 0.0; 730 if ($i > $m + 2) { 731 $this->H[$i][$i - 3] = 0.0; 732 } 733 } 734 735 // Double QR step involving rows l:n and columns m:n 736 for ($k = $m; $k <= $n - 1; ++$k) { 737 $notlast = ($k != $n - 1); 738 if ($k != $m) { 739 $p = $this->H[$k][$k - 1]; 740 $q = $this->H[$k + 1][$k - 1]; 741 $r = ($notlast ? $this->H[$k + 2][$k - 1] : 0.0); 742 $x = abs($p) + abs($q) + abs($r); 743 if ($x != 0.0) { 744 $p /= $x; 745 $q /= $x; 746 $r /= $x; 747 } 748 } 749 750 if ($x == 0.0) { 751 break; 752 } 753 754 $s = ($p * $p + $q * $q + $r * $r) ** .5; 755 if ($p < 0) { 756 $s = -$s; 757 } 758 759 if ($s != 0) { 760 if ($k != $m) { 761 $this->H[$k][$k - 1] = -$s * $x; 762 } elseif ($l != $m) { 763 $this->H[$k][$k - 1] = -$this->H[$k][$k - 1]; 764 } 765 766 $p += $s; 767 $x = $p / $s; 768 $y = $q / $s; 769 $z = $r / $s; 770 $q /= $p; 771 $r /= $p; 772 // Row modification 773 for ($j = $k; $j < $nn; ++$j) { 774 $p = $this->H[$k][$j] + $q * $this->H[$k + 1][$j]; 775 if ($notlast) { 776 $p += $r * $this->H[$k + 2][$j]; 777 $this->H[$k + 2][$j] -= $p * $z; 778 } 779 780 $this->H[$k][$j] -= $p * $x; 781 $this->H[$k + 1][$j] -= $p * $y; 782 } 783 784 // Column modification 785 for ($i = 0; $i <= min($n, $k + 3); ++$i) { 786 $p = $x * $this->H[$i][$k] + $y * $this->H[$i][$k + 1]; 787 if ($notlast) { 788 $p += $z * $this->H[$i][$k + 2]; 789 $this->H[$i][$k + 2] -= $p * $r; 790 } 791 792 $this->H[$i][$k] -= $p; 793 $this->H[$i][$k + 1] -= $p * $q; 794 } 795 796 // Accumulate transformations 797 for ($i = $low; $i <= $high; ++$i) { 798 $p = $x * $this->V[$i][$k] + $y * $this->V[$i][$k + 1]; 799 if ($notlast) { 800 $p += $z * $this->V[$i][$k + 2]; 801 $this->V[$i][$k + 2] -= $p * $r; 802 } 803 804 $this->V[$i][$k] -= $p; 805 $this->V[$i][$k + 1] -= $p * $q; 806 } 807 } // ($s != 0) 808 } // k loop 809 } // check convergence 810 } // while ($n >= $low) 811 812 // Backsubstitute to find vectors of upper triangular form 813 if ($norm == 0.0) { 814 return; 815 } 816 817 for ($n = $nn - 1; $n >= 0; --$n) { 818 $p = $this->d[$n]; 819 $q = $this->e[$n]; 820 // Real vector 821 if ($q == 0) { 822 $l = $n; 823 $this->H[$n][$n] = 1.0; 824 for ($i = $n - 1; $i >= 0; --$i) { 825 $w = $this->H[$i][$i] - $p; 826 $r = 0.0; 827 for ($j = $l; $j <= $n; ++$j) { 828 $r += $this->H[$i][$j] * $this->H[$j][$n]; 829 } 830 831 if ($this->e[$i] < 0.0) { 832 $z = $w; 833 $s = $r; 834 } else { 835 $l = $i; 836 if ($this->e[$i] == 0.0) { 837 if ($w != 0.0) { 838 $this->H[$i][$n] = -$r / $w; 839 } else { 840 $this->H[$i][$n] = -$r / ($eps * $norm); 841 } 842 843 // Solve real equations 844 } else { 845 $x = $this->H[$i][$i + 1]; 846 $y = $this->H[$i + 1][$i]; 847 $q = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i]; 848 $t = ($x * $s - $z * $r) / $q; 849 $this->H[$i][$n] = $t; 850 if (abs($x) > abs($z)) { 851 $this->H[$i + 1][$n] = (-$r - $w * $t) / $x; 852 } else { 853 $this->H[$i + 1][$n] = (-$s - $y * $t) / $z; 854 } 855 } 856 857 // Overflow control 858 $t = abs($this->H[$i][$n]); 859 if (($eps * $t) * $t > 1) { 860 for ($j = $i; $j <= $n; ++$j) { 861 $this->H[$j][$n] /= $t; 862 } 863 } 864 } 865 } 866 867 // Complex vector 868 } elseif ($q < 0) { 869 $l = $n - 1; 870 // Last vector component imaginary so matrix is triangular 871 if (abs($this->H[$n][$n - 1]) > abs($this->H[$n - 1][$n])) { 872 $this->H[$n - 1][$n - 1] = $q / $this->H[$n][$n - 1]; 873 $this->H[$n - 1][$n] = -($this->H[$n][$n] - $p) / $this->H[$n][$n - 1]; 874 } else { 875 $this->cdiv(0.0, -$this->H[$n - 1][$n], $this->H[$n - 1][$n - 1] - $p, $q); 876 $this->H[$n - 1][$n - 1] = $this->cdivr; 877 $this->H[$n - 1][$n] = $this->cdivi; 878 } 879 880 $this->H[$n][$n - 1] = 0.0; 881 $this->H[$n][$n] = 1.0; 882 for ($i = $n - 2; $i >= 0; --$i) { 883 // double ra,sa,vr,vi; 884 $ra = 0.0; 885 $sa = 0.0; 886 for ($j = $l; $j <= $n; ++$j) { 887 $ra += $this->H[$i][$j] * $this->H[$j][$n - 1]; 888 $sa += $this->H[$i][$j] * $this->H[$j][$n]; 889 } 890 891 $w = $this->H[$i][$i] - $p; 892 if ($this->e[$i] < 0.0) { 893 $z = $w; 894 $r = $ra; 895 $s = $sa; 896 } else { 897 $l = $i; 898 if ($this->e[$i] == 0) { 899 $this->cdiv(-$ra, -$sa, $w, $q); 900 $this->H[$i][$n - 1] = $this->cdivr; 901 $this->H[$i][$n] = $this->cdivi; 902 } else { 903 // Solve complex equations 904 $x = $this->H[$i][$i + 1]; 905 $y = $this->H[$i + 1][$i]; 906 $vr = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i] - $q * $q; 907 $vi = ($this->d[$i] - $p) * 2.0 * $q; 908 if ($vr == 0.0 && $vi == 0.0) { 909 $vr = $eps * $norm * (abs($w) + abs($q) + abs($x) + abs($y) + abs($z)); 910 } 911 912 $this->cdiv($x * $r - $z * $ra + $q * $sa, $x * $s - $z * $sa - $q * $ra, $vr, $vi); 913 $this->H[$i][$n - 1] = $this->cdivr; 914 $this->H[$i][$n] = $this->cdivi; 915 if (abs($x) > (abs($z) + abs($q))) { 916 $this->H[$i + 1][$n - 1] = (-$ra - $w * $this->H[$i][$n - 1] + $q * $this->H[$i][$n]) / $x; 917 $this->H[$i + 1][$n] = (-$sa - $w * $this->H[$i][$n] - $q * $this->H[$i][$n - 1]) / $x; 918 } else { 919 $this->cdiv(-$r - $y * $this->H[$i][$n - 1], -$s - $y * $this->H[$i][$n], $z, $q); 920 $this->H[$i + 1][$n - 1] = $this->cdivr; 921 $this->H[$i + 1][$n] = $this->cdivi; 922 } 923 } 924 925 // Overflow control 926 $t = max(abs($this->H[$i][$n - 1]), abs($this->H[$i][$n])); 927 if (($eps * $t) * $t > 1) { 928 for ($j = $i; $j <= $n; ++$j) { 929 $this->H[$j][$n - 1] /= $t; 930 $this->H[$j][$n] /= $t; 931 } 932 } 933 } // end else 934 } // end for 935 } // end else for complex case 936 } // end for 937 938 // Vectors of isolated roots 939 for ($i = 0; $i < $nn; ++$i) { 940 if ($i < $low || $i > $high) { 941 for ($j = $i; $j < $nn; ++$j) { 942 $this->V[$i][$j] = $this->H[$i][$j]; 943 } 944 } 945 } 946 947 // Back transformation to get eigenvectors of original matrix 948 for ($j = $nn - 1; $j >= $low; --$j) { 949 for ($i = $low; $i <= $high; ++$i) { 950 $z = 0.0; 951 for ($k = $low; $k <= min($j, $high); ++$k) { 952 $z += $this->V[$i][$k] * $this->H[$k][$j]; 953 } 954 955 $this->V[$i][$j] = $z; 956 } 957 } 958 } 959 }
title
Description
Body
title
Description
Body
title
Description
Body
title
Body
