ILIAS  Release_4_0_x_branch Revision 61816
 All Data Structures Namespaces Files Functions Variables Groups Pages
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.
 getV ()
 Return the eigenvector matrix.
 getRealEigenvalues ()
 Return the real parts of the eigenvalues.
 getImagEigenvalues ()
 Return the imaginary parts of the eigenvalues.
 getD ()
 Return the block diagonal eigenvalue matrix.

Private Member Functions

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

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

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 $issymmetric, $n, hqr2(), n, orthes(), tql2(), and tred2().

{
$this->A = $Arg->getArray();
$this->n = $Arg->getColumnDimension();
$issymmetric = true;
for ($j = 0; ($j < $this->n) & $issymmetric; ++$j) {
for ($i = 0; ($i < $this->n) & $issymmetric; ++$i) {
$issymmetric = ($this->A[$i][$j] == $this->A[$j][$i]);
}
}
if ($issymmetric) {
$this->V = $this->A;
// Tridiagonalize.
$this->tred2();
// Diagonalize.
$this->tql2();
} else {
$this->H = $this->A;
$this->ort = array();
// Reduce to Hessenberg form.
$this->orthes();
// Reduce Hessenberg to real Schur form.
$this->hqr2();
}
}

+ Here is the call graph for this function:

Member Function Documentation

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

Performs complex division.

private

Definition at line 373 of file EigenvalueDecomposition.php.

References $d.

Referenced by hqr2().

{
if (abs($yr) > abs($yi)) {
$r = $yi / $yr;
$d = $yr + $r * $yi;
$this->cdivr = ($xr + $r * $xi) / $d;
$this->cdivi = ($xi - $r * $xr) / $d;
} else {
$r = $yr / $yi;
$d = $yi + $r * $yr;
$this->cdivr = ($r * $xr + $xi) / $d;
$this->cdivi = ($r * $xi - $xr) / $d;
}
}

+ Here is the caller graph for this function:

EigenvalueDecomposition::getD ( )

Return the block diagonal eigenvalue matrix.

public

Returns
D

Definition at line 849 of file EigenvalueDecomposition.php.

References $n.

{
for ($i = 0; $i < $this->n; ++$i) {
$D[$i] = array_fill(0, $this->n, 0.0);
$D[$i][$i] = $this->d[$i];
if ($this->e[$i] == 0) {
continue;
}
$o = ($this->e[$i] > 0) ? $i + 1 : $i - 1;
$D[$i][$o] = $this->e[$i];
}
return new Matrix($D);
}
EigenvalueDecomposition::getImagEigenvalues ( )

Return the imaginary parts of the eigenvalues.

public

Returns
imag(diag(D))

Definition at line 838 of file EigenvalueDecomposition.php.

References $e.

{
return $this->e;
}
EigenvalueDecomposition::getRealEigenvalues ( )

Return the real parts of the eigenvalues.

public

Returns
real(diag(D))

Definition at line 827 of file EigenvalueDecomposition.php.

References $d.

{
return $this->d;
}
EigenvalueDecomposition::getV ( )

Return the eigenvector matrix.

public

Returns
V

Definition at line 816 of file EigenvalueDecomposition.php.

References n.

{
return new Matrix($this->V, $this->n, $this->n);
}
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, $l, $n, $q, $t, $x, $y, cdiv(), elseif(), and n.

Referenced by __construct().

