Subversion Repositories eFlore/Applications.cel

Compare Revisions

Ignore whitespace Rev 2388 → Rev 2910

/branches/v2.25-scarificateur/jrest/lib/PHPExcel/Classes/PHPExcel/Shared/JAMA/EigenvalueDecomposition.php
New file
0,0 → 1,862
<?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
} 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
 
 
/**
* 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);
}
 
} // class EigenvalueDecomposition
/branches/v2.25-scarificateur/jrest/lib/PHPExcel/Classes/PHPExcel/Shared/JAMA/Matrix.php
New file
0,0 → 1,1059
<?php
/**
* @package JAMA
*/
 
/** PHPExcel root directory */
if (!defined('PHPEXCEL_ROOT')) {
/**
* @ignore
*/
define('PHPEXCEL_ROOT', dirname(__FILE__) . '/../../../');
require(PHPEXCEL_ROOT . 'PHPExcel/Autoloader.php');
}
 
 
/*
* Matrix class
*
* @author Paul Meagher
* @author Michael Bommarito
* @author Lukasz Karapuda
* @author Bartek Matosiuk
* @version 1.8
* @license PHP v3.0
* @see http://math.nist.gov/javanumerics/jama/
*/
class PHPExcel_Shared_JAMA_Matrix {
 
 
const PolymorphicArgumentException = "Invalid argument pattern for polymorphic function.";
const ArgumentTypeException = "Invalid argument type.";
const ArgumentBoundsException = "Invalid argument range.";
const MatrixDimensionException = "Matrix dimensions are not equal.";
const ArrayLengthException = "Array length must be a multiple of m.";
 
/**
* Matrix storage
*
* @var array
* @access public
*/
public $A = array();
 
/**
* Matrix row dimension
*
* @var int
* @access private
*/
private $m;
 
/**
* Matrix column dimension
*
* @var int
* @access private
*/
private $n;
 
 
/**
* Polymorphic constructor
*
* As PHP has no support for polymorphic constructors, we hack our own sort of polymorphism using func_num_args, func_get_arg, and gettype. In essence, we're just implementing a simple RTTI filter and calling the appropriate constructor.
*/
public function __construct() {
if (func_num_args() > 0) {
$args = func_get_args();
$match = implode(",", array_map('gettype', $args));
 
switch($match) {
//Rectangular matrix - m x n initialized from 2D array
case 'array':
$this->m = count($args[0]);
$this->n = count($args[0][0]);
$this->A = $args[0];
break;
//Square matrix - n x n
case 'integer':
$this->m = $args[0];
$this->n = $args[0];
$this->A = array_fill(0, $this->m, array_fill(0, $this->n, 0));
break;
//Rectangular matrix - m x n
case 'integer,integer':
$this->m = $args[0];
$this->n = $args[1];
$this->A = array_fill(0, $this->m, array_fill(0, $this->n, 0));
break;
//Rectangular matrix - m x n initialized from packed array
case 'array,integer':
$this->m = $args[1];
if ($this->m != 0) {
$this->n = count($args[0]) / $this->m;
} else {
$this->n = 0;
}
if (($this->m * $this->n) == count($args[0])) {
for($i = 0; $i < $this->m; ++$i) {
for($j = 0; $j < $this->n; ++$j) {
$this->A[$i][$j] = $args[0][$i + $j * $this->m];
}
}
} else {
throw new PHPExcel_Calculation_Exception(self::ArrayLengthException);
}
break;
default:
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException);
break;
}
} else {
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException);
}
} // function __construct()
 
 
/**
* getArray
*
* @return array Matrix array
*/
public function getArray() {
return $this->A;
} // function getArray()
 
 
/**
* getRowDimension
*
* @return int Row dimension
*/
public function getRowDimension() {
return $this->m;
} // function getRowDimension()
 
 
/**
* getColumnDimension
*
* @return int Column dimension
*/
public function getColumnDimension() {
return $this->n;
} // function getColumnDimension()
 
 
/**
* get
*
* Get the i,j-th element of the matrix.
* @param int $i Row position
* @param int $j Column position
* @return mixed Element (int/float/double)
*/
public function get($i = null, $j = null) {
return $this->A[$i][$j];
} // function get()
 
 
/**
* getMatrix
*
* Get a submatrix
* @param int $i0 Initial row index
* @param int $iF Final row index
* @param int $j0 Initial column index
* @param int $jF Final column index
* @return Matrix Submatrix
*/
public function getMatrix() {
if (func_num_args() > 0) {
$args = func_get_args();
$match = implode(",", array_map('gettype', $args));
 
switch($match) {
//A($i0...; $j0...)
case 'integer,integer':
list($i0, $j0) = $args;
if ($i0 >= 0) { $m = $this->m - $i0; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentBoundsException); }
if ($j0 >= 0) { $n = $this->n - $j0; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentBoundsException); }
$R = new PHPExcel_Shared_JAMA_Matrix($m, $n);
for($i = $i0; $i < $this->m; ++$i) {
for($j = $j0; $j < $this->n; ++$j) {
$R->set($i, $j, $this->A[$i][$j]);
}
}
return $R;
break;
//A($i0...$iF; $j0...$jF)
case 'integer,integer,integer,integer':
list($i0, $iF, $j0, $jF) = $args;
if (($iF > $i0) && ($this->m >= $iF) && ($i0 >= 0)) { $m = $iF - $i0; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentBoundsException); }
if (($jF > $j0) && ($this->n >= $jF) && ($j0 >= 0)) { $n = $jF - $j0; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentBoundsException); }
$R = new PHPExcel_Shared_JAMA_Matrix($m+1, $n+1);
for($i = $i0; $i <= $iF; ++$i) {
for($j = $j0; $j <= $jF; ++$j) {
$R->set($i - $i0, $j - $j0, $this->A[$i][$j]);
}
}
return $R;
break;
//$R = array of row indices; $C = array of column indices
case 'array,array':
list($RL, $CL) = $args;
if (count($RL) > 0) { $m = count($RL); } else { throw new PHPExcel_Calculation_Exception(self::ArgumentBoundsException); }
if (count($CL) > 0) { $n = count($CL); } else { throw new PHPExcel_Calculation_Exception(self::ArgumentBoundsException); }
$R = new PHPExcel_Shared_JAMA_Matrix($m, $n);
for($i = 0; $i < $m; ++$i) {
for($j = 0; $j < $n; ++$j) {
$R->set($i - $i0, $j - $j0, $this->A[$RL[$i]][$CL[$j]]);
}
}
return $R;
break;
//$RL = array of row indices; $CL = array of column indices
case 'array,array':
list($RL, $CL) = $args;
if (count($RL) > 0) { $m = count($RL); } else { throw new PHPExcel_Calculation_Exception(self::ArgumentBoundsException); }
if (count($CL) > 0) { $n = count($CL); } else { throw new PHPExcel_Calculation_Exception(self::ArgumentBoundsException); }
$R = new PHPExcel_Shared_JAMA_Matrix($m, $n);
for($i = 0; $i < $m; ++$i) {
for($j = 0; $j < $n; ++$j) {
$R->set($i, $j, $this->A[$RL[$i]][$CL[$j]]);
}
}
return $R;
break;
//A($i0...$iF); $CL = array of column indices
case 'integer,integer,array':
list($i0, $iF, $CL) = $args;
if (($iF > $i0) && ($this->m >= $iF) && ($i0 >= 0)) { $m = $iF - $i0; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentBoundsException); }
if (count($CL) > 0) { $n = count($CL); } else { throw new PHPExcel_Calculation_Exception(self::ArgumentBoundsException); }
$R = new PHPExcel_Shared_JAMA_Matrix($m, $n);
for($i = $i0; $i < $iF; ++$i) {
for($j = 0; $j < $n; ++$j) {
$R->set($i - $i0, $j, $this->A[$RL[$i]][$j]);
}
}
return $R;
break;
//$RL = array of row indices
case 'array,integer,integer':
list($RL, $j0, $jF) = $args;
if (count($RL) > 0) { $m = count($RL); } else { throw new PHPExcel_Calculation_Exception(self::ArgumentBoundsException); }
if (($jF >= $j0) && ($this->n >= $jF) && ($j0 >= 0)) { $n = $jF - $j0; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentBoundsException); }
$R = new PHPExcel_Shared_JAMA_Matrix($m, $n+1);
for($i = 0; $i < $m; ++$i) {
for($j = $j0; $j <= $jF; ++$j) {
$R->set($i, $j - $j0, $this->A[$RL[$i]][$j]);
}
}
return $R;
break;
default:
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException);
break;
}
} else {
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException);
}
} // function getMatrix()
 
 
/**
* checkMatrixDimensions
*
* Is matrix B the same size?
* @param Matrix $B Matrix B
* @return boolean
*/
public function checkMatrixDimensions($B = null) {
if ($B instanceof PHPExcel_Shared_JAMA_Matrix) {
if (($this->m == $B->getRowDimension()) && ($this->n == $B->getColumnDimension())) {
return true;
} else {
throw new PHPExcel_Calculation_Exception(self::MatrixDimensionException);
}
} else {
throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException);
}
} // function checkMatrixDimensions()
 
 
 
