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
/*
 * Math-Matrix library
 *
 * @author Ashley Kitson <akitson@zf4.biz>
 * @copyright Ashley Kitson, UK, 2014
 * @licence GPL V3 or later : http://www.gnu.org/licenses/gpl.html
 * @link http://en.wikipedia.org/wiki/Matrix_(mathematics)
 */
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;

/**
 * This is lifted from the JAMA package and adapted for this library from its
 * PHP4 origin into PHP 5 and to fit the Chippyash\Matrix model
 * You can find the JAMA package at http://www.phpmath.com/build02/JAMA/downloads/
 * My thanks to the original authors mentioned below
 *
 * 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 Lu extends AbstractDecomposition
{
    use CreateCorrectMatrixType;
    use CreateCorrectScalarType;

    /**
     * Products of the decomposition
     * @var array [productName => Matrix,...]
     */
    protected $products = array(
        'LU' => null, //NumericMatrix: full decomposition
        'L' => null,  //NumericMatrix: lower triangle
        'U' => null,   //NumericMatrix: upper triangle
        'PivotVector' => null, //NumericMatrix: Pivot vector of the decomposition,
        'PermutationMatrix' => null, //NumericMatrix: permutation matrix of decomposition,
        'Det' => null //NumericTypeInterface|null: determinant of source matrix if matrix is square else null
    );

    /**
     * Internal storage of pivot vector.
     * @var array
     */
    protected $piv = [];

    /**
     * Pivot sign - used for Det
     * @var int
     */
    protected $pivsign;

    /**
     * Number of rows in original matrix
     * @var int
     */
    protected $rows;
    /**
     * Number of columns in original matrix
     * @var int
     */
    protected $cols;

    /**
     * Do the decomposition
     *
     * @param \Chippyash\Math\Matrix\NumericMatrix $mA
     * @param mixed $extra ignored
     * @return \Chippyash\Matrix\Transformation\Decomposition\Lu
     */
    public function decompose(NumericMatrix $mA, $extra = null)
    {
        $this->LUDecomposition($mA);

        $this->setOtherProducts($mA);

        return clone $this;
    }

    /**
     * LU Decomposition constructor.
     *
     * @param \Chippyash\Math\Matrix\NumericMatrix $mA
     */
    protected function LUDecomposition(NumericMatrix $mA)
    {
        // Use a "left-looking", dot-product, Crout/Doolittle algorithm.
        $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);

        // Outer loop.
        for ($j = 0; $j < $n; $j++) {
            // Make a copy of the j-th column to localize references.
            for ($i = 0; $i < $m; $i++) {
                $LUcolj[$i] = &$LU[$i][$j];
            }
            // Apply previous transformations.
            for ($i = 0; $i < $m; $i++) {
                $LUrowi = $LU[$i];
                // Most of the time is spent in the following dot product.
                $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];
            }
            // Find pivot and exchange if necessary.
            $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++) {
                    //swap
                    $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;
            }
            // Compute multipliers.
            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));
    }

    /**
     * Set other products of the decomposition
     */
    protected function setOtherProducts(NumericMatrix $mA)
    {
        $this->setLowerProduct($mA);
        $this->setUpperProduct($mA);
        $this->setPivotVector($mA);
        $this->setPermutationMatrix($mA);
        $this->setDeterminant($mA);
    }

    /**
     * Set lower triangular factor.
     */
    protected function setLowerProduct(NumericMatrix $mA)
    {
        $m = $this->LU->rows();
        $n = $this->LU->columns();
        $LU = $this->LU->toArray();
        $rcFactor = $this->rows -$this->cols;
        //set Lower
        $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);
                            }
                        }
                    }

                    //remove extra cols for non square matrices
                    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);
                    }
                }
        );
    }

    /**
     * Set upper triangular factor.
     */
    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);
                            }
                        }
                    }
                    //remove extra rows for non square matrices
                    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);
                    }
                }
        );
    }

    /**
     * Set pivot permutation vector.
     */
    protected function setPivotVector(NumericMatrix $mA)
    {
        $mB = $this->createCorrectMatrixType($mA, [$this->piv]);
        $this->set(
                'PivotVector',
                $mB("Add\\Scalar", 1));
    }

    /**
     * Set permutation matrix
     */
    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);
                }
                );
    }

    /**
     * Set determinant of original matrix if it is square
     */
    protected function setDeterminant(NumericMatrix $mA)
    {
        if (!$mA->is('square')) {
            //determinant undefined for non square matrix
            $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;
                }
                );
    }
}