ILIAS  release_5-2 Revision v5.2.25-18-g3f80b828510
EigenvalueDecomposition.php
Go to the documentation of this file.
1<?php
25
30 private $n;
31
36 private $issymmetric;
37
42 private $d = array();
43 private $e = array();
44
49 private $V = array();
50
55 private $H = array();
56
61 private $ort;
62
67 private $cdivr;
68 private $cdivi;
69
70
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
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
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
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
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
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
816 public function getV() {
817 return new Matrix($this->V, $this->n, $this->n);
818 }
819
820
827 public function getRealEigenvalues() {
828 return $this->d;
829 }
830
831
838 public function getImagEigenvalues() {
839 return $this->e;
840 }
841
842
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
hypo($a, $b)
Definition: Maths.php:14
global $l
Definition: afr.php:30
An exception for terminatinating execution or to throw for unit testing.
getRealEigenvalues()
Return the real parts of the eigenvalues.
cdiv($xr, $xi, $yr, $yi)
Performs complex division.
getV()
Return the eigenvector matrix.
orthes()
Nonsymmetric reduction to Hessenberg form.
hqr2()
Nonsymmetric reduction from Hessenberg to real Schur form.
getImagEigenvalues()
Return the imaginary parts of the eigenvalues.
tred2()
Symmetric Householder reduction to tridiagonal form.
__construct($Arg)
Constructor: Check for symmetry, then construct the eigenvalue decomposition.
tql2()
Symmetric tridiagonal QL algorithm.
getD()
Return the block diagonal eigenvalue matrix.
$y
Definition: example_007.php:83
$x
Definition: example_009.php:98
$h
$w
$r
Definition: example_031.php:79
e($cmd)
Definition: flush.php:14