1<?php 2 3namespace JAMA; 4 5/** 6 * @package JAMA 7 * 8 * Class to obtain eigenvalues and eigenvectors of a real matrix. 9 * 10 * If A is symmetric, then A = V*D*V' where the eigenvalue matrix D 11 * is diagonal and the eigenvector matrix V is orthogonal (i.e. 12 * A = V.times(D.times(V.transpose())) and V.times(V.transpose()) 13 * equals the identity matrix). 14 * 15 * If A is not symmetric, then the eigenvalue matrix D is block diagonal 16 * with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues, 17 * lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda]. The 18 * columns of V represent the eigenvectors in the sense that A*V = V*D, 19 * i.e. A.times(V) equals V.times(D). The matrix V may be badly 20 * conditioned, or even singular, so the validity of the equation 21 * A = V*D*inverse(V) depends upon V.cond(). 22 * 23 * @author Paul Meagher 24 * @license PHP v3.0 25 * @version 1.1 26 */ 27class EigenvalueDecomposition 28{ 29 /** 30 * Row and column dimension (square matrix). 31 * @var int 32 */ 33 private $n; 34 35 /** 36 * Internal symmetry flag. 37 * @var int 38 */ 39 private $issymmetric; 40 41 /** 42 * Arrays for internal storage of eigenvalues. 43 * @var array 44 */ 45 private $d = []; 46 private $e = []; 47 48 /** 49 * Array for internal storage of eigenvectors. 50 * @var array 51 */ 52 private $V = []; 53 54 /** 55 * Array for internal storage of nonsymmetric Hessenberg form. 56 * @var array 57 */ 58 private $H = []; 59 60 /** 61 * Working storage for nonsymmetric algorithm. 62 * @var array 63 */ 64 private $ort; 65 66 /** 67 * Used for complex scalar division. 68 * @var float 69 */ 70 private $cdivr; 71 private $cdivi; 72 73 /** 74 * Symmetric Householder reduction to tridiagonal form. 75 * 76 * @access private 77 */ 78 private function tred2() 79 { 80 // This is derived from the Algol procedures tred2 by 81 // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for 82 // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding 83 // Fortran subroutine in EISPACK. 84 $this->d = $this->V[$this->n - 1]; 85 // Householder reduction to tridiagonal form. 86 for ($i = $this->n - 1; $i > 0; --$i) { 87 $i_ = $i - 1; 88 // Scale to avoid under/overflow. 89 $h = $scale = 0.0; 90 $scale += array_sum(array_map('abs', $this->d)); 91 if ($scale == 0.0) { 92 $this->e[$i] = $this->d[$i_]; 93 $this->d = array_slice($this->V[$i_], 0, $i_); 94 for ($j = 0; $j < $i; ++$j) { 95 $this->V[$j][$i] = $this->V[$i][$j] = 0.0; 96 } 97 } else { 98 // Generate Householder vector. 99 for ($k = 0; $k < $i; ++$k) { 100 $this->d[$k] /= $scale; 101 $h += pow($this->d[$k], 2); 102 } 103 $f = $this->d[$i_]; 104 $g = sqrt($h); 105 if ($f > 0) { 106 $g = -$g; 107 } 108 $this->e[$i] = $scale * $g; 109 $h = $h - $f * $g; 110 $this->d[$i_] = $f - $g; 111 for ($j = 0; $j < $i; ++$j) { 112 $this->e[$j] = 0.0; 113 } 114 // Apply similarity transformation to remaining columns. 115 for ($j = 0; $j < $i; ++$j) { 116 $f = $this->d[$j]; 117 $this->V[$j][$i] = $f; 118 $g = $this->e[$j] + $this->V[$j][$j] * $f; 119 for ($k = $j + 1; $k <= $i_; ++$k) { 120 $g += $this->V[$k][$j] * $this->d[$k]; 121 $this->e[$k] += $this->V[$k][$j] * $f; 122 } 123 $this->e[$j] = $g; 124 } 125 $f = 0.0; 126 for ($j = 0; $j < $i; ++$j) { 127 $this->e[$j] /= $h; 128 $f += $this->e[$j] * $this->d[$j]; 129 } 130 $hh = $f / (2 * $h); 131 for ($j = 0; $j < $i; ++$j) { 132 $this->e[$j] -= $hh * $this->d[$j]; 133 } 134 for ($j = 0; $j < $i; ++$j) { 135 $f = $this->d[$j]; 136 $g = $this->e[$j]; 137 for ($k = $j; $k <= $i_; ++$k) { 138 $this->V[$k][$j] -= ($f * $this->e[$k] + $g * $this->d[$k]); 139 } 140 $this->d[$j] = $this->V[$i - 1][$j]; 141 $this->V[$i][$j] = 0.0; 142 } 143 } 144 $this->d[$i] = $h; 145 } 146 147 // Accumulate transformations. 148 for ($i = 0; $i < $this->n - 1; ++$i) { 149 $this->V[$this->n - 1][$i] = $this->V[$i][$i]; 150 $this->V[$i][$i] = 1.0; 151 $h = $this->d[$i + 1]; 152 if ($h != 0.0) { 153 for ($k = 0; $k <= $i; ++$k) { 154 $this->d[$k] = $this->V[$k][$i + 1] / $h; 155 } 156 for ($j = 0; $j <= $i; ++$j) { 157 $g = 0.0; 158 for ($k = 0; $k <= $i; ++$k) { 159 $g += $this->V[$k][$i + 1] * $this->V[$k][$j]; 160 } 161 for ($k = 0; $k <= $i; ++$k) { 162 $this->V[$k][$j] -= $g * $this->d[$k]; 163 } 164 } 165 } 166 for ($k = 0; $k <= $i; ++$k) { 167 $this->V[$k][$i + 1] = 0.0; 168 } 169 } 170 171 $this->d = $this->V[$this->n - 1]; 172 $this->V[$this->n - 1] = array_fill(0, $j, 0.0); 173 $this->V[$this->n - 1][$this->n - 1] = 1.0; 174 $this->e[0] = 0.0; 175 } 176 177 /** 178 * Symmetric tridiagonal QL algorithm. 179 * 180 * This is derived from the Algol procedures tql2, by 181 * Bowdler, Martin, Reinsch, and Wilkinson, Handbook for 182 * Auto. Comp., Vol.ii-Linear Algebra, and the corresponding 183 * Fortran subroutine in EISPACK. 184 * 185 * @access private 186 */ 187 private function tql2() 188 { 189 for ($i = 1; $i < $this->n; ++$i) { 190 $this->e[$i - 1] = $this->e[$i]; 191 } 192 $this->e[$this->n - 1] = 0.0; 193 $f = 0.0; 194 $tst1 = 0.0; 195 $eps = pow(2.0, -52.0); 196 197 for ($l = 0; $l < $this->n; ++$l) { 198 // Find small subdiagonal element 199 $tst1 = max($tst1, abs($this->d[$l]) + abs($this->e[$l])); 200 $m = $l; 201 while ($m < $this->n) { 202 if (abs($this->e[$m]) <= $eps * $tst1) { 203 break; 204 } 205 ++$m; 206 } 207 // If m == l, $this->d[l] is an eigenvalue, 208 // otherwise, iterate. 209 if ($m > $l) { 210 $iter = 0; 211 do { 212 // Could check iteration count here. 213 $iter += 1; 214 // Compute implicit shift 215 $g = $this->d[$l]; 216 $p = ($this->d[$l + 1] - $g) / (2.0 * $this->e[$l]); 217 $r = hypo($p, 1.0); 218 if ($p < 0) { 219 $r *= -1; 220 } 221 $this->d[$l] = $this->e[$l] / ($p + $r); 222 $this->d[$l + 1] = $this->e[$l] * ($p + $r); 223 $dl1 = $this->d[$l + 1]; 224 $h = $g - $this->d[$l]; 225 for ($i = $l + 2; $i < $this->n; ++$i) { 226 $this->d[$i] -= $h; 227 } 228 $f += $h; 229 // Implicit QL transformation. 230 $p = $this->d[$m]; 231 $c = 1.0; 232 $c2 = $c3 = $c; 233 $el1 = $this->e[$l + 1]; 234 $s = $s2 = 0.0; 235 for ($i = $m - 1; $i >= $l; --$i) { 236 $c3 = $c2; 237 $c2 = $c; 238 $s2 = $s; 239 $g = $c * $this->e[$i]; 240 $h = $c * $p; 241 $r = hypo($p, $this->e[$i]); 242 $this->e[$i + 1] = $s * $r; 243 $s = $this->e[$i] / $r; 244 $c = $p / $r; 245 $p = $c * $this->d[$i] - $s * $g; 246 $this->d[$i + 1] = $h + $s * ($c * $g + $s * $this->d[$i]); 247 // Accumulate transformation. 248 for ($k = 0; $k < $this->n; ++$k) { 249 $h = $this->V[$k][$i + 1]; 250 $this->V[$k][$i + 1] = $s * $this->V[$k][$i] + $c * $h; 251 $this->V[$k][$i] = $c * $this->V[$k][$i] - $s * $h; 252 } 253 } 254 $p = -$s * $s2 * $c3 * $el1 * $this->e[$l] / $dl1; 255 $this->e[$l] = $s * $p; 256 $this->d[$l] = $c * $p; 257 // Check for convergence. 258 } while (abs($this->e[$l]) > $eps * $tst1); 259 } 260 $this->d[$l] = $this->d[$l] + $f; 261 $this->e[$l] = 0.0; 262 } 263 264 // Sort eigenvalues and corresponding vectors. 265 for ($i = 0; $i < $this->n - 1; ++$i) { 266 $k = $i; 267 $p = $this->d[$i]; 268 for ($j = $i + 1; $j < $this->n; ++$j) { 269 if ($this->d[$j] < $p) { 270 $k = $j; 271 $p = $this->d[$j]; 272 } 273 } 274 if ($k != $i) { 275 $this->d[$k] = $this->d[$i]; 276 $this->d[$i] = $p; 277 for ($j = 0; $j < $this->n; ++$j) { 278 $p = $this->V[$j][$i]; 279 $this->V[$j][$i] = $this->V[$j][$k]; 280 $this->V[$j][$k] = $p; 281 } 282 } 283 } 284 } 285 286 /** 287 * Nonsymmetric reduction to Hessenberg form. 288 * 289 * This is derived from the Algol procedures orthes and ortran, 290 * by Martin and Wilkinson, Handbook for Auto. Comp., 291 * Vol.ii-Linear Algebra, and the corresponding 292 * Fortran subroutines in EISPACK. 293 * 294 * @access private 295 */ 296 private function orthes() 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 * @access private 377 * @param mixed $xr 378 * @param mixed $xi 379 * @param mixed $yr 380 * @param mixed $yi 381 */ 382 private function cdiv($xr, $xi, $yr, $yi) 383 { 384 if (abs($yr) > abs($yi)) { 385 $r = $yi / $yr; 386 $d = $yr + $r * $yi; 387 $this->cdivr = ($xr + $r * $xi) / $d; 388 $this->cdivi = ($xi - $r * $xr) / $d; 389 } else { 390 $r = $yr / $yi; 391 $d = $yi + $r * $yr; 392 $this->cdivr = ($r * $xr + $xi) / $d; 393 $this->cdivi = ($r * $xi - $xr) / $d; 394 } 395 } 396 397 /** 398 * Nonsymmetric reduction from Hessenberg to real Schur form. 399 * 400 * Code is derived from the Algol procedure hqr2, 401 * by Martin and Wilkinson, Handbook for Auto. Comp., 402 * Vol.ii-Linear Algebra, and the corresponding 403 * Fortran subroutine in EISPACK. 404 * 405 * @access private 406 */ 407 private function hqr2() 408 { 409 // Initialize 410 $nn = $this->n; 411 $n = $nn - 1; 412 $low = 0; 413 $high = $nn - 1; 414 $eps = pow(2.0, -52.0); 415 $exshift = 0.0; 416 $p = $q = $r = $s = $z = 0; 417 // Store roots isolated by balanc and compute matrix norm 418 $norm = 0.0; 419 420 for ($i = 0; $i < $nn; ++$i) { 421 if (($i < $low) or ($i > $high)) { 422 $this->d[$i] = $this->H[$i][$i]; 423 $this->e[$i] = 0.0; 424 } 425 for ($j = max($i - 1, 0); $j < $nn; ++$j) { 426 $norm = $norm + abs($this->H[$i][$j]); 427 } 428 } 429 430 // Outer loop over eigenvalue index 431 $iter = 0; 432 while ($n >= $low) { 433 // Look for single small sub-diagonal element 434 $l = $n; 435 while ($l > $low) { 436 $s = abs($this->H[$l - 1][$l - 1]) + abs($this->H[$l][$l]); 437 if ($s == 0.0) { 438 $s = $norm; 439 } 440 if (abs($this->H[$l][$l - 1]) < $eps * $s) { 441 break; 442 } 443 --$l; 444 } 445 // Check for convergence 446 // One root found 447 if ($l == $n) { 448 $this->H[$n][$n] = $this->H[$n][$n] + $exshift; 449 $this->d[$n] = $this->H[$n][$n]; 450 $this->e[$n] = 0.0; 451 --$n; 452 $iter = 0; 453 // Two roots found 454 } elseif ($l == $n - 1) { 455 $w = $this->H[$n][$n - 1] * $this->H[$n - 1][$n]; 456 $p = ($this->H[$n - 1][$n - 1] - $this->H[$n][$n]) / 2.0; 457 $q = $p * $p + $w; 458 $z = sqrt(abs($q)); 459 $this->H[$n][$n] = $this->H[$n][$n] + $exshift; 460 $this->H[$n - 1][$n - 1] = $this->H[$n - 1][$n - 1] + $exshift; 461 $x = $this->H[$n][$n]; 462 // Real pair 463 if ($q >= 0) { 464 if ($p >= 0) { 465 $z = $p + $z; 466 } else { 467 $z = $p - $z; 468 } 469 $this->d[$n - 1] = $x + $z; 470 $this->d[$n] = $this->d[$n - 1]; 471 if ($z != 0.0) { 472 $this->d[$n] = $x - $w / $z; 473 } 474 $this->e[$n - 1] = 0.0; 475 $this->e[$n] = 0.0; 476 $x = $this->H[$n][$n - 1]; 477 $s = abs($x) + abs($z); 478 $p = $x / $s; 479 $q = $z / $s; 480 $r = sqrt($p * $p + $q * $q); 481 $p = $p / $r; 482 $q = $q / $r; 483 // Row modification 484 for ($j = $n - 1; $j < $nn; ++$j) { 485 $z = $this->H[$n - 1][$j]; 486 $this->H[$n - 1][$j] = $q * $z + $p * $this->H[$n][$j]; 487 $this->H[$n][$j] = $q * $this->H[$n][$j] - $p * $z; 488 } 489 // Column modification 490 for ($i = 0; $i <= $n; ++$i) { 491 $z = $this->H[$i][$n - 1]; 492 $this->H[$i][$n - 1] = $q * $z + $p * $this->H[$i][$n]; 493 $this->H[$i][$n] = $q * $this->H[$i][$n] - $p * $z; 494 } 495 // Accumulate transformations 496 for ($i = $low; $i <= $high; ++$i) { 497 $z = $this->V[$i][$n - 1]; 498 $this->V[$i][$n - 1] = $q * $z + $p * $this->V[$i][$n]; 499 $this->V[$i][$n] = $q * $this->V[$i][$n] - $p * $z; 500 } 501 // Complex pair 502 } else { 503 $this->d[$n - 1] = $x + $p; 504 $this->d[$n] = $x + $p; 505 $this->e[$n - 1] = $z; 506 $this->e[$n] = -$z; 507 } 508 $n = $n - 2; 509 $iter = 0; 510 // No convergence yet 511 } else { 512 // Form shift 513 $x = $this->H[$n][$n]; 514 $y = 0.0; 515 $w = 0.0; 516 if ($l < $n) { 517 $y = $this->H[$n - 1][$n - 1]; 518 $w = $this->H[$n][$n - 1] * $this->H[$n - 1][$n]; 519 } 520 // Wilkinson's original ad hoc shift 521 if ($iter == 10) { 522 $exshift += $x; 523 for ($i = $low; $i <= $n; ++$i) { 524 $this->H[$i][$i] -= $x; 525 } 526 $s = abs($this->H[$n][$n - 1]) + abs($this->H[$n - 1][$n - 2]); 527 $x = $y = 0.75 * $s; 528 $w = -0.4375 * $s * $s; 529 } 530 // MATLAB's new ad hoc shift 531 if ($iter == 30) { 532 $s = ($y - $x) / 2.0; 533 $s = $s * $s + $w; 534 if ($s > 0) { 535 $s = sqrt($s); 536 if ($y < $x) { 537 $s = -$s; 538 } 539 $s = $x - $w / (($y - $x) / 2.0 + $s); 540 for ($i = $low; $i <= $n; ++$i) { 541 $this->H[$i][$i] -= $s; 542 } 543 $exshift += $s; 544 $x = $y = $w = 0.964; 545 } 546 } 547 // Could check iteration count here. 548 $iter = $iter + 1; 549 // Look for two consecutive small sub-diagonal elements 550 $m = $n - 2; 551 while ($m >= $l) { 552 $z = $this->H[$m][$m]; 553 $r = $x - $z; 554 $s = $y - $z; 555 $p = ($r * $s - $w) / $this->H[$m + 1][$m] + $this->H[$m][$m + 1]; 556 $q = $this->H[$m + 1][$m + 1] - $z - $r - $s; 557 $r = $this->H[$m + 2][$m + 1]; 558 $s = abs($p) + abs($q) + abs($r); 559 $p = $p / $s; 560 $q = $q / $s; 561 $r = $r / $s; 562 if ($m == $l) { 563 break; 564 } 565 if (abs($this->H[$m][$m - 1]) * (abs($q) + abs($r)) < 566 $eps * (abs($p) * (abs($this->H[$m - 1][$m - 1]) + abs($z) + abs($this->H[$m + 1][$m + 1])))) { 567 break; 568 } 569 --$m; 570 } 571 for ($i = $m + 2; $i <= $n; ++$i) { 572 $this->H[$i][$i - 2] = 0.0; 573 if ($i > $m + 2) { 574 $this->H[$i][$i - 3] = 0.0; 575 } 576 } 577 // Double QR step involving rows l:n and columns m:n 578 for ($k = $m; $k <= $n - 1; ++$k) { 579 $notlast = ($k != $n - 1); 580 if ($k != $m) { 581 $p = $this->H[$k][$k - 1]; 582 $q = $this->H[$k + 1][$k - 1]; 583 $r = ($notlast ? $this->H[$k + 2][$k - 1] : 0.0); 584 $x = abs($p) + abs($q) + abs($r); 585 if ($x != 0.0) { 586 $p = $p / $x; 587 $q = $q / $x; 588 $r = $r / $x; 589 } 590 } 591 if ($x == 0.0) { 592 break; 593 } 594 $s = sqrt($p * $p + $q * $q + $r * $r); 595 if ($p < 0) { 596 $s = -$s; 597 } 598 if ($s != 0) { 599 if ($k != $m) { 600 $this->H[$k][$k - 1] = -$s * $x; 601 } elseif ($l != $m) { 602 $this->H[$k][$k - 1] = -$this->H[$k][$k - 1]; 603 } 604 $p = $p + $s; 605 $x = $p / $s; 606 $y = $q / $s; 607 $z = $r / $s; 608 $q = $q / $p; 609 $r = $r / $p; 610 // Row modification 611 for ($j = $k; $j < $nn; ++$j) { 612 $p = $this->H[$k][$j] + $q * $this->H[$k + 1][$j]; 613 if ($notlast) { 614 $p = $p + $r * $this->H[$k + 2][$j]; 615 $this->H[$k + 2][$j] = $this->H[$k + 2][$j] - $p * $z; 616 } 617 $this->H[$k][$j] = $this->H[$k][$j] - $p * $x; 618 $this->H[$k + 1][$j] = $this->H[$k + 1][$j] - $p * $y; 619 } 620 // Column modification 621 for ($i = 0; $i <= min($n, $k + 3); ++$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 < $low | $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 for ($k = $low; $k <= min($j, $high); ++$k) { 777 $z = $z + $this->V[$i][$k] * $this->H[$k][$j]; 778 } 779 $this->V[$i][$j] = $z; 780 } 781 } 782 } // end hqr2 783 784 /** 785 * Constructor: Check for symmetry, then construct the eigenvalue decomposition 786 * 787 * @access public 788 * @param A Square matrix 789 * @param mixed $Arg 790 * @return Structure to access D and V. 791 */ 792 public function __construct($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 * @access public 824 * @return Matrix 825 */ 826 public function getV() 827 { 828 return new Matrix($this->V, $this->n, $this->n); 829 } 830 831 /** 832 * Return the real parts of the eigenvalues 833 * 834 * @access public 835 * @return array 836 */ 837 public function getRealEigenvalues() 838 { 839 return $this->d; 840 } 841 842 /** 843 * Return the imaginary parts of the eigenvalues 844 * 845 * @access public 846 * @return array 847 */ 848 public function getImagEigenvalues() 849 { 850 return $this->e; 851 } 852 853 /** 854 * Return the block diagonal eigenvalue matrix 855 * 856 * @access public 857 * @return Matrix 858 */ 859 public function getD() 860 { 861 for ($i = 0; $i < $this->n; ++$i) { 862 $D[$i] = array_fill(0, $this->n, 0.0); 863 $D[$i][$i] = $this->d[$i]; 864 if ($this->e[$i] == 0) { 865 continue; 866 } 867 $o = ($this->e[$i] > 0) ? $i + 1 : $i - 1; 868 $D[$i][$o] = $this->e[$i]; 869 } 870 return new Matrix($D); 871 } 872} // class EigenvalueDecomposition 873