Blame | Last modification | View Log | RSS feed
<?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