Init Repo

extend/phpexcel/PHPExcel/Shared/JAMA/CholeskyDecomposition.php

@@ -0,0 +1,149 @@
+<?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