/**
* set
*
* Set the i,j-th element of the matrix.
* @param int $i Row position
* @param int $j Column position
* @param mixed $c Int/float/double value
* @return mixed Element (int/float/double)
*/
public function set($i = null, $j = null, $c = null) {
// Optimized set version just has this
$this->A[$i][$j] = $c;
} // function set()
 
 
/**
* identity
*
* Generate an identity matrix.
* @param int $m Row dimension
* @param int $n Column dimension
* @return Matrix Identity matrix
*/
public function identity($m = null, $n = null) {
return $this->diagonal($m, $n, 1);
} // function identity()
 
 
/**
* diagonal
*
* Generate a diagonal matrix
* @param int $m Row dimension
* @param int $n Column dimension
* @param mixed $c Diagonal value
* @return Matrix Diagonal matrix
*/
public function diagonal($m = null, $n = null, $c = 1) {
$R = new PHPExcel_Shared_JAMA_Matrix($m, $n);
for($i = 0; $i < $m; ++$i) {
$R->set($i, $i, $c);
}
return $R;
} // function diagonal()
 
 
/**
* getMatrixByRow
*
* Get a submatrix by row index/range
* @param int $i0 Initial row index
* @param int $iF Final row index
* @return Matrix Submatrix
*/
public function getMatrixByRow($i0 = null, $iF = null) {
if (is_int($i0)) {
if (is_int($iF)) {
return $this->getMatrix($i0, 0, $iF + 1, $this->n);
} else {
return $this->getMatrix($i0, 0, $i0 + 1, $this->n);
}
} else {
throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException);
}
} // function getMatrixByRow()
 
 
/**
* getMatrixByCol
*
* Get a submatrix by column index/range
* @param int $i0 Initial column index
* @param int $iF Final column index
* @return Matrix Submatrix
*/
public function getMatrixByCol($j0 = null, $jF = null) {
if (is_int($j0)) {
if (is_int($jF)) {
return $this->getMatrix(0, $j0, $this->m, $jF + 1);
} else {
return $this->getMatrix(0, $j0, $this->m, $j0 + 1);
}
} else {
throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException);
}
} // function getMatrixByCol()
 
 
/**
* transpose
*
* Tranpose matrix
* @return Matrix Transposed matrix
*/
public function transpose() {
$R = new PHPExcel_Shared_JAMA_Matrix($this->n, $this->m);
for($i = 0; $i < $this->m; ++$i) {
for($j = 0; $j < $this->n; ++$j) {
$R->set($j, $i, $this->A[$i][$j]);
}
}
return $R;
} // function transpose()
 
 
/**
* trace
*
* Sum of diagonal elements
* @return float Sum of diagonal elements
*/
public function trace() {
$s = 0;
$n = min($this->m, $this->n);
for($i = 0; $i < $n; ++$i) {
$s += $this->A[$i][$i];
}
return $s;
} // function trace()
 
 
/**
* uminus
*
* Unary minus matrix -A
* @return Matrix Unary minus matrix
*/
public function uminus() {
} // function uminus()
 
 
/**
* plus
*
* A + B
* @param mixed $B Matrix/Array
* @return Matrix Sum
*/
public function plus() {
if (func_num_args() > 0) {
$args = func_get_args();
$match = implode(",", array_map('gettype', $args));
 
switch($match) {
case 'object':
if ($args[0] instanceof PHPExcel_Shared_JAMA_Matrix) { $M = $args[0]; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException); }
break;
case 'array':
$M = new PHPExcel_Shared_JAMA_Matrix($args[0]);
break;
default:
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException);
break;
}
$this->checkMatrixDimensions($M);
for($i = 0; $i < $this->m; ++$i) {
for($j = 0; $j < $this->n; ++$j) {
$M->set($i, $j, $M->get($i, $j) + $this->A[$i][$j]);
}
}
return $M;
} else {
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException);
}
} // function plus()
 
 
/**
* plusEquals
*
* A = A + B
* @param mixed $B Matrix/Array
* @return Matrix Sum
*/
public function plusEquals() {
if (func_num_args() > 0) {
$args = func_get_args();
$match = implode(",", array_map('gettype', $args));
 
switch($match) {
case 'object':
if ($args[0] instanceof PHPExcel_Shared_JAMA_Matrix) { $M = $args[0]; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException); }
break;
case 'array':
$M = new PHPExcel_Shared_JAMA_Matrix($args[0]);
break;
default:
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException);
break;
}
$this->checkMatrixDimensions($M);
for($i = 0; $i < $this->m; ++$i) {
for($j = 0; $j < $this->n; ++$j) {
$validValues = True;
$value = $M->get($i, $j);
if ((is_string($this->A[$i][$j])) && (strlen($this->A[$i][$j]) > 0) && (!is_numeric($this->A[$i][$j]))) {
$this->A[$i][$j] = trim($this->A[$i][$j],'"');
$validValues &= PHPExcel_Shared_String::convertToNumberIfFraction($this->A[$i][$j]);
}
if ((is_string($value)) && (strlen($value) > 0) && (!is_numeric($value))) {
$value = trim($value,'"');
$validValues &= PHPExcel_Shared_String::convertToNumberIfFraction($value);
}
if ($validValues) {
$this->A[$i][$j] += $value;
} else {
$this->A[$i][$j] = PHPExcel_Calculation_Functions::NaN();
}
}
}
return $this;
} else {
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException);
}
} // function plusEquals()
 
 
/**
* minus
*
* A - B
* @param mixed $B Matrix/Array
* @return Matrix Sum
*/
public function minus() {
if (func_num_args() > 0) {
$args = func_get_args();
$match = implode(",", array_map('gettype', $args));
 
switch($match) {
case 'object':
if ($args[0] instanceof PHPExcel_Shared_JAMA_Matrix) { $M = $args[0]; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException); }
break;
case 'array':
$M = new PHPExcel_Shared_JAMA_Matrix($args[0]);
break;
default:
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException);
break;
}
$this->checkMatrixDimensions($M);
for($i = 0; $i < $this->m; ++$i) {
for($j = 0; $j < $this->n; ++$j) {
$M->set($i, $j, $M->get($i, $j) - $this->A[$i][$j]);
}
}
return $M;
} else {
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException);
}
} // function minus()
 
 
/**
* minusEquals
*
* A = A - B
* @param mixed $B Matrix/Array
* @return Matrix Sum
*/
public function minusEquals() {
if (func_num_args() > 0) {
$args = func_get_args();
$match = implode(",", array_map('gettype', $args));
 
switch($match) {
case 'object':
if ($args[0] instanceof PHPExcel_Shared_JAMA_Matrix) { $M = $args[0]; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException); }
break;
case 'array':
$M = new PHPExcel_Shared_JAMA_Matrix($args[0]);
break;
default:
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException);
break;
}
$this->checkMatrixDimensions($M);
for($i = 0; $i < $this->m; ++$i) {
for($j = 0; $j < $this->n; ++$j) {
$validValues = True;
$value = $M->get($i, $j);
if ((is_string($this->A[$i][$j])) && (strlen($this->A[$i][$j]) > 0) && (!is_numeric($this->A[$i][$j]))) {
$this->A[$i][$j] = trim($this->A[$i][$j],'"');
$validValues &= PHPExcel_Shared_String::convertToNumberIfFraction($this->A[$i][$j]);
}
if ((is_string($value)) && (strlen($value) > 0) && (!is_numeric($value))) {
$value = trim($value,'"');
$validValues &= PHPExcel_Shared_String::convertToNumberIfFraction($value);
}
if ($validValues) {
$this->A[$i][$j] -= $value;
} else {
$this->A[$i][$j] = PHPExcel_Calculation_Functions::NaN();
}
}
}
return $this;
} else {
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException);
}
} // function minusEquals()
 
 
/**
* arrayTimes
*
* Element-by-element multiplication
* Cij = Aij * Bij
* @param mixed $B Matrix/Array
* @return Matrix Matrix Cij
*/
public function arrayTimes() {
if (func_num_args() > 0) {
$args = func_get_args();
$match = implode(",", array_map('gettype', $args));
 
switch($match) {
case 'object':
if ($args[0] instanceof PHPExcel_Shared_JAMA_Matrix) { $M = $args[0]; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException); }
break;
case 'array':
$M = new PHPExcel_Shared_JAMA_Matrix($args[0]);
break;
default:
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException);
break;
}
$this->checkMatrixDimensions($M);
for($i = 0; $i < $this->m; ++$i) {
for($j = 0; $j < $this->n; ++$j) {
$M->set($i, $j, $M->get($i, $j) * $this->A[$i][$j]);
}
}
return $M;
} else {
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException);
}
} // function arrayTimes()
 
 
/**
* arrayTimesEquals
*
* Element-by-element multiplication
* Aij = Aij * Bij
* @param mixed $B Matrix/Array
* @return Matrix Matrix Aij
*/
public function arrayTimesEquals() {
if (func_num_args() > 0) {
$args = func_get_args();
$match = implode(",", array_map('gettype', $args));
 
switch($match) {
case 'object':
if ($args[0] instanceof PHPExcel_Shared_JAMA_Matrix) { $M = $args[0]; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException); }
break;
case 'array':
$M = new PHPExcel_Shared_JAMA_Matrix($args[0]);
break;
default:
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException);
break;
}
$this->checkMatrixDimensions($M);
for($i = 0; $i < $this->m; ++$i) {
for($j = 0; $j < $this->n; ++$j) {
$validValues = True;
$value = $M->get($i, $j);
if ((is_string($this->A[$i][$j])) && (strlen($this->A[$i][$j]) > 0) && (!is_numeric($this->A[$i][$j]))) {
$this->A[$i][$j] = trim($this->A[$i][$j],'"');
$validValues &= PHPExcel_Shared_String::convertToNumberIfFraction($this->A[$i][$j]);
}
if ((is_string($value)) && (strlen($value) > 0) && (!is_numeric($value))) {
$value = trim($value,'"');
$validValues &= PHPExcel_Shared_String::convertToNumberIfFraction($value);
}
if ($validValues) {
$this->A[$i][$j] *= $value;
} else {
$this->A[$i][$j] = PHPExcel_Calculation_Functions::NaN();
}
}
}
return $this;
} else {
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException);
}
} // function arrayTimesEquals()
 
 
/**
* arrayRightDivide
*
* Element-by-element right division
* A / B
* @param Matrix $B Matrix B
* @return Matrix Division result
*/
public function arrayRightDivide() {
if (func_num_args() > 0) {
$args = func_get_args();
$match = implode(",", array_map('gettype', $args));
 
switch($match) {
case 'object':
if ($args[0] instanceof PHPExcel_Shared_JAMA_Matrix) { $M = $args[0]; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException); }
break;
case 'array':
$M = new PHPExcel_Shared_JAMA_Matrix($args[0]);
break;
default:
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException);
break;
}
$this->checkMatrixDimensions($M);
for($i = 0; $i < $this->m; ++$i) {
for($j = 0; $j < $this->n; ++$j) {
$validValues = True;
$value = $M->get($i, $j);
if ((is_string($this->A[$i][$j])) && (strlen($this->A[$i][$j]) > 0) && (!is_numeric($this->A[$i][$j]))) {
$this->A[$i][$j] = trim($this->A[$i][$j],'"');
$validValues &= PHPExcel_Shared_String::convertToNumberIfFraction($this->A[$i][$j]);
}
if ((is_string($value)) && (strlen($value) > 0) && (!is_numeric($value))) {
$value = trim($value,'"');
$validValues &= PHPExcel_Shared_String::convertToNumberIfFraction($value);
}
if ($validValues) {
if ($value == 0) {
// Trap for Divide by Zero error
$M->set($i, $j, '#DIV/0!');
} else {
$M->set($i, $j, $this->A[$i][$j] / $value);
}
} else {
$M->set($i, $j, PHPExcel_Calculation_Functions::NaN());
}
}
}
return $M;
} else {
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException);
}
} // function arrayRightDivide()
 
 
/**
* arrayRightDivideEquals
*
* Element-by-element right division
* Aij = Aij / Bij
* @param mixed $B Matrix/Array
* @return Matrix Matrix Aij
*/
public function arrayRightDivideEquals() {
if (func_num_args() > 0) {
$args = func_get_args();
$match = implode(",", array_map('gettype', $args));
 
switch($match) {
case 'object':
if ($args[0] instanceof PHPExcel_Shared_JAMA_Matrix) { $M = $args[0]; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException); }
break;
case 'array':
$M = new PHPExcel_Shared_JAMA_Matrix($args[0]);
break;
default:
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException);
break;
}
$this->checkMatrixDimensions($M);
for($i = 0; $i < $this->m; ++$i) {
for($j = 0; $j < $this->n; ++$j) {
$this->A[$i][$j] = $this->A[$i][$j] / $M->get($i, $j);
}
}
return $M;
} else {
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException);
}
} // function arrayRightDivideEquals()
 
 
/**
* arrayLeftDivide
*
* Element-by-element Left division
* A / B
* @param Matrix $B Matrix B
* @return Matrix Division result
*/
public function arrayLeftDivide() {
if (func_num_args() > 0) {
$args = func_get_args();
$match = implode(",", array_map('gettype', $args));
 
switch($match) {
case 'object':
if ($args[0] instanceof PHPExcel_Shared_JAMA_Matrix) { $M = $args[0]; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException); }
break;
case 'array':
$M = new PHPExcel_Shared_JAMA_Matrix($args[0]);
break;
default:
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException);
break;
}
$this->checkMatrixDimensions($M);
for($i = 0; $i < $this->m; ++$i) {
for($j = 0; $j < $this->n; ++$j) {
$M->set($i, $j, $M->get($i, $j) / $this->A[$i][$j]);
}
}
return $M;
} else {
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException);
}
} // function arrayLeftDivide()
 
 
/**
* arrayLeftDivideEquals
*
* Element-by-element Left division
* Aij = Aij / Bij
* @param mixed $B Matrix/Array
* @return Matrix Matrix Aij
*/
public function arrayLeftDivideEquals() {
if (func_num_args() > 0) {
$args = func_get_args();
$match = implode(",", array_map('gettype', $args));
 
switch($match) {
case 'object':
if ($args[0] instanceof PHPExcel_Shared_JAMA_Matrix) { $M = $args[0]; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException); }
break;
case 'array':
$M = new PHPExcel_Shared_JAMA_Matrix($args[0]);
break;
default:
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException);
break;
}
$this->checkMatrixDimensions($M);
for($i = 0; $i < $this->m; ++$i) {
for($j = 0; $j < $this->n; ++$j) {
$this->A[$i][$j] = $M->get($i, $j) / $this->A[$i][$j];
}
}
return $M;
} else {
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException);
}
} // function arrayLeftDivideEquals()
 
 
/**
* times
*
* Matrix multiplication
* @param mixed $n Matrix/Array/Scalar
* @return Matrix Product
*/
public function times() {
if (func_num_args() > 0) {
$args = func_get_args();
$match = implode(",", array_map('gettype', $args));
 
switch($match) {
case 'object':
if ($args[0] instanceof PHPExcel_Shared_JAMA_Matrix) { $B = $args[0]; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException); }
if ($this->n == $B->m) {
$C = new PHPExcel_Shared_JAMA_Matrix($this->m, $B->n);
for($j = 0; $j < $B->n; ++$j) {
for ($k = 0; $k < $this->n; ++$k) {
$Bcolj[$k] = $B->A[$k][$j];
}
for($i = 0; $i < $this->m; ++$i) {
$Arowi = $this->A[$i];
$s = 0;
for($k = 0; $k < $this->n; ++$k) {
$s += $Arowi[$k] * $Bcolj[$k];
}
$C->A[$i][$j] = $s;
}
}
return $C;
} else {
throw new PHPExcel_Calculation_Exception(JAMAError(MatrixDimensionMismatch));
}
break;
case 'array':
$B = new PHPExcel_Shared_JAMA_Matrix($args[0]);
if ($this->n == $B->m) {
$C = new PHPExcel_Shared_JAMA_Matrix($this->m, $B->n);
for($i = 0; $i < $C->m; ++$i) {
for($j = 0; $j < $C->n; ++$j) {
$s = "0";
for($k = 0; $k < $C->n; ++$k) {
$s += $this->A[$i][$k] * $B->A[$k][$j];
}
$C->A[$i][$j] = $s;
}
}
return $C;
} else {
throw new PHPExcel_Calculation_Exception(JAMAError(MatrixDimensionMismatch));
}
return $M;
break;
case 'integer':
$C = new PHPExcel_Shared_JAMA_Matrix($this->A);
for($i = 0; $i < $C->m; ++$i) {
for($j = 0; $j < $C->n; ++$j) {
$C->A[$i][$j] *= $args[0];
}
}
return $C;
break;
case 'double':
$C = new PHPExcel_Shared_JAMA_Matrix($this->m, $this->n);
for($i = 0; $i < $C->m; ++$i) {
for($j = 0; $j < $C->n; ++$j) {
$C->A[$i][$j] = $args[0] * $this->A[$i][$j];
}
}
return $C;
break;
case 'float':
$C = new PHPExcel_Shared_JAMA_Matrix($this->A);
for($i = 0; $i < $C->m; ++$i) {
for($j = 0; $j < $C->n; ++$j) {
$C->A[$i][$j] *= $args[0];
}
}
return $C;
break;
default:
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException);
break;
}
} else {
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException);
}
} // function times()
 
 
/**
* power
*
* A = A ^ B
* @param mixed $B Matrix/Array
* @return Matrix Sum
*/
public function power() {
if (func_num_args() > 0) {
$args = func_get_args();
$match = implode(",", array_map('gettype', $args));
 
switch($match) {
case 'object':
if ($args[0] instanceof PHPExcel_Shared_JAMA_Matrix) { $M = $args[0]; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException); }
break;
case 'array':
$M = new PHPExcel_Shared_JAMA_Matrix($args[0]);
break;
default:
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException);
break;
}
$this->checkMatrixDimensions($M);
for($i = 0; $i < $this->m; ++$i) {
for($j = 0; $j < $this->n; ++$j) {
$validValues = True;
$value = $M->get($i, $j);
if ((is_string($this->A[$i][$j])) && (strlen($this->A[$i][$j]) > 0) && (!is_numeric($this->A[$i][$j]))) {
$this->A[$i][$j] = trim($this->A[$i][$j],'"');
$validValues &= PHPExcel_Shared_String::convertToNumberIfFraction($this->A[$i][$j]);
}
if ((is_string($value)) && (strlen($value) > 0) && (!is_numeric($value))) {
$value = trim($value,'"');
$validValues &= PHPExcel_Shared_String::convertToNumberIfFraction($value);
}
if ($validValues) {
$this->A[$i][$j] = pow($this->A[$i][$j],$value);
} else {
$this->A[$i][$j] = PHPExcel_Calculation_Functions::NaN();
}
}
}
return $this;
} else {
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException);
}
} // function power()
 
 
/**
* concat
*
* A = A & B
* @param mixed $B Matrix/Array
* @return Matrix Sum
*/
public function concat() {
if (func_num_args() > 0) {
$args = func_get_args();
$match = implode(",", array_map('gettype', $args));
 
switch($match) {
case 'object':
if ($args[0] instanceof PHPExcel_Shared_JAMA_Matrix) { $M = $args[0]; } else { throw new PHPExcel_Calculation_Exception(self::ArgumentTypeException); }
case 'array':
$M = new PHPExcel_Shared_JAMA_Matrix($args[0]);
break;
default:
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException);
break;
}
$this->checkMatrixDimensions($M);
for($i = 0; $i < $this->m; ++$i) {
for($j = 0; $j < $this->n; ++$j) {
$this->A[$i][$j] = trim($this->A[$i][$j],'"').trim($M->get($i, $j),'"');
}
}
return $this;
} else {
throw new PHPExcel_Calculation_Exception(self::PolymorphicArgumentException);
}
} // function concat()
 
 
/**
* Solve A*X = B.
*
* @param Matrix $B Right hand side
* @return Matrix ... Solution if A is square, least squares solution otherwise
*/
public function solve($B) {
if ($this->m == $this->n) {
$LU = new PHPExcel_Shared_JAMA_LUDecomposition($this);
return $LU->solve($B);
} else {
$QR = new QRDecomposition($this);
return $QR->solve($B);
}
} // function solve()
 
 
/**
* Matrix inverse or pseudoinverse.
*
* @return Matrix ... Inverse(A) if A is square, pseudoinverse otherwise.
*/
public function inverse() {
return $this->solve($this->identity($this->m, $this->m));
} // function inverse()
 
 
/**
* det
*
* Calculate determinant
* @return float Determinant
*/
public function det() {
$L = new PHPExcel_Shared_JAMA_LUDecomposition($this);
return $L->det();
} // function det()
 
 
} // class PHPExcel_Shared_JAMA_Matrix
/branches/v2.25-scarificateur/jrest/lib/PHPExcel/Classes/PHPExcel/Shared/JAMA/utils/Error.php
New file
0,0 → 1,82
<?php
/**
* @package JAMA
*
* Error handling
* @author Michael Bommarito
* @version 01292005
*/
 
