ILIAS  release_5-4 Revision v5.4.26-12-gabc799a52e6
EigenvalueDecomposition.php
Go to the documentation of this file.
1<?php
2
4
26{
32 private $n;
33
39 private $d = [];
40
41 private $e = [];
42
48 private $V = [];
49
55 private $H = [];
56
62 private $ort;
63
69 private $cdivr;
70
71 private $cdivi;
72
76 private $A;
77
81 private function tred2(): void
82 {
83 // This is derived from the Algol procedures tred2 by
84 // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
85 // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
86 // Fortran subroutine in EISPACK.
87 $this->d = $this->V[$this->n - 1];
88 $j = 0;
89 // Householder reduction to tridiagonal form.
90 for ($i = $this->n - 1; $i > 0; --$i) {
91 $i_ = $i - 1;
92 // Scale to avoid under/overflow.
93 $h = $scale = 0.0;
94 $scale += array_sum(array_map('abs', $this->d));
95 if ($scale == 0.0) {
96 $this->e[$i] = $this->d[$i_];
97 $this->d = array_slice($this->V[$i_], 0, $i_);
98 for ($j = 0; $j < $i; ++$j) {
99 $this->V[$j][$i] = $this->V[$i][$j] = 0.0;
100 }
101 } else {
102 // Generate Householder vector.
103 for ($k = 0; $k < $i; ++$k) {
104 $this->d[$k] /= $scale;
105 $h += $this->d[$k] ** 2;
106 }
107 $f = $this->d[$i_];
108 $g = sqrt($h);
109 if ($f > 0) {
110 $g = -$g;
111 }
112 $this->e[$i] = $scale * $g;
113 $h = $h - $f * $g;
114 $this->d[$i_] = $f - $g;
115 for ($j = 0; $j < $i; ++$j) {
116 $this->e[$j] = 0.0;
117 }
118 // Apply similarity transformation to remaining columns.
119 for ($j = 0; $j < $i; ++$j) {
120 $f = $this->d[$j];
121 $this->V[$j][$i] = $f;
122 $g = $this->e[$j] + $this->V[$j][$j] * $f;
123 for ($k = $j + 1; $k <= $i_; ++$k) {
124 $g += $this->V[$k][$j] * $this->d[$k];
125 $this->e[$k] += $this->V[$k][$j] * $f;
126 }
127 $this->e[$j] = $g;
128 }
129 $f = 0.0;
130 for ($j = 0; $j < $i; ++$j) {
131 $this->e[$j] /= $h;
132 $f += $this->e[$j] * $this->d[$j];
133 }
134 $hh = $f / (2 * $h);
135 for ($j = 0; $j < $i; ++$j) {
136 $this->e[$j] -= $hh * $this->d[$j];
137 }
138 for ($j = 0; $j < $i; ++$j) {
139 $f = $this->d[$j];
140 $g = $this->e[$j];
141 for ($k = $j; $k <= $i_; ++$k) {
142 $this->V[$k][$j] -= ($f * $this->e[$k] + $g * $this->d[$k]);
143 }
144 $this->d[$j] = $this->V[$i - 1][$j];
145 $this->V[$i][$j] = 0.0;
146 }
147 }
148 $this->d[$i] = $h;
149 }
150
151 // Accumulate transformations.
152 for ($i = 0; $i < $this->n - 1; ++$i) {
153 $this->V[$this->n - 1][$i] = $this->V[$i][$i];
154 $this->V[$i][$i] = 1.0;
155 $h = $this->d[$i + 1];
156 if ($h != 0.0) {
157 for ($k = 0; $k <= $i; ++$k) {
158 $this->d[$k] = $this->V[$k][$i + 1] / $h;
159 }
160 for ($j = 0; $j <= $i; ++$j) {
161 $g = 0.0;
162 for ($k = 0; $k <= $i; ++$k) {
163 $g += $this->V[$k][$i + 1] * $this->V[$k][$j];
164 }
165 for ($k = 0; $k <= $i; ++$k) {
166 $this->V[$k][$j] -= $g * $this->d[$k];
167 }
168 }
169 }
170 for ($k = 0; $k <= $i; ++$k) {
171 $this->V[$k][$i + 1] = 0.0;
172 }
173 }
174
175 $this->d = $this->V[$this->n - 1];
176 $this->V[$this->n - 1] = array_fill(0, $j, 0.0);
177 $this->V[$this->n - 1][$this->n - 1] = 1.0;
178 $this->e[0] = 0.0;
179 }
180
189 private function tql2(): void
190 {
191 for ($i = 1; $i < $this->n; ++$i) {
192 $this->e[$i - 1] = $this->e[$i];
193 }
194 $this->e[$this->n - 1] = 0.0;
195 $f = 0.0;
196 $tst1 = 0.0;
197 $eps = 2.0 ** (-52.0);
198
199 for ($l = 0; $l < $this->n; ++$l) {
200 // Find small subdiagonal element
201 $tst1 = max($tst1, abs($this->d[$l]) + abs($this->e[$l]));
202 $m = $l;
203 while ($m < $this->n) {
204 if (abs($this->e[$m]) <= $eps * $tst1) {
205 break;
206 }
207 ++$m;
208 }
209 // If m == l, $this->d[l] is an eigenvalue,
210 // otherwise, iterate.
211 if ($m > $l) {
212 $iter = 0;
213 do {
214 // Could check iteration count here.
215 ++$iter;
216 // Compute implicit shift
217 $g = $this->d[$l];
218 $p = ($this->d[$l + 1] - $g) / (2.0 * $this->e[$l]);
219 $r = hypo($p, 1.0);
220 if ($p < 0) {
221 $r *= -1;
222 }
223 $this->d[$l] = $this->e[$l] / ($p + $r);
224 $this->d[$l + 1] = $this->e[$l] * ($p + $r);
225 $dl1 = $this->d[$l + 1];
226 $h = $g - $this->d[$l];
227 for ($i = $l + 2; $i < $this->n; ++$i) {
228 $this->d[$i] -= $h;
229 }
230 $f += $h;
231 // Implicit QL transformation.
232 $p = $this->d[$m];
233 $c = 1.0;
234 $c2 = $c3 = $c;
235 $el1 = $this->e[$l + 1];
236 $s = $s2 = 0.0;
237 for ($i = $m - 1; $i >= $l; --$i) {
238 $c3 = $c2;
239 $c2 = $c;
240 $s2 = $s;
241 $g = $c * $this->e[$i];
242 $h = $c * $p;
243 $r = hypo($p, $this->e[$i]);
244 $this->e[$i + 1] = $s * $r;
245 $s = $this->e[$i] / $r;
246 $c = $p / $r;
247 $p = $c * $this->d[$i] - $s * $g;
248 $this->d[$i + 1] = $h + $s * ($c * $g + $s * $this->d[$i]);
249 // Accumulate transformation.
250 for ($k = 0; $k < $this->n; ++$k) {
251 $h = $this->V[$k][$i + 1];
252 $this->V[$k][$i + 1] = $s * $this->V[$k][$i] + $c * $h;
253 $this->V[$k][$i] = $c * $this->V[$k][$i] - $s * $h;
254 }
255 }
256 $p = -$s * $s2 * $c3 * $el1 * $this->e[$l] / $dl1;
257 $this->e[$l] = $s * $p;
258 $this->d[$l] = $c * $p;
259 // Check for convergence.
260 } while (abs($this->e[$l]) > $eps * $tst1);
261 }
262 $this->d[$l] = $this->d[$l] + $f;
263 $this->e[$l] = 0.0;
264 }
265
266 // Sort eigenvalues and corresponding vectors.
267 for ($i = 0; $i < $this->n - 1; ++$i) {
268 $k = $i;
269 $p = $this->d[$i];
270 for ($j = $i + 1; $j < $this->n; ++$j) {
271 if ($this->d[$j] < $p) {
272 $k = $j;
273 $p = $this->d[$j];
274 }
275 }
276 if ($k != $i) {
277 $this->d[$k] = $this->d[$i];
278 $this->d[$i] = $p;
279 for ($j = 0; $j < $this->n; ++$j) {
280 $p = $this->V[$j][$i];
281 $this->V[$j][$i] = $this->V[$j][$k];
282 $this->V[$j][$k] = $p;
283 }
284 }
285 }
286 }
287
296 private function orthes(): void
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
381 private function cdiv($xr, $xi, $yr, $yi): void
382 {
383 if (abs($yr) > abs($yi)) {
384 $r = $yi / $yr;
385 $d = $yr + $r * $yi;
386 $this->cdivr = ($xr + $r * $xi) / $d;
387 $this->cdivi = ($xi - $r * $xr) / $d;
388 } else {
389 $r = $yr / $yi;
390 $d = $yi + $r * $yr;
391 $this->cdivr = ($r * $xr + $xi) / $d;
392 $this->cdivi = ($r * $xi - $xr) / $d;
393 }
394 }
395
404 private function hqr2(): void
405 {
406 // Initialize
407 $nn = $this->n;
408 $n = $nn - 1;
409 $low = 0;
410 $high = $nn - 1;
411 $eps = 2.0 ** (-52.0);
412 $exshift = 0.0;
413 $p = $q = $r = $s = $z = 0;
414 // Store roots isolated by balanc and compute matrix norm
415 $norm = 0.0;
416
417 for ($i = 0; $i < $nn; ++$i) {
418 if (($i < $low) || ($i > $high)) {
419 $this->d[$i] = $this->H[$i][$i];
420 $this->e[$i] = 0.0;
421 }
422 for ($j = max($i - 1, 0); $j < $nn; ++$j) {
423 $norm = $norm + abs($this->H[$i][$j]);
424 }
425 }
426
427 // Outer loop over eigenvalue index
428 $iter = 0;
429 while ($n >= $low) {
430 // Look for single small sub-diagonal element
431 $l = $n;
432 while ($l > $low) {
433 $s = abs($this->H[$l - 1][$l - 1]) + abs($this->H[$l][$l]);
434 if ($s == 0.0) {
435 $s = $norm;
436 }
437 if (abs($this->H[$l][$l - 1]) < $eps * $s) {
438 break;
439 }
440 --$l;
441 }
442 // Check for convergence
443 // One root found
444 if ($l == $n) {
445 $this->H[$n][$n] = $this->H[$n][$n] + $exshift;
446 $this->d[$n] = $this->H[$n][$n];
447 $this->e[$n] = 0.0;
448 --$n;
449 $iter = 0;
450 // Two roots found
451 } elseif ($l == $n - 1) {
452 $w = $this->H[$n][$n - 1] * $this->H[$n - 1][$n];
453 $p = ($this->H[$n - 1][$n - 1] - $this->H[$n][$n]) / 2.0;
454 $q = $p * $p + $w;
455 $z = sqrt(abs($q));
456 $this->H[$n][$n] = $this->H[$n][$n] + $exshift;
457 $this->H[$n - 1][$n - 1] = $this->H[$n - 1][$n - 1] + $exshift;
458 $x = $this->H[$n][$n];
459 // Real pair
460 if ($q >= 0) {
461 if ($p >= 0) {
462 $z = $p + $z;
463 } else {
464 $z = $p - $z;
465 }
466 $this->d[$n - 1] = $x + $z;
467 $this->d[$n] = $this->d[$n - 1];
468 if ($z != 0.0) {
469 $this->d[$n] = $x - $w / $z;
470 }
471 $this->e[$n - 1] = 0.0;
472 $this->e[$n] = 0.0;
473 $x = $this->H[$n][$n - 1];
474 $s = abs($x) + abs($z);
475 $p = $x / $s;
476 $q = $z / $s;
477 $r = sqrt($p * $p + $q * $q);
478 $p = $p / $r;
479 $q = $q / $r;
480 // Row modification
481 for ($j = $n - 1; $j < $nn; ++$j) {
482 $z = $this->H[$n - 1][$j];
483 $this->H[$n - 1][$j] = $q * $z + $p * $this->H[$n][$j];
484 $this->H[$n][$j] = $q * $this->H[$n][$j] - $p * $z;
485 }
486 // Column modification
487 for ($i = 0; $i <= $n; ++$i) {
488 $z = $this->H[$i][$n - 1];
489 $this->H[$i][$n - 1] = $q * $z + $p * $this->H[$i][$n];
490 $this->H[$i][$n] = $q * $this->H[$i][$n] - $p * $z;
491 }
492 // Accumulate transformations
493 for ($i = $low; $i <= $high; ++$i) {
494 $z = $this->V[$i][$n - 1];
495 $this->V[$i][$n - 1] = $q * $z + $p * $this->V[$i][$n];
496 $this->V[$i][$n] = $q * $this->V[$i][$n] - $p * $z;
497 }
498 // Complex pair
499 } else {
500 $this->d[$n - 1] = $x + $p;
501 $this->d[$n] = $x + $p;
502 $this->e[$n - 1] = $z;
503 $this->e[$n] = -$z;
504 }
505 $n = $n - 2;
506 $iter = 0;
507 // No convergence yet
508 } else {
509 // Form shift
510 $x = $this->H[$n][$n];
511 $y = 0.0;
512 $w = 0.0;
513 if ($l < $n) {
514 $y = $this->H[$n - 1][$n - 1];
515 $w = $this->H[$n][$n - 1] * $this->H[$n - 1][$n];
516 }
517 // Wilkinson's original ad hoc shift
518 if ($iter == 10) {
519 $exshift += $x;
520 for ($i = $low; $i <= $n; ++$i) {
521 $this->H[$i][$i] -= $x;
522 }
523 $s = abs($this->H[$n][$n - 1]) + abs($this->H[$n - 1][$n - 2]);
524 $x = $y = 0.75 * $s;
525 $w = -0.4375 * $s * $s;
526 }
527 // MATLAB's new ad hoc shift
528 if ($iter == 30) {
529 $s = ($y - $x) / 2.0;
530 $s = $s * $s + $w;
531 if ($s > 0) {
532 $s = sqrt($s);
533 if ($y < $x) {
534 $s = -$s;
535 }
536 $s = $x - $w / (($y - $x) / 2.0 + $s);
537 for ($i = $low; $i <= $n; ++$i) {
538 $this->H[$i][$i] -= $s;
539 }
540 $exshift += $s;
541 $x = $y = $w = 0.964;
542 }
543 }
544 // Could check iteration count here.
545 $iter = $iter + 1;
546 // Look for two consecutive small sub-diagonal elements
547 $m = $n - 2;
548 while ($m >= $l) {
549 $z = $this->H[$m][$m];
550 $r = $x - $z;
551 $s = $y - $z;
552 $p = ($r * $s - $w) / $this->H[$m + 1][$m] + $this->H[$m][$m + 1];
553 $q = $this->H[$m + 1][$m + 1] - $z - $r - $s;
554 $r = $this->H[$m + 2][$m + 1];
555 $s = abs($p) + abs($q) + abs($r);
556 $p = $p / $s;
557 $q = $q / $s;
558 $r = $r / $s;
559 if ($m == $l) {
560 break;
561 }
562 if (
563 abs($this->H[$m][$m - 1]) * (abs($q) + abs($r)) <
564 $eps * (abs($p) * (abs($this->H[$m - 1][$m - 1]) + abs($z) + abs($this->H[$m + 1][$m + 1])))
565 ) {
566 break;
567 }
568 --$m;
569 }
570 for ($i = $m + 2; $i <= $n; ++$i) {
571 $this->H[$i][$i - 2] = 0.0;
572 if ($i > $m + 2) {
573 $this->H[$i][$i - 3] = 0.0;
574 }
575 }
576 // Double QR step involving rows l:n and columns m:n
577 for ($k = $m; $k <= $n - 1; ++$k) {
578 $notlast = ($k != $n - 1);
579 if ($k != $m) {
580 $p = $this->H[$k][$k - 1];
581 $q = $this->H[$k + 1][$k - 1];
582 $r = ($notlast ? $this->H[$k + 2][$k - 1] : 0.0);
583 $x = abs($p) + abs($q) + abs($r);
584 if ($x != 0.0) {
585 $p = $p / $x;
586 $q = $q / $x;
587 $r = $r / $x;
588 }
589 }
590 if ($x == 0.0) {
591 break;
592 }
593 $s = sqrt($p * $p + $q * $q + $r * $r);
594 if ($p < 0) {
595 $s = -$s;
596 }
597 if ($s != 0) {
598 if ($k != $m) {
599 $this->H[$k][$k - 1] = -$s * $x;
600 } elseif ($l != $m) {
601 $this->H[$k][$k - 1] = -$this->H[$k][$k - 1];
602 }
603 $p = $p + $s;
604 $x = $p / $s;
605 $y = $q / $s;
606 $z = $r / $s;
607 $q = $q / $p;
608 $r = $r / $p;
609 // Row modification
610 for ($j = $k; $j < $nn; ++$j) {
611 $p = $this->H[$k][$j] + $q * $this->H[$k + 1][$j];
612 if ($notlast) {
613 $p = $p + $r * $this->H[$k + 2][$j];
614 $this->H[$k + 2][$j] = $this->H[$k + 2][$j] - $p * $z;
615 }
616 $this->H[$k][$j] = $this->H[$k][$j] - $p * $x;
617 $this->H[$k + 1][$j] = $this->H[$k + 1][$j] - $p * $y;
618 }
619 // Column modification
620 $iMax = min($n, $k + 3);
621 for ($i = 0; $i <= $iMax; ++$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 $kMax = min($j, $high);
777 for ($k = $low; $k <= $kMax; ++$k) {
778 $z = $z + $this->V[$i][$k] * $this->H[$k][$j];
779 }
780 $this->V[$i][$j] = $z;
781 }
782 }
783 }
784
785 // end hqr2
786
792 public function __construct(Matrix $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
825 public function getV()
826 {
827 return new Matrix($this->V, $this->n, $this->n);
828 }
829
835 public function getRealEigenvalues()
836 {
837 return $this->d;
838 }
839
845 public function getImagEigenvalues()
846 {
847 return $this->e;
848 }
849
855 public function getD()
856 {
857 $D = [];
858 for ($i = 0; $i < $this->n; ++$i) {
859 $D[$i] = array_fill(0, $this->n, 0.0);
860 $D[$i][$i] = $this->d[$i];
861 if ($this->e[$i] == 0) {
862 continue;
863 }
864 $o = ($this->e[$i] > 0) ? $i + 1 : $i - 1;
865 $D[$i][$o] = $this->e[$i];
866 }
867
868 return new Matrix($D);
869 }
870}
hypo($a, $b)
Pythagorean Theorem:.
Definition: Maths.php:18
global $l
Definition: afr.php:30
An exception for terminatinating execution or to throw for unit testing.
Class to obtain eigenvalues and eigenvectors of a real matrix.
tred2()
Symmetric Householder reduction to tridiagonal form.
getImagEigenvalues()
Return the imaginary parts of the eigenvalues.
hqr2()
Nonsymmetric reduction from Hessenberg to real Schur form.
__construct(Matrix $Arg)
Constructor: Check for symmetry, then construct the eigenvalue decomposition.
getColumnDimension()
getColumnDimension.
Definition: Matrix.php:137
$x
Definition: complexTest.php:9
$i
Definition: disco.tpl.php:19
$y
Definition: example_007.php:83
$h
$w
$r
Definition: example_031.php:79
$eps
Definition: metadata.php:61
Class for the creating "special" Matrices.
Definition: Builder.php:11
$s
Definition: pwgen.php:45