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