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