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.
257 lines
7.9 KiB
257 lines
7.9 KiB
<?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 MATRIX_SINGULAR_EXCEPTION = "Can only perform operation on singular matrix."; |
|
const MATRIX_SQUARE_EXCEPTION = "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::ARGUMENT_TYPE_EXCEPTION); |
|
} |
|
} // 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::MATRIX_DIMENSION_EXCEPTION); |
|
} |
|
} // 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::MATRIX_SINGULAR_EXCEPTION); |
|
} |
|
} else { |
|
throw new PHPExcel_Calculation_Exception(self::MATRIX_SQUARE_EXCEPTION); |
|
} |
|
} |
|
}
|
|
|