LUDecomposition.php

Go to the documentation of this file.

<?php
class LUDecomposition {

        private $LU = array();

        private $m;

        private $n;

        private $pivsign;

        private $piv = array();


        public function __construct($A) {
                if ($A instanceof Matrix) {
                        // Use a "left-looking", dot-product, Crout/Doolittle algorithm.
                        $this->LU = $A->getArrayCopy();
                        $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 Exception(JAMAError(ArgumentTypeException));
                }
        }       //      function __construct()


        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 Matrix($L);
        }       //      function getL()


        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 Matrix($U);
        }       //      function getU()


        public function getPivot() {
                return $this->piv;
        }       //      function getPivot()


        public function getDoublePivot() {
                return $this->getPivot();
        }       //      function getDoublePivot()


        public function isNonsingular() {
                for ($j = 0; $j < $this->n; ++$j) {
                        if ($this->LU[$j][$j] == 0) {
                                return false;
                        }
                }
                return true;
        }       //      function isNonsingular()


        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 Exception(JAMAError(MatrixDimensionException));
                }
        }       //      function det()


        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 Exception(JAMAError(MatrixSingularException));
                        }
                } else {
                        throw new Exception(JAMAError(MatrixSquareException));
                }
        }       //      function solve()

}       //      class LUDecomposition