{
// Initialize
$nn = $this->n;
$n = $nn - 1;
$low = 0;
$high = $nn - 1;
$eps = pow(2.0, -52.0);
$exshift = 0.0;
$p = $q = $r = $s = $z = 0;
// Store roots isolated by balanc and compute matrix norm
$norm = 0.0;
for ($i = 0; $i < $nn; ++$i) {
if (($i < $low) OR ($i > $high)) {
$this->d[$i] = $this->H[$i][$i];
$this->e[$i] = 0.0;
}
for ($j = max($i-1, 0); $j < $nn; ++$j) {
$norm = $norm + abs($this->H[$i][$j]);
}
}
// Outer loop over eigenvalue index
$iter = 0;
while ($n >= $low) {
// Look for single small sub-diagonal element
$l = $n;
while ($l > $low) {
$s = abs($this->H[$l-1][$l-1]) + abs($this->H[$l][$l]);
if ($s == 0.0) {
$s = $norm;
}
if (abs($this->H[$l][$l-1]) < $eps * $s) {
break;
}
--$l;
}
// Check for convergence
// One root found
if ($l == $n) {
$this->H[$n][$n] = $this->H[$n][$n] + $exshift;
$this->d[$n] = $this->H[$n][$n];
$this->e[$n] = 0.0;
--$n;
$iter = 0;
// Two roots found
} else if ($l == $n-1) {
$w = $this->H[$n][$n-1] * $this->H[$n-1][$n];
$p = ($this->H[$n-1][$n-1] - $this->H[$n][$n]) / 2.0;
$q = $p * $p + $w;
$z = sqrt(abs($q));
$this->H[$n][$n] = $this->H[$n][$n] + $exshift;
$this->H[$n-1][$n-1] = $this->H[$n-1][$n-1] + $exshift;
$x = $this->H[$n][$n];
// Real pair
if ($q >= 0) {
if ($p >= 0) {
$z = $p + $z;
} else {
$z = $p - $z;
}
$this->d[$n-1] = $x + $z;
$this->d[$n] = $this->d[$n-1];
if ($z != 0.0) {
$this->d[$n] = $x - $w / $z;
}
$this->e[$n-1] = 0.0;
$this->e[$n] = 0.0;
$x = $this->H[$n][$n-1];
$s = abs($x) + abs($z);
$p = $x / $s;
$q = $z / $s;
$r = sqrt($p * $p + $q * $q);
$p = $p / $r;
$q = $q / $r;
// Row modification
for ($j = $n-1; $j < $nn; ++$j) {
$z = $this->H[$n-1][$j];
$this->H[$n-1][$j] = $q * $z + $p * $this->H[$n][$j];
$this->H[$n][$j] = $q * $this->H[$n][$j] - $p * $z;
}
// Column modification
for ($i = 0; $i <= n; ++$i) {
$z = $this->H[$i][$n-1];
$this->H[$i][$n-1] = $q * $z + $p * $this->H[$i][$n];
$this->H[$i][$n] = $q * $this->H[$i][$n] - $p * $z;
}
// Accumulate transformations
for ($i = $low; $i <= $high; ++$i) {
$z = $this->V[$i][$n-1];
$this->V[$i][$n-1] = $q * $z + $p * $this->V[$i][$n];
$this->V[$i][$n] = $q * $this->V[$i][$n] - $p * $z;
}
// Complex pair
} else {
$this->d[$n-1] = $x + $p;
$this->d[$n] = $x + $p;
$this->e[$n-1] = $z;
$this->e[$n] = -$z;
}
$n = $n - 2;
$iter = 0;
// No convergence yet
} else {
// Form shift
$x = $this->H[$n][$n];
$y = 0.0;
$w = 0.0;
if ($l < $n) {
$y = $this->H[$n-1][$n-1];
$w = $this->H[$n][$n-1] * $this->H[$n-1][$n];
}
// Wilkinson's original ad hoc shift
if ($iter == 10) {
$exshift += $x;
for ($i = $low; $i <= $n; ++$i) {
$this->H[$i][$i] -= $x;
}
$s = abs($this->H[$n][$n-1]) + abs($this->H[$n-1][$n-2]);
$x = $y = 0.75 * $s;
$w = -0.4375 * $s * $s;
}
// MATLAB's new ad hoc shift
if ($iter == 30) {
$s = ($y - $x) / 2.0;
$s = $s * $s + $w;
if ($s > 0) {
$s = sqrt($s);
if ($y < $x) {
$s = -$s;
}
$s = $x - $w / (($y - $x) / 2.0 + $s);
for ($i = $low; $i <= $n; ++$i) {
$this->H[$i][$i] -= $s;
}
$exshift += $s;
$x = $y = $w = 0.964;
}
}
// Could check iteration count here.
$iter = $iter + 1;
// Look for two consecutive small sub-diagonal elements
$m = $n - 2;
while ($m >= $l) {
$z = $this->H[$m][$m];
$r = $x - $z;
$s = $y - $z;
$p = ($r * $s - $w) / $this->H[$m+1][$m] + $this->H[$m][$m+1];
$q = $this->H[$m+1][$m+1] - $z - $r - $s;
$r = $this->H[$m+2][$m+1];
$s = abs($p) + abs($q) + abs($r);
$p = $p / $s;
$q = $q / $s;
$r = $r / $s;
if ($m == $l) {
break;
}
if (abs($this->H[$m][$m-1]) * (abs($q) + abs($r)) <
$eps * (abs($p) * (abs($this->H[$m-1][$m-1]) + abs($z) + abs($this->H[$m+1][$m+1])))) {
break;
}
--$m;
}
for ($i = $m + 2; $i <= $n; ++$i) {
$this->H[$i][$i-2] = 0.0;
if ($i > $m+2) {
$this->H[$i][$i-3] = 0.0;
}
}
// Double QR step involving rows l:n and columns m:n
for ($k = $m; $k <= $n-1; ++$k) {
$notlast = ($k != $n-1);
if ($k != $m) {
$p = $this->H[$k][$k-1];
$q = $this->H[$k+1][$k-1];
$r = ($notlast ? $this->H[$k+2][$k-1] : 0.0);
$x = abs($p) + abs($q) + abs($r);
if ($x != 0.0) {
$p = $p / $x;
$q = $q / $x;
$r = $r / $x;
}
}
if ($x == 0.0) {
break;
}
$s = sqrt($p * $p + $q * $q + $r * $r);
if ($p < 0) {
$s = -$s;
}
if ($s != 0) {
if ($k != $m) {
$this->H[$k][$k-1] = -$s * $x;
} elseif ($l != $m) {
$this->H[$k][$k-1] = -$this->H[$k][$k-1];
}
$p = $p + $s;
$x = $p / $s;
$y = $q / $s;
$z = $r / $s;
$q = $q / $p;
$r = $r / $p;
// Row modification
for ($j = $k; $j < $nn; ++$j) {
$p = $this->H[$k][$j] + $q * $this->H[$k+1][$j];
if ($notlast) {
$p = $p + $r * $this->H[$k+2][$j];
$this->H[$k+2][$j] = $this->H[$k+2][$j] - $p * $z;
}
$this->H[$k][$j] = $this->H[$k][$j] - $p * $x;
$this->H[$k+1][$j] = $this->H[$k+1][$j] - $p * $y;
}
// Column modification
for ($i = 0; $i <= min($n, $k+3); ++$i) {
$p = $x * $this->H[$i][$k] + $y * $this->H[$i][$k+1];
if ($notlast) {
$p = $p + $z * $this->H[$i][$k+2];
$this->H[$i][$k+2] = $this->H[$i][$k+2] - $p * $r;
}
$this->H[$i][$k] = $this->H[$i][$k] - $p;
$this->H[$i][$k+1] = $this->H[$i][$k+1] - $p * $q;
}
// Accumulate transformations
for ($i = $low; $i <= $high; ++$i) {
$p = $x * $this->V[$i][$k] + $y * $this->V[$i][$k+1];
if ($notlast) {
$p = $p + $z * $this->V[$i][$k+2];
$this->V[$i][$k+2] = $this->V[$i][$k+2] - $p * $r;
}
$this->V[$i][$k] = $this->V[$i][$k] - $p;
$this->V[$i][$k+1] = $this->V[$i][$k+1] - $p * $q;
}
} // ($s != 0)
} // k loop
} // check convergence
} // while ($n >= $low)
// Backsubstitute to find vectors of upper triangular form
if ($norm == 0.0) {
return;
}
for ($n = $nn-1; $n >= 0; --$n) {
$p = $this->d[$n];
$q = $this->e[$n];
// Real vector
if ($q == 0) {
$l = $n;
$this->H[$n][$n] = 1.0;
for ($i = $n-1; $i >= 0; --$i) {
$w = $this->H[$i][$i] - $p;
$r = 0.0;
for ($j = $l; $j <= $n; ++$j) {
$r = $r + $this->H[$i][$j] * $this->H[$j][$n];
}
if ($this->e[$i] < 0.0) {
$z = $w;
$s = $r;
} else {
$l = $i;
if ($this->e[$i] == 0.0) {
if ($w != 0.0) {
$this->H[$i][$n] = -$r / $w;
} else {
$this->H[$i][$n] = -$r / ($eps * $norm);
}
// Solve real equations
} else {
$x = $this->H[$i][$i+1];
$y = $this->H[$i+1][$i];
$q = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i];
$t = ($x * $s - $z * $r) / $q;
$this->H[$i][$n] = $t;
if (abs($x) > abs($z)) {
$this->H[$i+1][$n] = (-$r - $w * $t) / $x;
} else {
$this->H[$i+1][$n] = (-$s - $y * $t) / $z;
}
}
// Overflow control
$t = abs($this->H[$i][$n]);
if (($eps * $t) * $t > 1) {
for ($j = $i; $j <= $n; ++$j) {
$this->H[$j][$n] = $this->H[$j][$n] / $t;
}
}
}
}
// Complex vector
} else if ($q < 0) {
$l = $n-1;
// Last vector component imaginary so matrix is triangular
if (abs($this->H[$n][$n-1]) > abs($this->H[$n-1][$n])) {
$this->H[$n-1][$n-1] = $q / $this->H[$n][$n-1];
$this->H[$n-1][$n] = -($this->H[$n][$n] - $p) / $this->H[$n][$n-1];
} else {
$this->cdiv(0.0, -$this->H[$n-1][$n], $this->H[$n-1][$n-1] - $p, $q);
$this->H[$n-1][$n-1] = $this->cdivr;
$this->H[$n-1][$n] = $this->cdivi;
}
$this->H[$n][$n-1] = 0.0;
$this->H[$n][$n] = 1.0;
for ($i = $n-2; $i >= 0; --$i) {
// double ra,sa,vr,vi;
$ra = 0.0;
$sa = 0.0;
for ($j = $l; $j <= $n; ++$j) {
$ra = $ra + $this->H[$i][$j] * $this->H[$j][$n-1];
$sa = $sa + $this->H[$i][$j] * $this->H[$j][$n];
}
$w = $this->H[$i][$i] - $p;
if ($this->e[$i] < 0.0) {
$z = $w;
$r = $ra;
$s = $sa;
} else {
$l = $i;
if ($this->e[$i] == 0) {
$this->cdiv(-$ra, -$sa, $w, $q);
$this->H[$i][$n-1] = $this->cdivr;
$this->H[$i][$n] = $this->cdivi;
} else {
// Solve complex equations
$x = $this->H[$i][$i+1];
$y = $this->H[$i+1][$i];
$vr = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i] - $q * $q;
$vi = ($this->d[$i] - $p) * 2.0 * $q;
if ($vr == 0.0 & $vi == 0.0) {
$vr = $eps * $norm * (abs($w) + abs($q) + abs($x) + abs($y) + abs($z));
}
$this->cdiv($x * $r - $z * $ra + $q * $sa, $x * $s - $z * $sa - $q * $ra, $vr, $vi);
$this->H[$i][$n-1] = $this->cdivr;
$this->H[$i][$n] = $this->cdivi;
if (abs($x) > (abs($z) + abs($q))) {
$this->H[$i+1][$n-1] = (-$ra - $w * $this->H[$i][$n-1] + $q * $this->H[$i][$n]) / $x;
$this->H[$i+1][$n] = (-$sa - $w * $this->H[$i][$n] - $q * $this->H[$i][$n-1]) / $x;
} else {
$this->cdiv(-$r - $y * $this->H[$i][$n-1], -$s - $y * $this->H[$i][$n], $z, $q);
$this->H[$i+1][$n-1] = $this->cdivr;
$this->H[$i+1][$n] = $this->cdivi;
}
}
// Overflow control
$t = max(abs($this->H[$i][$n-1]),abs($this->H[$i][$n]));
if (($eps * $t) * $t > 1) {
for ($j = $i; $j <= $n; ++$j) {
$this->H[$j][$n-1] = $this->H[$j][$n-1] / $t;
$this->H[$j][$n] = $this->H[$j][$n] / $t;
}
}
} // end else
} // end for
} // end else for complex case
} // end for
// Vectors of isolated roots
for ($i = 0; $i < $nn; ++$i) {
if ($i < $low | $i > $high) {
for ($j = $i; $j < $nn; ++$j) {
$this->V[$i][$j] = $this->H[$i][$j];
}
}
}
// Back transformation to get eigenvectors of original matrix
for ($j = $nn-1; $j >= $low; --$j) {
for ($i = $low; $i <= $high; ++$i) {
$z = 0.0;
for ($k = $low; $k <= min($j,$high); ++$k) {
$z = $z + $this->V[$i][$k] * $this->H[$k][$j];
}
$this->V[$i][$j] = $z;
}
}
} // end hqr2

