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