ILIAS  eassessment Revision 61809
 All Data Structures Namespaces Files Functions Variables Groups Pages
EigenvalueDecomposition.php
Go to the documentation of this file.
1 <?php
24 class EigenvalueDecomposition {
25 
30  private $n;
31 
36  private $issymmetric;
37 
42  private $d = array();
43  private $e = array();
44 
49  private $V = array();
50 
55  private $H = array();
56 
61  private $ort;
62 
67  private $cdivr;
68  private $cdivi;
69 
70 
76  private function tred2 () {
77  // This is derived from the Algol procedures tred2 by
78  // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
79  // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
80  // Fortran subroutine in EISPACK.
81  $this->d = $this->V[$this->n-1];
82  // Householder reduction to tridiagonal form.
83  for ($i = $this->n-1; $i > 0; --$i) {
84  $i_ = $i -1;
85  // Scale to avoid under/overflow.
86  $h = $scale = 0.0;
87  $scale += array_sum(array_map(abs, $this->d));
88  if ($scale == 0.0) {
89  $this->e[$i] = $this->d[$i_];
90  $this->d = array_slice($this->V[$i_], 0, $i_);
91  for ($j = 0; $j < $i; ++$j) {
92  $this->V[$j][$i] = $this->V[$i][$j] = 0.0;
93  }
94  } else {
95  // Generate Householder vector.
96  for ($k = 0; $k < $i; ++$k) {
97  $this->d[$k] /= $scale;
98  $h += pow($this->d[$k], 2);
99  }
100  $f = $this->d[$i_];
101  $g = sqrt($h);
102  if ($f > 0) {
103  $g = -$g;
104  }
105  $this->e[$i] = $scale * $g;
106  $h = $h - $f * $g;
107  $this->d[$i_] = $f - $g;
108  for ($j = 0; $j < $i; ++$j) {
109  $this->e[$j] = 0.0;
110  }
111  // Apply similarity transformation to remaining columns.
112  for ($j = 0; $j < $i; ++$j) {
113  $f = $this->d[$j];
114  $this->V[$j][$i] = $f;
115  $g = $this->e[$j] + $this->V[$j][$j] * $f;
116  for ($k = $j+1; $k <= $i_; ++$k) {
117  $g += $this->V[$k][$j] * $this->d[$k];
118  $this->e[$k] += $this->V[$k][$j] * $f;
119  }
120  $this->e[$j] = $g;
121  }
122  $f = 0.0;
123  for ($j = 0; $j < $i; ++$j) {
124  $this->e[$j] /= $h;
125  $f += $this->e[$j] * $this->d[$j];
126  }
127  $hh = $f / (2 * $h);
128  for ($j=0; $j < $i; ++$j) {
129  $this->e[$j] -= $hh * $this->d[$j];
130  }
131  for ($j = 0; $j < $i; ++$j) {
132  $f = $this->d[$j];
133  $g = $this->e[$j];
134  for ($k = $j; $k <= $i_; ++$k) {
135  $this->V[$k][$j] -= ($f * $this->e[$k] + $g * $this->d[$k]);
136  }
137  $this->d[$j] = $this->V[$i-1][$j];
138  $this->V[$i][$j] = 0.0;
139  }
140  }
141  $this->d[$i] = $h;
142  }
143 
144  // Accumulate transformations.
145  for ($i = 0; $i < $this->n-1; ++$i) {
146  $this->V[$this->n-1][$i] = $this->V[$i][$i];
147  $this->V[$i][$i] = 1.0;
148  $h = $this->d[$i+1];
149  if ($h != 0.0) {
150  for ($k = 0; $k <= $i; ++$k) {
151  $this->d[$k] = $this->V[$k][$i+1] / $h;
152  }
153  for ($j = 0; $j <= $i; ++$j) {
154  $g = 0.0;
155  for ($k = 0; $k <= $i; ++$k) {
156  $g += $this->V[$k][$i+1] * $this->V[$k][$j];
157  }
158  for ($k = 0; $k <= $i; ++$k) {
159  $this->V[$k][$j] -= $g * $this->d[$k];
160  }
161  }
162  }
163  for ($k = 0; $k <= $i; ++$k) {
164  $this->V[$k][$i+1] = 0.0;
165  }
166  }
167 
168  $this->d = $this->V[$this->n-1];
169  $this->V[$this->n-1] = array_fill(0, $j, 0.0);
170  $this->V[$this->n-1][$this->n-1] = 1.0;
171  $this->e[0] = 0.0;
172  }
173 
174 
185  private function tql2() {
186  for ($i = 1; $i < $this->n; ++$i) {
187  $this->e[$i-1] = $this->e[$i];
188  }
189  $this->e[$this->n-1] = 0.0;
190  $f = 0.0;
191  $tst1 = 0.0;
192  $eps = pow(2.0,-52.0);
193 
194  for ($l = 0; $l < $this->n; ++$l) {
195  // Find small subdiagonal element
196  $tst1 = max($tst1, abs($this->d[$l]) + abs($this->e[$l]));
197  $m = $l;
198  while ($m < $this->n) {
199  if (abs($this->e[$m]) <= $eps * $tst1)
200  break;
201  ++$m;
202  }
203  // If m == l, $this->d[l] is an eigenvalue,
204  // otherwise, iterate.
205  if ($m > $l) {
206  $iter = 0;
207  do {
208  // Could check iteration count here.
209  $iter += 1;
210  // Compute implicit shift
211  $g = $this->d[$l];
212  $p = ($this->d[$l+1] - $g) / (2.0 * $this->e[$l]);
213  $r = hypo($p, 1.0);
214  if ($p < 0)
215  $r *= -1;
216  $this->d[$l] = $this->e[$l] / ($p + $r);
217  $this->d[$l+1] = $this->e[$l] * ($p + $r);
218  $dl1 = $this->d[$l+1];
219  $h = $g - $this->d[$l];
220  for ($i = $l + 2; $i < $this->n; ++$i)
221  $this->d[$i] -= $h;
222  $f += $h;
223  // Implicit QL transformation.
224  $p = $this->d[$m];
225  $c = 1.0;
226  $c2 = $c3 = $c;
227  $el1 = $this->e[$l + 1];
228  $s = $s2 = 0.0;
229  for ($i = $m-1; $i >= $l; --$i) {
230  $c3 = $c2;
231  $c2 = $c;
232  $s2 = $s;
233  $g = $c * $this->e[$i];
234  $h = $c * $p;
235  $r = hypo($p, $this->e[$i]);
236  $this->e[$i+1] = $s * $r;
237  $s = $this->e[$i] / $r;
238  $c = $p / $r;
239  $p = $c * $this->d[$i] - $s * $g;
240  $this->d[$i+1] = $h + $s * ($c * $g + $s * $this->d[$i]);
241  // Accumulate transformation.
242  for ($k = 0; $k < $this->n; ++$k) {
243  $h = $this->V[$k][$i+1];
244  $this->V[$k][$i+1] = $s * $this->V[$k][$i] + $c * $h;
245  $this->V[$k][$i] = $c * $this->V[$k][$i] - $s * $h;
246  }
247  }
248  $p = -$s * $s2 * $c3 * $el1 * $this->e[$l] / $dl1;
249  $this->e[$l] = $s * $p;
250  $this->d[$l] = $c * $p;
251  // Check for convergence.
252  } while (abs($this->e[$l]) > $eps * $tst1);
253  }
254  $this->d[$l] = $this->d[$l] + $f;
255  $this->e[$l] = 0.0;
256  }
257 
258  // Sort eigenvalues and corresponding vectors.
259  for ($i = 0; $i < $this->n - 1; ++$i) {
260  $k = $i;
261  $p = $this->d[$i];
262  for ($j = $i+1; $j < $this->n; ++$j) {
263  if ($this->d[$j] < $p) {
264  $k = $j;
265  $p = $this->d[$j];
266  }
267  }
268  if ($k != $i) {
269  $this->d[$k] = $this->d[$i];
270  $this->d[$i] = $p;
271  for ($j = 0; $j < $this->n; ++$j) {
272  $p = $this->V[$j][$i];
273  $this->V[$j][$i] = $this->V[$j][$k];
274  $this->V[$j][$k] = $p;
275  }
276  }
277  }
278  }
279 
280 
291  private function orthes () {
292  $low = 0;
293  $high = $this->n-1;
294 
295  for ($m = $low+1; $m <= $high-1; ++$m) {
296  // Scale column.
297  $scale = 0.0;
298  for ($i = $m; $i <= $high; ++$i) {
299  $scale = $scale + abs($this->H[$i][$m-1]);
300  }
301  if ($scale != 0.0) {
302  // Compute Householder transformation.
303  $h = 0.0;
304  for ($i = $high; $i >= $m; --$i) {
305  $this->ort[$i] = $this->H[$i][$m-1] / $scale;
306  $h += $this->ort[$i] * $this->ort[$i];
307  }
308  $g = sqrt($h);
309  if ($this->ort[$m] > 0) {
310  $g *= -1;
311  }
312  $h -= $this->ort[$m] * $g;
313  $this->ort[$m] -= $g;
314  // Apply Householder similarity transformation
315  // H = (I -u * u' / h) * H * (I -u * u') / h)
316  for ($j = $m; $j < $this->n; ++$j) {
317  $f = 0.0;
318  for ($i = $high; $i >= $m; --$i) {
319  $f += $this->ort[$i] * $this->H[$i][$j];
320  }
321  $f /= $h;
322  for ($i = $m; $i <= $high; ++$i) {
323  $this->H[$i][$j] -= $f * $this->ort[$i];
324  }
325  }
326  for ($i = 0; $i <= $high; ++$i) {
327  $f = 0.0;
328  for ($j = $high; $j >= $m; --$j) {
329  $f += $this->ort[$j] * $this->H[$i][$j];
330  }
331  $f = $f / $h;
332  for ($j = $m; $j <= $high; ++$j) {
333  $this->H[$i][$j] -= $f * $this->ort[$j];
334  }
335  }
336  $this->ort[$m] = $scale * $this->ort[$m];
337  $this->H[$m][$m-1] = $scale * $g;
338  }
339  }
340 
341  // Accumulate transformations (Algol's ortran).
342  for ($i = 0; $i < $this->n; ++$i) {
343  for ($j = 0; $j < $this->n; ++$j) {
344  $this->V[$i][$j] = ($i == $j ? 1.0 : 0.0);
345  }
346  }
347  for ($m = $high-1; $m >= $low+1; --$m) {
348  if ($this->H[$m][$m-1] != 0.0) {
349  for ($i = $m+1; $i <= $high; ++$i) {
350  $this->ort[$i] = $this->H[$i][$m-1];
351  }
352  for ($j = $m; $j <= $high; ++$j) {
353  $g = 0.0;
354  for ($i = $m; $i <= $high; ++$i) {
355  $g += $this->ort[$i] * $this->V[$i][$j];
356  }
357  // Double division avoids possible underflow
358  $g = ($g / $this->ort[$m]) / $this->H[$m][$m-1];
359  for ($i = $m; $i <= $high; ++$i) {
360  $this->V[$i][$j] += $g * $this->ort[$i];
361  }
362  }
363  }
364  }
365  }
366 
367 
373  private function cdiv($xr, $xi, $yr, $yi) {
374  if (abs($yr) > abs($yi)) {
375  $r = $yi / $yr;
376  $d = $yr + $r * $yi;
377  $this->cdivr = ($xr + $r * $xi) / $d;
378  $this->cdivi = ($xi - $r * $xr) / $d;
379  } else {
380  $r = $yr / $yi;
381  $d = $yi + $r * $yr;
382  $this->cdivr = ($r * $xr + $xi) / $d;
383  $this->cdivi = ($r * $xi - $xr) / $d;
384  }
385  }
386 
387 
398  private function hqr2 () {
399  // Initialize
400  $nn = $this->n;
401  $n = $nn - 1;
402  $low = 0;
403  $high = $nn - 1;
404  $eps = pow(2.0, -52.0);
405  $exshift = 0.0;
406  $p = $q = $r = $s = $z = 0;
407  // Store roots isolated by balanc and compute matrix norm
408  $norm = 0.0;
409 
410  for ($i = 0; $i < $nn; ++$i) {
411  if (($i < $low) OR ($i > $high)) {
412  $this->d[$i] = $this->H[$i][$i];
413  $this->e[$i] = 0.0;
414  }
415  for ($j = max($i-1, 0); $j < $nn; ++$j) {
416  $norm = $norm + abs($this->H[$i][$j]);
417  }
418  }
419 
420  // Outer loop over eigenvalue index
421  $iter = 0;
422  while ($n >= $low) {
423  // Look for single small sub-diagonal element
424  $l = $n;
425  while ($l > $low) {
426  $s = abs($this->H[$l-1][$l-1]) + abs($this->H[$l][$l]);
427  if ($s == 0.0) {
428  $s = $norm;
429  }
430  if (abs($this->H[$l][$l-1]) < $eps * $s) {
431  break;
432  }
433  --$l;
434  }
435  // Check for convergence
436  // One root found
437  if ($l == $n) {
438  $this->H[$n][$n] = $this->H[$n][$n] + $exshift;
439  $this->d[$n] = $this->H[$n][$n];
440  $this->e[$n] = 0.0;
441  --$n;
442  $iter = 0;
443  // Two roots found
444  } else if ($l == $n-1) {
445  $w = $this->H[$n][$n-1] * $this->H[$n-1][$n];
446  $p = ($this->H[$n-1][$n-1] - $this->H[$n][$n]) / 2.0;
447  $q = $p * $p + $w;
448  $z = sqrt(abs($q));
449  $this->H[$n][$n] = $this->H[$n][$n] + $exshift;
450  $this->H[$n-1][$n-1] = $this->H[$n-1][$n-1] + $exshift;
451  $x = $this->H[$n][$n];
452  // Real pair
453  if ($q >= 0) {
454  if ($p >= 0) {
455  $z = $p + $z;
456  } else {
457  $z = $p - $z;
458  }
459  $this->d[$n-1] = $x + $z;
460  $this->d[$n] = $this->d[$n-1];
461  if ($z != 0.0) {
462  $this->d[$n] = $x - $w / $z;
463  }
464  $this->e[$n-1] = 0.0;
465  $this->e[$n] = 0.0;
466  $x = $this->H[$n][$n-1];
467  $s = abs($x) + abs($z);
468  $p = $x / $s;
469  $q = $z / $s;
470  $r = sqrt($p * $p + $q * $q);
471  $p = $p / $r;
472  $q = $q / $r;
473  // Row modification
474  for ($j = $n-1; $j < $nn; ++$j) {
475  $z = $this->H[$n-1][$j];
476  $this->H[$n-1][$j] = $q * $z + $p * $this->H[$n][$j];
477  $this->H[$n][$j] = $q * $this->H[$n][$j] - $p * $z;
478  }
479  // Column modification
480  for ($i = 0; $i <= n; ++$i) {
481  $z = $this->H[$i][$n-1];
482  $this->H[$i][$n-1] = $q * $z + $p * $this->H[$i][$n];
483  $this->H[$i][$n] = $q * $this->H[$i][$n] - $p * $z;
484  }
485  // Accumulate transformations
486  for ($i = $low; $i <= $high; ++$i) {
487  $z = $this->V[$i][$n-1];
488  $this->V[$i][$n-1] = $q * $z + $p * $this->V[$i][$n];
489  $this->V[$i][$n] = $q * $this->V[$i][$n] - $p * $z;
490  }
491  // Complex pair
492  } else {
493  $this->d[$n-1] = $x + $p;
494  $this->d[$n] = $x + $p;
495  $this->e[$n-1] = $z;
496  $this->e[$n] = -$z;
497  }
498  $n = $n - 2;
499  $iter = 0;
500  // No convergence yet
501  } else {
502  // Form shift
503  $x = $this->H[$n][$n];
504  $y = 0.0;
505  $w = 0.0;
506  if ($l < $n) {
507  $y = $this->H[$n-1][$n-1];
508  $w = $this->H[$n][$n-1] * $this->H[$n-1][$n];
509  }
510  // Wilkinson's original ad hoc shift
511  if ($iter == 10) {
512  $exshift += $x;
513  for ($i = $low; $i <= $n; ++$i) {
514  $this->H[$i][$i] -= $x;
515  }
516  $s = abs($this->H[$n][$n-1]) + abs($this->H[$n-1][$n-2]);
517  $x = $y = 0.75 * $s;
518  $w = -0.4375 * $s * $s;
519  }
520  // MATLAB's new ad hoc shift
521  if ($iter == 30) {
522  $s = ($y - $x) / 2.0;
523  $s = $s * $s + $w;
524  if ($s > 0) {
525  $s = sqrt($s);
526  if ($y < $x) {
527  $s = -$s;
528  }
529  $s = $x - $w / (($y - $x) / 2.0 + $s);
530  for ($i = $low; $i <= $n; ++$i) {
531  $this->H[$i][$i] -= $s;
532  }
533  $exshift += $s;
534  $x = $y = $w = 0.964;
535  }
536  }
537  // Could check iteration count here.
538  $iter = $iter + 1;
539  // Look for two consecutive small sub-diagonal elements
540  $m = $n - 2;
541  while ($m >= $l) {
542  $z = $this->H[$m][$m];
543  $r = $x - $z;
544  $s = $y - $z;
545  $p = ($r * $s - $w) / $this->H[$m+1][$m] + $this->H[$m][$m+1];
546  $q = $this->H[$m+1][$m+1] - $z - $r - $s;
547  $r = $this->H[$m+2][$m+1];
548  $s = abs($p) + abs($q) + abs($r);
549  $p = $p / $s;
550  $q = $q / $s;
551  $r = $r / $s;
552  if ($m == $l) {
553  break;
554  }
555  if (abs($this->H[$m][$m-1]) * (abs($q) + abs($r)) <
556  $eps * (abs($p) * (abs($this->H[$m-1][$m-1]) + abs($z) + abs($this->H[$m+1][$m+1])))) {
557  break;
558  }
559  --$m;
560  }
561  for ($i = $m + 2; $i <= $n; ++$i) {
562  $this->H[$i][$i-2] = 0.0;
563  if ($i > $m+2) {
564  $this->H[$i][$i-3] = 0.0;
565  }
566  }
567  // Double QR step involving rows l:n and columns m:n
568  for ($k = $m; $k <= $n-1; ++$k) {
569  $notlast = ($k != $n-1);
570  if ($k != $m) {
571  $p = $this->H[$k][$k-1];
572  $q = $this->H[$k+1][$k-1];
573  $r = ($notlast ? $this->H[$k+2][$k-1] : 0.0);
574  $x = abs($p) + abs($q) + abs($r);
575  if ($x != 0.0) {
576  $p = $p / $x;
577  $q = $q / $x;
578  $r = $r / $x;
579  }
580  }
581  if ($x == 0.0) {
582  break;
583  }
584  $s = sqrt($p * $p + $q * $q + $r * $r);
585  if ($p < 0) {
586  $s = -$s;
587  }
588  if ($s != 0) {
589  if ($k != $m) {
590  $this->H[$k][$k-1] = -$s * $x;
591  } elseif ($l != $m) {
592  $this->H[$k][$k-1] = -$this->H[$k][$k-1];
593  }
594  $p = $p + $s;
595  $x = $p / $s;
596  $y = $q / $s;
597  $z = $r / $s;
598  $q = $q / $p;
599  $r = $r / $p;
600  // Row modification
601  for ($j = $k; $j < $nn; ++$j) {
602  $p = $this->H[$k][$j] + $q * $this->H[$k+1][$j];
603  if ($notlast) {
604  $p = $p + $r * $this->H[$k+2][$j];
605  $this->H[$k+2][$j] = $this->H[$k+2][$j] - $p * $z;
606  }
607  $this->H[$k][$j] = $this->H[$k][$j] - $p * $x;
608  $this->H[$k+1][$j] = $this->H[$k+1][$j] - $p * $y;
609  }
610  // Column modification
611  for ($i = 0; $i <= min($n, $k+3); ++$i) {
612  $p = $x * $this->H[$i][$k] + $y * $this->H[$i][$k+1];
613  if ($notlast) {
614  $p = $p + $z * $this->H[$i][$k+2];
615  $this->H[$i][$k+2] = $this->H[$i][$k+2] - $p * $r;
616  }
617  $this->H[$i][$k] = $this->H[$i][$k] - $p;
618  $this->H[$i][$k+1] = $this->H[$i][$k+1] - $p * $q;
619  }
620  // Accumulate transformations
621  for ($i = $low; $i <= $high; ++$i) {
622  $p = $x * $this->V[$i][$k] + $y * $this->V[$i][$k+1];
623  if ($notlast) {
624  $p = $p + $z * $this->V[$i][$k+2];
625  $this->V[$i][$k+2] = $this->V[$i][$k+2] - $p * $r;
626  }
627  $this->V[$i][$k] = $this->V[$i][$k] - $p;
628  $this->V[$i][$k+1] = $this->V[$i][$k+1] - $p * $q;
629  }
630  } // ($s != 0)
631  } // k loop
632  } // check convergence
633  } // while ($n >= $low)
634 
635  // Backsubstitute to find vectors of upper triangular form
636  if ($norm == 0.0) {
637  return;
638  }
639 
640  for ($n = $nn-1; $n >= 0; --$n) {
641  $p = $this->d[$n];
642  $q = $this->e[$n];
643  // Real vector
644  if ($q == 0) {
645  $l = $n;
646  $this->H[$n][$n] = 1.0;
647  for ($i = $n-1; $i >= 0; --$i) {
648  $w = $this->H[$i][$i] - $p;
649  $r = 0.0;
650  for ($j = $l; $j <= $n; ++$j) {
651  $r = $r + $this->H[$i][$j] * $this->H[$j][$n];
652  }
653  if ($this->e[$i] < 0.0) {
654  $z = $w;
655  $s = $r;
656  } else {
657  $l = $i;
658  if ($this->e[$i] == 0.0) {
659  if ($w != 0.0) {
660  $this->H[$i][$n] = -$r / $w;
661  } else {
662  $this->H[$i][$n] = -$r / ($eps * $norm);
663  }
664  // Solve real equations
665  } else {
666  $x = $this->H[$i][$i+1];
667  $y = $this->H[$i+1][$i];
668  $q = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i];
669  $t = ($x * $s - $z * $r) / $q;
670  $this->H[$i][$n] = $t;
671  if (abs($x) > abs($z)) {
672  $this->H[$i+1][$n] = (-$r - $w * $t) / $x;
673  } else {
674  $this->H[$i+1][$n] = (-$s - $y * $t) / $z;
675  }
676  }
677  // Overflow control
678  $t = abs($this->H[$i][$n]);
679  if (($eps * $t) * $t > 1) {
680  for ($j = $i; $j <= $n; ++$j) {
681  $this->H[$j][$n] = $this->H[$j][$n] / $t;
682  }
683  }
684  }
685  }
686  // Complex vector
687  } else if ($q < 0) {
688  $l = $n-1;
689  // Last vector component imaginary so matrix is triangular
690  if (abs($this->H[$n][$n-1]) > abs($this->H[$n-1][$n])) {
691  $this->H[$n-1][$n-1] = $q / $this->H[$n][$n-1];
692  $this->H[$n-1][$n] = -($this->H[$n][$n] - $p) / $this->H[$n][$n-1];
693  } else {
694  $this->cdiv(0.0, -$this->H[$n-1][$n], $this->H[$n-1][$n-1] - $p, $q);
695  $this->H[$n-1][$n-1] = $this->cdivr;
696  $this->H[$n-1][$n] = $this->cdivi;
697  }
698  $this->H[$n][$n-1] = 0.0;
699  $this->H[$n][$n] = 1.0;
700  for ($i = $n-2; $i >= 0; --$i) {
701  // double ra,sa,vr,vi;
702  $ra = 0.0;
703  $sa = 0.0;
704  for ($j = $l; $j <= $n; ++$j) {
705  $ra = $ra + $this->H[$i][$j] * $this->H[$j][$n-1];
706  $sa = $sa + $this->H[$i][$j] * $this->H[$j][$n];
707  }
708  $w = $this->H[$i][$i] - $p;
709  if ($this->e[$i] < 0.0) {
710  $z = $w;
711  $r = $ra;
712  $s = $sa;
713  } else {
714  $l = $i;
715  if ($this->e[$i] == 0) {
716  $this->cdiv(-$ra, -$sa, $w, $q);
717  $this->H[$i][$n-1] = $this->cdivr;
718  $this->H[$i][$n] = $this->cdivi;
719  } else {
720  // Solve complex equations
721  $x = $this->H[$i][$i+1];
722  $y = $this->H[$i+1][$i];
723  $vr = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i] - $q * $q;
724  $vi = ($this->d[$i] - $p) * 2.0 * $q;
725  if ($vr == 0.0 & $vi == 0.0) {
726  $vr = $eps * $norm * (abs($w) + abs($q) + abs($x) + abs($y) + abs($z));
727  }
728  $this->cdiv($x * $r - $z * $ra + $q * $sa, $x * $s - $z * $sa - $q * $ra, $vr, $vi);
729  $this->H[$i][$n-1] = $this->cdivr;
730  $this->H[$i][$n] = $this->cdivi;
731  if (abs($x) > (abs($z) + abs($q))) {
732  $this->H[$i+1][$n-1] = (-$ra - $w * $this->H[$i][$n-1] + $q * $this->H[$i][$n]) / $x;
733  $this->H[$i+1][$n] = (-$sa - $w * $this->H[$i][$n] - $q * $this->H[$i][$n-1]) / $x;
734  } else {
735  $this->cdiv(-$r - $y * $this->H[$i][$n-1], -$s - $y * $this->H[$i][$n], $z, $q);
736  $this->H[$i+1][$n-1] = $this->cdivr;
737  $this->H[$i+1][$n] = $this->cdivi;
738  }
739  }
740  // Overflow control
741  $t = max(abs($this->H[$i][$n-1]),abs($this->H[$i][$n]));
742  if (($eps * $t) * $t > 1) {
743  for ($j = $i; $j <= $n; ++$j) {
744  $this->H[$j][$n-1] = $this->H[$j][$n-1] / $t;
745  $this->H[$j][$n] = $this->H[$j][$n] / $t;
746  }
747  }
748  } // end else
749  } // end for
750  } // end else for complex case
751  } // end for
752 
753  // Vectors of isolated roots
754  for ($i = 0; $i < $nn; ++$i) {
755  if ($i < $low | $i > $high) {
756  for ($j = $i; $j < $nn; ++$j) {
757  $this->V[$i][$j] = $this->H[$i][$j];
758  }
759  }
760  }
761 
762  // Back transformation to get eigenvectors of original matrix
763  for ($j = $nn-1; $j >= $low; --$j) {
764  for ($i = $low; $i <= $high; ++$i) {
765  $z = 0.0;
766  for ($k = $low; $k <= min($j,$high); ++$k) {
767  $z = $z + $this->V[$i][$k] * $this->H[$k][$j];
768  }
769  $this->V[$i][$j] = $z;
770  }
771  }
772  } // end hqr2
773 
774 
782  public function __construct($Arg) {
783  $this->A = $Arg->getArray();
784  $this->n = $Arg->getColumnDimension();
785 
786  $issymmetric = true;
787  for ($j = 0; ($j < $this->n) & $issymmetric; ++$j) {
788  for ($i = 0; ($i < $this->n) & $issymmetric; ++$i) {
789  $issymmetric = ($this->A[$i][$j] == $this->A[$j][$i]);
790  }
791  }
792 
793  if ($issymmetric) {
794  $this->V = $this->A;
795  // Tridiagonalize.
796  $this->tred2();
797  // Diagonalize.
798  $this->tql2();
799  } else {
800  $this->H = $this->A;
801  $this->ort = array();
802  // Reduce to Hessenberg form.
803  $this->orthes();
804  // Reduce Hessenberg to real Schur form.
805  $this->hqr2();
806  }
807  }
808 
809 
816  public function getV() {
817  return new Matrix($this->V, $this->n, $this->n);
818  }
819 
820 
827  public function getRealEigenvalues() {
828  return $this->d;
829  }
830 
831 
838  public function getImagEigenvalues() {
839  return $this->e;
840  }
841 
842 
849  public function getD() {
850  for ($i = 0; $i < $this->n; ++$i) {
851  $D[$i] = array_fill(0, $this->n, 0.0);
852  $D[$i][$i] = $this->d[$i];
853  if ($this->e[$i] == 0) {
854  continue;
855  }
856  $o = ($this->e[$i] > 0) ? $i + 1 : $i - 1;
857  $D[$i][$o] = $this->e[$i];
858  }
859  return new Matrix($D);
860  }
861 
862 } // class EigenvalueDecomposition