ILIAS  release_5-3 Revision v5.3.23-19-g915713cf615
Public Member Functions | Private Member Functions | Private Attributes
EigenvalueDecomposition Class Reference
+ Collaboration diagram for EigenvalueDecomposition:

Public Member Functions

 __construct ($Arg)
 Constructor: Check for symmetry, then construct the eigenvalue decomposition. More...
 
 getV ()
 Return the eigenvector matrix. More...
 
 getRealEigenvalues ()
 Return the real parts of the eigenvalues. More...
 
 getImagEigenvalues ()
 Return the imaginary parts of the eigenvalues. More...
 
 getD ()
 Return the block diagonal eigenvalue matrix. More...
 

Private Member Functions

 tred2 ()
 Symmetric Householder reduction to tridiagonal form. More...
 
 tql2 ()
 Symmetric tridiagonal QL algorithm. More...
 
 orthes ()
 Nonsymmetric reduction to Hessenberg form. More...
 
 cdiv ($xr, $xi, $yr, $yi)
 Performs complex division. More...
 
 hqr2 ()
 Nonsymmetric reduction from Hessenberg to real Schur form. More...
 

Private Attributes

 $n
 
 $issymmetric
 
 $d = array()
 
 $e = array()
 
 $V = array()
 
 $H = array()
 
 $ort
 
 $cdivr
 
 $cdivi
 

Detailed Description

Definition at line 24 of file EigenvalueDecomposition.php.

Constructor & Destructor Documentation

◆ __construct()

EigenvalueDecomposition::__construct (   $Arg)

Constructor: Check for symmetry, then construct the eigenvalue decomposition.

public

Parameters
ASquare matrix
Returns
Structure to access D and V.

Definition at line 782 of file EigenvalueDecomposition.php.

References $i, $issymmetric, $n, array, hqr2(), n, orthes(), tql2(), and tred2().

782  {
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  }
EigenvalueDecomposition\tred2
tred2()
Symmetric Householder reduction to tridiagonal form.
Definition: EigenvalueDecomposition.php:76
EigenvalueDecomposition\tql2
tql2()
Symmetric tridiagonal QL algorithm.
Definition: EigenvalueDecomposition.php:185
EigenvalueDecomposition\hqr2
hqr2()
Nonsymmetric reduction from Hessenberg to real Schur form.
Definition: EigenvalueDecomposition.php:398
EigenvalueDecomposition\orthes
orthes()
Nonsymmetric reduction to Hessenberg form.
Definition: EigenvalueDecomposition.php:291
n
if(! $in) print Initializing normalization quick check tables n
Definition: UtfNormalGenerate.php:37
array
Create styles array
The data for the language used.
Definition: 40duplicateStyle.php:19
$i
$i
Definition: disco.tpl.php:19
+ Here is the call graph for this function:

Member Function Documentation

◆ cdiv()

EigenvalueDecomposition::cdiv (   $xr,
  $xi,
  $yr,
  $yi 
)
private

Performs complex division.

private

Definition at line 373 of file EigenvalueDecomposition.php.

References $d, and $r.

Referenced by hqr2().

373  {
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  }
$r
$r
Definition: example_031.php:79
+ Here is the caller graph for this function:

◆ getD()

EigenvalueDecomposition::getD ( )

Return the block diagonal eigenvalue matrix.

public

Returns
D

Definition at line 849 of file EigenvalueDecomposition.php.

References $i, $n, and e().

849  {
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  }
e
e($cmd)
Definition: flush.php:14
$i
$i
Definition: disco.tpl.php:19
+ Here is the call graph for this function:

◆ getImagEigenvalues()

EigenvalueDecomposition::getImagEigenvalues ( )

Return the imaginary parts of the eigenvalues.

public

Returns
imag(diag(D))

Definition at line 838 of file EigenvalueDecomposition.php.

References $e.

838  {
839  return $this->e;
840  }

◆ getRealEigenvalues()

EigenvalueDecomposition::getRealEigenvalues ( )

Return the real parts of the eigenvalues.

public

Returns
real(diag(D))

Definition at line 827 of file EigenvalueDecomposition.php.

References $d.

827  {
828  return $this->d;
829  }

◆ getV()

EigenvalueDecomposition::getV ( )

Return the eigenvector matrix.

public

Returns
V

Definition at line 816 of file EigenvalueDecomposition.php.

References n.

816  {
817  return new Matrix($this->V, $this->n, $this->n);
818  }
n
if(! $in) print Initializing normalization quick check tables n
Definition: UtfNormalGenerate.php:37

◆ hqr2()

EigenvalueDecomposition::hqr2 ( )
private

Nonsymmetric reduction from Hessenberg to real Schur form.

Code is derived from the Algol procedure hqr2, by Martin and Wilkinson, Handbook for Auto. Comp., Vol.ii-Linear Algebra, and the corresponding Fortran subroutine in EISPACK.

private

Definition at line 398 of file EigenvalueDecomposition.php.

