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