1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19: 20: 21: 22: 23: 24: 25: 26: 27: 28: 29: 30: 31: 32: 33: 34: 35: 36: 37: 38: 39: 40: 41: 42: 43: 44: 45: 46: 47: 48: 49: 50: 51: 52: 53: 54: 55: 56: 57: 58: 59: 60: 61: 62: 63: 64: 65: 66: 67: 68: 69: 70: 71: 72: 73: 74: 75: 76: 77: 78: 79: 80: 81: 82: 83: 84: 85: 86: 87: 88: 89: 90: 91: 92: 93: 94: 95: 96: 97: 98: 99: 100: 101: 102: 103: 104: 105: 106: 107: 108: 109: 110: 111: 112: 113: 114: 115: 116: 117: 118: 119: 120: 121: 122: 123: 124: 125: 126: 127: 128: 129: 130: 131: 132: 133: 134: 135: 136: 137: 138: 139: 140: 141: 142: 143: 144: 145: 146: 147: 148: 149: 150: 151: 152: 153: 154: 155: 156: 157: 158: 159: 160: 161: 162: 163: 164: 165: 166: 167: 168: 169: 170: 171: 172: 173: 174: 175: 176: 177: 178: 179: 180: 181: 182: 183: 184: 185: 186: 187: 188: 189: 190: 191: 192: 193: 194: 195: 196: 197: 198: 199: 200: 201: 202: 203: 204: 205: 206: 207: 208: 209: 210: 211: 212: 213: 214: 215: 216: 217: 218: 219: 220: 221: 222: 223: 224: 225: 226: 227: 228: 229: 230: 231: 232: 233: 234: 235: 236: 237: 238: 239: 240: 241: 242: 243: 244: 245: 246: 247: 248: 249: 250: 251: 252: 253: 254: 255: 256: 257: 258: 259: 260: 261: 262: 263: 264: 265: 266: 267: 268: 269: 270: 271: 272: 273: 274: 275: 276: 277: 278: 279: 280: 281: 282: 283: 284: 285: 286: 287: 288: 289: 290: 291: 292: 293: 294: 295: 296: 297: 298: 299: 300: 301: 302: 303: 304: 305: 306: 307: 308: 309: 310: 311: 312: 313: 314:
<?php
namespace Chippyash\Math\Matrix\Decomposition;
use Chippyash\Math\Matrix\NumericMatrix;
use Chippyash\Math\Type\Calculator;
use Chippyash\Math\Type\Comparator;
use Chippyash\Math\Matrix\Traits\CreateCorrectMatrixType;
use Chippyash\Math\Matrix\Traits\CreateCorrectScalarType;
class Lu extends AbstractDecomposition
{
use CreateCorrectMatrixType;
use CreateCorrectScalarType;
protected $products = array(
'LU' => null,
'L' => null,
'U' => null,
'PivotVector' => null,
'PermutationMatrix' => null,
'Det' => null
);
protected $piv = [];
protected $pivsign;
protected $rows;
protected $cols;
public function decompose(NumericMatrix $mA, $extra = null)
{
$this->LUDecomposition($mA);
$this->setOtherProducts($mA);
return clone $this;
}
protected function LUDecomposition(NumericMatrix $mA)
{
$LU = $mA->toArray();
$m = $this->rows = $mA->rows();
$n = $this->cols = $mA->columns();
for ($i = 0; $i < $m; $i++) {
$this->piv[$i] = $i;
}
$this->pivsign = 1;
$LUrowi = [];
$LUcolj = [];
$calc = new Calculator();
$comp = new Comparator();
$zeroInt = $this->createCorrectScalarType($mA, 0);
for ($j = 0; $j < $n; $j++) {
for ($i = 0; $i < $m; $i++) {
$LUcolj[$i] = &$LU[$i][$j];
}
for ($i = 0; $i < $m; $i++) {
$LUrowi = $LU[$i];
$kmax = min($i, $j);
$s = clone $zeroInt;
for ($k = 0; $k < $kmax; $k++) {
$s = $calc->add($s, $calc->mul($LUrowi[$k], $LUcolj[$k]));
}
$LUcolj[$i] = $calc->sub($LUcolj[$i], $s);
$LUrowi[$j] = $LUcolj[$i];
}
$p = $j;
for ($i = $j + 1; $i < $m; $i++) {
if ($comp->gt($LUcolj[$i]->abs(), $LUcolj[$p]->abs())) {
$p = $i;
}
}
if ($p != $j) {
for ($k = 0; $k < $n; $k++) {
$t = $LU[$p][$k];
$LU[$p][$k] = $LU[$j][$k];
$LU[$j][$k] = $t;
}
$k = $this->piv[$p];
$this->piv[$p] = $this->piv[$j];
$this->piv[$j] = $k;
$this->pivsign = $this->pivsign * -1;
}
if (($j < $m) && $comp->neq($LU[$j][$j], $zeroInt)) {
for ($i = $j + 1; $i < $m; $i++) {
$LU[$i][$j] = $calc->div($LU[$i][$j],$LU[$j][$j]);
}
}
}
$this->set('LU', $this->createCorrectMatrixType($mA, $LU));
}
protected function setOtherProducts(NumericMatrix $mA)
{
$this->setLowerProduct($mA);
$this->setUpperProduct($mA);
$this->setPivotVector($mA);
$this->setPermutationMatrix($mA);
$this->setDeterminant($mA);
}
protected function setLowerProduct(NumericMatrix $mA)
{
$m = $this->LU->rows();
$n = $this->LU->columns();
$LU = $this->LU->toArray();
$rcFactor = $this->rows -$this->cols;
$this->set(
'L',
function() use ($m, $n, $LU, $rcFactor, $mA) {
$L = [];
for ($i = 0; $i < $m; $i++) {
for ($j = 0; $j < $n; $j++) {
if ($i > $j) {
$L[$i][$j] = $LU[$i][$j];
} elseif ($i == $j) {
$L[$i][$j] = $this->createCorrectScalarType($mA, 1);
} else {
$L[$i][$j] = $this->createCorrectScalarType($mA, 0);
}
}
}
if ($rcFactor < 0) {
$mLL = new NumericMatrix($L);
return $this->createCorrectMatrixType(
$mA,
$mLL('Colslice', array(1,$mLL->columns()+$rcFactor))->toArray());
} else {
return $this->createCorrectMatrixType($mA, $L);
}
}
);
}
protected function setUpperProduct(NumericMatrix $mA)
{
$n = $this->LU->columns();
$LU = $this->LU->toArray();
$rcFactor = $this->cols -$this->rows;
$this->set(
'U',
function() use ($n, $LU, $rcFactor, $mA) {
$U = [];
for ($i = 0; $i < $n; $i++) {
for ($j = 0; $j < $n; $j++) {
if ($i <= $j) {
$U[$i][$j] = (isset($LU[$i][$j]) ? $LU[$i][$j] : $this->createCorrectScalarType($mA, 0));
} else {
$U[$i][$j] = $this->createCorrectScalarType($mA, 0);
}
}
}
if ($rcFactor > 0) {
$mUU = new NumericMatrix($U);
return $this->createCorrectMatrixType(
$mA,
$mUU('Rowslice', array(1,$mUU->rows()-$rcFactor))->toArray());
} else {
return $this->createCorrectMatrixType($mA, $U);
}
}
);
}
protected function setPivotVector(NumericMatrix $mA)
{
$mB = $this->createCorrectMatrixType($mA, [$this->piv]);
$this->set(
'PivotVector',
$mB("Add\\Scalar", 1));
}
protected function setPermutationMatrix(NumericMatrix $mA)
{
$p = $this->piv;
$this->set(
'PermutationMatrix',
function() use ($p, $mA) {
$size = count($p);
$perm = array_fill(0, $size, array_fill(0, $size, 0));
for ($j=0; $j<$size; $j++) {
$perm[array_shift($p)][$j] = 1;
}
return $this->createCorrectMatrixType($mA,$perm);
}
);
}
protected function setDeterminant(NumericMatrix $mA)
{
if (!$mA->is('square')) {
$this->set('Det', null);
return;
}
if ($mA->is('empty')) {
$this->set('Det', $this->createCorrectScalarType($mA, 1));
return;
}
$det = $this->createCorrectScalarType($mA, $this->pivsign);
$LU = $this->LU->toArray();
$calc = new Calculator();
$this->set(
'Det',
function() use ($det, $LU, $calc) {
$c = count($LU);
for ($j = 0; $j < $c; $j++) {
$det = $calc->mul($det, $LU[$j][$j]);
}
return $det;
}
);
}
}