<?php

/**

 *	@package JAMA

 *

 *	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 PHPExcel_Shared_JAMA_LUDecomposition {



	const MatrixSingularException	= "Can only perform operation on singular matrix.";

	const MatrixSquareException		= "Mismatched Row dimension";



	/**

	 *	Decomposition storage

	 *	@var array

	 */

	private $LU = array();



	/**

	 *	Row dimension.

	 *	@var int

	 */

	private $m;



	/**

	 *	Column dimension.

	 *	@var int

	 */

	private $n;



	/**

	 *	Pivot sign.

	 *	@var int

	 */

	private $pivsign;



	/**

	 *	Internal storage of pivot vector.

	 *	@var array

	 */

	private $piv = array();





	/**

	 *	LU Decomposition constructor.

	 *

	 *	@param $A Rectangular matrix

	 *	@return Structure to access L, U and piv.

	 */

	public function __construct($A) {

		if ($A instanceof PHPExcel_Shared_JAMA_Matrix) {

			// Use a "left-looking", dot-product, Crout/Doolittle algorithm.

			$this->LU = $A->getArray();

			$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 PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::ArgumentTypeException);

		}

	}	//	function __construct()





	/**

	 *	Get lower triangular factor.

	 *

	 *	@return array Lower triangular factor

	 */

	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 PHPExcel_Shared_JAMA_Matrix($L);

	}	//	function getL()





	/**

	 *	Get upper triangular factor.

	 *

	 *	@return array Upper triangular factor

	 */

	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 PHPExcel_Shared_JAMA_Matrix($U);

	}	//	function getU()





	/**

	 *	Return pivot permutation vector.

	 *

	 *	@return array Pivot vector

	 */

	public function getPivot() {

		return $this->piv;

	}	//	function getPivot()





	/**

	 *	Alias for getPivot

	 *

	 *	@see getPivot

	 */

	public function getDoublePivot() {

		return $this->getPivot();

	}	//	function getDoublePivot()





	/**

	 *	Is the matrix nonsingular?

	 *

	 *	@return true if U, and hence A, is nonsingular.

	 */

	public function isNonsingular() {

		for ($j = 0; $j < $this->n; ++$j) {

			if ($this->LU[$j][$j] == 0) {

				return false;

			}

		}

		return true;

	}	//	function isNonsingular()





	/**

	 *	Count determinants

	 *

	 *	@return array d matrix deterninat

	 */

	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 PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::MatrixDimensionException);

		}

	}	//	function det()





	/**

	 *	Solve A*X = B

	 *

	 *	@param  $B  A Matrix with as many rows as A and any number of columns.

	 *	@return  X so that L*U*X = B(piv,:)

	 *	@PHPExcel_Calculation_Exception  IllegalArgumentException Matrix row dimensions must agree.

	 *	@PHPExcel_Calculation_Exception  RuntimeException  Matrix is singular.

	 */

	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 PHPExcel_Calculation_Exception(self::MatrixSingularException);

			}

		} else {

			throw new PHPExcel_Calculation_Exception(self::MatrixSquareException);

		}

	}	//	function solve()



}	//	class PHPExcel_Shared_JAMA_LUDecomposition