//Language constant
define('JAMALANG', 'EN');
 
 
//All errors may be defined by the following format:
//define('ExceptionName', N);
//$error['lang'][ExceptionName] = 'Error message';
$error = array();
 
/*
I've used Babelfish and a little poor knowledge of Romance/Germanic languages for the translations here.
Feel free to correct anything that looks amiss to you.
*/
 
define('PolymorphicArgumentException', -1);
$error['EN'][PolymorphicArgumentException] = "Invalid argument pattern for polymorphic function.";
$error['FR'][PolymorphicArgumentException] = "Modèle inadmissible d'argument pour la fonction polymorphe.".
$error['DE'][PolymorphicArgumentException] = "Unzulässiges Argumentmuster für polymorphe Funktion.";
 
define('ArgumentTypeException', -2);
$error['EN'][ArgumentTypeException] = "Invalid argument type.";
$error['FR'][ArgumentTypeException] = "Type inadmissible d'argument.";
$error['DE'][ArgumentTypeException] = "Unzulässige Argumentart.";
 
define('ArgumentBoundsException', -3);
$error['EN'][ArgumentBoundsException] = "Invalid argument range.";
$error['FR'][ArgumentBoundsException] = "Gamme inadmissible d'argument.";
$error['DE'][ArgumentBoundsException] = "Unzulässige Argumentstrecke.";
 
