| [ Index ] |
PHP Cross Reference of moodle-2.8 |
[Summary view] [Print] [Text view]
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 */ 24 class 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 /** 72 * Symmetric Householder reduction to tridiagonal form. 73 * 74 * @access private 75 */ 76 private function tred2 () { 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 /** 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 * @access private 184 */ 185 private function tql2() { 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 ++$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 $this->d[$l] = $this->e[$l] / ($p + $r); 217 $this->d[$l+1] = $this->e[$l] * ($p + $r); 218 $dl1 = $this->d[$l+1]; 219 $h = $g - $this->d[$l]; 220 for ($i = $l + 2; $i < $this->n; ++$i) 221 $this->d[$i] -= $h; 222 $f += $h; 223 // Implicit QL transformation. 224 $p = $this->d[$m]; 225 $c = 1.0; 226 $c2 = $c3 = $c; 227 $el1 = $this->e[$l + 1]; 228 $s = $s2 = 0.0; 229 for ($i = $m-1; $i >= $l; --$i) { 230 $c3 = $c2; 231 $c2 = $c; 232 $s2 = $s; 233 $g = $c * $this->e[$i]; 234 $h = $c * $p; 235 $r = hypo($p, $this->e[$i]); 236 $this->e[$i+1] = $s * $r; 237 $s = $this->e[$i] / $r; 238 $c = $p / $r; 239 $p = $c * $this->d[$i] - $s * $g; 240 $this->d[$i+1] = $h + $s * ($c * $g + $s * $this->d[$i]); 241 // Accumulate transformation. 242 for ($k = 0; $k < $this->n; ++$k) { 243 $h = $this->V[$k][$i+1]; 244 $this->V[$k][$i+1] = $s * $this->V[$k][$i] + $c * $h; 245 $this->V[$k][$i] = $c * $this->V[$k][$i] - $s * $h; 246 } 247 } 248 $p = -$s * $s2 * $c3 * $el1 * $this->e[$l] / $dl1; 249 $this->e[$l] = $s * $p; 250 $this->d[$l] = $c * $p; 251 // Check for convergence. 252 } while (abs($this->e[$l]) > $eps * $tst1); 253 } 254 $this->d[$l] = $this->d[$l] + $f; 255 $this->e[$l] = 0.0; 256 } 257 258 // Sort eigenvalues and corresponding vectors. 259 for ($i = 0; $i < $this->n - 1; ++$i) { 260 $k = $i; 261 $p = $this->d[$i]; 262 for ($j = $i+1; $j < $this->n; ++$j) { 263 if ($this->d[$j] < $p) { 264 $k = $j; 265 $p = $this->d[$j]; 266 } 267 } 268 if ($k != $i) { 269 $this->d[$k] = $this->d[$i]; 270 $this->d[$i] = $p; 271 for ($j = 0; $j < $this->n; ++$j) { 272 $p = $this->V[$j][$i]; 273 $this->V[$j][$i] = $this->V[$j][$k]; 274 $this->V[$j][$k] = $p; 275 } 276 } 277 } 278 } 279 280 281 /** 282 * Nonsymmetric reduction to Hessenberg form. 283 * 284 * This is derived from the Algol procedures orthes and ortran, 285 * by Martin and Wilkinson, Handbook for Auto. Comp., 286 * Vol.ii-Linear Algebra, and the corresponding 287 * Fortran subroutines in EISPACK. 288 * 289 * @access private 290 */ 291 private function orthes () { 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 /** 369 * Performs complex division. 370 * 371 * @access private 372 */ 373 private function cdiv($xr, $xi, $yr, $yi) { 374 if (abs($yr) > abs($yi)) { 375 $r = $yi / $yr; 376 $d = $yr + $r * $yi; 377 $this->cdivr = ($xr + $r * $xi) / $d; 378 $this->cdivi = ($xi - $r * $xr) / $d; 379 } else { 380 $r = $yr / $yi; 381 $d = $yi + $r * $yr; 382 $this->cdivr = ($r * $xr + $xi) / $d; 383 $this->cdivi = ($r * $xi - $xr) / $d; 384 } 385 } 386 387 388 /** 389 * Nonsymmetric reduction from Hessenberg to real Schur form. 390 * 391 * Code is derived from the Algol procedure hqr2, 392 * by Martin and Wilkinson, Handbook for Auto. Comp., 393 * Vol.ii-Linear Algebra, and the corresponding 394 * Fortran subroutine in EISPACK. 395 * 396 * @access private 397 */ 398 private function hqr2 () { 399 // Initialize 400 $nn = $this->n; 401 $n = $nn - 1; 402 $low = 0; 403 $high = $nn - 1; 404 $eps = pow(2.0, -52.0); 405 $exshift = 0.0; 406 $p = $q = $r = $s = $z = 0; 407 // Store roots isolated by balanc and compute matrix norm 408 $norm = 0.0; 409 410 for ($i = 0; $i < $nn; ++$i) { 411 if (($i < $low) OR ($i > $high)) { 412 $this->d[$i] = $this->H[$i][$i]; 413 $this->e[$i] = 0.0; 414 } 415 for ($j = max($i-1, 0); $j < $nn; ++$j) { 416 $norm = $norm + abs($this->H[$i][$j]); 417 } 418 } 419 420 // Outer loop over eigenvalue index 421 $iter = 0; 422 while ($n >= $low) { 423 // Look for single small sub-diagonal element 424 $l = $n; 425 while ($l > $low) { 426 $s = abs($this->H[$l-1][$l-1]) + abs($this->H[$l][$l]); 427 if ($s == 0.0) { 428 $s = $norm; 429 } 430 if (abs($this->H[$l][$l-1]) < $eps * $s) { 431 break; 432 } 433 --$l; 434 } 435 // Check for convergence 436 // One root found 437 if ($l == $n) { 438 $this->H[$n][$n] = $this->H[$n][$n] + $exshift; 439 $this->d[$n] = $this->H[$n][$n]; 440 $this->e[$n] = 0.0; 441 --$n; 442 $iter = 0; 443 // Two roots found 444 } else if ($l == $n-1) { 445 $w = $this->H[$n][$n-1] * $this->H[$n-1][$n]; 446 $p = ($this->H[$n-1][$n-1] - $this->H[$n][$n]) / 2.0; 447 $q = $p * $p + $w; 448 $z = sqrt(abs($q)); 449 $this->H[$n][$n] = $this->H[$n][$n] + $exshift; 450 $this->H[$n-1][$n-1] = $this->H[$n-1][$n-1] + $exshift; 451 $x = $this->H[$n][$n]; 452 // Real pair 453 if ($q >= 0) { 454 if ($p >= 0) { 455 $z = $p + $z; 456 } else { 457 $z = $p - $z; 458 } 459 $this->d[$n-1] = $x + $z; 460 $this->d[$n] = $this->d[$n-1]; 461 if ($z != 0.0) { 462 $this->d[$n] = $x - $w / $z; 463 } 464 $this->e[$n-1] = 0.0; 465 $this->e[$n] = 0.0; 466 $x = $this->H[$n][$n-1]; 467 $s = abs($x) + abs($z); 468 $p = $x / $s; 469 $q = $z / $s; 470 $r = sqrt($p * $p + $q * $q); 471 $p = $p / $r; 472 $q = $q / $r; 473 // Row modification 474 for ($j = $n-1; $j < $nn; ++$j) { 475 $z = $this->H[$n-1][$j]; 476 $this->H[$n-1][$j] = $q * $z + $p * $this->H[$n][$j]; 477 $this->H[$n][$j] = $q * $this->H[$n][$j] - $p * $z; 478 } 479 // Column modification 480 for ($i = 0; $i <= n; ++$i) { 481 $z = $this->H[$i][$n-1]; 482 $this->H[$i][$n-1] = $q * $z + $p * $this->H[$i][$n]; 483 $this->H[$i][$n] = $q * $this->H[$i][$n] - $p * $z; 484 } 485 // Accumulate transformations 486 for ($i = $low; $i <= $high; ++$i) { 487 $z = $this->V[$i][$n-1]; 488 $this->V[$i][$n-1] = $q * $z + $p * $this->V[$i][$n]; 489 $this->V[$i][$n] = $q * $this->V[$i][$n] - $p * $z; 490 } 491 // Complex pair 492 } else { 493 $this->d[$n-1] = $x + $p; 494 $this->d[$n] = $x + $p; 495 $this->e[$n-1] = $z; 496 $this->e[$n] = -$z; 497 } 498 $n = $n - 2; 499 $iter = 0; 500 // No convergence yet 501 } else { 502 // Form shift 503 $x = $this->H[$n][$n]; 504 $y = 0.0; 505 $w = 0.0; 506 if ($l < $n) { 507 $y = $this->H[$n-1][$n-1]; 508 $w = $this->H[$n][$n-1] * $this->H[$n-1][$n]; 509 } 510 // Wilkinson's original ad hoc shift 511 if ($iter == 10) { 512 $exshift += $x; 513 for ($i = $low; $i <= $n; ++$i) { 514 $this->H[$i][$i] -= $x; 515 } 516 $s = abs($this->H[$n][$n-1]) + abs($this->H[$n-1][$n-2]); 517 $x = $y = 0.75 * $s; 518 $w = -0.4375 * $s * $s; 519 } 520 // MATLAB's new ad hoc shift 521 if ($iter == 30) { 522 $s = ($y - $x) / 2.0; 523 $s = $s * $s + $w; 524 if ($s > 0) { 525 $s = sqrt($s); 526 if ($y < $x) { 527 $s = -$s; 528 } 529 $s = $x - $w / (($y - $x) / 2.0 + $s); 530 for ($i = $low; $i <= $n; ++$i) { 531 $this->H[$i][$i] -= $s; 532 } 533 $exshift += $s; 534 $x = $y = $w = 0.964; 535 } 536 } 537 // Could check iteration count here. 538 $iter = $iter + 1; 539 // Look for two consecutive small sub-diagonal elements 540 $m = $n - 2; 541 while ($m >= $l) { 542 $z = $this->H[$m][$m]; 543 $r = $x - $z; 544 $s = $y - $z; 545 $p = ($r * $s - $w) / $this->H[$m+1][$m] + $this->H[$m][$m+1]; 546 $q = $this->H[$m+1][$m+1] - $z - $r - $s; 547 $r = $this->H[$m+2][$m+1]; 548 $s = abs($p) + abs($q) + abs($r); 549 $p = $p / $s; 550 $q = $q / $s; 551 $r = $r / $s; 552 if ($m == $l) { 553 break; 554 } 555 if (abs($this->H[$m][$m-1]) * (abs($q) + abs($r)) < 556 $eps * (abs($p) * (abs($this->H[$m-1][$m-1]) + abs($z) + abs($this->H[$m+1][$m+1])))) { 557 break; 558 } 559 --$m; 560 } 561 for ($i = $m + 2; $i <= $n; ++$i) { 562 $this->H[$i][$i-2] = 0.0; 563 if ($i > $m+2) { 564 $this->H[$i][$i-3] = 0.0; 565 } 566 } 567 // Double QR step involving rows l:n and columns m:n 568 for ($k = $m; $k <= $n-1; ++$k) { 569 $notlast = ($k != $n-1); 570 if ($k != $m) { 571 $p = $this->H[$k][$k-1]; 572 $q = $this->H[$k+1][$k-1]; 573 $r = ($notlast ? $this->H[$k+2][$k-1] : 0.0); 574 $x = abs($p) + abs($q) + abs($r); 575 if ($x != 0.0) { 576 $p = $p / $x; 577 $q = $q / $x; 578 $r = $r / $x; 579 } 580 } 581 if ($x == 0.0) { 582 break; 583 } 584 $s = sqrt($p * $p + $q * $q + $r * $r); 585 if ($p < 0) { 586 $s = -$s; 587 } 588 if ($s != 0) { 589 if ($k != $m) { 590 $this->H[$k][$k-1] = -$s * $x; 591 } elseif ($l != $m) { 592 $this->H[$k][$k-1] = -$this->H[$k][$k-1]; 593 } 594 $p = $p + $s; 595 $x = $p / $s; 596 $y = $q / $s; 597 $z = $r / $s; 598 $q = $q / $p; 599 $r = $r / $p; 600 // Row modification 601 for ($j = $k; $j < $nn; ++$j) { 602 $p = $this->H[$k][$j] + $q * $this->H[$k+1][$j]; 603 if ($notlast) { 604 $p = $p + $r * $this->H[$k+2][$j]; 605 $this->H[$k+2][$j] = $this->H[$k+2][$j] - $p * $z; 606 } 607 $this->H[$k][$j] = $this->H[$k][$j] - $p * $x; 608 $this->H[$k+1][$j] = $this->H[$k+1][$j] - $p * $y; 609 } 610 // Column modification 611 for ($i = 0; $i <= min($n, $k+3); ++$i) { 612 $p = $x * $this->H[$i][$k] + $y * $this->H[$i][$k+1]; 613 if ($notlast) { 614 $p = $p + $z * $this->H[$i][$k+2]; 615 $this->H[$i][$k+2] = $this->H[$i][$k+2] - $p * $r; 616 } 617 $this->H[$i][$k] = $this->H[$i][$k] - $p; 618 $this->H[$i][$k+1] = $this->H[$i][$k+1] - $p * $q; 619 } 620 // Accumulate transformations 621 for ($i = $low; $i <= $high; ++$i) { 622 $p = $x * $this->V[$i][$k] + $y * $this->V[$i][$k+1]; 623 if ($notlast) { 624 $p = $p + $z * $this->V[$i][$k+2]; 625 $this->V[$i][$k+2] = $this->V[$i][$k+2] - $p * $r; 626 } 627 $this->V[$i][$k] = $this->V[$i][$k] - $p; 628 $this->V[$i][$k+1] = $this->V[$i][$k+1] - $p * $q; 629 } 630 } // ($s != 0) 631 } // k loop 632 } // check convergence 633 } // while ($n >= $low) 634 635 // Backsubstitute to find vectors of upper triangular form 636 if ($norm == 0.0) { 637 return; 638 } 639 640 for ($n = $nn-1; $n >= 0; --$n) { 641 $p = $this->d[$n]; 642 $q = $this->e[$n]; 643 // Real vector 644 if ($q == 0) { 645 $l = $n; 646 $this->H[$n][$n] = 1.0; 647 for ($i = $n-1; $i >= 0; --$i) { 648 $w = $this->H[$i][$i] - $p; 649 $r = 0.0; 650 for ($j = $l; $j <= $n; ++$j) { 651 $r = $r + $this->H[$i][$j] * $this->H[$j][$n]; 652 } 653 if ($this->e[$i] < 0.0) { 654 $z = $w; 655 $s = $r; 656 } else { 657 $l = $i; 658 if ($this->e[$i] == 0.0) { 659 if ($w != 0.0) { 660 $this->H[$i][$n] = -$r / $w; 661 } else { 662 $this->H[$i][$n] = -$r / ($eps * $norm); 663 } 664 // Solve real equations 665 } else { 666 $x = $this->H[$i][$i+1]; 667 $y = $this->H[$i+1][$i]; 668 $q = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i]; 669 $t = ($x * $s - $z * $r) / $q; 670 $this->H[$i][$n] = $t; 671 if (abs($x) > abs($z)) { 672 $this->H[$i+1][$n] = (-$r - $w * $t) / $x; 673 } else { 674 $this->H[$i+1][$n] = (-$s - $y * $t) / $z; 675 } 676 } 677 // Overflow control 678 $t = abs($this->H[$i][$n]); 679 if (($eps * $t) * $t > 1) { 680 for ($j = $i; $j <= $n; ++$j) { 681 $this->H[$j][$n] = $this->H[$j][$n] / $t; 682 } 683 } 684 } 685 } 686 // Complex vector 687 } else if ($q < 0) { 688 $l = $n-1; 689 // Last vector component imaginary so matrix is triangular 690 if (abs($this->H[$n][$n-1]) > abs($this->H[$n-1][$n])) { 691 $this->H[$n-1][$n-1] = $q / $this->H[$n][$n-1]; 692 $this->H[$n-1][$n] = -($this->H[$n][$n] - $p) / $this->H[$n][$n-1]; 693 } else { 694 $this->cdiv(0.0, -$this->H[$n-1][$n], $this->H[$n-1][$n-1] - $p, $q); 695 $this->H[$n-1][$n-1] = $this->cdivr; 696 $this->H[$n-1][$n] = $this->cdivi; 697 } 698 $this->H[$n][$n-1] = 0.0; 699 $this->H[$n][$n] = 1.0; 700 for ($i = $n-2; $i >= 0; --$i) { 701 // double ra,sa,vr,vi; 702 $ra = 0.0; 703 $sa = 0.0; 704 for ($j = $l; $j <= $n; ++$j) { 705 $ra = $ra + $this->H[$i][$j] * $this->H[$j][$n-1]; 706 $sa = $sa + $this->H[$i][$j] * $this->H[$j][$n]; 707 } 708 $w = $this->H[$i][$i] - $p; 709 if ($this->e[$i] < 0.0) { 710 $z = $w; 711 $r = $ra; 712 $s = $sa; 713 } else { 714 $l = $i; 715 if ($this->e[$i] == 0) { 716 $this->cdiv(-$ra, -$sa, $w, $q); 717 $this->H[$i][$n-1] = $this->cdivr; 718 $this->H[$i][$n] = $this->cdivi; 719 } else { 720 // Solve complex equations 721 $x = $this->H[$i][$i+1]; 722 $y = $this->H[$i+1][$i]; 723 $vr = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i] - $q * $q; 724 $vi = ($this->d[$i] - $p) * 2.0 * $q; 725 if ($vr == 0.0 & $vi == 0.0) { 726 $vr = $eps * $norm * (abs($w) + abs($q) + abs($x) + abs($y) + abs($z)); 727 } 728 $this->cdiv($x * $r - $z * $ra + $q * $sa, $x * $s - $z * $sa - $q * $ra, $vr, $vi); 729 $this->H[$i][$n-1] = $this->cdivr; 730 $this->H[$i][$n] = $this->cdivi; 731 if (abs($x) > (abs($z) + abs($q))) { 732 $this->H[$i+1][$n-1] = (-$ra - $w * $this->H[$i][$n-1] + $q * $this->H[$i][$n]) / $x; 733 $this->H[$i+1][$n] = (-$sa - $w * $this->H[$i][$n] - $q * $this->H[$i][$n-1]) / $x; 734 } else { 735 $this->cdiv(-$r - $y * $this->H[$i][$n-1], -$s - $y * $this->H[$i][$n], $z, $q); 736 $this->H[$i+1][$n-1] = $this->cdivr; 737 $this->H[$i+1][$n] = $this->cdivi; 738 } 739 } 740 // Overflow control 741 $t = max(abs($this->H[$i][$n-1]),abs($this->H[$i][$n])); 742 if (($eps * $t) * $t > 1) { 743 for ($j = $i; $j <= $n; ++$j) { 744 $this->H[$j][$n-1] = $this->H[$j][$n-1] / $t; 745 $this->H[$j][$n] = $this->H[$j][$n] / $t; 746 } 747 } 748 } // end else 749 } // end for 750 } // end else for complex case 751 } // end for 752 753 // Vectors of isolated roots 754 for ($i = 0; $i < $nn; ++$i) { 755 if ($i < $low | $i > $high) { 756 for ($j = $i; $j < $nn; ++$j) { 757 $this->V[$i][$j] = $this->H[$i][$j]; 758 } 759 } 760 } 761 762 // Back transformation to get eigenvectors of original matrix 763 for ($j = $nn-1; $j >= $low; --$j) { 764 for ($i = $low; $i <= $high; ++$i) { 765 $z = 0.0; 766 for ($k = $low; $k <= min($j,$high); ++$k) { 767 $z = $z + $this->V[$i][$k] * $this->H[$k][$j]; 768 } 769 $this->V[$i][$j] = $z; 770 } 771 } 772 } // end hqr2 773 774 775 /** 776 * Constructor: Check for symmetry, then construct the eigenvalue decomposition 777 * 778 * @access public 779 * @param A Square matrix 780 * @return Structure to access D and V. 781 */ 782 public function __construct($Arg) { 783 $this->A = $Arg->getArray(); 784 $this->n = $Arg->getColumnDimension(); 785 786 $issymmetric = true; 787 for ($j = 0; ($j < $this->n) & $issymmetric; ++$j) { 788 for ($i = 0; ($i < $this->n) & $issymmetric; ++$i) { 789 $issymmetric = ($this->A[$i][$j] == $this->A[$j][$i]); 790 } 791 } 792 793 if ($issymmetric) { 794 $this->V = $this->A; 795 // Tridiagonalize. 796 $this->tred2(); 797 // Diagonalize. 798 $this->tql2(); 799 } else { 800 $this->H = $this->A; 801 $this->ort = array(); 802 // Reduce to Hessenberg form. 803 $this->orthes(); 804 // Reduce Hessenberg to real Schur form. 805 $this->hqr2(); 806 } 807 } 808 809 810 /** 811 * Return the eigenvector matrix 812 * 813 * @access public 814 * @return V 815 */ 816 public function getV() { 817 return new Matrix($this->V, $this->n, $this->n); 818 } 819 820 821 /** 822 * Return the real parts of the eigenvalues 823 * 824 * @access public 825 * @return real(diag(D)) 826 */ 827 public function getRealEigenvalues() { 828 return $this->d; 829 } 830 831 832 /** 833 * Return the imaginary parts of the eigenvalues 834 * 835 * @access public 836 * @return imag(diag(D)) 837 */ 838 public function getImagEigenvalues() { 839 return $this->e; 840 } 841 842 843 /** 844 * Return the block diagonal eigenvalue matrix 845 * 846 * @access public 847 * @return D 848 */ 849 public function getD() { 850 for ($i = 0; $i < $this->n; ++$i) { 851 $D[$i] = array_fill(0, $this->n, 0.0); 852 $D[$i][$i] = $this->d[$i]; 853 if ($this->e[$i] == 0) { 854 continue; 855 } 856 $o = ($this->e[$i] > 0) ? $i + 1 : $i - 1; 857 $D[$i][$o] = $this->e[$i]; 858 } 859 return new Matrix($D); 860 } 861 862 } // class EigenvalueDecomposition
title
Description
Body
title
Description
Body
title
Description
Body
title
Body
| Generated: Fri Nov 28 20:29:05 2014 | Cross-referenced by PHPXref 0.7.1 |