Passed
Push — master ( e83f7b...d953ef )
by Arkadiusz
03:28
created

Math/LinearAlgebra/EigenvalueDecomposition.php (2 issues)

Upgrade to new PHP Analysis Engine

These results are based on our legacy PHP analysis, consider migrating to our new PHP analysis engine instead. Learn more

1
<?php
2
3
declare(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
32
namespace Phpml\Math\LinearAlgebra;
33
34
use Phpml\Math\Matrix;
35
36
class EigenvalueDecomposition
37
{
38
    /**
39
     *	Row and column dimension (square matrix).
40
     *
41
     *	@var int
42
     */
43
    private $n;
44
45
    /**
46
     *	Internal symmetry flag.
47
     *
48
     *	@var bool
49
     */
50
    private $symmetric;
51
52
    /**
53
     *	Arrays for internal storage of eigenvalues.
54
     *
55
     *	@var array
56
     */
57
    private $d = [];
58
59
    private $e = [];
60
61
    /**
62
     *	Array for internal storage of eigenvectors.
63
     *
64
     *	@var array
65
     */
66
    private $V = [];
67
68
    /**
69
     *	Array for internal storage of nonsymmetric Hessenberg form.
70
     *
71
     *	@var array
72
     */
73
    private $H = [];
74
75
    /**
76
     *	Working storage for nonsymmetric algorithm.
77
     *
78
     *	@var array
79
     */
80
    private $ort = [];
81
82
    /**
83
     *	Used for complex scalar division.
84
     *
85
     *	@var float
86
     */
87
    private $cdivr;
88
89
    private $cdivi;
90
91
    private $A;
92
93
    /**
94
     *	Constructor: Check for symmetry, then construct the eigenvalue decomposition
95
     */
96
    public function __construct(array $Arg)
97
    {
98
        $this->A = $Arg;
99
        $this->n = count($Arg[0]);
100
        $this->symmetric = true;
101
102
        for ($j = 0; ($j < $this->n) && $this->symmetric; ++$j) {
103
            for ($i = 0; ($i < $this->n) & $this->symmetric; ++$i) {
104
                $this->symmetric = ($this->A[$i][$j] == $this->A[$j][$i]);
105
            }
106
        }
107
108
        if ($this->symmetric) {
109
            $this->V = $this->A;
110
            // Tridiagonalize.
111
            $this->tred2();
112
            // Diagonalize.
113
            $this->tql2();
114
        } else {
115
            $this->H = $this->A;
116
            $this->ort = [];
117
            // Reduce to Hessenberg form.
118
            $this->orthes();
119
            // Reduce Hessenberg to real Schur form.
120
            $this->hqr2();
121
        }
122
    }
123
124
    /**
125
     * Return the eigenvector matrix
126
     */
127
    public function getEigenvectors(): array
128
    {
129
        $vectors = $this->V;
130
131
        // Always return the eigenvectors of length 1.0
132
        $vectors = new Matrix($vectors);
133
        $vectors = array_map(function ($vect) {
134
            $sum = 0;
135 View Code Duplication
            for ($i = 0; $i < count($vect); ++$i) {
0 ignored issues
show
Performance Best Practice introduced by
It seems like you are calling the size function count() as part of the test condition. You might want to compute the size beforehand, and not on each iteration.

If the size of the collection does not change during the iteration, it is generally a good practice to compute it beforehand, and not on each iteration:

for ($i=0; $i<count($array); $i++) { // calls count() on each iteration
}

// Better
for ($i=0, $c=count($array); $i<$c; $i++) { // calls count() just once
}
Loading history...
136
                $sum += $vect[$i] ** 2;
137
            }
138
139
            $sum = sqrt($sum);
140 View Code Duplication
            for ($i = 0; $i < count($vect); ++$i) {
0 ignored issues
show
Performance Best Practice introduced by
It seems like you are calling the size function count() as part of the test condition. You might want to compute the size beforehand, and not on each iteration.

If the size of the collection does not change during the iteration, it is generally a good practice to compute it beforehand, and not on each iteration:

for ($i=0; $i<count($array); $i++) { // calls count() on each iteration
}

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