Subversion Repositories eFlore/Applications.cel

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
2388 jpm 1 
<?php
2 
/**
3 
 *	@package JAMA
4 
 *
5 
 *	Cholesky decomposition class
6 
 *
7 
 *	For a symmetric, positive definite matrix A, the Cholesky decomposition
8 
 *	is an lower triangular matrix L so that A = L*L'.
9 
 *
10 
 *	If the matrix is not symmetric or positive definite, the constructor
11 
 *	returns a partial decomposition and sets an internal flag that may
12 
 *	be queried by the isSPD() method.
13 
 *
14 
 *	@author Paul Meagher
15 
 *	@author Michael Bommarito
16 
 *	@version 1.2
17 
 */
18 
class CholeskyDecomposition {
19 
 
20 
	/**
21 
	 *	Decomposition storage
22 
	 *	@var array
23 
	 *	@access private
24 
	 */
25 
	private $L = array();
26 
 
27 
	/**
28 
	 *	Matrix row and column dimension
29 
	 *	@var int
30 
	 *	@access private
31 
	 */
32 
	private $m;
33 
 
34 
	/**
35 
	 *	Symmetric positive definite flag
36 
	 *	@var boolean
37 
	 *	@access private
38 
	 */
39 
	private $isspd = true;
40 
 
41 
 
42 
	/**
43 
	 *	CholeskyDecomposition
44 
	 *
45 
	 *	Class constructor - decomposes symmetric positive definite matrix
46 
	 *	@param mixed Matrix square symmetric positive definite matrix
47 
	 */
48 
	public function __construct($A = null) {
49 
		if ($A instanceof Matrix) {
50 
			$this->L = $A->getArray();
51 
			$this->m = $A->getRowDimension();
52 
 
53 
			for($i = 0; $i < $this->m; ++$i) {
54 
				for($j = $i; $j < $this->m; ++$j) {
55 
					for($sum = $this->L[$i][$j], $k = $i - 1; $k >= 0; --$k) {
56 
						$sum -= $this->L[$i][$k] * $this->L[$j][$k];
57 
					}
58 
					if ($i == $j) {
59 
						if ($sum >= 0) {
60 
							$this->L[$i][$i] = sqrt($sum);
61 
						} else {
62 
							$this->isspd = false;
63 
						}
64 
					} else {
65 
						if ($this->L[$i][$i] != 0) {
66 
							$this->L[$j][$i] = $sum / $this->L[$i][$i];
67 
						}
68 
					}
69 
				}
70 
 
71 
				for ($k = $i+1; $k < $this->m; ++$k) {
72 
					$this->L[$i][$k] = 0.0;
73 
				}
74 
			}
75 
		} else {
76 
			throw new PHPExcel_Calculation_Exception(JAMAError(ArgumentTypeException));
77 
		}
78 
	}	//	function __construct()
79 
 
80 
 
81 
	/**
82 
	 *	Is the matrix symmetric and positive definite?
83 
	 *
84 
	 *	@return boolean
85 
	 */
86 
	public function isSPD() {
87 
		return $this->isspd;
88 
	}	//	function isSPD()
89 
 
90 
 
91 
	/**
92 
	 *	getL
93 
	 *
94 
	 *	Return triangular factor.
95 
	 *	@return Matrix Lower triangular matrix
96 
	 */
97 
	public function getL() {
98 
		return new Matrix($this->L);
99 
	}	//	function getL()
100 
 
101 
 
102 
	/**
103 
	 *	Solve A*X = B
104 
	 *
105 
	 *	@param $B Row-equal matrix
106 
	 *	@return Matrix L * L' * X = B
107 
	 */
108 
	public function solve($B = null) {
109 
		if ($B instanceof Matrix) {
110 
			if ($B->getRowDimension() == $this->m) {
111 
				if ($this->isspd) {
112 
					$X  = $B->getArrayCopy();
113 
					$nx = $B->getColumnDimension();
114 
 
115 
					for ($k = 0; $k < $this->m; ++$k) {
116 
						for ($i = $k + 1; $i < $this->m; ++$i) {
117 
							for ($j = 0; $j < $nx; ++$j) {
118 
								$X[$i][$j] -= $X[$k][$j] * $this->L[$i][$k];
119 
							}
120 
						}
121 
						for ($j = 0; $j < $nx; ++$j) {
122 
							$X[$k][$j] /= $this->L[$k][$k];
123 
						}
124 
					}
125 
 
126 
					for ($k = $this->m - 1; $k >= 0; --$k) {
127 
						for ($j = 0; $j < $nx; ++$j) {
128 
							$X[$k][$j] /= $this->L[$k][$k];
129 
						}
130 
						for ($i = 0; $i < $k; ++$i) {
131 
							for ($j = 0; $j < $nx; ++$j) {
132 
								$X[$i][$j] -= $X[$k][$j] * $this->L[$k][$i];
133 
							}
134 
						}
135 
					}
136 
 
137 
					return new Matrix($X, $this->m, $nx);
138 
				} else {
139 
					throw new PHPExcel_Calculation_Exception(JAMAError(MatrixSPDException));
140 
				}
141 
			} else {
142 
				throw new PHPExcel_Calculation_Exception(JAMAError(MatrixDimensionException));
143 
			}
144 
		} else {
145 
			throw new PHPExcel_Calculation_Exception(JAMAError(ArgumentTypeException));
146 
		}
147 
	}	//	function solve()
148 
 
149 
}	//	class CholeskyDecomposition