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:
<?php
namespace Chippyash\Math\Matrix\Decomposition;
use Chippyash\Math\Matrix\NumericMatrix;
use Chippyash\Math\Matrix\Traits\CreateCorrectMatrixType;
use Chippyash\Math\Matrix\Traits\AssertMatrixIsNumeric;
use Chippyash\Math\Matrix\Exceptions\SingularMatrixException;
use Chippyash\Matrix\Traits\AssertParameterIsMatrix;
use Chippyash\Matrix\Traits\AssertMatrixIsSquare;
use Chippyash\Matrix\Traits\AssertMatrixRowsAreEqual;
use Chippyash\Math\Type\Calculator;
use Chippyash\Math\Type\Comparator;
use Chippyash\Type\Number\Rational\RationalTypeFactory;
class GaussJordanElimination extends AbstractDecomposition
{
use AssertParameterIsMatrix;
use AssertMatrixIsNumeric;
use AssertMatrixIsSquare;
use AssertMatrixRowsAreEqual;
use CreateCorrectMatrixType;
protected $products = array(
'left' => null,
'right' => null
);
public function decompose(NumericMatrix $mA, $extra = null)
{
$this->assertParameterIsMatrix($extra, 'Parameter extra is not a matrix')
->assertMatrixIsNumeric($extra, 'Parameter extra is not a numeric matrix')
->assertMatrixIsSquare($mA,
'Parameter mA is not a square matrix')
->assertMatrixRowsAreEqual($mA, $extra, 'mA->rows != extra->rows');
$rows = $mA->rows();
$dA = $mA->toArray();
$dB = $extra->toArray();
$zero = function(){return RationalTypeFactory::create(0, 1);};
$one = function(){return RationalTypeFactory::create(1, 1);};
$calc = new Calculator();
$comp = new Comparator();
$ipiv = array_fill(0, $rows, $zero());
$indxr = array_fill(0, $rows, 0);
$indxc = array_fill(0, $rows, 0);
$irow = $icol = 0;
for ($i = 0; $i < $rows; $i++) {
$big = $zero();
for ($j = 0; $j < $rows; $j++) {
if ($comp->neq($ipiv[$j], $one())) {
for ($k = 0; $k < $rows; $k++) {
if ($comp->eq($ipiv[$k], $zero())) {
$absVal = $dA[$j][$k]->abs();
if ($comp->gt($absVal, $big)) {
unset($big);
$big = clone $absVal;
$irow = $j;
$icol = $k;
}
} elseif ($comp->gt($ipiv[$k], $one())) {
throw new SingularMatrixException('GaussJordanElimination');
}
}
}
}
$ipiv[$icol] = $calc->add($ipiv[$icol], $one());
if ($irow != $icol) {
$this->swapRows($dA, $irow, $icol);
$this->swapRows($dB, $irow, $icol);
}
$indxr[$i] = $irow;
$indxc[$i] = $icol;
if ($comp->eq($dA[$icol][$icol], $zero())) {
throw new SingularMatrixException('GaussJordanElimination');
}
$pivinv = $calc->reciprocal($dA[$icol][$icol]);
$this->multRow($dA, $icol, $pivinv, $calc);
$this->multRow($dB, $icol, $pivinv, $calc);
for ($ll = 0; $ll < $rows; $ll++) {
if ($ll != $icol) {
$multiplier = clone $dA[$ll][$icol];
$multiplier->negate();
$this->addMultipleOfOtherRowToRow($dA, $multiplier, $icol, $ll, $calc);
$this->addMultipleOfOtherRowToRow($dB, $multiplier, $icol, $ll, $calc);
}
}
}
$this->set('left', $this->createCorrectMatrixType($mA, $dA));
$this->set('right', $this->createCorrectMatrixType($extra, $dB));
return clone $this;
}
protected function swapRows(array &$a, $r1, $r2)
{
$tmp = $a[$r1];
$a[$r1] = $a[$r2];
$a[$r2] = $tmp;
}
protected function multRow(array &$a, $row, $num, Calculator $calc)
{
foreach ($a[$row] as &$value) {
$value = $calc->mul($value, $num);
}
}
protected function addMultipleOfOtherRowToRow(array &$a, $multiple, $rowToMultiplyWith,
$rowToAddTo, Calculator $calc)
{
$numberOfColumns = count($a[0]);
for ($l = 0; $l < $numberOfColumns; $l++) {
$a[$rowToAddTo][$l] = $calc->add(
$a[$rowToAddTo][$l],
$calc->mul(
$a[$rowToMultiplyWith][$l],
$multiple));
}
}
}