<?php

/**

 *	@package JAMA

 *

 *	Cholesky decomposition class

 *

 *	For a symmetric, positive definite matrix A, the Cholesky decomposition

 *	is an lower triangular matrix L so that A = L*L'.

 *

 *	If the matrix is not symmetric or positive definite, the constructor

 *	returns a partial decomposition and sets an internal flag that may

 *	be queried by the isSPD() method.

 *

 *	@author Paul Meagher

 *	@author Michael Bommarito

 *	@version 1.2

 */

class CholeskyDecomposition {



	/**

	 *	Decomposition storage

	 *	@var array

	 *	@access private

	 */

	private $L = array();



	/**

	 *	Matrix row and column dimension

	 *	@var int

	 *	@access private

	 */

	private $m;



	/**

	 *	Symmetric positive definite flag

	 *	@var boolean

	 *	@access private

	 */

	private $isspd = true;





	/**

	 *	CholeskyDecomposition

	 *

	 *	Class constructor - decomposes symmetric positive definite matrix

	 *	@param mixed Matrix square symmetric positive definite matrix

	 */

	public function __construct($A = null) {

		if ($A instanceof Matrix) {

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

			$this->m = $A->getRowDimension();



			for($i = 0; $i < $this->m; ++$i) {

				for($j = $i; $j < $this->m; ++$j) {

					for($sum = $this->L[$i][$j], $k = $i - 1; $k >= 0; --$k) {

						$sum -= $this->L[$i][$k] * $this->L[$j][$k];

					}

					if ($i == $j) {

						if ($sum >= 0) {

							$this->L[$i][$i] = sqrt($sum);

						} else {

							$this->isspd = false;

						}

					} else {

						if ($this->L[$i][$i] != 0) {

							$this->L[$j][$i] = $sum / $this->L[$i][$i];

						}

					}

				}



				for ($k = $i+1; $k < $this->m; ++$k) {

					$this->L[$i][$k] = 0.0;

				}

			}

		} else {

			throw new PHPExcel_Calculation_Exception(JAMAError(ArgumentTypeException));

		}

	}	//	function __construct()





	/**

	 *	Is the matrix symmetric and positive definite?

	 *

	 *	@return boolean

	 */

	public function isSPD() {

		return $this->isspd;

	}	//	function isSPD()





	/**

	 *	getL

	 *

	 *	Return triangular factor.

	 *	@return Matrix Lower triangular matrix

	 */

	public function getL() {

		return new Matrix($this->L);

	}	//	function getL()





	/**

	 *	Solve A*X = B

	 *

	 *	@param $B Row-equal matrix

	 *	@return Matrix L * L' * X = B

	 */

	public function solve($B = null) {

		if ($B instanceof Matrix) {

			if ($B->getRowDimension() == $this->m) {

				if ($this->isspd) {

					$X  = $B->getArrayCopy();

					$nx = $B->getColumnDimension();



					for ($k = 0; $k < $this->m; ++$k) {

						for ($i = $k + 1; $i < $this->m; ++$i) {

							for ($j = 0; $j < $nx; ++$j) {

								$X[$i][$j] -= $X[$k][$j] * $this->L[$i][$k];

							}

						}

						for ($j = 0; $j < $nx; ++$j) {

							$X[$k][$j] /= $this->L[$k][$k];

						}

					}



					for ($k = $this->m - 1; $k >= 0; --$k) {

						for ($j = 0; $j < $nx; ++$j) {

							$X[$k][$j] /= $this->L[$k][$k];

						}

						for ($i = 0; $i < $k; ++$i) {

							for ($j = 0; $j < $nx; ++$j) {

								$X[$i][$j] -= $X[$k][$j] * $this->L[$k][$i];

							}

						}

					}



					return new Matrix($X, $this->m, $nx);

				} else {

					throw new PHPExcel_Calculation_Exception(JAMAError(MatrixSPDException));

				}

			} else {

				throw new PHPExcel_Calculation_Exception(JAMAError(MatrixDimensionException));

			}

		} else {

			throw new PHPExcel_Calculation_Exception(JAMAError(ArgumentTypeException));

		}

	}	//	function solve()



}	//	class CholeskyDecomposition

