forked from SimplesIP/pabx-app
You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
864 lines
32 KiB
864 lines
32 KiB
<?php |
|
/** |
|
* @package JAMA |
|
* |
|
* Class to obtain eigenvalues and eigenvectors of a real matrix. |
|
* |
|
* If A is symmetric, then A = V*D*V' where the eigenvalue matrix D |
|
* is diagonal and the eigenvector matrix V is orthogonal (i.e. |
|
* A = V.times(D.times(V.transpose())) and V.times(V.transpose()) |
|
* equals the identity matrix). |
|
* |
|
* If A is not symmetric, then the eigenvalue matrix D is block diagonal |
|
* with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues, |
|
* lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda]. The |
|
* columns of V represent the eigenvectors in the sense that A*V = V*D, |
|
* i.e. A.times(V) equals V.times(D). The matrix V may be badly |
|
* conditioned, or even singular, so the validity of the equation |
|
* A = V*D*inverse(V) depends upon V.cond(). |
|
* |
|
* @author Paul Meagher |
|
* @license PHP v3.0 |
|
* @version 1.1 |
|
*/ |
|
class EigenvalueDecomposition |
|
{ |
|
/** |
|
* Row and column dimension (square matrix). |
|
* @var int |
|
*/ |
|
private $n; |
|
|
|
/** |
|
* Internal symmetry flag. |
|
* @var int |
|
*/ |
|
private $issymmetric; |
|
|
|
/** |
|
* Arrays for internal storage of eigenvalues. |
|
* @var array |
|
*/ |
|
private $d = array(); |
|
private $e = array(); |
|
|
|
/** |
|
* Array for internal storage of eigenvectors. |
|
* @var array |
|
*/ |
|
private $V = array(); |
|
|
|
/** |
|
* Array for internal storage of nonsymmetric Hessenberg form. |
|
* @var array |
|
*/ |
|
private $H = array(); |
|
|
|
/** |
|
* Working storage for nonsymmetric algorithm. |
|
* @var array |
|
*/ |
|
private $ort; |
|
|
|
/** |
|
* Used for complex scalar division. |
|
* @var float |
|
*/ |
|
private $cdivr; |
|
private $cdivi; |
|
|
|
/** |
|
* Symmetric Householder reduction to tridiagonal form. |
|
* |
|
* @access private |
|
*/ |
|
private function tred2() |
|
{ |
|
// 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; |
|
} |
|
|
|
/** |
|
* 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. |
|
* |
|
* @access private |
|
*/ |
|
private function tql2() |
|
{ |
|
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; |
|
} |
|
} |
|
} |
|
} |
|
|
|
/** |
|
* 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. |
|
* |
|
* @access private |
|
*/ |
|
private function orthes() |
|
{ |
|
$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]; |
|
} |
|
} |
|
} |
|
} |
|
} |
|
|
|
/** |
|
* Performs complex division. |
|
* |
|
* @access private |
|
*/ |
|
private function cdiv($xr, $xi, $yr, $yi) |
|
{ |
|
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; |
|
} |
|
} |
|
|
|
/** |
|
* 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. |
|
* |
|
* @access private |
|
*/ |
|
private function hqr2() |
|
{ |
|
// 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 |
|
} elseif ($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 |
|
} elseif ($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 |
|
|
|
/** |
|
* Constructor: Check for symmetry, then construct the eigenvalue decomposition |
|
* |
|
* @access public |
|
* @param A Square matrix |
|
* @return Structure to access D and V. |
|
*/ |
|
public function __construct($Arg) |
|
{ |
|
$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(); |
|
} |
|
} |
|
|
|
/** |
|
* Return the eigenvector matrix |
|
* |
|
* @access public |
|
* @return V |
|
*/ |
|
public function getV() |
|
{ |
|
return new Matrix($this->V, $this->n, $this->n); |
|
} |
|
|
|
/** |
|
* Return the real parts of the eigenvalues |
|
* |
|
* @access public |
|
* @return real(diag(D)) |
|
*/ |
|
public function getRealEigenvalues() |
|
{ |
|
return $this->d; |
|
} |
|
|
|
/** |
|
* Return the imaginary parts of the eigenvalues |
|
* |
|
* @access public |
|
* @return imag(diag(D)) |
|
*/ |
|
public function getImagEigenvalues() |
|
{ |
|
return $this->e; |
|
} |
|
|
|
/** |
|
* Return the block diagonal eigenvalue matrix |
|
* |
|
* @access public |
|
* @return D |
|
*/ |
|
public function getD() |
|
{ |
|
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); |
|
} |
|
}
|
|
|