ILIAS  release_5-3 Revision v5.3.23-19-g915713cf615
EigenvalueDecomposition.php
Go to the documentation of this file.
1<?php
24class EigenvalueDecomposition {
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
hypo($a, $b)
Definition: Maths.php:14
$l
global $l
Definition: afr.php:30
php
An exception for terminatinating execution or to throw for unit testing.
EigenvalueDecomposition\getRealEigenvalues
getRealEigenvalues()
Return the real parts of the eigenvalues.
Definition: EigenvalueDecomposition.php:827
EigenvalueDecomposition\cdiv
cdiv($xr, $xi, $yr, $yi)
Performs complex division.
Definition: EigenvalueDecomposition.php:373
EigenvalueDecomposition\getV
getV()
Return the eigenvector matrix.
Definition: EigenvalueDecomposition.php:816
EigenvalueDecomposition\orthes
orthes()
Nonsymmetric reduction to Hessenberg form.
Definition: EigenvalueDecomposition.php:291
EigenvalueDecomposition\hqr2
hqr2()
Nonsymmetric reduction from Hessenberg to real Schur form.
Definition: EigenvalueDecomposition.php:398
EigenvalueDecomposition\getImagEigenvalues
getImagEigenvalues()
Return the imaginary parts of the eigenvalues.
Definition: EigenvalueDecomposition.php:838
EigenvalueDecomposition\tred2
tred2()
Symmetric Householder reduction to tridiagonal form.
Definition: EigenvalueDecomposition.php:76
EigenvalueDecomposition\__construct
__construct($Arg)
Constructor: Check for symmetry, then construct the eigenvalue decomposition.
Definition: EigenvalueDecomposition.php:782
EigenvalueDecomposition\tql2
tql2()
Symmetric tridiagonal QL algorithm.
Definition: EigenvalueDecomposition.php:185
EigenvalueDecomposition\getD
getD()
Return the block diagonal eigenvalue matrix.
Definition: EigenvalueDecomposition.php:849
$i
$i
Definition: disco.tpl.php:19
$y
$y
Definition: example_007.php:83
$x
$x
Definition: example_009.php:98
$h
$h
Definition: example_009.php:101
$w
$w
Definition: example_009.php:100
$r
$r
Definition: example_031.php:79
e
e($cmd)
Definition: flush.php:14
$eps
$eps
Definition: metadata.php:61
$s
$s
Definition: pwgen.php:45
$m
$m
Definition: show_metadata.php:25