define('MatrixDimensionException', -4);
$error['EN'][MatrixDimensionException] = "Matrix dimensions are not equal.";
$error['FR'][MatrixDimensionException] = "Les dimensions de Matrix ne sont pas égales.";
$error['DE'][MatrixDimensionException] = "Matrixmaße sind nicht gleich.";
 
define('PrecisionLossException', -5);
$error['EN'][PrecisionLossException] = "Significant precision loss detected.";
$error['FR'][PrecisionLossException] = "Perte significative de précision détectée.";
$error['DE'][PrecisionLossException] = "Bedeutender Präzision Verlust ermittelte.";
 
define('MatrixSPDException', -6);
$error['EN'][MatrixSPDException] = "Can only perform operation on symmetric positive definite matrix.";
$error['FR'][MatrixSPDException] = "Perte significative de précision détectée.";
$error['DE'][MatrixSPDException] = "Bedeutender Präzision Verlust ermittelte.";
 
define('MatrixSingularException', -7);
$error['EN'][MatrixSingularException] = "Can only perform operation on singular matrix.";
 
define('MatrixRankException', -8);
$error['EN'][MatrixRankException] = "Can only perform operation on full-rank matrix.";
 
define('ArrayLengthException', -9);
$error['EN'][ArrayLengthException] = "Array length must be a multiple of m.";
 
define('RowLengthException', -10);
$error['EN'][RowLengthException] = "All rows must have the same length.";
 