References $cdivi, $cdivr, $eps, $i, $l, $m, $n, $r, $s, $t, $w, $x, $y, cdiv(), e(), and n.

Referenced by __construct().

398  {
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
$x
$x
Definition: example_009.php:98
n
if(! $in) print Initializing normalization quick check tables n
Definition: UtfNormalGenerate.php:37
$s
$s
Definition: pwgen.php:45
$m
$m
Definition: show_metadata.php:25
$w
$w
Definition: example_009.php:100
$r
$r
Definition: example_031.php:79
$y
$y
Definition: example_007.php:83
e
e($cmd)
Definition: flush.php:14
$eps
$eps
Definition: metadata.php:61
$l
global $l
Definition: afr.php:30
$i
$i
Definition: disco.tpl.php:19
EigenvalueDecomposition\cdiv
cdiv($xr, $xi, $yr, $yi)
Performs complex division.
Definition: EigenvalueDecomposition.php:373
+ Here is the call graph for this function:
+ Here is the caller graph for this function:

◆ orthes()

EigenvalueDecomposition::orthes ( )
private

Nonsymmetric reduction to Hessenberg form.

This is derived from the Algol procedures orthes and ortran, by Martin and Wilkinson, Handbook for Auto. Comp., Vol.ii-Linear Algebra, and the corresponding Fortran subroutines in EISPACK.

private

Definition at line 291 of file EigenvalueDecomposition.php.

References $h, $i, $m, $n, and n.

Referenced by __construct().

291  {
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  }
$h
$h
Definition: example_009.php:101
n
if(! $in) print Initializing normalization quick check tables n
Definition: UtfNormalGenerate.php:37
$m
$m
Definition: show_metadata.php:25
$i
$i
Definition: disco.tpl.php:19
+ Here is the caller graph for this function:

◆ tql2()

EigenvalueDecomposition::tql2 ( )
private

Symmetric tridiagonal QL algorithm.

This is derived from the Algol procedures tql2, by Bowdler, Martin, Reinsch, and Wilkinson, Handbook for Auto. Comp., Vol.ii-Linear Algebra, and the corresponding Fortran subroutine in EISPACK.

private

Definition at line 185 of file EigenvalueDecomposition.php.

References $eps, $h, $i, $l, $m, $n, $r, $s, e(), and hypo().

Referenced by __construct().

185  {
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  }
$h
$h
Definition: example_009.php:101
$s
$s
Definition: pwgen.php:45
$m
$m
Definition: show_metadata.php:25
$r
$r
Definition: example_031.php:79
e
e($cmd)
Definition: flush.php:14
$eps
$eps
Definition: metadata.php:61
hypo
hypo($a, $b)
Definition: Maths.php:14
$l
global $l
Definition: afr.php:30
$i
$i
Definition: disco.tpl.php:19
+ Here is the call graph for this function:
+ Here is the caller graph for this function:

◆ tred2()

EigenvalueDecomposition::tred2 ( )
private

Symmetric Householder reduction to tridiagonal form.

private

Definition at line 76 of file EigenvalueDecomposition.php.

References $h, $i, e(), and n.

Referenced by __construct().

76  {
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  }
$h
$h
Definition: example_009.php:101
n
if(! $in) print Initializing normalization quick check tables n
Definition: UtfNormalGenerate.php:37
e
e($cmd)
Definition: flush.php:14
$i
$i
Definition: disco.tpl.php:19
+ Here is the call graph for this function:
+ Here is the caller graph for this function:

Field Documentation

◆ $cdivi

EigenvalueDecomposition::$cdivi
private

Definition at line 68 of file EigenvalueDecomposition.php.

Referenced by hqr2().

◆ $cdivr

EigenvalueDecomposition::$cdivr
private

Definition at line 67 of file EigenvalueDecomposition.php.

Referenced by hqr2().

◆ $d

EigenvalueDecomposition::$d = array()
private

Definition at line 42 of file EigenvalueDecomposition.php.

Referenced by cdiv(), and getRealEigenvalues().

◆ $e

EigenvalueDecomposition::$e = array()
private

Definition at line 43 of file EigenvalueDecomposition.php.

Referenced by getImagEigenvalues().

◆ $H

EigenvalueDecomposition::$H = array()
private

Definition at line 55 of file EigenvalueDecomposition.php.

◆ $issymmetric

EigenvalueDecomposition::$issymmetric
private

Definition at line 36 of file EigenvalueDecomposition.php.

Referenced by __construct().

◆ $n

EigenvalueDecomposition::$n
private

Definition at line 30 of file EigenvalueDecomposition.php.

Referenced by __construct(), getD(), hqr2(), orthes(), and tql2().

◆ $ort

EigenvalueDecomposition::$ort
private

Definition at line 61 of file EigenvalueDecomposition.php.

◆ $V

EigenvalueDecomposition::$V = array()
private

Definition at line 49 of file EigenvalueDecomposition.php.

The documentation for this class was generated from the following file: