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