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

Math/LinearAlgebra/EigenvalueDecomposition.php (1 issue)

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) {
136
                $sum += $vect[$i] ** 2;
137
            }
138
139
            $sum = sqrt($sum);
140 View Code Duplication
            for ($i = 0; $i < count($vect); ++$i) {
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);
0 ignored issues
show
The variable $j does not seem to be defined for all execution paths leading up to this point.

If you define a variable conditionally, it can happen that it is not defined for all execution paths.

Let’s take a look at an example:

function myFunction($a) {
    switch ($a) {
        case 'foo':
            $x = 1;
            break;

        case 'bar':
            $x = 2;
            break;
    }

    // $x is potentially undefined here.
    echo $x;
}

In the above example, the variable $x is defined if you pass “foo” or “bar” as argument for $a. However, since the switch statement has no default case statement, if you pass any other value, the variable $x would be undefined.

Available Fixes

  1. Check for existence of the variable explicitly:

    function myFunction($a) {
        switch ($a) {
            case 'foo':
                $x = 1;
                break;
    
            case 'bar':
                $x = 2;
                break;
        }
    
        if (isset($x)) { // Make sure it's always set.
            echo $x;
        }
    }
    
  2. Define a default value for the variable:

    function myFunction($a) {
        $x = ''; // Set a default which gets overridden for certain paths.
        switch ($a) {
            case 'foo':
                $x = 1;
                break;
    
            case 'bar':
                $x = 2;
                break;
        }
    
        echo $x;
    }
    
  3. Add a value for the missing path:

    function myFunction($a) {
        switch ($a) {
            case 'foo':
                $x = 1;
                break;
    
            case 'bar':
                $x = 2;
                break;
    
            // We add support for the missing case.
            default:
                $x = '';
                break;
        }
    
        echo $x;
    }
    
Loading history...
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