30 private $LU = array();
66 $this->LU = $A->getArray();
67 $this->m = $A->getRowDimension();
68 $this->n = $A->getColumnDimension();
73 $LUrowi = $LUcolj = array();
79 $LUcolj[$i] = &$this->LU[$i][$j];
83 $LUrowi = $this->LU[$i];
87 for ($k = 0; $k < $kmax; ++$k) {
88 $s += $LUrowi[$k] * $LUcolj[$k];
90 $LUrowi[$j] = $LUcolj[$i] -= $s;
94 for ($i = $j+1; $i <
$this->m; ++$i) {
95 if (abs($LUcolj[$i]) > abs($LUcolj[$p])) {
101 $t = $this->LU[$p][$k];
102 $this->LU[$p][$k] = $this->LU[$j][$k];
103 $this->LU[$j][$k] =
$t;
106 $this->piv[$p] = $this->piv[$j];
108 $this->pivsign = $this->pivsign * -1;
111 if (($j < $this->m) && ($this->LU[$j][$j] != 0.0)) {
112 for ($i = $j+1; $i <
$this->m; ++$i) {
113 $this->LU[$i][$j] /= $this->LU[$j][$j];
132 $L[$i][$j] = $this->LU[$i][$j];
133 } elseif ($i == $j) {
153 $U[$i][$j] = $this->LU[$i][$j];
190 if ($this->LU[$j][$j] == 0) {
204 if ($this->m == $this->n) {
207 $d *= $this->LU[$j][$j];
225 if ($B->getRowDimension() == $this->m) {
228 $nx = $B->getColumnDimension();
229 $X = $B->getMatrix($this->piv, 0, $nx-1);
232 for ($i = $k+1; $i <
$this->n; ++$i) {
233 for ($j = 0; $j < $nx; ++$j) {
234 $X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k];
239 for ($k = $this->n-1; $k >= 0; --$k) {
240 for ($j = 0; $j < $nx; ++$j) {
241 $X->A[$k][$j] /= $this->LU[$k][$k];
243 for ($i = 0; $i < $k; ++$i) {
244 for ($j = 0; $j < $nx; ++$j) {
245 $X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k];
for($col=0; $col< 50; $col++) $d
const MatrixSingularException
An exception for terminatinating execution or to throw for unit testing.
__construct($A)
LU Decomposition constructor.
getPivot()
Return pivot permutation vector.
getL()
Get lower triangular factor.
getDoublePivot()
Alias for getPivot.
const MatrixSingularException
const MatrixSquareException
getU()
Get upper triangular factor.
isNonsingular()
Is the matrix nonsingular?
const ArgumentTypeException
const MatrixDimensionException