/**
* Custom error handler
* @param int $num Error number
*/
function JAMAError($errorNumber = null) {
global $error;
 
if (isset($errorNumber)) {
if (isset($error[JAMALANG][$errorNumber])) {
return $error[JAMALANG][$errorNumber];
} else {
return $error['EN'][$errorNumber];
}
} else {
return ("Invalid argument to JAMAError()");
}
}
/branches/v2.25-scarificateur/jrest/lib/PHPExcel/Classes/PHPExcel/Shared/JAMA/utils/Maths.php
New file
0,0 → 1,43
<?php
/**
* @package JAMA
*
* Pythagorean Theorem:
*
* a = 3
* b = 4
* r = sqrt(square(a) + square(b))
* r = 5
*
* r = sqrt(a^2 + b^2) without under/overflow.
*/
function hypo($a, $b) {
if (abs($a) > abs($b)) {
$r = $b / $a;
$r = abs($a) * sqrt(1 + $r * $r);
} elseif ($b != 0) {
$r = $a / $b;
$r = abs($b) * sqrt(1 + $r * $r);
} else {
$r = 0.0;
}
return $r;
} // function hypo()
 
 
/**
* Mike Bommarito's version.
* Compute n-dimensional hyotheneuse.
*
function hypot() {
$s = 0;
foreach (func_get_args() as $d) {
if (is_numeric($d)) {
$s += pow($d, 2);
} else {
throw new PHPExcel_Calculation_Exception(JAMAError(ArgumentTypeException));
}
}
return sqrt($s);
}
*/
/branches/v2.25-scarificateur/jrest/lib/PHPExcel/Classes/PHPExcel/Shared/JAMA/LUDecomposition.php
New file
0,0 → 1,258
<?php
/**
* @package JAMA
*
* For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n
* unit lower triangular matrix L, an n-by-n upper triangular matrix U,
* and a permutation vector piv of length m so that A(piv,:) = L*U.
* If m < n, then L is m-by-m and U is m-by-n.
*
* The LU decompostion with pivoting always exists, even if the matrix is
* singular, so the constructor will never fail. The primary use of the
* LU decomposition is in the solution of square systems of simultaneous
* linear equations. This will fail if isNonsingular() returns false.
*
* @author Paul Meagher
* @author Bartosz Matosiuk
* @author Michael Bommarito
* @version 1.1
* @license PHP v3.0
*/
class PHPExcel_Shared_JAMA_LUDecomposition {
 
const MatrixSingularException = "Can only perform operation on singular matrix.";
const MatrixSquareException = "Mismatched Row dimension";
 
/**
* Decomposition storage
* @var array
*/
private $LU = array();
 
/**
* Row dimension.
* @var int
*/
private $m;
 
/**
* Column dimension.
* @var int
*/
private $n;
 
/**
* Pivot sign.
* @var int
*/
private $pivsign;
 
/**
* Internal storage of pivot vector.
* @var array
*/
private $piv = array();
 
 
/**
* LU Decomposition constructor.
*
* @param $A Rectangular matrix
* @return Structure to access L, U and piv.
*/
public function __construct($A) {
if ($A instanceof PHPExcel_Shared_JAMA_Matrix) {
// Use a "left-looking", dot-product, Crout/Doolittle algorithm.
$this->LU = $A->getArray();
$this->m = $A->getRowDimension();
$this->n = $A->getColumnDimension();
for ($i = 0; $i < $this->m; ++$i) {
$this->piv[$i] = $i;
}
$this->pivsign = 1;
$LUrowi = $LUcolj = array();
 
// Outer loop.
for ($j = 0; $j < $this->n; ++$j) {
// Make a copy of the j-th column to localize references.
for ($i = 0; $i < $this->m; ++$i) {
$LUcolj[$i] = &$this->LU[$i][$j];
}
// Apply previous transformations.
for ($i = 0; $i < $this->m; ++$i) {
$LUrowi = $this->LU[$i];
// Most of the time is spent in the following dot product.
$kmax = min($i,$j);
$s = 0.0;
for ($k = 0; $k < $kmax; ++$k) {
$s += $LUrowi[$k] * $LUcolj[$k];
}
$LUrowi[$j] = $LUcolj[$i] -= $s;
}
// Find pivot and exchange if necessary.
$p = $j;
for ($i = $j+1; $i < $this->m; ++$i) {
if (abs($LUcolj[$i]) > abs($LUcolj[$p])) {
$p = $i;
}
}
if ($p != $j) {
for ($k = 0; $k < $this->n; ++$k) {
$t = $this->LU[$p][$k];
$this->LU[$p][$k] = $this->LU[$j][$k];
$this->LU[$j][$k] = $t;
}
$k = $this->piv[$p];
$this->piv[$p] = $this->piv[$j];
$this->piv[$j] = $k;
$this->pivsign = $this->pivsign * -1;
}
// Compute multipliers.
if (($j < $this->m) && ($this->LU[$j][$j] != 0.0)) {
for ($i = $j+1; $i < $this->m; ++$i) {
$this->LU[$i][$j] /= $this->LU[$j][$j];
}
}
}
} else {
throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::ArgumentTypeException);
}
} // function __construct()
 
 
/**
* Get lower triangular factor.
*
* @return array Lower triangular factor
*/
public function getL() {
for ($i = 0; $i < $this->m; ++$i) {
for ($j = 0; $j < $this->n; ++$j) {
if ($i > $j) {
$L[$i][$j] = $this->LU[$i][$j];
} elseif ($i == $j) {
$L[$i][$j] = 1.0;
} else {
$L[$i][$j] = 0.0;
}
}
}
return new PHPExcel_Shared_JAMA_Matrix($L);
} // function getL()
 
 
/**
* Get upper triangular factor.
*
* @return array Upper triangular factor
*/
public function getU() {
for ($i = 0; $i < $this->n; ++$i) {
for ($j = 0; $j < $this->n; ++$j) {
if ($i <= $j) {
$U[$i][$j] = $this->LU[$i][$j];
} else {
$U[$i][$j] = 0.0;
}
}
}
return new PHPExcel_Shared_JAMA_Matrix($U);
} // function getU()
 
 
/**
* Return pivot permutation vector.
*
* @return array Pivot vector
*/
public function getPivot() {
return $this->piv;
} // function getPivot()
 
 
/**
* Alias for getPivot
*
* @see getPivot
*/
public function getDoublePivot() {
return $this->getPivot();
} // function getDoublePivot()
 
 
/**
* Is the matrix nonsingular?
*
* @return true if U, and hence A, is nonsingular.
*/
public function isNonsingular() {
for ($j = 0; $j < $this->n; ++$j) {
if ($this->LU[$j][$j] == 0) {
return false;
}
}
return true;
} // function isNonsingular()
 
 
/**
* Count determinants
*
* @return array d matrix deterninat
*/
public function det() {
if ($this->m == $this->n) {
$d = $this->pivsign;
for ($j = 0; $j < $this->n; ++$j) {
$d *= $this->LU[$j][$j];
}
return $d;
} else {
throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::MatrixDimensionException);
}
} // function det()
 
 
/**
* Solve A*X = B
*
* @param $B A Matrix with as many rows as A and any number of columns.
* @return X so that L*U*X = B(piv,:)
* @PHPExcel_Calculation_Exception IllegalArgumentException Matrix row dimensions must agree.
* @PHPExcel_Calculation_Exception RuntimeException Matrix is singular.
*/
public function solve($B) {
if ($B->getRowDimension() == $this->m) {
if ($this->isNonsingular()) {
// Copy right hand side with pivoting
$nx = $B->getColumnDimension();
$X = $B->getMatrix($this->piv, 0, $nx-1);
// Solve L*Y = B(piv,:)
for ($k = 0; $k < $this->n; ++$k) {
for ($i = $k+1; $i < $this->n; ++$i) {
for ($j = 0; $j < $nx; ++$j) {
$X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k];
}
}
}
// Solve U*X = Y;
for ($k = $this->n-1; $k >= 0; --$k) {
for ($j = 0; $j < $nx; ++$j) {
$X->A[$k][$j] /= $this->LU[$k][$k];
}
for ($i = 0; $i < $k; ++$i) {
for ($j = 0; $j < $nx; ++$j) {
$X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k];
}
}
}
return $X;
} else {
throw new PHPExcel_Calculation_Exception(self::MatrixSingularException);
}
} else {
throw new PHPExcel_Calculation_Exception(self::MatrixSquareException);
}
} // function solve()
 
} // class PHPExcel_Shared_JAMA_LUDecomposition
/branches/v2.25-scarificateur/jrest/lib/PHPExcel/Classes/PHPExcel/Shared/JAMA/CholeskyDecomposition.php
New file
0,0 → 1,149
<?php
/**
* @package JAMA
*
* Cholesky decomposition class
*
* For a symmetric, positive definite matrix A, the Cholesky decomposition
* is an lower triangular matrix L so that A = L*L'.
*
* If the matrix is not symmetric or positive definite, the constructor
* returns a partial decomposition and sets an internal flag that may
* be queried by the isSPD() method.
*
* @author Paul Meagher
* @author Michael Bommarito
* @version 1.2
*/
class CholeskyDecomposition {
 
/**
* Decomposition storage
* @var array
* @access private
*/
private $L = array();
 
/**
* Matrix row and column dimension
* @var int
* @access private
*/
private $m;
 
/**
* Symmetric positive definite flag
* @var boolean
* @access private
*/
private $isspd = true;
 
 
/**
* CholeskyDecomposition
*
* Class constructor - decomposes symmetric positive definite matrix
* @param mixed Matrix square symmetric positive definite matrix
*/
public function __construct($A = null) {
if ($A instanceof Matrix) {
$this->L = $A->getArray();
$this->m = $A->getRowDimension();
 
for($i = 0; $i < $this->m; ++$i) {
for($j = $i; $j < $this->m; ++$j) {
for($sum = $this->L[$i][$j], $k = $i - 1; $k >= 0; --$k) {
$sum -= $this->L[$i][$k] * $this->L[$j][$k];
}
if ($i == $j) {
if ($sum >= 0) {
$this->L[$i][$i] = sqrt($sum);
} else {
$this->isspd = false;
}
} else {
if ($this->L[$i][$i] != 0) {
$this->L[$j][$i] = $sum / $this->L[$i][$i];
}
}
}
 
for ($k = $i+1; $k < $this->m; ++$k) {
$this->L[$i][$k] = 0.0;
}
}
} else {
throw new PHPExcel_Calculation_Exception(JAMAError(ArgumentTypeException));
}
} // function __construct()
 
 
/**
* Is the matrix symmetric and positive definite?
*
* @return boolean
*/
public function isSPD() {
return $this->isspd;
} // function isSPD()
 
 
/**
* getL
*
* Return triangular factor.
* @return Matrix Lower triangular matrix
*/
public function getL() {
return new Matrix($this->L);
} // function getL()
 
 
/**
* Solve A*X = B
*
* @param $B Row-equal matrix
* @return Matrix L * L' * X = B
*/
public function solve($B = null) {
if ($B instanceof Matrix) {
if ($B->getRowDimension() == $this->m) {
if ($this->isspd) {
$X = $B->getArrayCopy();
$nx = $B->getColumnDimension();
 
for ($k = 0; $k < $this->m; ++$k) {
for ($i = $k + 1; $i < $this->m; ++$i) {
for ($j = 0; $j < $nx; ++$j) {
$X[$i][$j] -= $X[$k][$j] * $this->L[$i][$k];
}
}
for ($j = 0; $j < $nx; ++$j) {
$X[$k][$j] /= $this->L[$k][$k];
}
}
 
for ($k = $this->m - 1; $k >= 0; --$k) {
for ($j = 0; $j < $nx; ++$j) {
$X[$k][$j] /= $this->L[$k][$k];
}
for ($i = 0; $i < $k; ++$i) {
for ($j = 0; $j < $nx; ++$j) {
$X[$i][$j] -= $X[$k][$j] * $this->L[$k][$i];
}
}
}
 
return new Matrix($X, $this->m, $nx);
} else {
throw new PHPExcel_Calculation_Exception(JAMAError(MatrixSPDException));
}
} else {
throw new PHPExcel_Calculation_Exception(JAMAError(MatrixDimensionException));
}
} else {
throw new PHPExcel_Calculation_Exception(JAMAError(ArgumentTypeException));
}
} // function solve()
 
} // class CholeskyDecomposition
/branches/v2.25-scarificateur/jrest/lib/PHPExcel/Classes/PHPExcel/Shared/JAMA/SingularValueDecomposition.php
New file
0,0 → 1,526
<?php
/**
* @package JAMA
*
* For an m-by-n matrix A with m >= n, the singular value decomposition is
* an m-by-n orthogonal matrix U, an n-by-n diagonal matrix S, and
* an n-by-n orthogonal matrix V so that A = U*S*V'.
*
* The singular values, sigma[$k] = S[$k][$k], are ordered so that
* sigma[0] >= sigma[1] >= ... >= sigma[n-1].
*
* The singular value decompostion always exists, so the constructor will
* never fail. The matrix condition number and the effective numerical
* rank can be computed from this decomposition.
*
* @author Paul Meagher
* @license PHP v3.0
* @version 1.1
*/
class SingularValueDecomposition {
 
/**
* Internal storage of U.
* @var array
*/
private $U = array();
 
/**
* Internal storage of V.
* @var array
*/
private $V = array();
 
/**
* Internal storage of singular values.
* @var array
*/
private $s = array();
 
/**
* Row dimension.
* @var int
*/
private $m;
 
/**
* Column dimension.
* @var int
*/
private $n;
 
 
/**
* Construct the singular value decomposition
*
* Derived from LINPACK code.
*
* @param $A Rectangular matrix
* @return Structure to access U, S and V.
*/
public function __construct($Arg) {
 
// Initialize.
$A = $Arg->getArrayCopy();
$this->m = $Arg->getRowDimension();
$this->n = $Arg->getColumnDimension();
$nu = min($this->m, $this->n);
$e = array();
$work = array();
$wantu = true;
$wantv = true;
$nct = min($this->m - 1, $this->n);
$nrt = max(0, min($this->n - 2, $this->m));
 
// Reduce A to bidiagonal form, storing the diagonal elements
// in s and the super-diagonal elements in e.
for ($k = 0; $k < max($nct,$nrt); ++$k) {
 
if ($k < $nct) {
// Compute the transformation for the k-th column and
// place the k-th diagonal in s[$k].
// Compute 2-norm of k-th column without under/overflow.
$this->s[$k] = 0;
for ($i = $k; $i < $this->m; ++$i) {
$this->s[$k] = hypo($this->s[$k], $A[$i][$k]);
}
if ($this->s[$k] != 0.0) {
if ($A[$k][$k] < 0.0) {
$this->s[$k] = -$this->s[$k];
}
for ($i = $k; $i < $this->m; ++$i) {
$A[$i][$k] /= $this->s[$k];
}
$A[$k][$k] += 1.0;
}
$this->s[$k] = -$this->s[$k];
}
 
for ($j = $k + 1; $j < $this->n; ++$j) {
if (($k < $nct) & ($this->s[$k] != 0.0)) {
// Apply the transformation.
$t = 0;
for ($i = $k; $i < $this->m; ++$i) {
$t += $A[$i][$k] * $A[$i][$j];
}
$t = -$t / $A[$k][$k];
for ($i = $k; $i < $this->m; ++$i) {
$A[$i][$j] += $t * $A[$i][$k];
}
// Place the k-th row of A into e for the
// subsequent calculation of the row transformation.
$e[$j] = $A[$k][$j];
}
}
 
if ($wantu AND ($k < $nct)) {
// Place the transformation in U for subsequent back
// multiplication.
for ($i = $k; $i < $this->m; ++$i) {
$this->U[$i][$k] = $A[$i][$k];
}
}
 
if ($k < $nrt) {
// Compute the k-th row transformation and place the
// k-th super-diagonal in e[$k].
// Compute 2-norm without under/overflow.
$e[$k] = 0;
for ($i = $k + 1; $i < $this->n; ++$i) {
$e[$k] = hypo($e[$k], $e[$i]);
}
if ($e[$k] != 0.0) {
if ($e[$k+1] < 0.0) {
$e[$k] = -$e[$k];
}
for ($i = $k + 1; $i < $this->n; ++$i) {
$e[$i] /= $e[$k];
}
$e[$k+1] += 1.0;
}
$e[$k] = -$e[$k];
if (($k+1 < $this->m) AND ($e[$k] != 0.0)) {
// Apply the transformation.
for ($i = $k+1; $i < $this->m; ++$i) {
$work[$i] = 0.0;
}
for ($j = $k+1; $j < $this->n; ++$j) {
for ($i = $k+1; $i < $this->m; ++$i) {
$work[$i] += $e[$j] * $A[$i][$j];
}
}
for ($j = $k + 1; $j < $this->n; ++$j) {
$t = -$e[$j] / $e[$k+1];
for ($i = $k + 1; $i < $this->m; ++$i) {
$A[$i][$j] += $t * $work[$i];
}
}
}
if ($wantv) {
// Place the transformation in V for subsequent
// back multiplication.
for ($i = $k + 1; $i < $this->n; ++$i) {
$this->V[$i][$k] = $e[$i];
}
}
}
}
 
// Set up the final bidiagonal matrix or order p.
$p = min($this->n, $this->m + 1);
if ($nct < $this->n) {
$this->s[$nct] = $A[$nct][$nct];
}
if ($this->m < $p) {
$this->s[$p-1] = 0.0;
}
if ($nrt + 1 < $p) {
$e[$nrt] = $A[$nrt][$p-1];
}
$e[$p-1] = 0.0;
// If required, generate U.
if ($wantu) {
for ($j = $nct; $j < $nu; ++$j) {
for ($i = 0; $i < $this->m; ++$i) {
$this->U[$i][$j] = 0.0;
}
$this->U[$j][$j] = 1.0;
}
for ($k = $nct - 1; $k >= 0; --$k) {
if ($this->s[$k] != 0.0) {
for ($j = $k + 1; $j < $nu; ++$j) {
$t = 0;
for ($i = $k; $i < $this->m; ++$i) {
$t += $this->U[$i][$k] * $this->U[$i][$j];
}
$t = -$t / $this->U[$k][$k];
for ($i = $k; $i < $this->m; ++$i) {
$this->U[$i][$j] += $t * $this->U[$i][$k];
}
}
for ($i = $k; $i < $this->m; ++$i ) {
$this->U[$i][$k] = -$this->U[$i][$k];
}
$this->U[$k][$k] = 1.0 + $this->U[$k][$k];
for ($i = 0; $i < $k - 1; ++$i) {
$this->U[$i][$k] = 0.0;
}
} else {
for ($i = 0; $i < $this->m; ++$i) {
$this->U[$i][$k] = 0.0;
}
$this->U[$k][$k] = 1.0;
}
}
}
 
// If required, generate V.
if ($wantv) {
for ($k = $this->n - 1; $k >= 0; --$k) {
if (($k < $nrt) AND ($e[$k] != 0.0)) {
for ($j = $k + 1; $j < $nu; ++$j) {
$t = 0;
for ($i = $k + 1; $i < $this->n; ++$i) {
$t += $this->V[$i][$k]* $this->V[$i][$j];
}
$t = -$t / $this->V[$k+1][$k];
for ($i = $k + 1; $i < $this->n; ++$i) {
$this->V[$i][$j] += $t * $this->V[$i][$k];
}
}
}
for ($i = 0; $i < $this->n; ++$i) {
$this->V[$i][$k] = 0.0;
}
$this->V[$k][$k] = 1.0;
}
}
 
// Main iteration loop for the singular values.
$pp = $p - 1;
$iter = 0;
$eps = pow(2.0, -52.0);
 
while ($p > 0) {
// Here is where a test for too many iterations would go.
// This section of the program inspects for negligible
// elements in the s and e arrays. On completion the
// variables kase and k are set as follows:
// kase = 1 if s(p) and e[k-1] are negligible and k<p
// kase = 2 if s(k) is negligible and k<p
// kase = 3 if e[k-1] is negligible, k<p, and
// s(k), ..., s(p) are not negligible (qr step).
// kase = 4 if e(p-1) is negligible (convergence).
for ($k = $p - 2; $k >= -1; --$k) {
if ($k == -1) {
break;
}
if (abs($e[$k]) <= $eps * (abs($this->s[$k]) + abs($this->s[$k+1]))) {
$e[$k] = 0.0;
break;
}
}
if ($k == $p - 2) {
$kase = 4;
} else {
for ($ks = $p - 1; $ks >= $k; --$ks) {
if ($ks == $k) {
break;
}
$t = ($ks != $p ? abs($e[$ks]) : 0.) + ($ks != $k + 1 ? abs($e[$ks-1]) : 0.);
if (abs($this->s[$ks]) <= $eps * $t) {
$this->s[$ks] = 0.0;
break;
}
}
if ($ks == $k) {
$kase = 3;
} else if ($ks == $p-1) {
$kase = 1;
} else {
$kase = 2;
$k = $ks;
}
}
++$k;
 
// Perform the task indicated by kase.
switch ($kase) {
// Deflate negligible s(p).
case 1:
$f = $e[$p-2];
$e[$p-2] = 0.0;
for ($j = $p - 2; $j >= $k; --$j) {
$t = hypo($this->s[$j],$f);
$cs = $this->s[$j] / $t;
$sn = $f / $t;
$this->s[$j] = $t;
if ($j != $k) {
$f = -$sn * $e[$j-1];
$e[$j-1] = $cs * $e[$j-1];
}
if ($wantv) {
for ($i = 0; $i < $this->n; ++$i) {
$t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$p-1];
$this->V[$i][$p-1] = -$sn * $this->V[$i][$j] + $cs * $this->V[$i][$p-1];
$this->V[$i][$j] = $t;
}
}
}
break;
// Split at negligible s(k).
case 2:
$f = $e[$k-1];
$e[$k-1] = 0.0;
for ($j = $k; $j < $p; ++$j) {
$t = hypo($this->s[$j], $f);
$cs = $this->s[$j] / $t;
$sn = $f / $t;
$this->s[$j] = $t;
$f = -$sn * $e[$j];
$e[$j] = $cs * $e[$j];
if ($wantu) {
for ($i = 0; $i < $this->m; ++$i) {
$t = $cs * $this->U[$i][$j] + $sn * $this->U[$i][$k-1];
$this->U[$i][$k-1] = -$sn * $this->U[$i][$j] + $cs * $this->U[$i][$k-1];
$this->U[$i][$j] = $t;
}
}
}
break;
// Perform one qr step.
case 3:
// Calculate the shift.
$scale = max(max(max(max(
abs($this->s[$p-1]),abs($this->s[$p-2])),abs($e[$p-2])),
abs($this->s[$k])), abs($e[$k]));
$sp = $this->s[$p-1] / $scale;
$spm1 = $this->s[$p-2] / $scale;
$epm1 = $e[$p-2] / $scale;
$sk = $this->s[$k] / $scale;
$ek = $e[$k] / $scale;
$b = (($spm1 + $sp) * ($spm1 - $sp) + $epm1 * $epm1) / 2.0;
$c = ($sp * $epm1) * ($sp * $epm1);
$shift = 0.0;
if (($b != 0.0) || ($c != 0.0)) {
$shift = sqrt($b * $b + $c);
if ($b < 0.0) {
$shift = -$shift;
}
$shift = $c / ($b + $shift);
}
$f = ($sk + $sp) * ($sk - $sp) + $shift;
$g = $sk * $ek;
// Chase zeros.
for ($j = $k; $j < $p-1; ++$j) {
$t = hypo($f,$g);
$cs = $f/$t;
$sn = $g/$t;
if ($j != $k) {
$e[$j-1] = $t;
}
$f = $cs * $this->s[$j] + $sn * $e[$j];
$e[$j] = $cs * $e[$j] - $sn * $this->s[$j];
$g = $sn * $this->s[$j+1];
$this->s[$j+1] = $cs * $this->s[$j+1];
if ($wantv) {
for ($i = 0; $i < $this->n; ++$i) {
$t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$j+1];
$this->V[$i][$j+1] = -$sn * $this->V[$i][$j] + $cs * $this->V[$i][$j+1];
$this->V[$i][$j] = $t;
}
}
$t = hypo($f,$g);
$cs = $f/$t;
$sn = $g/$t;
$this->s[$j] = $t;
$f = $cs * $e[$j] + $sn * $this->s[$j+1];
$this->s[$j+1] = -$sn * $e[$j] + $cs * $this->s[$j+1];
$g = $sn * $e[$j+1];
$e[$j+1] = $cs * $e[$j+1];
if ($wantu && ($j < $this->m - 1)) {
for ($i = 0; $i < $this->m; ++$i) {
$t = $cs * $this->U[$i][$j] + $sn * $this->U[$i][$j+1];
$this->U[$i][$j+1] = -$sn * $this->U[$i][$j] + $cs * $this->U[$i][$j+1];
$this->U[$i][$j] = $t;
}
}
}
$e[$p-2] = $f;
$iter = $iter + 1;
break;
// Convergence.
case 4:
// Make the singular values positive.
if ($this->s[$k] <= 0.0) {
$this->s[$k] = ($this->s[$k] < 0.0 ? -$this->s[$k] : 0.0);
if ($wantv) {
for ($i = 0; $i <= $pp; ++$i) {
$this->V[$i][$k] = -$this->V[$i][$k];
}
}
}
// Order the singular values.
while ($k < $pp) {
if ($this->s[$k] >= $this->s[$k+1]) {
break;
}
$t = $this->s[$k];
$this->s[$k] = $this->s[$k+1];
$this->s[$k+1] = $t;
if ($wantv AND ($k < $this->n - 1)) {
for ($i = 0; $i < $this->n; ++$i) {
$t = $this->V[$i][$k+1];
$this->V[$i][$k+1] = $this->V[$i][$k];
$this->V[$i][$k] = $t;
}
}
if ($wantu AND ($k < $this->m-1)) {
for ($i = 0; $i < $this->m; ++$i) {
$t = $this->U[$i][$k+1];
$this->U[$i][$k+1] = $this->U[$i][$k];
$this->U[$i][$k] = $t;
}
}
++$k;
}
$iter = 0;
--$p;
break;
} // end switch
} // end while
 
} // end constructor
 
 
/**
* Return the left singular vectors
*
* @access public
* @return U
*/
public function getU() {
return new Matrix($this->U, $this->m, min($this->m + 1, $this->n));
}
 
 
/**
* Return the right singular vectors
*
* @access public
* @return V
*/
public function getV() {
return new Matrix($this->V);
}
 
 
/**
* Return the one-dimensional array of singular values
*
* @access public
* @return diagonal of S.
*/
public function getSingularValues() {
return $this->s;
}
 
 
/**
* Return the diagonal matrix of singular values
*
* @access public
* @return S
*/
public function getS() {
for ($i = 0; $i < $this->n; ++$i) {
for ($j = 0; $j < $this->n; ++$j) {
$S[$i][$j] = 0.0;
}
$S[$i][$i] = $this->s[$i];
}
return new Matrix($S);
}
 
 
/**
* Two norm
*
* @access public
* @return max(S)
*/
public function norm2() {
return $this->s[0];
}
 
 
/**
* Two norm condition number
*
* @access public
* @return max(S)/min(S)
*/
public function cond() {
return $this->s[0] / $this->s[min($this->m, $this->n) - 1];
}
 
 
/**
* Effective numerical matrix rank
*
* @access public
* @return Number of nonnegligible singular values.
*/
public function rank() {
$eps = pow(2.0, -52.0);
$tol = max($this->m, $this->n) * $this->s[0] * $eps;
$r = 0;
for ($i = 0; $i < count($this->s); ++$i) {
if ($this->s[$i] > $tol) {
++$r;
}
}
return $r;
}
 
} // class SingularValueDecomposition
/branches/v2.25-scarificateur/jrest/lib/PHPExcel/Classes/PHPExcel/Shared/JAMA/CHANGELOG.TXT
New file
0,0 → 1,16
Mar 1, 2005 11:15 AST by PM
 