+ Here is the call graph for this function:

+ Here is the caller graph for this function:

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 $f, $n, and n.

Referenced by __construct().

{
$low = 0;
$high = $this->n-1;
for ($m = $low+1; $m <= $high-1; ++$m) {
// Scale column.
$scale = 0.0;
for ($i = $m; $i <= $high; ++$i) {
$scale = $scale + abs($this->H[$i][$m-1]);
}
if ($scale != 0.0) {
// Compute Householder transformation.
$h = 0.0;
for ($i = $high; $i >= $m; --$i) {
$this->ort[$i] = $this->H[$i][$m-1] / $scale;
$h += $this->ort[$i] * $this->ort[$i];
}
$g = sqrt($h);
if ($this->ort[$m] > 0) {
$g *= -1;
}
$h -= $this->ort[$m] * $g;
$this->ort[$m] -= $g;
// Apply Householder similarity transformation
// H = (I -u * u' / h) * H * (I -u * u') / h)
for ($j = $m; $j < $this->n; ++$j) {
$f = 0.0;
for ($i = $high; $i >= $m; --$i) {
$f += $this->ort[$i] * $this->H[$i][$j];
}
$f /= $h;
for ($i = $m; $i <= $high; ++$i) {
$this->H[$i][$j] -= $f * $this->ort[$i];
}
}
for ($i = 0; $i <= $high; ++$i) {
$f = 0.0;
for ($j = $high; $j >= $m; --$j) {
$f += $this->ort[$j] * $this->H[$i][$j];
}
$f = $f / $h;
for ($j = $m; $j <= $high; ++$j) {
$this->H[$i][$j] -= $f * $this->ort[$j];
}
}
$this->ort[$m] = $scale * $this->ort[$m];
$this->H[$m][$m-1] = $scale * $g;
}
}
// Accumulate transformations (Algol's ortran).
for ($i = 0; $i < $this->n; ++$i) {
for ($j = 0; $j < $this->n; ++$j) {
$this->V[$i][$j] = ($i == $j ? 1.0 : 0.0);
}
}
for ($m = $high-1; $m >= $low+1; --$m) {
if ($this->H[$m][$m-1] != 0.0) {
for ($i = $m+1; $i <= $high; ++$i) {
$this->ort[$i] = $this->H[$i][$m-1];
}
for ($j = $m; $j <= $high; ++$j) {
$g = 0.0;
for ($i = $m; $i <= $high; ++$i) {
$g += $this->ort[$i] * $this->V[$i][$j];
}
// Double division avoids possible underflow
$g = ($g / $this->ort[$m]) / $this->H[$m][$m-1];
for ($i = $m; $i <= $high; ++$i) {
$this->V[$i][$j] += $g * $this->ort[$i];
}
}
}
}
}

