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); |
||
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) { |
||
0 ignored issues
–
show
|
|||
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 |
When comparing the result of a bit operation, we suggest to add explicit parenthesis and not to rely on PHP’s built-in operator precedence to ensure the code behaves as intended and to make it more readable.
Let’s take a look at these examples: