1<?php
2
3declare(strict_types=1);
4
5/**
6 * Class to obtain eigenvalues and eigenvectors of a real matrix.
7 *
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).
12 *
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().
20 *
21 * @author Paul Meagher
22 * @license PHP v3.0
23 *
24 * @version 1.1
25 *
26 *  Slightly changed to adapt the original code to PHP-ML library
27 *  @date 2017/04/11
28 *
29 *  @author Mustafa Karabulut
30 */
31
32namespace Phpml\Math\LinearAlgebra;
33
34use Phpml\Math\Matrix;
35
36class EigenvalueDecomposition
37{
38    /**
39     * Row and column dimension (square matrix).
40     *
41     * @var int
42     */
43    private $n;
44
45    /**
46     * Arrays for internal storage of eigenvalues.
47     *
48     * @var array
49     */
50    private $d = [];
51
52    /**
53     * @var array
54     */
55    private $e = [];
56
57    /**
58     * Array for internal storage of eigenvectors.
59     *
60     * @var array
61     */
62    private $V = [];
63
64    /**
65     * Array for internal storage of nonsymmetric Hessenberg form.
66     *
67     * @var array
68     */
69    private $H = [];
70
71    /**
72     * Working storage for nonsymmetric algorithm.
73     *
74     * @var array
75     */
76    private $ort = [];
77
78    /**
79     * Used for complex scalar division.
80     *
81     * @var float
82     */
83    private $cdivr;
84
85    /**
86     * @var float
87     */
88    private $cdivi;
89
90    /**
91     * Constructor: Check for symmetry, then construct the eigenvalue decomposition
92     */
93    public function __construct(array $arg)
94    {
95        $this->n = count($arg[0]);
96        $symmetric = true;
97
98        for ($j = 0; ($j < $this->n) & $symmetric; ++$j) {
99            for ($i = 0; ($i < $this->n) & $symmetric; ++$i) {
100                $symmetric = $arg[$i][$j] == $arg[$j][$i];
101            }
102        }
103
104        if ($symmetric) {
105            $this->V = $arg;
106            // Tridiagonalize.
107            $this->tred2();
108            // Diagonalize.
109            $this->tql2();
110        } else {
111            $this->H = $arg;
112            $this->ort = [];
113            // Reduce to Hessenberg form.
114            $this->orthes();
115            // Reduce Hessenberg to real Schur form.
116            $this->hqr2();
117        }
118    }
119
120    /**
121     * Return the eigenvector matrix
122     */
123    public function getEigenvectors(): array
124    {
125        $vectors = $this->V;
126
127        // Always return the eigenvectors of length 1.0
128        $vectors = new Matrix($vectors);
129        $vectors = array_map(function ($vect) {
130            $sum = 0;
131            $count = count($vect);
132            for ($i = 0; $i < $count; ++$i) {
133                $sum += $vect[$i] ** 2;
134            }
135
136            $sum **= .5;
137            for ($i = 0; $i < $count; ++$i) {
138                $vect[$i] /= $sum;
139            }
140
141            return $vect;
142        }, $vectors->transpose()->toArray());
143
144        return $vectors;
145    }
146
147    /**
148     * Return the real parts of the eigenvalues<br>
149     *  d = real(diag(D));
150     */
151    public function getRealEigenvalues(): array
152    {
153        return $this->d;
154    }
155
156    /**
157     * Return the imaginary parts of the eigenvalues <br>
158     *  d = imag(diag(D))
159     */
160    public function getImagEigenvalues(): array
161    {
162        return $this->e;
163    }
164
165    /**
166     * Return the block diagonal eigenvalue matrix
167     */
168    public function getDiagonalEigenvalues(): array
169    {
170        $D = [];
171
172        for ($i = 0; $i < $this->n; ++$i) {
173            $D[$i] = array_fill(0, $this->n, 0.0);
174            $D[$i][$i] = $this->d[$i];
175            if ($this->e[$i] == 0) {
176                continue;
177            }
178
179            $o = $this->e[$i] > 0 ? $i + 1 : $i - 1;
180            $D[$i][$o] = $this->e[$i];
181        }
182
183        return $D;
184    }
185
186    /**
187     * Symmetric Householder reduction to tridiagonal form.
188     */
189    private function tred2(): void
190    {
191        //  This is derived from the Algol procedures tred2 by
192        //  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
193        //  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
194        //  Fortran subroutine in EISPACK.
195        $this->d = $this->V[$this->n - 1];
196        // Householder reduction to tridiagonal form.
197        for ($i = $this->n - 1; $i > 0; --$i) {
198            $i_ = $i - 1;
199            // Scale to avoid under/overflow.
200            $h = $scale = 0.0;
201            $scale += array_sum(array_map('abs', $this->d));
202            if ($scale == 0.0) {
203                $this->e[$i] = $this->d[$i_];
204                $this->d = array_slice($this->V[$i_], 0, $this->n - 1);
205                for ($j = 0; $j < $i; ++$j) {
206                    $this->V[$j][$i] = $this->V[$i][$j] = 0.0;
207                }
208            } else {
209                // Generate Householder vector.
210                for ($k = 0; $k < $i; ++$k) {
211                    $this->d[$k] /= $scale;
212                    $h += $this->d[$k] ** 2;
213                }
214
215                $f = $this->d[$i_];
216                $g = $h ** .5;
217                if ($f > 0) {
218                    $g = -$g;
219                }
220
221                $this->e[$i] = $scale * $g;
222                $h -= $f * $g;
223                $this->d[$i_] = $f - $g;
224
225                for ($j = 0; $j < $i; ++$j) {
226                    $this->e[$j] = 0.0;
227                }
228
229                // Apply similarity transformation to remaining columns.
230                for ($j = 0; $j < $i; ++$j) {
231                    $f = $this->d[$j];
232                    $this->V[$j][$i] = $f;
233                    $g = $this->e[$j] + $this->V[$j][$j] * $f;
234
235                    for ($k = $j + 1; $k <= $i_; ++$k) {
236                        $g += $this->V[$k][$j] * $this->d[$k];
237                        $this->e[$k] += $this->V[$k][$j] * $f;
238                    }
239
240                    $this->e[$j] = $g;
241                }
242
243                $f = 0.0;
244
245                if ($h == 0.0) {
246                    $h = 1e-32;
247                }
248
249                for ($j = 0; $j < $i; ++$j) {
250                    $this->e[$j] /= $h;
251                    $f += $this->e[$j] * $this->d[$j];
252                }
253
254                $hh = $f / (2 * $h);
255                for ($j = 0; $j < $i; ++$j) {
256                    $this->e[$j] -= $hh * $this->d[$j];
257                }
258
259                for ($j = 0; $j < $i; ++$j) {
260                    $f = $this->d[$j];
261                    $g = $this->e[$j];
262                    for ($k = $j; $k <= $i_; ++$k) {
263                        $this->V[$k][$j] -= ($f * $this->e[$k] + $g * $this->d[$k]);
264                    }
265
266                    $this->d[$j] = $this->V[$i - 1][$j];
267                    $this->V[$i][$j] = 0.0;
268                }
269            }
270
271            $this->d[$i] = $h;
272        }
273
274        // Accumulate transformations.
275        for ($i = 0; $i < $this->n - 1; ++$i) {
276            $this->V[$this->n - 1][$i] = $this->V[$i][$i];
277            $this->V[$i][$i] = 1.0;
278            $h = $this->d[$i + 1];
279            if ($h != 0.0) {
280                for ($k = 0; $k <= $i; ++$k) {
281                    $this->d[$k] = $this->V[$k][$i + 1] / $h;
282                }
283
284                for ($j = 0; $j <= $i; ++$j) {
285                    $g = 0.0;
286                    for ($k = 0; $k <= $i; ++$k) {
287                        $g += $this->V[$k][$i + 1] * $this->V[$k][$j];
288                    }
289
290                    for ($k = 0; $k <= $i; ++$k) {
291                        $this->V[$k][$j] -= $g * $this->d[$k];
292                    }
293                }
294            }
295
296            for ($k = 0; $k <= $i; ++$k) {
297                $this->V[$k][$i + 1] = 0.0;
298            }
299        }
300
301        $this->d = $this->V[$this->n - 1];
302        $this->V[$this->n - 1] = array_fill(0, $this->n, 0.0);
303        $this->V[$this->n - 1][$this->n - 1] = 1.0;
304        $this->e[0] = 0.0;
305    }
306
307    /**
308     * Symmetric tridiagonal QL algorithm.
309     *
310     * This is derived from the Algol procedures tql2, by
311     * Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
312     * Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
313     * Fortran subroutine in EISPACK.
314     */
315    private function tql2(): void
316    {
317        for ($i = 1; $i < $this->n; ++$i) {
318            $this->e[$i - 1] = $this->e[$i];
319        }
320
321        $this->e[$this->n - 1] = 0.0;
322        $f = 0.0;
323        $tst1 = 0.0;
324        $eps = 2.0 ** -52.0;
325
326        for ($l = 0; $l < $this->n; ++$l) {
327            // Find small subdiagonal element
328            $tst1 = max($tst1, abs($this->d[$l]) + abs($this->e[$l]));
329            $m = $l;
330            while ($m < $this->n) {
331                if (abs($this->e[$m]) <= $eps * $tst1) {
332                    break;
333                }
334
335                ++$m;
336            }
337
338            // If m == l, $this->d[l] is an eigenvalue,
339            // otherwise, iterate.
340            if ($m > $l) {
341                do {
342                    // Compute implicit shift
343                    $g = $this->d[$l];
344                    $p = ($this->d[$l + 1] - $g) / (2.0 * $this->e[$l]);
345                    $r = hypot($p, 1.0);
346                    if ($p < 0) {
347                        $r *= -1;
348                    }
349
350                    $this->d[$l] = $this->e[$l] / ($p + $r);
351                    $this->d[$l + 1] = $this->e[$l] * ($p + $r);
352                    $dl1 = $this->d[$l + 1];
353                    $h = $g - $this->d[$l];
354                    for ($i = $l + 2; $i < $this->n; ++$i) {
355                        $this->d[$i] -= $h;
356                    }
357
358                    $f += $h;
359                    // Implicit QL transformation.
360                    $p = $this->d[$m];
361                    $c = 1.0;
362                    $c2 = $c3 = $c;
363                    $el1 = $this->e[$l + 1];
364                    $s = $s2 = 0.0;
365                    for ($i = $m - 1; $i >= $l; --$i) {
366                        $c3 = $c2;
367                        $c2 = $c;
368                        $s2 = $s;
369                        $g = $c * $this->e[$i];
370                        $h = $c * $p;
371                        $r = hypot($p, $this->e[$i]);
372                        $this->e[$i + 1] = $s * $r;
373                        $s = $this->e[$i] / $r;
374                        $c = $p / $r;
375                        $p = $c * $this->d[$i] - $s * $g;
376                        $this->d[$i + 1] = $h + $s * ($c * $g + $s * $this->d[$i]);
377                        // Accumulate transformation.
378                        for ($k = 0; $k < $this->n; ++$k) {
379                            $h = $this->V[$k][$i + 1];
380                            $this->V[$k][$i + 1] = $s * $this->V[$k][$i] + $c * $h;
381                            $this->V[$k][$i] = $c * $this->V[$k][$i] - $s * $h;
382                        }
383                    }
384
385                    $p = -$s * $s2 * $c3 * $el1 * $this->e[$l] / $dl1;
386                    $this->e[$l] = $s * $p;
387                    $this->d[$l] = $c * $p;
388                    // Check for convergence.
389                } while (abs($this->e[$l]) > $eps * $tst1);
390            }
391
392            $this->d[$l] += $f;
393            $this->e[$l] = 0.0;
394        }
395
396        // Sort eigenvalues and corresponding vectors.
397        for ($i = 0; $i < $this->n - 1; ++$i) {
398            $k = $i;
399            $p = $this->d[$i];
400            for ($j = $i + 1; $j < $this->n; ++$j) {
401                if ($this->d[$j] < $p) {
402                    $k = $j;
403                    $p = $this->d[$j];
404                }
405            }
406
407            if ($k != $i) {
408                $this->d[$k] = $this->d[$i];
409                $this->d[$i] = $p;
410                for ($j = 0; $j < $this->n; ++$j) {
411                    $p = $this->V[$j][$i];
412                    $this->V[$j][$i] = $this->V[$j][$k];
413                    $this->V[$j][$k] = $p;
414                }
415            }
416        }
417    }
418
419    /**
420     * Nonsymmetric reduction to Hessenberg form.
421     *
422     * This is derived from the Algol procedures orthes and ortran,
423     * by Martin and Wilkinson, Handbook for Auto. Comp.,
424     * Vol.ii-Linear Algebra, and the corresponding
425     * Fortran subroutines in EISPACK.
426     */
427    private function orthes(): void
428    {
429        $low = 0;
430        $high = $this->n - 1;
431
432        for ($m = $low + 1; $m <= $high - 1; ++$m) {
433            // Scale column.
434            $scale = 0.0;
435            for ($i = $m; $i <= $high; ++$i) {
436                $scale += abs($this->H[$i][$m - 1]);
437            }
438
439            if ($scale != 0.0) {
440                // Compute Householder transformation.
441                $h = 0.0;
442                for ($i = $high; $i >= $m; --$i) {
443                    $this->ort[$i] = $this->H[$i][$m - 1] / $scale;
444                    $h += $this->ort[$i] * $this->ort[$i];
445                }
446
447                $g = $h ** .5;
448                if ($this->ort[$m] > 0) {
449                    $g *= -1;
450                }
451
452                $h -= $this->ort[$m] * $g;
453                $this->ort[$m] -= $g;
454                // Apply Householder similarity transformation
455                // H = (I -u * u' / h) * H * (I -u * u') / h)
456                for ($j = $m; $j < $this->n; ++$j) {
457                    $f = 0.0;
458                    for ($i = $high; $i >= $m; --$i) {
459                        $f += $this->ort[$i] * $this->H[$i][$j];
460                    }
461
462                    $f /= $h;
463                    for ($i = $m; $i <= $high; ++$i) {
464                        $this->H[$i][$j] -= $f * $this->ort[$i];
465                    }
466                }
467
468                for ($i = 0; $i <= $high; ++$i) {
469                    $f = 0.0;
470                    for ($j = $high; $j >= $m; --$j) {
471                        $f += $this->ort[$j] * $this->H[$i][$j];
472                    }
473
474                    $f /= $h;
475                    for ($j = $m; $j <= $high; ++$j) {
476                        $this->H[$i][$j] -= $f * $this->ort[$j];
477                    }
478                }
479
480                $this->ort[$m] = $scale * $this->ort[$m];
481                $this->H[$m][$m - 1] = $scale * $g;
482            }
483        }
484
485        // Accumulate transformations (Algol's ortran).
486        for ($i = 0; $i < $this->n; ++$i) {
487            for ($j = 0; $j < $this->n; ++$j) {
488                $this->V[$i][$j] = ($i == $j ? 1.0 : 0.0);
489            }
490        }
491
492        for ($m = $high - 1; $m >= $low + 1; --$m) {
493            if ($this->H[$m][$m - 1] != 0.0) {
494                for ($i = $m + 1; $i <= $high; ++$i) {
495                    $this->ort[$i] = $this->H[$i][$m - 1];
496                }
497
498                for ($j = $m; $j <= $high; ++$j) {
499                    $g = 0.0;
500                    for ($i = $m; $i <= $high; ++$i) {
501                        $g += $this->ort[$i] * $this->V[$i][$j];
502                    }
503
504                    // Double division avoids possible underflow
505                    $g = ($g / $this->ort[$m]) / $this->H[$m][$m - 1];
506                    for ($i = $m; $i <= $high; ++$i) {
507                        $this->V[$i][$j] += $g * $this->ort[$i];
508                    }
509                }
510            }
511        }
512    }
513
514    /**
515     * Performs complex division.
516     *
517     * @param int|float $xr
518     * @param int|float $xi
519     * @param int|float $yr
520     * @param int|float $yi
521     */
522    private function cdiv($xr, $xi, $yr, $yi): void
523    {
524        if (abs($yr) > abs($yi)) {
525            $r = $yi / $yr;
526            $d = $yr + $r * $yi;
527            $this->cdivr = ($xr + $r * $xi) / $d;
528            $this->cdivi = ($xi - $r * $xr) / $d;
529        } else {
530            $r = $yr / $yi;
531            $d = $yi + $r * $yr;
532            $this->cdivr = ($r * $xr + $xi) / $d;
533            $this->cdivi = ($r * $xi - $xr) / $d;
534        }
535    }
536
537    /**
538     * Nonsymmetric reduction from Hessenberg to real Schur form.
539     *
540     * Code is derived from the Algol procedure hqr2,
541     * by Martin and Wilkinson, Handbook for Auto. Comp.,
542     * Vol.ii-Linear Algebra, and the corresponding
543     * Fortran subroutine in EISPACK.
544     */
545    private function hqr2(): void
546    {
547        //  Initialize
548        $nn = $this->n;
549        $n = $nn - 1;
550        $low = 0;
551        $high = $nn - 1;
552        $eps = 2.0 ** -52.0;
553        $exshift = 0.0;
554        $p = $q = $r = $s = $z = 0;
555        // Store roots isolated by balanc and compute matrix norm
556        $norm = 0.0;
557
558        for ($i = 0; $i < $nn; ++$i) {
559            if (($i < $low) or ($i > $high)) {
560                $this->d[$i] = $this->H[$i][$i];
561                $this->e[$i] = 0.0;
562            }
563
564            for ($j = max($i - 1, 0); $j < $nn; ++$j) {
565                $norm += abs($this->H[$i][$j]);
566            }
567        }
568
569        // Outer loop over eigenvalue index
570        $iter = 0;
571        while ($n >= $low) {
572            // Look for single small sub-diagonal element
573            $l = $n;
574            while ($l > $low) {
575                $s = abs($this->H[$l - 1][$l - 1]) + abs($this->H[$l][$l]);
576                if ($s == 0.0) {
577                    $s = $norm;
578                }
579
580                if (abs($this->H[$l][$l - 1]) < $eps * $s) {
581                    break;
582                }
583
584                --$l;
585            }
586
587            // Check for convergence
588            // One root found
589            if ($l == $n) {
590                $this->H[$n][$n] += $exshift;
591                $this->d[$n] = $this->H[$n][$n];
592                $this->e[$n] = 0.0;
593                --$n;
594                $iter = 0;
595            // Two roots found
596            } elseif ($l == $n - 1) {
597                $w = $this->H[$n][$n - 1] * $this->H[$n - 1][$n];
598                $p = ($this->H[$n - 1][$n - 1] - $this->H[$n][$n]) / 2.0;
599                $q = $p * $p + $w;
600                $z = abs($q) ** .5;
601                $this->H[$n][$n] += $exshift;
602                $this->H[$n - 1][$n - 1] += $exshift;
603                $x = $this->H[$n][$n];
604                // Real pair
605                if ($q >= 0) {
606                    if ($p >= 0) {
607                        $z = $p + $z;
608                    } else {
609                        $z = $p - $z;
610                    }
611
612                    $this->d[$n - 1] = $x + $z;
613                    $this->d[$n] = $this->d[$n - 1];
614                    if ($z != 0.0) {
615                        $this->d[$n] = $x - $w / $z;
616                    }
617
618                    $this->e[$n - 1] = 0.0;
619                    $this->e[$n] = 0.0;
620                    $x = $this->H[$n][$n - 1];
621                    $s = abs($x) + abs($z);
622                    $p = $x / $s;
623                    $q = $z / $s;
624                    $r = ($p * $p + $q * $q) ** .5;
625                    $p /= $r;
626                    $q /= $r;
627                    // Row modification
628                    for ($j = $n - 1; $j < $nn; ++$j) {
629                        $z = $this->H[$n - 1][$j];
630                        $this->H[$n - 1][$j] = $q * $z + $p * $this->H[$n][$j];
631                        $this->H[$n][$j] = $q * $this->H[$n][$j] - $p * $z;
632                    }
633
634                    // Column modification
635                    for ($i = 0; $i <= $n; ++$i) {
636                        $z = $this->H[$i][$n - 1];
637                        $this->H[$i][$n - 1] = $q * $z + $p * $this->H[$i][$n];
638                        $this->H[$i][$n] = $q * $this->H[$i][$n] - $p * $z;
639                    }
640
641                    // Accumulate transformations
642                    for ($i = $low; $i <= $high; ++$i) {
643                        $z = $this->V[$i][$n - 1];
644                        $this->V[$i][$n - 1] = $q * $z + $p * $this->V[$i][$n];
645                        $this->V[$i][$n] = $q * $this->V[$i][$n] - $p * $z;
646                    }
647
648                    // Complex pair
649                } else {
650                    $this->d[$n - 1] = $x + $p;
651                    $this->d[$n] = $x + $p;
652                    $this->e[$n - 1] = $z;
653                    $this->e[$n] = -$z;
654                }
655
656                $n -= 2;
657                $iter = 0;
658            // No convergence yet
659            } else {
660                // Form shift
661                $x = $this->H[$n][$n];
662                $y = 0.0;
663                $w = 0.0;
664                if ($l < $n) {
665                    $y = $this->H[$n - 1][$n - 1];
666                    $w = $this->H[$n][$n - 1] * $this->H[$n - 1][$n];
667                }
668
669                // Wilkinson's original ad hoc shift
670                if ($iter == 10) {
671                    $exshift += $x;
672                    for ($i = $low; $i <= $n; ++$i) {
673                        $this->H[$i][$i] -= $x;
674                    }
675
676                    $s = abs($this->H[$n][$n - 1]) + abs($this->H[$n - 1][$n - 2]);
677                    $x = $y = 0.75 * $s;
678                    $w = -0.4375 * $s * $s;
679                }
680
681                // MATLAB's new ad hoc shift
682                if ($iter == 30) {
683                    $s = ($y - $x) / 2.0;
684                    $s *= $s + $w;
685                    if ($s > 0) {
686                        $s **= .5;
687                        if ($y < $x) {
688                            $s = -$s;
689                        }
690
691                        $s = $x - $w / (($y - $x) / 2.0 + $s);
692                        for ($i = $low; $i <= $n; ++$i) {
693                            $this->H[$i][$i] -= $s;
694                        }
695
696                        $exshift += $s;
697                        $x = $y = $w = 0.964;
698                    }
699                }
700
701                // Could check iteration count here.
702                ++$iter;
703                // Look for two consecutive small sub-diagonal elements
704                $m = $n - 2;
705                while ($m >= $l) {
706                    $z = $this->H[$m][$m];
707                    $r = $x - $z;
708                    $s = $y - $z;
709                    $p = ($r * $s - $w) / $this->H[$m + 1][$m] + $this->H[$m][$m + 1];
710                    $q = $this->H[$m + 1][$m + 1] - $z - $r - $s;
711                    $r = $this->H[$m + 2][$m + 1];
712                    $s = abs($p) + abs($q) + abs($r);
713                    $p /= $s;
714                    $q /= $s;
715                    $r /= $s;
716                    if ($m == $l) {
717                        break;
718                    }
719
720                    if (abs($this->H[$m][$m - 1]) * (abs($q) + abs($r)) <
721                        $eps * (abs($p) * (abs($this->H[$m - 1][$m - 1]) + abs($z) + abs($this->H[$m + 1][$m + 1])))) {
722                        break;
723                    }
724
725                    --$m;
726                }
727
728                for ($i = $m + 2; $i <= $n; ++$i) {
729                    $this->H[$i][$i - 2] = 0.0;
730                    if ($i > $m + 2) {
731                        $this->H[$i][$i - 3] = 0.0;
732                    }
733                }
734
735                // Double QR step involving rows l:n and columns m:n
736                for ($k = $m; $k <= $n - 1; ++$k) {
737                    $notlast = ($k != $n - 1);
738                    if ($k != $m) {
739                        $p = $this->H[$k][$k - 1];
740                        $q = $this->H[$k + 1][$k - 1];
741                        $r = ($notlast ? $this->H[$k + 2][$k - 1] : 0.0);
742                        $x = abs($p) + abs($q) + abs($r);
743                        if ($x != 0.0) {
744                            $p /= $x;
745                            $q /= $x;
746                            $r /= $x;
747                        }
748                    }
749
750                    if ($x == 0.0) {
751                        break;
752                    }
753
754                    $s = ($p * $p + $q * $q + $r * $r) ** .5;
755                    if ($p < 0) {
756                        $s = -$s;
757                    }
758
759                    if ($s != 0) {
760                        if ($k != $m) {
761                            $this->H[$k][$k - 1] = -$s * $x;
762                        } elseif ($l != $m) {
763                            $this->H[$k][$k - 1] = -$this->H[$k][$k - 1];
764                        }
765
766                        $p += $s;
767                        $x = $p / $s;
768                        $y = $q / $s;
769                        $z = $r / $s;
770                        $q /= $p;
771                        $r /= $p;
772                        // Row modification
773                        for ($j = $k; $j < $nn; ++$j) {
774                            $p = $this->H[$k][$j] + $q * $this->H[$k + 1][$j];
775                            if ($notlast) {
776                                $p += $r * $this->H[$k + 2][$j];
777                                $this->H[$k + 2][$j] -= $p * $z;
778                            }
779
780                            $this->H[$k][$j] -= $p * $x;
781                            $this->H[$k + 1][$j] -= $p * $y;
782                        }
783
784                        // Column modification
785                        for ($i = 0; $i <= min($n, $k + 3); ++$i) {
786                            $p = $x * $this->H[$i][$k] + $y * $this->H[$i][$k + 1];
787                            if ($notlast) {
788                                $p += $z * $this->H[$i][$k + 2];
789                                $this->H[$i][$k + 2] -= $p * $r;
790                            }
791
792                            $this->H[$i][$k] -= $p;
793                            $this->H[$i][$k + 1] -= $p * $q;
794                        }
795
796                        // Accumulate transformations
797                        for ($i = $low; $i <= $high; ++$i) {
798                            $p = $x * $this->V[$i][$k] + $y * $this->V[$i][$k + 1];
799                            if ($notlast) {
800                                $p += $z * $this->V[$i][$k + 2];
801                                $this->V[$i][$k + 2] -= $p * $r;
802                            }
803
804                            $this->V[$i][$k] -= $p;
805                            $this->V[$i][$k + 1] -= $p * $q;
806                        }
807                    }  // ($s != 0)
808                }  // k loop
809            }  // check convergence
810        }  // while ($n >= $low)
811
812        // Backsubstitute to find vectors of upper triangular form
813        if ($norm == 0.0) {
814            return;
815        }
816
817        for ($n = $nn - 1; $n >= 0; --$n) {
818            $p = $this->d[$n];
819            $q = $this->e[$n];
820            // Real vector
821            if ($q == 0) {
822                $l = $n;
823                $this->H[$n][$n] = 1.0;
824                for ($i = $n - 1; $i >= 0; --$i) {
825                    $w = $this->H[$i][$i] - $p;
826                    $r = 0.0;
827                    for ($j = $l; $j <= $n; ++$j) {
828                        $r += $this->H[$i][$j] * $this->H[$j][$n];
829                    }
830
831                    if ($this->e[$i] < 0.0) {
832                        $z = $w;
833                        $s = $r;
834                    } else {
835                        $l = $i;
836                        if ($this->e[$i] == 0.0) {
837                            if ($w != 0.0) {
838                                $this->H[$i][$n] = -$r / $w;
839                            } else {
840                                $this->H[$i][$n] = -$r / ($eps * $norm);
841                            }
842
843                            // Solve real equations
844                        } else {
845                            $x = $this->H[$i][$i + 1];
846                            $y = $this->H[$i + 1][$i];
847                            $q = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i];
848                            $t = ($x * $s - $z * $r) / $q;
849                            $this->H[$i][$n] = $t;
850                            if (abs($x) > abs($z)) {
851                                $this->H[$i + 1][$n] = (-$r - $w * $t) / $x;
852                            } else {
853                                $this->H[$i + 1][$n] = (-$s - $y * $t) / $z;
854                            }
855                        }
856
857                        // Overflow control
858                        $t = abs($this->H[$i][$n]);
859                        if (($eps * $t) * $t > 1) {
860                            for ($j = $i; $j <= $n; ++$j) {
861                                $this->H[$j][$n] /= $t;
862                            }
863                        }
864                    }
865                }
866
867                // Complex vector
868            } elseif ($q < 0) {
869                $l = $n - 1;
870                // Last vector component imaginary so matrix is triangular
871                if (abs($this->H[$n][$n - 1]) > abs($this->H[$n - 1][$n])) {
872                    $this->H[$n - 1][$n - 1] = $q / $this->H[$n][$n - 1];
873                    $this->H[$n - 1][$n] = -($this->H[$n][$n] - $p) / $this->H[$n][$n - 1];
874                } else {
875                    $this->cdiv(0.0, -$this->H[$n - 1][$n], $this->H[$n - 1][$n - 1] - $p, $q);
876                    $this->H[$n - 1][$n - 1] = $this->cdivr;
877                    $this->H[$n - 1][$n] = $this->cdivi;
878                }
879
880                $this->H[$n][$n - 1] = 0.0;
881                $this->H[$n][$n] = 1.0;
882                for ($i = $n - 2; $i >= 0; --$i) {
883                    // double ra,sa,vr,vi;
884                    $ra = 0.0;
885                    $sa = 0.0;
886                    for ($j = $l; $j <= $n; ++$j) {
887                        $ra += $this->H[$i][$j] * $this->H[$j][$n - 1];
888                        $sa += $this->H[$i][$j] * $this->H[$j][$n];
889                    }
890
891                    $w = $this->H[$i][$i] - $p;
892                    if ($this->e[$i] < 0.0) {
893                        $z = $w;
894                        $r = $ra;
895                        $s = $sa;
896                    } else {
897                        $l = $i;
898                        if ($this->e[$i] == 0) {
899                            $this->cdiv(-$ra, -$sa, $w, $q);
900                            $this->H[$i][$n - 1] = $this->cdivr;
901                            $this->H[$i][$n] = $this->cdivi;
902                        } else {
903                            // Solve complex equations
904                            $x = $this->H[$i][$i + 1];
905                            $y = $this->H[$i + 1][$i];
906                            $vr = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i] - $q * $q;
907                            $vi = ($this->d[$i] - $p) * 2.0 * $q;
908                            if ($vr == 0.0 && $vi == 0.0) {
909                                $vr = $eps * $norm * (abs($w) + abs($q) + abs($x) + abs($y) + abs($z));
910                            }
911
912                            $this->cdiv($x * $r - $z * $ra + $q * $sa, $x * $s - $z * $sa - $q * $ra, $vr, $vi);
913                            $this->H[$i][$n - 1] = $this->cdivr;
914                            $this->H[$i][$n] = $this->cdivi;
915                            if (abs($x) > (abs($z) + abs($q))) {
916                                $this->H[$i + 1][$n - 1] = (-$ra - $w * $this->H[$i][$n - 1] + $q * $this->H[$i][$n]) / $x;
917                                $this->H[$i + 1][$n] = (-$sa - $w * $this->H[$i][$n] - $q * $this->H[$i][$n - 1]) / $x;
918                            } else {
919                                $this->cdiv(-$r - $y * $this->H[$i][$n - 1], -$s - $y * $this->H[$i][$n], $z, $q);
920                                $this->H[$i + 1][$n - 1] = $this->cdivr;
921                                $this->H[$i + 1][$n] = $this->cdivi;
922                            }
923                        }
924
925                        // Overflow control
926                        $t = max(abs($this->H[$i][$n - 1]), abs($this->H[$i][$n]));
927                        if (($eps * $t) * $t > 1) {
928                            for ($j = $i; $j <= $n; ++$j) {
929                                $this->H[$j][$n - 1] /= $t;
930                                $this->H[$j][$n] /= $t;
931                            }
932                        }
933                    } // end else
934                } // end for
935            } // end else for complex case
936        } // end for
937
938        // Vectors of isolated roots
939        for ($i = 0; $i < $nn; ++$i) {
940            if ($i < $low || $i > $high) {
941                for ($j = $i; $j < $nn; ++$j) {
942                    $this->V[$i][$j] = $this->H[$i][$j];
943                }
944            }
945        }
946
947        // Back transformation to get eigenvectors of original matrix
948        for ($j = $nn - 1; $j >= $low; --$j) {
949            for ($i = $low; $i <= $high; ++$i) {
950                $z = 0.0;
951                for ($k = $low; $k <= min($j, $high); ++$k) {
952                    $z += $this->V[$i][$k] * $this->H[$k][$j];
953                }
954
955                $this->V[$i][$j] = $z;
956            }
957        }
958    }
959}
960