+ For consistency, renamed Math.php to Maths.java, utils to util,
tests to test, docs to doc -
 
+ Removed conditional logic from top of Matrix class.
 
+ Switched to using hypo function in Maths.php for all php-hypot calls.
NOTE TO SELF: Need to make sure that all decompositions have been
switched over to using the bundled hypo.
 
Feb 25, 2005 at 10:00 AST by PM
 
+ Recommend using simpler Error.php instead of JAMA_Error.php but
can be persuaded otherwise.
 
/branches/v2.25-scarificateur/jrest/lib/PHPExcel/Classes/PHPExcel/Shared/JAMA/QRDecomposition.php
New file
0,0 → 1,234
<?php
/**
* @package JAMA
*
* For an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n
* orthogonal matrix Q and an n-by-n upper triangular matrix R so that
* A = Q*R.
*
* The QR decompostion always exists, even if the matrix does not have
* full rank, so the constructor will never fail. The primary use of the
* QR decomposition is in the least squares solution of nonsquare systems
* of simultaneous linear equations. This will fail if isFullRank()
* returns false.
*
* @author Paul Meagher
* @license PHP v3.0
* @version 1.1
*/
class PHPExcel_Shared_JAMA_QRDecomposition {
 
const MatrixRankException = "Can only perform operation on full-rank matrix.";
 
/**
* Array for internal storage of decomposition.
* @var array
*/
private $QR = array();
 
/**
* Row dimension.
* @var integer
*/
private $m;
 
/**
* Column dimension.
* @var integer
*/
private $n;
 
/**
* Array for internal storage of diagonal of R.
* @var array
*/
private $Rdiag = array();
 
 
/**
* QR Decomposition computed by Householder reflections.
*
* @param matrix $A Rectangular matrix
* @return Structure to access R and the Householder vectors and compute Q.
*/
public function __construct($A) {
if($A instanceof PHPExcel_Shared_JAMA_Matrix) {
// Initialize.
$this->QR = $A->getArrayCopy();
$this->m = $A->getRowDimension();
$this->n = $A->getColumnDimension();
// Main loop.
for ($k = 0; $k < $this->n; ++$k) {
// Compute 2-norm of k-th column without under/overflow.
$nrm = 0.0;
for ($i = $k; $i < $this->m; ++$i) {
$nrm = hypo($nrm, $this->QR[$i][$k]);
}
if ($nrm != 0.0) {
// Form k-th Householder vector.
if ($this->QR[$k][$k] < 0) {
$nrm = -$nrm;
}
for ($i = $k; $i < $this->m; ++$i) {
$this->QR[$i][$k] /= $nrm;
}
$this->QR[$k][$k] += 1.0;
// Apply transformation to remaining columns.
for ($j = $k+1; $j < $this->n; ++$j) {
$s = 0.0;
for ($i = $k; $i < $this->m; ++$i) {
$s += $this->QR[$i][$k] * $this->QR[$i][$j];
}
$s = -$s/$this->QR[$k][$k];
for ($i = $k; $i < $this->m; ++$i) {
$this->QR[$i][$j] += $s * $this->QR[$i][$k];
}
}
}
$this->Rdiag[$k] = -$nrm;
}
} else {
throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::ArgumentTypeException);
}
} // function __construct()
 
 
/**
* Is the matrix full rank?
*
* @return boolean true if R, and hence A, has full rank, else false.
*/
public function isFullRank() {
for ($j = 0; $j < $this->n; ++$j) {
if ($this->Rdiag[$j] == 0) {
return false;
}
}
return true;
} // function isFullRank()
 
 
/**
* Return the Householder vectors
*
* @return Matrix Lower trapezoidal matrix whose columns define the reflections
*/
public function getH() {
for ($i = 0; $i < $this->m; ++$i) {
for ($j = 0; $j < $this->n; ++$j) {
if ($i >= $j) {
$H[$i][$j] = $this->QR[$i][$j];
} else {
$H[$i][$j] = 0.0;
}
}
}
return new PHPExcel_Shared_JAMA_Matrix($H);
} // function getH()
 
 
/**
* Return the upper triangular factor
*
* @return Matrix upper triangular factor
*/
public function getR() {
for ($i = 0; $i < $this->n; ++$i) {
for ($j = 0; $j < $this->n; ++$j) {
if ($i < $j) {
$R[$i][$j] = $this->QR[$i][$j];
} elseif ($i == $j) {
$R[$i][$j] = $this->Rdiag[$i];
} else {
$R[$i][$j] = 0.0;
}
}
}
return new PHPExcel_Shared_JAMA_Matrix($R);
} // function getR()
 
 
/**
* Generate and return the (economy-sized) orthogonal factor
*
* @return Matrix orthogonal factor
*/
public function getQ() {
for ($k = $this->n-1; $k >= 0; --$k) {
for ($i = 0; $i < $this->m; ++$i) {
$Q[$i][$k] = 0.0;
}
$Q[$k][$k] = 1.0;
for ($j = $k; $j < $this->n; ++$j) {
if ($this->QR[$k][$k] != 0) {
$s = 0.0;
for ($i = $k; $i < $this->m; ++$i) {
$s += $this->QR[$i][$k] * $Q[$i][$j];
}
$s = -$s/$this->QR[$k][$k];
for ($i = $k; $i < $this->m; ++$i) {
$Q[$i][$j] += $s * $this->QR[$i][$k];
}
}
}
}
/*
for($i = 0; $i < count($Q); ++$i) {
for($j = 0; $j < count($Q); ++$j) {
if(! isset($Q[$i][$j]) ) {
$Q[$i][$j] = 0;
}
}
}
*/
return new PHPExcel_Shared_JAMA_Matrix($Q);
} // function getQ()
 
 
/**
* Least squares solution of A*X = B
*
* @param Matrix $B A Matrix with as many rows as A and any number of columns.
* @return Matrix Matrix that minimizes the two norm of Q*R*X-B.
*/
public function solve($B) {
if ($B->getRowDimension() == $this->m) {
if ($this->isFullRank()) {
// Copy right hand side
$nx = $B->getColumnDimension();
$X = $B->getArrayCopy();
// Compute Y = transpose(Q)*B
for ($k = 0; $k < $this->n; ++$k) {
for ($j = 0; $j < $nx; ++$j) {
$s = 0.0;
for ($i = $k; $i < $this->m; ++$i) {
$s += $this->QR[$i][$k] * $X[$i][$j];
}
$s = -$s/$this->QR[$k][$k];
for ($i = $k; $i < $this->m; ++$i) {
$X[$i][$j] += $s * $this->QR[$i][$k];
}
}
}
// Solve R*X = Y;
for ($k = $this->n-1; $k >= 0; --$k) {
for ($j = 0; $j < $nx; ++$j) {
$X[$k][$j] /= $this->Rdiag[$k];
}
for ($i = 0; $i < $k; ++$i) {
for ($j = 0; $j < $nx; ++$j) {
$X[$i][$j] -= $X[$k][$j]* $this->QR[$i][$k];
}
}
}
$X = new PHPExcel_Shared_JAMA_Matrix($X);
return ($X->getMatrix(0, $this->n-1, 0, $nx));
} else {
throw new PHPExcel_Calculation_Exception(self::MatrixRankException);
}
} else {
throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::MatrixDimensionException);
}
} // function solve()
 
} // PHPExcel_Shared_JAMA_class QRDecomposition