EigenvalueDecomposition::hqr2() - Code Metrics - Inspection of "Update CHANGELOG with #238 fix" - php-ai/php-ml - Measure and Improve Code Quality continuously with Scrutinizer

Passed

Branch — master (83f3e8)

by Arkadiusz

created 2018-02-23 22:04 UTC

EigenvalueDecomposition::hqr2() F

↳ Parent: Project

Complexity

Conditions	71
Paths	> 20000

Size

Total Lines	414
Code Lines	270

Duplication

Lines	48
Ratio	11.59 %

Importance

Changes

Metric	Value
dl	48
loc	414
rs	2
c	0
b	0
f	0
cc	71
eloc	270
nc	1766255
nop	0

How to fix Long Method Complexity

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

php-ai / php-ml

Branch — master (83f3e8)

EigenvalueDecomposition::hqr2() F

Complexity

Size

Duplication

Importance

How to fix Long Method Complexity

Long Method

Duplication Side-by-Side

Filter issues like