1<?php
2
3namespace PhpOffice\PhpSpreadsheet\Shared\JAMA;
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 *
23 *    @version 1.1
24 */
25class EigenvalueDecomposition
26{
27    /**
28     * Row and column dimension (square matrix).
29     *
30     * @var int
31     */
32    private $n;
33
34    /**
35     * Arrays for internal storage of eigenvalues.
36     *
37     * @var array
38     */
39    private $d = [];
40
41    private $e = [];
42
43    /**
44     * Array for internal storage of eigenvectors.
45     *
46     * @var array
47     */
48    private $V = [];
49
50    /**
51     * Array for internal storage of nonsymmetric Hessenberg form.
52     *
53     * @var array
54     */
55    private $H = [];
56
57    /**
58     * Working storage for nonsymmetric algorithm.
59     *
60     * @var array
61     */
62    private $ort;
63
64    /**
65     * Used for complex scalar division.
66     *
67     * @var float
68     */
69    private $cdivr;
70
71    private $cdivi;
72
73    /**
74     * Symmetric Householder reduction to tridiagonal form.
75     */
76    private function tred2()
77    {
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) {
85            $i_ = $i - 1;
86            // Scale to avoid under/overflow.
87            $h = $scale = 0.0;
88            $scale += array_sum(array_map('abs', $this->d));
89            if ($scale == 0.0) {
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;
94                }
95            } else {
96                // Generate Householder vector.
97                for ($k = 0; $k < $i; ++$k) {
98                    $this->d[$k] /= $scale;
99                    $h += pow($this->d[$k], 2);
100                }
101                $f = $this->d[$i_];
102                $g = sqrt($h);
103                if ($f > 0) {
104                    $g = -$g;
105                }
106                $this->e[$i] = $scale * $g;
107                $h = $h - $f * $g;
108                $this->d[$i_] = $f - $g;
109                for ($j = 0; $j < $i; ++$j) {
110                    $this->e[$j] = 0.0;
111                }
112                // Apply similarity transformation to remaining columns.
113                for ($j = 0; $j < $i; ++$j) {
114                    $f = $this->d[$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;
120                    }
121                    $this->e[$j] = $g;
122                }
123                $f = 0.0;
124                for ($j = 0; $j < $i; ++$j) {
125                    $this->e[$j] /= $h;
126                    $f += $this->e[$j] * $this->d[$j];
127                }
128                $hh = $f / (2 * $h);
129                for ($j = 0; $j < $i; ++$j) {
130                    $this->e[$j] -= $hh * $this->d[$j];
131                }
132                for ($j = 0; $j < $i; ++$j) {
133                    $f = $this->d[$j];
134                    $g = $this->e[$j];
135                    for ($k = $j; $k <= $i_; ++$k) {
136                        $this->V[$k][$j] -= ($f * $this->e[$k] + $g * $this->d[$k]);
137                    }
138                    $this->d[$j] = $this->V[$i - 1][$j];
139                    $this->V[$i][$j] = 0.0;
140                }
141            }
142            $this->d[$i] = $h;
143        }
144
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];
150            if ($h != 0.0) {
151                for ($k = 0; $k <= $i; ++$k) {
152                    $this->d[$k] = $this->V[$k][$i + 1] / $h;
153                }
154                for ($j = 0; $j <= $i; ++$j) {
155                    $g = 0.0;
156                    for ($k = 0; $k <= $i; ++$k) {
157                        $g += $this->V[$k][$i + 1] * $this->V[$k][$j];
158                    }
159                    for ($k = 0; $k <= $i; ++$k) {
160                        $this->V[$k][$j] -= $g * $this->d[$k];
161                    }
162                }
163            }
164            for ($k = 0; $k <= $i; ++$k) {
165                $this->V[$k][$i + 1] = 0.0;
166            }
167        }
168
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;
172        $this->e[0] = 0.0;
173    }
174
175    /**
176     * Symmetric tridiagonal QL algorithm.
177     *
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.
182     */
183    private function tql2()
184    {
185        for ($i = 1; $i < $this->n; ++$i) {
186            $this->e[$i - 1] = $this->e[$i];
187        }
188        $this->e[$this->n - 1] = 0.0;
189        $f = 0.0;
190        $tst1 = 0.0;
191        $eps = pow(2.0, -52.0);
192
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]));
196            $m = $l;
197            while ($m < $this->n) {
198                if (abs($this->e[$m]) <= $eps * $tst1) {
199                    break;
200                }
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                    }
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) {
222                        $this->d[$i] -= $h;
223                    }
224                    $f += $h;
225                    // Implicit QL transformation.
226                    $p = $this->d[$m];
227                    $c = 1.0;
228                    $c2 = $c3 = $c;
229                    $el1 = $this->e[$l + 1];
230                    $s = $s2 = 0.0;
231                    for ($i = $m - 1; $i >= $l; --$i) {
232                        $c3 = $c2;
233                        $c2 = $c;
234                        $s2 = $s;
235                        $g = $c * $this->e[$i];
236                        $h = $c * $p;
237                        $r = hypo($p, $this->e[$i]);
238                        $this->e[$i + 1] = $s * $r;
239                        $s = $this->e[$i] / $r;
240                        $c = $p / $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;
248                        }
249                    }
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);
255            }
256            $this->d[$l] = $this->d[$l] + $f;
257            $this->e[$l] = 0.0;
258        }
259
260        // Sort eigenvalues and corresponding vectors.
261        for ($i = 0; $i < $this->n - 1; ++$i) {
262            $k = $i;
263            $p = $this->d[$i];
264            for ($j = $i + 1; $j < $this->n; ++$j) {
265                if ($this->d[$j] < $p) {
266                    $k = $j;
267                    $p = $this->d[$j];
268                }
269            }
270            if ($k != $i) {
271                $this->d[$k] = $this->d[$i];
272                $this->d[$i] = $p;
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;
277                }
278            }
279        }
280    }
281
282    /**
283     * Nonsymmetric reduction to Hessenberg form.
284     *
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.
289     */
290    private function orthes()
291    {
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    /**
368     * Performs complex division.
369     *
370     * @param mixed $xr
371     * @param mixed $xi
372     * @param mixed $yr
373     * @param mixed $yi
374     */
375    private function cdiv($xr, $xi, $yr, $yi)
376    {
377        if (abs($yr) > abs($yi)) {
378            $r = $yi / $yr;
379            $d = $yr + $r * $yi;
380            $this->cdivr = ($xr + $r * $xi) / $d;
381            $this->cdivi = ($xi - $r * $xr) / $d;
382        } else {
383            $r = $yr / $yi;
384            $d = $yi + $r * $yr;
385            $this->cdivr = ($r * $xr + $xi) / $d;
386            $this->cdivi = ($r * $xi - $xr) / $d;
387        }
388    }
389
390    /**
391     * Nonsymmetric reduction from Hessenberg to real Schur form.
392     *
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.
397     */
398    private function hqr2()
399    {
400        //  Initialize
401        $nn = $this->n;
402        $n = $nn - 1;
403        $low = 0;
404        $high = $nn - 1;
405        $eps = pow(2.0, -52.0);
406        $exshift = 0.0;
407        $p = $q = $r = $s = $z = 0;
408        // Store roots isolated by balanc and compute matrix norm
409        $norm = 0.0;
410
411        for ($i = 0; $i < $nn; ++$i) {
412            if (($i < $low) or ($i > $high)) {
413                $this->d[$i] = $this->H[$i][$i];
414                $this->e[$i] = 0.0;
415            }
416            for ($j = max($i - 1, 0); $j < $nn; ++$j) {
417                $norm = $norm + abs($this->H[$i][$j]);
418            }
419        }
420
421        // Outer loop over eigenvalue index
422        $iter = 0;
423        while ($n >= $low) {
424            // Look for single small sub-diagonal element
425            $l = $n;
426            while ($l > $low) {
427                $s = abs($this->H[$l - 1][$l - 1]) + abs($this->H[$l][$l]);
428                if ($s == 0.0) {
429                    $s = $norm;
430                }
431                if (abs($this->H[$l][$l - 1]) < $eps * $s) {
432                    break;
433                }
434                --$l;
435            }
436            // Check for convergence
437            // One root found
438            if ($l == $n) {
439                $this->H[$n][$n] = $this->H[$n][$n] + $exshift;
440                $this->d[$n] = $this->H[$n][$n];
441                $this->e[$n] = 0.0;
442                --$n;
443                $iter = 0;
444            // Two roots found
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;
448                $q = $p * $p + $w;
449                $z = sqrt(abs($q));
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];
453                // Real pair
454                if ($q >= 0) {
455                    if ($p >= 0) {
456                        $z = $p + $z;
457                    } else {
458                        $z = $p - $z;
459                    }
460                    $this->d[$n - 1] = $x + $z;
461                    $this->d[$n] = $this->d[$n - 1];
462                    if ($z != 0.0) {
463                        $this->d[$n] = $x - $w / $z;
464                    }
465                    $this->e[$n - 1] = 0.0;
466                    $this->e[$n] = 0.0;
467                    $x = $this->H[$n][$n - 1];
468                    $s = abs($x) + abs($z);
469                    $p = $x / $s;
470                    $q = $z / $s;
471                    $r = sqrt($p * $p + $q * $q);
472                    $p = $p / $r;
473                    $q = $q / $r;
474                    // Row modification
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;
479                    }
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;
485                    }
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;
491                    }
492                    // Complex pair
493                } else {
494                    $this->d[$n - 1] = $x + $p;
495                    $this->d[$n] = $x + $p;
496                    $this->e[$n - 1] = $z;
497                    $this->e[$n] = -$z;
498                }
499                $n = $n - 2;
500                $iter = 0;
501            // No convergence yet
502            } else {
503                // Form shift
504                $x = $this->H[$n][$n];
505                $y = 0.0;
506                $w = 0.0;
507                if ($l < $n) {
508                    $y = $this->H[$n - 1][$n - 1];
509                    $w = $this->H[$n][$n - 1] * $this->H[$n - 1][$n];
510                }
511                // Wilkinson's original ad hoc shift
512                if ($iter == 10) {
513                    $exshift += $x;
514                    for ($i = $low; $i <= $n; ++$i) {
515                        $this->H[$i][$i] -= $x;
516                    }
517                    $s = abs($this->H[$n][$n - 1]) + abs($this->H[$n - 1][$n - 2]);
518                    $x = $y = 0.75 * $s;
519                    $w = -0.4375 * $s * $s;
520                }
521                // MATLAB's new ad hoc shift
522                if ($iter == 30) {
523                    $s = ($y - $x) / 2.0;
524                    $s = $s * $s + $w;
525                    if ($s > 0) {
526                        $s = sqrt($s);
527                        if ($y < $x) {
528                            $s = -$s;
529                        }
530                        $s = $x - $w / (($y - $x) / 2.0 + $s);
531                        for ($i = $low; $i <= $n; ++$i) {
532                            $this->H[$i][$i] -= $s;
533                        }
534                        $exshift += $s;
535                        $x = $y = $w = 0.964;
536                    }
537                }
538                // Could check iteration count here.
539                $iter = $iter + 1;
540                // Look for two consecutive small sub-diagonal elements
541                $m = $n - 2;
542                while ($m >= $l) {
543                    $z = $this->H[$m][$m];
544                    $r = $x - $z;
545                    $s = $y - $z;
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);
550                    $p = $p / $s;
551                    $q = $q / $s;
552                    $r = $r / $s;
553                    if ($m == $l) {
554                        break;
555                    }
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])))) {
558                        break;
559                    }
560                    --$m;
561                }
562                for ($i = $m + 2; $i <= $n; ++$i) {
563                    $this->H[$i][$i - 2] = 0.0;
564                    if ($i > $m + 2) {
565                        $this->H[$i][$i - 3] = 0.0;
566                    }
567                }
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);
571                    if ($k != $m) {
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);
576                        if ($x != 0.0) {
577                            $p = $p / $x;
578                            $q = $q / $x;
579                            $r = $r / $x;
580                        }
581                    }
582                    if ($x == 0.0) {
583                        break;
584                    }
585                    $s = sqrt($p * $p + $q * $q + $r * $r);
586                    if ($p < 0) {
587                        $s = -$s;
588                    }
589                    if ($s != 0) {
590                        if ($k != $m) {
591                            $this->H[$k][$k - 1] = -$s * $x;
592                        } elseif ($l != $m) {
593                            $this->H[$k][$k - 1] = -$this->H[$k][$k - 1];
594                        }
595                        $p = $p + $s;
596                        $x = $p / $s;
597                        $y = $q / $s;
598                        $z = $r / $s;
599                        $q = $q / $p;
600                        $r = $r / $p;
601                        // Row modification
602                        for ($j = $k; $j < $nn; ++$j) {
603                            $p = $this->H[$k][$j] + $q * $this->H[$k + 1][$j];
604                            if ($notlast) {
605                                $p = $p + $r * $this->H[$k + 2][$j];
606                                $this->H[$k + 2][$j] = $this->H[$k + 2][$j] - $p * $z;
607                            }
608                            $this->H[$k][$j] = $this->H[$k][$j] - $p * $x;
609                            $this->H[$k + 1][$j] = $this->H[$k + 1][$j] - $p * $y;
610                        }
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];
615                            if ($notlast) {
616                                $p = $p + $z * $this->H[$i][$k + 2];
617                                $this->H[$i][$k + 2] = $this->H[$i][$k + 2] - $p * $r;
618                            }
619                            $this->H[$i][$k] = $this->H[$i][$k] - $p;
620                            $this->H[$i][$k + 1] = $this->H[$i][$k + 1] - $p * $q;
621                        }
622                        // Accumulate transformations
623                        for ($i = $low; $i <= $high; ++$i) {
624                            $p = $x * $this->V[$i][$k] + $y * $this->V[$i][$k + 1];
625                            if ($notlast) {
626                                $p = $p + $z * $this->V[$i][$k + 2];
627                                $this->V[$i][$k + 2] = $this->V[$i][$k + 2] - $p * $r;
628                            }
629                            $this->V[$i][$k] = $this->V[$i][$k] - $p;
630                            $this->V[$i][$k + 1] = $this->V[$i][$k + 1] - $p * $q;
631                        }
632                    }  // ($s != 0)
633                }  // k loop
634            }  // check convergence
635        }  // while ($n >= $low)
636
637        // Backsubstitute to find vectors of upper triangular form
638        if ($norm == 0.0) {
639            return;
640        }
641
642        for ($n = $nn - 1; $n >= 0; --$n) {
643            $p = $this->d[$n];
644            $q = $this->e[$n];
645            // Real vector
646            if ($q == 0) {
647                $l = $n;
648                $this->H[$n][$n] = 1.0;
649                for ($i = $n - 1; $i >= 0; --$i) {
650                    $w = $this->H[$i][$i] - $p;
651                    $r = 0.0;
652                    for ($j = $l; $j <= $n; ++$j) {
653                        $r = $r + $this->H[$i][$j] * $this->H[$j][$n];
654                    }
655                    if ($this->e[$i] < 0.0) {
656                        $z = $w;
657                        $s = $r;
658                    } else {
659                        $l = $i;
660                        if ($this->e[$i] == 0.0) {
661                            if ($w != 0.0) {
662                                $this->H[$i][$n] = -$r / $w;
663                            } else {
664                                $this->H[$i][$n] = -$r / ($eps * $norm);
665                            }
666                            // Solve real equations
667                        } else {
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;
675                            } else {
676                                $this->H[$i + 1][$n] = (-$s - $y * $t) / $z;
677                            }
678                        }
679                        // Overflow control
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;
684                            }
685                        }
686                    }
687                }
688                // Complex vector
689            } elseif ($q < 0) {
690                $l = $n - 1;
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];
695                } else {
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;
699                }
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;
704                    $ra = 0.0;
705                    $sa = 0.0;
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];
709                    }
710                    $w = $this->H[$i][$i] - $p;
711                    if ($this->e[$i] < 0.0) {
712                        $z = $w;
713                        $r = $ra;
714                        $s = $sa;
715                    } else {
716                        $l = $i;
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;
721                        } else {
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));
729                            }
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;
736                            } else {
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;
740                            }
741                        }
742                        // Overflow control
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;
748                            }
749                        }
750                    } // end else
751                } // end for
752            } // end else for complex case
753        } // end for
754
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];
760                }
761            }
762        }
763
764        // Back transformation to get eigenvectors of original matrix
765        for ($j = $nn - 1; $j >= $low; --$j) {
766            for ($i = $low; $i <= $high; ++$i) {
767                $z = 0.0;
768                $kMax = min($j, $high);
769                for ($k = $low; $k <= $kMax; ++$k) {
770                    $z = $z + $this->V[$i][$k] * $this->H[$k][$j];
771                }
772                $this->V[$i][$j] = $z;
773            }
774        }
775    }
776
777    // end hqr2
778
779    /**
780     * Constructor: Check for symmetry, then construct the eigenvalue decomposition.
781     *
782     * @param mixed $Arg A Square matrix
783     */
784    public function __construct($Arg)
785    {
786        $this->A = $Arg->getArray();
787        $this->n = $Arg->getColumnDimension();
788
789        $issymmetric = true;
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]);
793            }
794        }
795
796        if ($issymmetric) {
797            $this->V = $this->A;
798            // Tridiagonalize.
799            $this->tred2();
800            // Diagonalize.
801            $this->tql2();
802        } else {
803            $this->H = $this->A;
804            $this->ort = [];
805            // Reduce to Hessenberg form.
806            $this->orthes();
807            // Reduce Hessenberg to real Schur form.
808            $this->hqr2();
809        }
810    }
811
812    /**
813     * Return the eigenvector matrix.
814     *
815     * @return Matrix V
816     */
817    public function getV()
818    {
819        return new Matrix($this->V, $this->n, $this->n);
820    }
821
822    /**
823     * Return the real parts of the eigenvalues.
824     *
825     * @return array real(diag(D))
826     */
827    public function getRealEigenvalues()
828    {
829        return $this->d;
830    }
831
832    /**
833     * Return the imaginary parts of the eigenvalues.
834     *
835     * @return array imag(diag(D))
836     */
837    public function getImagEigenvalues()
838    {
839        return $this->e;
840    }
841
842    /**
843     * Return the block diagonal eigenvalue matrix.
844     *
845     * @return Matrix D
846     */
847    public function getD()
848    {
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) {
853                continue;
854            }
855            $o = ($this->e[$i] > 0) ? $i + 1 : $i - 1;
856            $D[$i][$o] = $this->e[$i];
857        }
858
859        return new Matrix($D);
860    }
861}
862