+ Here is the caller graph for this function:

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 $f, $l, $n, and hypo().

Referenced by __construct().

{
for ($i = 1; $i < $this->n; ++$i) {
$this->e[$i-1] = $this->e[$i];
}
$this->e[$this->n-1] = 0.0;
$f = 0.0;
$tst1 = 0.0;
$eps = pow(2.0,-52.0);
for ($l = 0; $l < $this->n; ++$l) {
// Find small subdiagonal element
$tst1 = max($tst1, abs($this->d[$l]) + abs($this->e[$l]));
$m = $l;
while ($m < $this->n) {
if (abs($this->e[$m]) <= $eps * $tst1)
break;
++$m;
}
// If m == l, $this->d[l] is an eigenvalue,
// otherwise, iterate.
if ($m > $l) {
$iter = 0;
do {
// Could check iteration count here.
$iter += 1;
// Compute implicit shift
$g = $this->d[$l];
$p = ($this->d[$l+1] - $g) / (2.0 * $this->e[$l]);
$r = hypo($p, 1.0);
if ($p < 0)
$r *= -1;
$this->d[$l] = $this->e[$l] / ($p + $r);
$this->d[$l+1] = $this->e[$l] * ($p + $r);
$dl1 = $this->d[$l+1];
$h = $g - $this->d[$l];
for ($i = $l + 2; $i < $this->n; ++$i)
$this->d[$i] -= $h;
$f += $h;
// Implicit QL transformation.
$p = $this->d[$m];
$c = 1.0;
$c2 = $c3 = $c;
$el1 = $this->e[$l + 1];
$s = $s2 = 0.0;
for ($i = $m-1; $i >= $l; --$i) {
$c3 = $c2;
$c2 = $c;
$s2 = $s;
$g = $c * $this->e[$i];
$h = $c * $p;
$r = hypo($p, $this->e[$i]);
$this->e[$i+1] = $s * $r;
$s = $this->e[$i] / $r;
$c = $p / $r;
$p = $c * $this->d[$i] - $s * $g;
$this->d[$i+1] = $h + $s * ($c * $g + $s * $this->d[$i]);
// Accumulate transformation.
for ($k = 0; $k < $this->n; ++$k) {
$h = $this->V[$k][$i+1];
$this->V[$k][$i+1] = $s * $this->V[$k][$i] + $c * $h;
$this->V[$k][$i] = $c * $this->V[$k][$i] - $s * $h;
}
}
$p = -$s * $s2 * $c3 * $el1 * $this->e[$l] / $dl1;
$this->e[$l] = $s * $p;
$this->d[$l] = $c * $p;
// Check for convergence.
} while (abs($this->e[$l]) > $eps * $tst1);
}
$this->d[$l] = $this->d[$l] + $f;
$this->e[$l] = 0.0;
}
// Sort eigenvalues and corresponding vectors.
for ($i = 0; $i < $this->n - 1; ++$i) {
$k = $i;
$p = $this->d[$i];
for ($j = $i+1; $j < $this->n; ++$j) {
if ($this->d[$j] < $p) {
$k = $j;
$p = $this->d[$j];
}
}
if ($k != $i) {
$this->d[$k] = $this->d[$i];
$this->d[$i] = $p;
for ($j = 0; $j < $this->n; ++$j) {
$p = $this->V[$j][$i];
$this->V[$j][$i] = $this->V[$j][$k];
$this->V[$j][$k] = $p;
}
}
}
}

