These results are based on our legacy PHP analysis, consider migrating to our new PHP analysis engine instead. Learn more
1 | <?php |
|||||||||||
2 | ||||||||||||
3 | declare(strict_types=1); |
|||||||||||
4 | ||||||||||||
5 | /** |
|||||||||||
6 | * Class to obtain eigenvalues and eigenvectors of a real matrix. |
|||||||||||
7 | * |
|||||||||||
8 | * If A is symmetric, then A = V*D*V' where the eigenvalue matrix D |
|||||||||||
9 | * is diagonal and the eigenvector matrix V is orthogonal (i.e. |
|||||||||||
10 | * A = V.times(D.times(V.transpose())) and V.times(V.transpose()) |
|||||||||||
11 | * equals the identity matrix). |
|||||||||||
12 | * |
|||||||||||
13 | * If A is not symmetric, then the eigenvalue matrix D is block diagonal |
|||||||||||
14 | * with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues, |
|||||||||||
15 | * lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda]. The |
|||||||||||
16 | * columns of V represent the eigenvectors in the sense that A*V = V*D, |
|||||||||||
17 | * i.e. A.times(V) equals V.times(D). The matrix V may be badly |
|||||||||||
18 | * conditioned, or even singular, so the validity of the equation |
|||||||||||
19 | * A = V*D*inverse(V) depends upon V.cond(). |
|||||||||||
20 | * |
|||||||||||
21 | * @author Paul Meagher |
|||||||||||
22 | * @license PHP v3.0 |
|||||||||||
23 | * |
|||||||||||
24 | * @version 1.1 |
|||||||||||
25 | * |
|||||||||||
26 | * Slightly changed to adapt the original code to PHP-ML library |
|||||||||||
27 | * @date 2017/04/11 |
|||||||||||
28 | * |
|||||||||||
29 | * @author Mustafa Karabulut |
|||||||||||
30 | */ |
|||||||||||
31 | ||||||||||||
32 | namespace Phpml\Math\LinearAlgebra; |
|||||||||||
33 | ||||||||||||
34 | use Phpml\Math\Matrix; |
|||||||||||
35 | ||||||||||||
36 | class EigenvalueDecomposition |
|||||||||||
37 | { |
|||||||||||
38 | /** |
|||||||||||
39 | * Row and column dimension (square matrix). |
|||||||||||
40 | * |
|||||||||||
41 | * @var int |
|||||||||||
42 | */ |
|||||||||||
43 | private $n; |
|||||||||||
44 | ||||||||||||
45 | /** |
|||||||||||
46 | * Internal symmetry flag. |
|||||||||||
47 | * |
|||||||||||
48 | * @var bool |
|||||||||||
49 | */ |
|||||||||||
50 | private $symmetric; |
|||||||||||
51 | ||||||||||||
52 | /** |
|||||||||||
53 | * Arrays for internal storage of eigenvalues. |
|||||||||||
54 | * |
|||||||||||
55 | * @var array |
|||||||||||
56 | */ |
|||||||||||
57 | private $d = []; |
|||||||||||
58 | ||||||||||||
59 | private $e = []; |
|||||||||||
60 | ||||||||||||
61 | /** |
|||||||||||
62 | * Array for internal storage of eigenvectors. |
|||||||||||
63 | * |
|||||||||||
64 | * @var array |
|||||||||||
65 | */ |
|||||||||||
66 | private $V = []; |
|||||||||||
67 | ||||||||||||
68 | /** |
|||||||||||
69 | * Array for internal storage of nonsymmetric Hessenberg form. |
|||||||||||
70 | * |
|||||||||||
71 | * @var array |
|||||||||||
72 | */ |
|||||||||||
73 | private $H = []; |
|||||||||||
74 | ||||||||||||
75 | /** |
|||||||||||
76 | * Working storage for nonsymmetric algorithm. |
|||||||||||
77 | * |
|||||||||||
78 | * @var array |
|||||||||||
79 | */ |
|||||||||||
80 | private $ort = []; |
|||||||||||
81 | ||||||||||||
82 | /** |
|||||||||||
83 | * Used for complex scalar division. |
|||||||||||
84 | * |
|||||||||||
85 | * @var float |
|||||||||||
86 | */ |
|||||||||||
87 | private $cdivr; |
|||||||||||
88 | ||||||||||||
89 | private $cdivi; |
|||||||||||
90 | ||||||||||||
91 | private $A; |
|||||||||||
92 | ||||||||||||
93 | /** |
|||||||||||
94 | * Constructor: Check for symmetry, then construct the eigenvalue decomposition |
|||||||||||
95 | */ |
|||||||||||
96 | public function __construct(array $Arg) |
|||||||||||
97 | { |
|||||||||||
98 | $this->A = $Arg; |
|||||||||||
99 | $this->n = count($Arg[0]); |
|||||||||||
100 | $this->symmetric = true; |
|||||||||||
101 | ||||||||||||
102 | for ($j = 0; ($j < $this->n) && $this->symmetric; ++$j) { |
|||||||||||
103 | for ($i = 0; ($i < $this->n) & $this->symmetric; ++$i) { |
|||||||||||
104 | $this->symmetric = ($this->A[$i][$j] == $this->A[$j][$i]); |
|||||||||||
105 | } |
|||||||||||
106 | } |
|||||||||||
107 | ||||||||||||
108 | if ($this->symmetric) { |
|||||||||||
109 | $this->V = $this->A; |
|||||||||||
110 | // Tridiagonalize. |
|||||||||||
111 | $this->tred2(); |
|||||||||||
112 | // Diagonalize. |
|||||||||||
113 | $this->tql2(); |
|||||||||||
114 | } else { |
|||||||||||
115 | $this->H = $this->A; |
|||||||||||
116 | $this->ort = []; |
|||||||||||
117 | // Reduce to Hessenberg form. |
|||||||||||
118 | $this->orthes(); |
|||||||||||
119 | // Reduce Hessenberg to real Schur form. |
|||||||||||
120 | $this->hqr2(); |
|||||||||||
121 | } |
|||||||||||
122 | } |
|||||||||||
123 | ||||||||||||
124 | /** |
|||||||||||
125 | * Return the eigenvector matrix |
|||||||||||
126 | */ |
|||||||||||
127 | public function getEigenvectors(): array |
|||||||||||
128 | { |
|||||||||||
129 | $vectors = $this->V; |
|||||||||||
130 | ||||||||||||
131 | // Always return the eigenvectors of length 1.0 |
|||||||||||
132 | $vectors = new Matrix($vectors); |
|||||||||||
133 | $vectors = array_map(function ($vect) { |
|||||||||||
134 | $sum = 0; |
|||||||||||
135 | View Code Duplication | for ($i = 0; $i < count($vect); ++$i) { |
||||||||||
0 ignored issues
–
show
This code seems to be duplicated across your project.
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository.
Loading history...
|
||||||||||||
136 | $sum += $vect[$i] ** 2; |
|||||||||||
137 | } |
|||||||||||
138 | ||||||||||||
139 | $sum = sqrt($sum); |
|||||||||||
140 | View Code Duplication | for ($i = 0; $i < count($vect); ++$i) { |
||||||||||
0 ignored issues
–
show
It seems like you are calling the size function
count() as part of the test condition. You might want to compute the size beforehand, and not on each iteration.
If the size of the collection does not change during the iteration, it is generally a good practice to compute it beforehand, and not on each iteration: for ($i=0; $i<count($array); $i++) { // calls count() on each iteration
}
// Better
for ($i=0, $c=count($array); $i<$c; $i++) { // calls count() just once
}
Loading history...
This code seems to be duplicated across your project.
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository.
Loading history...
|
||||||||||||
141 | $vect[$i] /= $sum; |
|||||||||||
142 | } |
|||||||||||
143 | ||||||||||||
144 | return $vect; |
|||||||||||
145 | }, $vectors->transpose()->toArray()); |
|||||||||||
146 | ||||||||||||
147 | return $vectors; |
|||||||||||
148 | } |
|||||||||||
149 | ||||||||||||
150 | /** |
|||||||||||
151 | * Return the real parts of the eigenvalues<br> |
|||||||||||
152 | * d = real(diag(D)); |
|||||||||||
153 | */ |
|||||||||||
154 | public function getRealEigenvalues(): array |
|||||||||||
155 | { |
|||||||||||
156 | return $this->d; |
|||||||||||
157 | } |
|||||||||||
158 | ||||||||||||
159 | /** |
|||||||||||
160 | * Return the imaginary parts of the eigenvalues <br> |
|||||||||||
161 | * d = imag(diag(D)) |
|||||||||||
162 | */ |
|||||||||||
163 | public function getImagEigenvalues(): array |
|||||||||||
164 | { |
|||||||||||
165 | return $this->e; |
|||||||||||
166 | } |
|||||||||||
167 | ||||||||||||
168 | /** |
|||||||||||
169 | * Return the block diagonal eigenvalue matrix |
|||||||||||
170 | */ |
|||||||||||
171 | public function getDiagonalEigenvalues(): array |
|||||||||||
172 | { |
|||||||||||
173 | $D = []; |
|||||||||||
174 | ||||||||||||
175 | for ($i = 0; $i < $this->n; ++$i) { |
|||||||||||
176 | $D[$i] = array_fill(0, $this->n, 0.0); |
|||||||||||
177 | $D[$i][$i] = $this->d[$i]; |
|||||||||||
178 | if ($this->e[$i] == 0) { |
|||||||||||
179 | continue; |
|||||||||||
180 | } |
|||||||||||
181 | ||||||||||||
182 | $o = ($this->e[$i] > 0) ? $i + 1 : $i - 1; |
|||||||||||
183 | $D[$i][$o] = $this->e[$i]; |
|||||||||||
184 | } |
|||||||||||
185 | ||||||||||||
186 | return $D; |
|||||||||||
187 | } |
|||||||||||
188 | ||||||||||||
189 | /** |
|||||||||||
190 | * Symmetric Householder reduction to tridiagonal form. |
|||||||||||
191 | */ |
|||||||||||
192 | private function tred2(): void |
|||||||||||
193 | { |
|||||||||||
194 | // This is derived from the Algol procedures tred2 by |
|||||||||||
195 | // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for |
|||||||||||
196 | // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding |
|||||||||||
197 | // Fortran subroutine in EISPACK. |
|||||||||||
198 | $this->d = $this->V[$this->n - 1]; |
|||||||||||
199 | // Householder reduction to tridiagonal form. |
|||||||||||
200 | for ($i = $this->n - 1; $i > 0; --$i) { |
|||||||||||
201 | $i_ = $i - 1; |
|||||||||||
202 | // Scale to avoid under/overflow. |
|||||||||||
203 | $h = $scale = 0.0; |
|||||||||||
204 | $scale += array_sum(array_map('abs', $this->d)); |
|||||||||||
205 | if ($scale == 0.0) { |
|||||||||||
206 | $this->e[$i] = $this->d[$i_]; |
|||||||||||
207 | $this->d = array_slice($this->V[$i_], 0, $i_); |
|||||||||||
208 | for ($j = 0; $j < $i; ++$j) { |
|||||||||||
209 | $this->V[$j][$i] = $this->V[$i][$j] = 0.0; |
|||||||||||
210 | } |
|||||||||||
211 | } else { |
|||||||||||
212 | // Generate Householder vector. |
|||||||||||
213 | for ($k = 0; $k < $i; ++$k) { |
|||||||||||
214 | $this->d[$k] /= $scale; |
|||||||||||
215 | $h += pow($this->d[$k], 2); |
|||||||||||
216 | } |
|||||||||||
217 | ||||||||||||
218 | $f = $this->d[$i_]; |
|||||||||||
219 | $g = sqrt($h); |
|||||||||||
220 | if ($f > 0) { |
|||||||||||
221 | $g = -$g; |
|||||||||||
222 | } |
|||||||||||
223 | ||||||||||||
224 | $this->e[$i] = $scale * $g; |
|||||||||||
225 | $h = $h - $f * $g; |
|||||||||||
226 | $this->d[$i_] = $f - $g; |
|||||||||||
227 | ||||||||||||
228 | for ($j = 0; $j < $i; ++$j) { |
|||||||||||
229 | $this->e[$j] = 0.0; |
|||||||||||
230 | } |
|||||||||||
231 | ||||||||||||
232 | // Apply similarity transformation to remaining columns. |
|||||||||||
233 | for ($j = 0; $j < $i; ++$j) { |
|||||||||||
234 | $f = $this->d[$j]; |
|||||||||||
235 | $this->V[$j][$i] = $f; |
|||||||||||
236 | $g = $this->e[$j] + $this->V[$j][$j] * $f; |
|||||||||||
237 | ||||||||||||
238 | for ($k = $j + 1; $k <= $i_; ++$k) { |
|||||||||||
239 | $g += $this->V[$k][$j] * $this->d[$k]; |
|||||||||||
240 | $this->e[$k] += $this->V[$k][$j] * $f; |
|||||||||||
241 | } |
|||||||||||
242 | ||||||||||||
243 | $this->e[$j] = $g; |
|||||||||||
244 | } |
|||||||||||
245 | ||||||||||||
246 | $f = 0.0; |
|||||||||||
247 | if ($h === 0 || $h < 1e-32) { |
|||||||||||
248 | $h = 1e-32; |
|||||||||||
249 | } |
|||||||||||
250 | ||||||||||||
251 | for ($j = 0; $j < $i; ++$j) { |
|||||||||||
252 | $this->e[$j] /= $h; |
|||||||||||
253 | $f += $this->e[$j] * $this->d[$j]; |
|||||||||||
254 | } |
|||||||||||
255 | ||||||||||||
256 | $hh = $f / (2 * $h); |
|||||||||||
257 | for ($j = 0; $j < $i; ++$j) { |
|||||||||||
258 | $this->e[$j] -= $hh * $this->d[$j]; |
|||||||||||
259 | } |
|||||||||||
260 | ||||||||||||
261 | for ($j = 0; $j < $i; ++$j) { |
|||||||||||
262 | $f = $this->d[$j]; |
|||||||||||
263 | $g = $this->e[$j]; |
|||||||||||
264 | for ($k = $j; $k <= $i_; ++$k) { |
|||||||||||
265 | $this->V[$k][$j] -= ($f * $this->e[$k] + $g * $this->d[$k]); |
|||||||||||
266 | } |
|||||||||||
267 | ||||||||||||
268 | $this->d[$j] = $this->V[$i - 1][$j]; |
|||||||||||
269 | $this->V[$i][$j] = 0.0; |
|||||||||||
270 | } |
|||||||||||
271 | } |
|||||||||||
272 | ||||||||||||
273 | $this->d[$i] = $h; |
|||||||||||
274 | } |
|||||||||||
275 | ||||||||||||
276 | // Accumulate transformations. |
|||||||||||
277 | for ($i = 0; $i < $this->n - 1; ++$i) { |
|||||||||||
278 | $this->V[$this->n - 1][$i] = $this->V[$i][$i]; |
|||||||||||
279 | $this->V[$i][$i] = 1.0; |
|||||||||||
280 | $h = $this->d[$i + 1]; |
|||||||||||
281 | if ($h != 0.0) { |
|||||||||||
282 | for ($k = 0; $k <= $i; ++$k) { |
|||||||||||
283 | $this->d[$k] = $this->V[$k][$i + 1] / $h; |
|||||||||||
284 | } |
|||||||||||
285 | ||||||||||||
286 | for ($j = 0; $j <= $i; ++$j) { |
|||||||||||
287 | $g = 0.0; |
|||||||||||
288 | for ($k = 0; $k <= $i; ++$k) { |
|||||||||||
289 | $g += $this->V[$k][$i + 1] * $this->V[$k][$j]; |
|||||||||||
290 | } |
|||||||||||
291 | ||||||||||||
292 | for ($k = 0; $k <= $i; ++$k) { |
|||||||||||
293 | $this->V[$k][$j] -= $g * $this->d[$k]; |
|||||||||||
294 | } |
|||||||||||
295 | } |
|||||||||||
296 | } |
|||||||||||
297 | ||||||||||||
298 | for ($k = 0; $k <= $i; ++$k) { |
|||||||||||
299 | $this->V[$k][$i + 1] = 0.0; |
|||||||||||
300 | } |
|||||||||||
301 | } |
|||||||||||
302 | ||||||||||||
303 | $this->d = $this->V[$this->n - 1]; |
|||||||||||
304 | $this->V[$this->n - 1] = array_fill(0, $j, 0.0); |
|||||||||||
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
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) { |
||||||||||
0 ignored issues
–
show
This code seems to be duplicated across your project.
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository.
Loading history...
|
||||||||||||
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) { |
||||||||||
0 ignored issues
–
show
This code seems to be duplicated across your project.
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository.
Loading history...
|
||||||||||||
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) { |
||||||||||
0 ignored issues
–
show
This code seems to be duplicated across your project.
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository.
Loading history...
|
||||||||||||
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) { |
||||||||||
0 ignored issues
–
show
This code seems to be duplicated across your project.
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository.
Loading history...
|
||||||||||||
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) { |
||||||||||
0 ignored issues
–
show
This code seems to be duplicated across your project.
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository.
Loading history...
|
||||||||||||
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)) { |
|||||||||||
0 ignored issues
–
show
Comprehensibility
Best Practice
introduced
by
Using logical operators such as
or instead of || is generally not recommended.
PHP has two types of connecting operators (logical operators, and boolean operators):
The difference between these is the order in which they are executed. In most cases,
you would want to use a boolean operator like Let’s take a look at a few examples: // Logical operators have lower precedence:
$f = false or true;
// is executed like this:
($f = false) or true;
// Boolean operators have higher precedence:
$f = false || true;
// is executed like this:
$f = (false || true);
Logical Operators are used for Control-FlowOne case where you explicitly want to use logical operators is for control-flow such as this: $x === 5
or die('$x must be 5.');
// Instead of
if ($x !== 5) {
die('$x must be 5.');
}
Since // The following is currently a parse error.
$x === 5
or throw new RuntimeException('$x must be 5.');
These limitations lead to logical operators rarely being of use in current PHP code.
Loading history...
|
||||||||||||
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) { |
||||||||||
0 ignored issues
–
show
This code seems to be duplicated across your project.
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository.
Loading history...
|
||||||||||||
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) { |
||||||||||
0 ignored issues
–
show
This code seems to be duplicated across your project.
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository.
Loading history...
|
||||||||||||
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) { |
||||||||||
0 ignored issues
–
show
This code seems to be duplicated across your project.
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository.
Loading history...
|
||||||||||||
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) { |
||||||||||
0 ignored issues
–
show
This code seems to be duplicated across your project.
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository.
Loading history...
|
||||||||||||
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) { |
||||||||||
0 ignored issues
–
show
This code seems to be duplicated across your project.
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository.
Loading history...
|
||||||||||||
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) { |
||||||||||
0 ignored issues
–
show
This code seems to be duplicated across your project.
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository.
Loading history...
|
||||||||||||
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) { |
||||||||||
0 ignored issues
–
show
This code seems to be duplicated across your project.
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository.
Loading history...
|
||||||||||||
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) { |
||||||||||
0 ignored issues
–
show
This code seems to be duplicated across your project.
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository.
Loading history...
|
||||||||||||
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 |
If the size of the collection does not change during the iteration, it is generally a good practice to compute it beforehand, and not on each iteration: