3 namespace PhpOffice\PhpSpreadsheet\Shared\JAMA
;
6 * Class to obtain eigenvalues and eigenvectors of a real matrix.
8 * If A is symmetric, then A = V*D*V' where the eigenvalue matrix D
9 * is diagonal and the eigenvector matrix V is orthogonal (i.e.
10 * A = V.times(D.times(V.transpose())) and V.times(V.transpose())
11 * equals the identity matrix).
13 * If A is not symmetric, then the eigenvalue matrix D is block diagonal
14 * with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues,
15 * lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda]. The
16 * columns of V represent the eigenvectors in the sense that A*V = V*D,
17 * i.e. A.times(V) equals V.times(D). The matrix V may be badly
18 * conditioned, or even singular, so the validity of the equation
19 * A = V*D*inverse(V) depends upon V.cond().
21 * @author Paul Meagher
25 class EigenvalueDecomposition
28 * Row and column dimension (square matrix).
35 * Arrays for internal storage of eigenvalues.
44 * Array for internal storage of eigenvectors.
51 * Array for internal storage of nonsymmetric Hessenberg form.
58 * Working storage for nonsymmetric algorithm.
65 * Used for complex scalar division.
74 * Symmetric Householder reduction to tridiagonal form.
76 private function tred2()
78 // This is derived from the Algol procedures tred2 by
79 // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
80 // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
81 // Fortran subroutine in EISPACK.
82 $this->d
= $this->V
[$this->n
- 1];
83 // Householder reduction to tridiagonal form.
84 for ($i = $this->n
- 1; $i > 0; --$i) {
86 // Scale to avoid under/overflow.
88 $scale +
= array_sum(array_map('abs', $this->d
));
90 $this->e
[$i] = $this->d
[$i_];
91 $this->d
= array_slice($this->V
[$i_], 0, $i_);
92 for ($j = 0; $j < $i; ++
$j) {
93 $this->V
[$j][$i] = $this->V
[$i][$j] = 0.0;
96 // Generate Householder vector.
97 for ($k = 0; $k < $i; ++
$k) {
98 $this->d
[$k] /= $scale;
99 $h +
= pow($this->d
[$k], 2);
106 $this->e
[$i] = $scale * $g;
108 $this->d
[$i_] = $f - $g;
109 for ($j = 0; $j < $i; ++
$j) {
112 // Apply similarity transformation to remaining columns.
113 for ($j = 0; $j < $i; ++
$j) {
115 $this->V
[$j][$i] = $f;
116 $g = $this->e
[$j] +
$this->V
[$j][$j] * $f;
117 for ($k = $j +
1; $k <= $i_; ++
$k) {
118 $g +
= $this->V
[$k][$j] * $this->d
[$k];
119 $this->e
[$k] +
= $this->V
[$k][$j] * $f;
124 for ($j = 0; $j < $i; ++
$j) {
126 $f +
= $this->e
[$j] * $this->d
[$j];
129 for ($j = 0; $j < $i; ++
$j) {
130 $this->e
[$j] -= $hh * $this->d
[$j];
132 for ($j = 0; $j < $i; ++
$j) {
135 for ($k = $j; $k <= $i_; ++
$k) {
136 $this->V
[$k][$j] -= ($f * $this->e
[$k] +
$g * $this->d
[$k]);
138 $this->d
[$j] = $this->V
[$i - 1][$j];
139 $this->V
[$i][$j] = 0.0;
145 // Accumulate transformations.
146 for ($i = 0; $i < $this->n
- 1; ++
$i) {
147 $this->V
[$this->n
- 1][$i] = $this->V
[$i][$i];
148 $this->V
[$i][$i] = 1.0;
149 $h = $this->d
[$i +
1];
151 for ($k = 0; $k <= $i; ++
$k) {
152 $this->d
[$k] = $this->V
[$k][$i +
1] / $h;
154 for ($j = 0; $j <= $i; ++
$j) {
156 for ($k = 0; $k <= $i; ++
$k) {
157 $g +
= $this->V
[$k][$i +
1] * $this->V
[$k][$j];
159 for ($k = 0; $k <= $i; ++
$k) {
160 $this->V
[$k][$j] -= $g * $this->d
[$k];
164 for ($k = 0; $k <= $i; ++
$k) {
165 $this->V
[$k][$i +
1] = 0.0;
169 $this->d
= $this->V
[$this->n
- 1];
170 $this->V
[$this->n
- 1] = array_fill(0, $j, 0.0);
171 $this->V
[$this->n
- 1][$this->n
- 1] = 1.0;
176 * Symmetric tridiagonal QL algorithm.
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.
183 private function tql2()
185 for ($i = 1; $i < $this->n
; ++
$i) {
186 $this->e
[$i - 1] = $this->e
[$i];
188 $this->e
[$this->n
- 1] = 0.0;
191 $eps = pow(2.0, -52.0);
193 for ($l = 0; $l < $this->n
; ++
$l) {
194 // Find small subdiagonal element
195 $tst1 = max($tst1, abs($this->d
[$l]) +
abs($this->e
[$l]));
197 while ($m < $this->n
) {
198 if (abs($this->e
[$m]) <= $eps * $tst1) {
203 // If m == l, $this->d[l] is an eigenvalue,
204 // otherwise, iterate.
208 // Could check iteration count here.
210 // Compute implicit shift
212 $p = ($this->d
[$l +
1] - $g) / (2.0 * $this->e
[$l]);
217 $this->d
[$l] = $this->e
[$l] / ($p +
$r);
218 $this->d
[$l +
1] = $this->e
[$l] * ($p +
$r);
219 $dl1 = $this->d
[$l +
1];
220 $h = $g - $this->d
[$l];
221 for ($i = $l +
2; $i < $this->n
; ++
$i) {
225 // Implicit QL transformation.
229 $el1 = $this->e
[$l +
1];
231 for ($i = $m - 1; $i >= $l; --$i) {
235 $g = $c * $this->e
[$i];
237 $r = hypo($p, $this->e
[$i]);
238 $this->e
[$i +
1] = $s * $r;
239 $s = $this->e
[$i] / $r;
241 $p = $c * $this->d
[$i] - $s * $g;
242 $this->d
[$i +
1] = $h +
$s * ($c * $g +
$s * $this->d
[$i]);
243 // Accumulate transformation.
244 for ($k = 0; $k < $this->n
; ++
$k) {
245 $h = $this->V
[$k][$i +
1];
246 $this->V
[$k][$i +
1] = $s * $this->V
[$k][$i] +
$c * $h;
247 $this->V
[$k][$i] = $c * $this->V
[$k][$i] - $s * $h;
250 $p = -$s * $s2 * $c3 * $el1 * $this->e
[$l] / $dl1;
251 $this->e
[$l] = $s * $p;
252 $this->d
[$l] = $c * $p;
253 // Check for convergence.
254 } while (abs($this->e
[$l]) > $eps * $tst1);
256 $this->d
[$l] = $this->d
[$l] +
$f;
260 // Sort eigenvalues and corresponding vectors.
261 for ($i = 0; $i < $this->n
- 1; ++
$i) {
264 for ($j = $i +
1; $j < $this->n
; ++
$j) {
265 if ($this->d
[$j] < $p) {
271 $this->d
[$k] = $this->d
[$i];
273 for ($j = 0; $j < $this->n
; ++
$j) {
274 $p = $this->V
[$j][$i];
275 $this->V
[$j][$i] = $this->V
[$j][$k];
276 $this->V
[$j][$k] = $p;
283 * Nonsymmetric reduction to Hessenberg form.
285 * This is derived from the Algol procedures orthes and ortran,
286 * by Martin and Wilkinson, Handbook for Auto. Comp.,
287 * Vol.ii-Linear Algebra, and the corresponding
288 * Fortran subroutines in EISPACK.
290 private function orthes()
293 $high = $this->n
- 1;
295 for ($m = $low +
1; $m <= $high - 1; ++
$m) {
298 for ($i = $m; $i <= $high; ++
$i) {
299 $scale = $scale +
abs($this->H
[$i][$m - 1]);
302 // Compute Householder transformation.
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];
309 if ($this->ort
[$m] > 0) {
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) {
318 for ($i = $high; $i >= $m; --$i) {
319 $f +
= $this->ort
[$i] * $this->H
[$i][$j];
322 for ($i = $m; $i <= $high; ++
$i) {
323 $this->H
[$i][$j] -= $f * $this->ort
[$i];
326 for ($i = 0; $i <= $high; ++
$i) {
328 for ($j = $high; $j >= $m; --$j) {
329 $f +
= $this->ort
[$j] * $this->H
[$i][$j];
332 for ($j = $m; $j <= $high; ++
$j) {
333 $this->H
[$i][$j] -= $f * $this->ort
[$j];
336 $this->ort
[$m] = $scale * $this->ort
[$m];
337 $this->H
[$m][$m - 1] = $scale * $g;
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);
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];
352 for ($j = $m; $j <= $high; ++
$j) {
354 for ($i = $m; $i <= $high; ++
$i) {
355 $g +
= $this->ort
[$i] * $this->V
[$i][$j];
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];
368 * Performs complex division.
375 private function cdiv($xr, $xi, $yr, $yi)
377 if (abs($yr) > abs($yi)) {
380 $this->cdivr
= ($xr +
$r * $xi) / $d;
381 $this->cdivi
= ($xi - $r * $xr) / $d;
385 $this->cdivr
= ($r * $xr +
$xi) / $d;
386 $this->cdivi
= ($r * $xi - $xr) / $d;
391 * Nonsymmetric reduction from Hessenberg to real Schur form.
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.
398 private function hqr2()
405 $eps = pow(2.0, -52.0);
407 $p = $q = $r = $s = $z = 0;
408 // Store roots isolated by balanc and compute matrix norm
411 for ($i = 0; $i < $nn; ++
$i) {
412 if (($i < $low) or ($i > $high)) {
413 $this->d
[$i] = $this->H
[$i][$i];
416 for ($j = max($i - 1, 0); $j < $nn; ++
$j) {
417 $norm = $norm +
abs($this->H
[$i][$j]);
421 // Outer loop over eigenvalue index
424 // Look for single small sub-diagonal element
427 $s = abs($this->H
[$l - 1][$l - 1]) +
abs($this->H
[$l][$l]);
431 if (abs($this->H
[$l][$l - 1]) < $eps * $s) {
436 // Check for convergence
439 $this->H
[$n][$n] = $this->H
[$n][$n] +
$exshift;
440 $this->d
[$n] = $this->H
[$n][$n];
445 } elseif ($l == $n - 1) {
446 $w = $this->H
[$n][$n - 1] * $this->H
[$n - 1][$n];
447 $p = ($this->H
[$n - 1][$n - 1] - $this->H
[$n][$n]) / 2.0;
450 $this->H
[$n][$n] = $this->H
[$n][$n] +
$exshift;
451 $this->H
[$n - 1][$n - 1] = $this->H
[$n - 1][$n - 1] +
$exshift;
452 $x = $this->H
[$n][$n];
460 $this->d
[$n - 1] = $x +
$z;
461 $this->d
[$n] = $this->d
[$n - 1];
463 $this->d
[$n] = $x - $w / $z;
465 $this->e
[$n - 1] = 0.0;
467 $x = $this->H
[$n][$n - 1];
468 $s = abs($x) +
abs($z);
471 $r = sqrt($p * $p +
$q * $q);
475 for ($j = $n - 1; $j < $nn; ++
$j) {
476 $z = $this->H
[$n - 1][$j];
477 $this->H
[$n - 1][$j] = $q * $z +
$p * $this->H
[$n][$j];
478 $this->H
[$n][$j] = $q * $this->H
[$n][$j] - $p * $z;
480 // Column modification
481 for ($i = 0; $i <= $n; ++
$i) {
482 $z = $this->H
[$i][$n - 1];
483 $this->H
[$i][$n - 1] = $q * $z +
$p * $this->H
[$i][$n];
484 $this->H
[$i][$n] = $q * $this->H
[$i][$n] - $p * $z;
486 // Accumulate transformations
487 for ($i = $low; $i <= $high; ++
$i) {
488 $z = $this->V
[$i][$n - 1];
489 $this->V
[$i][$n - 1] = $q * $z +
$p * $this->V
[$i][$n];
490 $this->V
[$i][$n] = $q * $this->V
[$i][$n] - $p * $z;
494 $this->d
[$n - 1] = $x +
$p;
495 $this->d
[$n] = $x +
$p;
496 $this->e
[$n - 1] = $z;
501 // No convergence yet
504 $x = $this->H
[$n][$n];
508 $y = $this->H
[$n - 1][$n - 1];
509 $w = $this->H
[$n][$n - 1] * $this->H
[$n - 1][$n];
511 // Wilkinson's original ad hoc shift
514 for ($i = $low; $i <= $n; ++
$i) {
515 $this->H
[$i][$i] -= $x;
517 $s = abs($this->H
[$n][$n - 1]) +
abs($this->H
[$n - 1][$n - 2]);
519 $w = -0.4375 * $s * $s;
521 // MATLAB's new ad hoc shift
523 $s = ($y - $x) / 2.0;
530 $s = $x - $w / (($y - $x) / 2.0 +
$s);
531 for ($i = $low; $i <= $n; ++
$i) {
532 $this->H
[$i][$i] -= $s;
535 $x = $y = $w = 0.964;
538 // Could check iteration count here.
540 // Look for two consecutive small sub-diagonal elements
543 $z = $this->H
[$m][$m];
546 $p = ($r * $s - $w) / $this->H
[$m +
1][$m] +
$this->H
[$m][$m +
1];
547 $q = $this->H
[$m +
1][$m +
1] - $z - $r - $s;
548 $r = $this->H
[$m +
2][$m +
1];
549 $s = abs($p) +
abs($q) +
abs($r);
556 if (abs($this->H
[$m][$m - 1]) * (abs($q) +
abs($r)) <
557 $eps * (abs($p) * (abs($this->H
[$m - 1][$m - 1]) +
abs($z) +
abs($this->H
[$m +
1][$m +
1])))) {
562 for ($i = $m +
2; $i <= $n; ++
$i) {
563 $this->H
[$i][$i - 2] = 0.0;
565 $this->H
[$i][$i - 3] = 0.0;
568 // Double QR step involving rows l:n and columns m:n
569 for ($k = $m; $k <= $n - 1; ++
$k) {
570 $notlast = ($k != $n - 1);
572 $p = $this->H
[$k][$k - 1];
573 $q = $this->H
[$k +
1][$k - 1];
574 $r = ($notlast ?
$this->H
[$k +
2][$k - 1] : 0.0);
575 $x = abs($p) +
abs($q) +
abs($r);
585 $s = sqrt($p * $p +
$q * $q +
$r * $r);
591 $this->H
[$k][$k - 1] = -$s * $x;
592 } elseif ($l != $m) {
593 $this->H
[$k][$k - 1] = -$this->H
[$k][$k - 1];
602 for ($j = $k; $j < $nn; ++
$j) {
603 $p = $this->H
[$k][$j] +
$q * $this->H
[$k +
1][$j];
605 $p = $p +
$r * $this->H
[$k +
2][$j];
606 $this->H
[$k +
2][$j] = $this->H
[$k +
2][$j] - $p * $z;
608 $this->H
[$k][$j] = $this->H
[$k][$j] - $p * $x;
609 $this->H
[$k +
1][$j] = $this->H
[$k +
1][$j] - $p * $y;
611 // Column modification
612 $iMax = min($n, $k +
3);
613 for ($i = 0; $i <= $iMax; ++
$i) {
614 $p = $x * $this->H
[$i][$k] +
$y * $this->H
[$i][$k +
1];
616 $p = $p +
$z * $this->H
[$i][$k +
2];
617 $this->H
[$i][$k +
2] = $this->H
[$i][$k +
2] - $p * $r;
619 $this->H
[$i][$k] = $this->H
[$i][$k] - $p;
620 $this->H
[$i][$k +
1] = $this->H
[$i][$k +
1] - $p * $q;
622 // Accumulate transformations
623 for ($i = $low; $i <= $high; ++
$i) {
624 $p = $x * $this->V
[$i][$k] +
$y * $this->V
[$i][$k +
1];
626 $p = $p +
$z * $this->V
[$i][$k +
2];
627 $this->V
[$i][$k +
2] = $this->V
[$i][$k +
2] - $p * $r;
629 $this->V
[$i][$k] = $this->V
[$i][$k] - $p;
630 $this->V
[$i][$k +
1] = $this->V
[$i][$k +
1] - $p * $q;
634 } // check convergence
635 } // while ($n >= $low)
637 // Backsubstitute to find vectors of upper triangular form
642 for ($n = $nn - 1; $n >= 0; --$n) {
648 $this->H
[$n][$n] = 1.0;
649 for ($i = $n - 1; $i >= 0; --$i) {
650 $w = $this->H
[$i][$i] - $p;
652 for ($j = $l; $j <= $n; ++
$j) {
653 $r = $r +
$this->H
[$i][$j] * $this->H
[$j][$n];
655 if ($this->e
[$i] < 0.0) {
660 if ($this->e
[$i] == 0.0) {
662 $this->H
[$i][$n] = -$r / $w;
664 $this->H
[$i][$n] = -$r / ($eps * $norm);
666 // Solve real equations
668 $x = $this->H
[$i][$i +
1];
669 $y = $this->H
[$i +
1][$i];
670 $q = ($this->d
[$i] - $p) * ($this->d
[$i] - $p) +
$this->e
[$i] * $this->e
[$i];
671 $t = ($x * $s - $z * $r) / $q;
672 $this->H
[$i][$n] = $t;
673 if (abs($x) > abs($z)) {
674 $this->H
[$i +
1][$n] = (-$r - $w * $t) / $x;
676 $this->H
[$i +
1][$n] = (-$s - $y * $t) / $z;
680 $t = abs($this->H
[$i][$n]);
681 if (($eps * $t) * $t > 1) {
682 for ($j = $i; $j <= $n; ++
$j) {
683 $this->H
[$j][$n] = $this->H
[$j][$n] / $t;
691 // Last vector component imaginary so matrix is triangular
692 if (abs($this->H
[$n][$n - 1]) > abs($this->H
[$n - 1][$n])) {
693 $this->H
[$n - 1][$n - 1] = $q / $this->H
[$n][$n - 1];
694 $this->H
[$n - 1][$n] = -($this->H
[$n][$n] - $p) / $this->H
[$n][$n - 1];
696 $this->cdiv(0.0, -$this->H
[$n - 1][$n], $this->H
[$n - 1][$n - 1] - $p, $q);
697 $this->H
[$n - 1][$n - 1] = $this->cdivr
;
698 $this->H
[$n - 1][$n] = $this->cdivi
;
700 $this->H
[$n][$n - 1] = 0.0;
701 $this->H
[$n][$n] = 1.0;
702 for ($i = $n - 2; $i >= 0; --$i) {
703 // double ra,sa,vr,vi;
706 for ($j = $l; $j <= $n; ++
$j) {
707 $ra = $ra +
$this->H
[$i][$j] * $this->H
[$j][$n - 1];
708 $sa = $sa +
$this->H
[$i][$j] * $this->H
[$j][$n];
710 $w = $this->H
[$i][$i] - $p;
711 if ($this->e
[$i] < 0.0) {
717 if ($this->e
[$i] == 0) {
718 $this->cdiv(-$ra, -$sa, $w, $q);
719 $this->H
[$i][$n - 1] = $this->cdivr
;
720 $this->H
[$i][$n] = $this->cdivi
;
722 // Solve complex equations
723 $x = $this->H
[$i][$i +
1];
724 $y = $this->H
[$i +
1][$i];
725 $vr = ($this->d
[$i] - $p) * ($this->d
[$i] - $p) +
$this->e
[$i] * $this->e
[$i] - $q * $q;
726 $vi = ($this->d
[$i] - $p) * 2.0 * $q;
727 if ($vr == 0.0 & $vi == 0.0) {
728 $vr = $eps * $norm * (abs($w) +
abs($q) +
abs($x) +
abs($y) +
abs($z));
730 $this->cdiv($x * $r - $z * $ra +
$q * $sa, $x * $s - $z * $sa - $q * $ra, $vr, $vi);
731 $this->H
[$i][$n - 1] = $this->cdivr
;
732 $this->H
[$i][$n] = $this->cdivi
;
733 if (abs($x) > (abs($z) +
abs($q))) {
734 $this->H
[$i +
1][$n - 1] = (-$ra - $w * $this->H
[$i][$n - 1] +
$q * $this->H
[$i][$n]) / $x;
735 $this->H
[$i +
1][$n] = (-$sa - $w * $this->H
[$i][$n] - $q * $this->H
[$i][$n - 1]) / $x;
737 $this->cdiv(-$r - $y * $this->H
[$i][$n - 1], -$s - $y * $this->H
[$i][$n], $z, $q);
738 $this->H
[$i +
1][$n - 1] = $this->cdivr
;
739 $this->H
[$i +
1][$n] = $this->cdivi
;
743 $t = max(abs($this->H
[$i][$n - 1]), abs($this->H
[$i][$n]));
744 if (($eps * $t) * $t > 1) {
745 for ($j = $i; $j <= $n; ++
$j) {
746 $this->H
[$j][$n - 1] = $this->H
[$j][$n - 1] / $t;
747 $this->H
[$j][$n] = $this->H
[$j][$n] / $t;
752 } // end else for complex case
755 // Vectors of isolated roots
756 for ($i = 0; $i < $nn; ++
$i) {
757 if ($i < $low |
$i > $high) {
758 for ($j = $i; $j < $nn; ++
$j) {
759 $this->V
[$i][$j] = $this->H
[$i][$j];
764 // Back transformation to get eigenvectors of original matrix
765 for ($j = $nn - 1; $j >= $low; --$j) {
766 for ($i = $low; $i <= $high; ++
$i) {
768 $kMax = min($j, $high);
769 for ($k = $low; $k <= $kMax; ++
$k) {
770 $z = $z +
$this->V
[$i][$k] * $this->H
[$k][$j];
772 $this->V
[$i][$j] = $z;
780 * Constructor: Check for symmetry, then construct the eigenvalue decomposition.
782 * @param mixed $Arg A Square matrix
784 public function __construct($Arg)
786 $this->A
= $Arg->getArray();
787 $this->n
= $Arg->getColumnDimension();
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]);
805 // Reduce to Hessenberg form.
807 // Reduce Hessenberg to real Schur form.
813 * Return the eigenvector matrix.
817 public function getV()
819 return new Matrix($this->V
, $this->n
, $this->n
);
823 * Return the real parts of the eigenvalues.
825 * @return array real(diag(D))
827 public function getRealEigenvalues()
833 * Return the imaginary parts of the eigenvalues.
835 * @return array imag(diag(D))
837 public function getImagEigenvalues()
843 * Return the block diagonal eigenvalue matrix.
847 public function getD()
849 for ($i = 0; $i < $this->n
; ++
$i) {
850 $D[$i] = array_fill(0, $this->n
, 0.0);
851 $D[$i][$i] = $this->d
[$i];
852 if ($this->e
[$i] == 0) {
855 $o = ($this->e
[$i] > 0) ?
$i +
1 : $i - 1;
856 $D[$i][$o] = $this->e
[$i];
859 return new Matrix($D);