+ Here is the call graph for this function:

+ Here is the caller graph for this function:

EigenvalueDecomposition::tred2 ( )
private

Symmetric Householder reduction to tridiagonal form.

private

Definition at line 76 of file EigenvalueDecomposition.php.

References $f, and n.

Referenced by __construct().

{
// This is derived from the Algol procedures tred2 by
// Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
// Fortran subroutine in EISPACK.
$this->d = $this->V[$this->n-1];
// Householder reduction to tridiagonal form.
for ($i = $this->n-1; $i > 0; --$i) {
$i_ = $i -1;
// Scale to avoid under/overflow.
$h = $scale = 0.0;
$scale += array_sum(array_map(abs, $this->d));
if ($scale == 0.0) {
$this->e[$i] = $this->d[$i_];
$this->d = array_slice($this->V[$i_], 0, $i_);
for ($j = 0; $j < $i; ++$j) {
$this->V[$j][$i] = $this->V[$i][$j] = 0.0;
}
} else {
// Generate Householder vector.
for ($k = 0; $k < $i; ++$k) {
$this->d[$k] /= $scale;
$h += pow($this->d[$k], 2);
}
$f = $this->d[$i_];
$g = sqrt($h);
if ($f > 0) {
$g = -$g;
}
$this->e[$i] = $scale * $g;
$h = $h - $f * $g;
$this->d[$i_] = $f - $g;
for ($j = 0; $j < $i; ++$j) {
$this->e[$j] = 0.0;
}
// Apply similarity transformation to remaining columns.
for ($j = 0; $j < $i; ++$j) {
$f = $this->d[$j];
$this->V[$j][$i] = $f;
$g = $this->e[$j] + $this->V[$j][$j] * $f;
for ($k = $j+1; $k <= $i_; ++$k) {
$g += $this->V[$k][$j] * $this->d[$k];
$this->e[$k] += $this->V[$k][$j] * $f;
}
$this->e[$j] = $g;
}
$f = 0.0;
for ($j = 0; $j < $i; ++$j) {
$this->e[$j] /= $h;
$f += $this->e[$j] * $this->d[$j];
}
$hh = $f / (2 * $h);
for ($j=0; $j < $i; ++$j) {
$this->e[$j] -= $hh * $this->d[$j];
}
for ($j = 0; $j < $i; ++$j) {
$f = $this->d[$j];
$g = $this->e[$j];
for ($k = $j; $k <= $i_; ++$k) {
$this->V[$k][$j] -= ($f * $this->e[$k] + $g * $this->d[$k]);
}
$this->d[$j] = $this->V[$i-1][$j];
$this->V[$i][$j] = 0.0;
}
}
$this->d[$i] = $h;
}
// Accumulate transformations.
for ($i = 0; $i < $this->n-1; ++$i) {
$this->V[$this->n-1][$i] = $this->V[$i][$i];
$this->V[$i][$i] = 1.0;
$h = $this->d[$i+1];
if ($h != 0.0) {
for ($k = 0; $k <= $i; ++$k) {
$this->d[$k] = $this->V[$k][$i+1] / $h;
}
for ($j = 0; $j <= $i; ++$j) {
$g = 0.0;
for ($k = 0; $k <= $i; ++$k) {
$g += $this->V[$k][$i+1] * $this->V[$k][$j];
}
for ($k = 0; $k <= $i; ++$k) {
$this->V[$k][$j] -= $g * $this->d[$k];
}
}
}
for ($k = 0; $k <= $i; ++$k) {
$this->V[$k][$i+1] = 0.0;
}
}
$this->d = $this->V[$this->n-1];
$this->V[$this->n-1] = array_fill(0, $j, 0.0);
$this->V[$this->n-1][$this->n-1] = 1.0;
$this->e[0] = 0.0;
}

+ Here is the caller graph for this function:

Field Documentation

EigenvalueDecomposition::$cdivi
private

Definition at line 68 of file EigenvalueDecomposition.php.

Referenced by hqr2().

EigenvalueDecomposition::$cdivr
private

Definition at line 67 of file EigenvalueDecomposition.php.

Referenced by hqr2().

EigenvalueDecomposition::$d = array()
private

Definition at line 42 of file EigenvalueDecomposition.php.

Referenced by cdiv(), and getRealEigenvalues().

EigenvalueDecomposition::$e = array()
private

Definition at line 43 of file EigenvalueDecomposition.php.

Referenced by getImagEigenvalues().

EigenvalueDecomposition::$H = array()
private

Definition at line 55 of file EigenvalueDecomposition.php.

EigenvalueDecomposition::$issymmetric
private

Definition at line 36 of file EigenvalueDecomposition.php.

Referenced by __construct().

EigenvalueDecomposition::$n
private

Definition at line 30 of file EigenvalueDecomposition.php.

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

EigenvalueDecomposition::$ort
private

Definition at line 61 of file EigenvalueDecomposition.php.

EigenvalueDecomposition::$V = array()
private

Definition at line 49 of file EigenvalueDecomposition.php.

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