Subversion Repositories eFlore/Applications.cel

Rev

Go to most recent revision | Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
2388 jpm 1
<?php
2
/**
3
 *	@package JAMA
4
 *
5
 *	For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n
6
 *	unit lower triangular matrix L, an n-by-n upper triangular matrix U,
7
 *	and a permutation vector piv of length m so that A(piv,:) = L*U.
8
 *	If m < n, then L is m-by-m and U is m-by-n.
9
 *
10
 *	The LU decompostion with pivoting always exists, even if the matrix is
11
 *	singular, so the constructor will never fail. The primary use of the
12
 *	LU decomposition is in the solution of square systems of simultaneous
13
 *	linear equations. This will fail if isNonsingular() returns false.
14
 *
15
 *	@author Paul Meagher
16
 *	@author Bartosz Matosiuk
17
 *	@author Michael Bommarito
18
 *	@version 1.1
19
 *	@license PHP v3.0
20
 */
21
class PHPExcel_Shared_JAMA_LUDecomposition {
22
 
23
	const MatrixSingularException	= "Can only perform operation on singular matrix.";
24
	const MatrixSquareException		= "Mismatched Row dimension";
25
 
26
	/**
27
	 *	Decomposition storage
28
	 *	@var array
29
	 */
30
	private $LU = array();
31
 
32
	/**
33
	 *	Row dimension.
34
	 *	@var int
35
	 */
36
	private $m;
37
 
38
	/**
39
	 *	Column dimension.
40
	 *	@var int
41
	 */
42
	private $n;
43
 
44
	/**
45
	 *	Pivot sign.
46
	 *	@var int
47
	 */
48
	private $pivsign;
49
 
50
	/**
51
	 *	Internal storage of pivot vector.
52
	 *	@var array
53
	 */
54
	private $piv = array();
55
 
56
 
57
	/**
58
	 *	LU Decomposition constructor.
59
	 *
60
	 *	@param $A Rectangular matrix
61
	 *	@return Structure to access L, U and piv.
62
	 */
63
	public function __construct($A) {
64
		if ($A instanceof PHPExcel_Shared_JAMA_Matrix) {
65
			// Use a "left-looking", dot-product, Crout/Doolittle algorithm.
66
			$this->LU = $A->getArray();
67
			$this->m  = $A->getRowDimension();
68
			$this->n  = $A->getColumnDimension();
69
			for ($i = 0; $i < $this->m; ++$i) {
70
				$this->piv[$i] = $i;
71
			}
72
			$this->pivsign = 1;
73
			$LUrowi = $LUcolj = array();
74
 
75
			// Outer loop.
76
			for ($j = 0; $j < $this->n; ++$j) {
77
				// Make a copy of the j-th column to localize references.
78
				for ($i = 0; $i < $this->m; ++$i) {
79
					$LUcolj[$i] = &$this->LU[$i][$j];
80
				}
81
				// Apply previous transformations.
82
				for ($i = 0; $i < $this->m; ++$i) {
83
					$LUrowi = $this->LU[$i];
84
					// Most of the time is spent in the following dot product.
85
					$kmax = min($i,$j);
86
					$s = 0.0;
87
					for ($k = 0; $k < $kmax; ++$k) {
88
						$s += $LUrowi[$k] * $LUcolj[$k];
89
					}
90
					$LUrowi[$j] = $LUcolj[$i] -= $s;
91
				}
92
				// Find pivot and exchange if necessary.
93
				$p = $j;
94
				for ($i = $j+1; $i < $this->m; ++$i) {
95
					if (abs($LUcolj[$i]) > abs($LUcolj[$p])) {
96
						$p = $i;
97
					}
98
				}
99
				if ($p != $j) {
100
					for ($k = 0; $k < $this->n; ++$k) {
101
						$t = $this->LU[$p][$k];
102
						$this->LU[$p][$k] = $this->LU[$j][$k];
103
						$this->LU[$j][$k] = $t;
104
					}
105
					$k = $this->piv[$p];
106
					$this->piv[$p] = $this->piv[$j];
107
					$this->piv[$j] = $k;
108
					$this->pivsign = $this->pivsign * -1;
109
				}
110
				// Compute multipliers.
111
				if (($j < $this->m) && ($this->LU[$j][$j] != 0.0)) {
112
					for ($i = $j+1; $i < $this->m; ++$i) {
113
						$this->LU[$i][$j] /= $this->LU[$j][$j];
114
					}
115
				}
116
			}
117
		} else {
118
			throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::ArgumentTypeException);
119
		}
120
	}	//	function __construct()
121
 
122
 
123
	/**
124
	 *	Get lower triangular factor.
125
	 *
126
	 *	@return array Lower triangular factor
127
	 */
128
	public function getL() {
129
		for ($i = 0; $i < $this->m; ++$i) {
130
			for ($j = 0; $j < $this->n; ++$j) {
131
				if ($i > $j) {
132
					$L[$i][$j] = $this->LU[$i][$j];
133
				} elseif ($i == $j) {
134
					$L[$i][$j] = 1.0;
135
				} else {
136
					$L[$i][$j] = 0.0;
137
				}
138
			}
139
		}
140
		return new PHPExcel_Shared_JAMA_Matrix($L);
141
	}	//	function getL()
142
 
143
 
144
	/**
145
	 *	Get upper triangular factor.
146
	 *
147
	 *	@return array Upper triangular factor
148
	 */
149
	public function getU() {
150
		for ($i = 0; $i < $this->n; ++$i) {
151
			for ($j = 0; $j < $this->n; ++$j) {
152
				if ($i <= $j) {
153
					$U[$i][$j] = $this->LU[$i][$j];
154
				} else {
155
					$U[$i][$j] = 0.0;
156
				}
157
			}
158
		}
159
		return new PHPExcel_Shared_JAMA_Matrix($U);
160
	}	//	function getU()
161
 
162
 
163
	/**
164
	 *	Return pivot permutation vector.
165
	 *
166
	 *	@return array Pivot vector
167
	 */
168
	public function getPivot() {
169
		return $this->piv;
170
	}	//	function getPivot()
171
 
172
 
173
	/**
174
	 *	Alias for getPivot
175
	 *
176
	 *	@see getPivot
177
	 */
178
	public function getDoublePivot() {
179
		return $this->getPivot();
180
	}	//	function getDoublePivot()
181
 
182
 
183
	/**
184
	 *	Is the matrix nonsingular?
185
	 *
186
	 *	@return true if U, and hence A, is nonsingular.
187
	 */
188
	public function isNonsingular() {
189
		for ($j = 0; $j < $this->n; ++$j) {
190
			if ($this->LU[$j][$j] == 0) {
191
				return false;
192
			}
193
		}
194
		return true;
195
	}	//	function isNonsingular()
196
 
197
 
198
	/**
199
	 *	Count determinants
200
	 *
201
	 *	@return array d matrix deterninat
202
	 */
203
	public function det() {
204
		if ($this->m == $this->n) {
205
			$d = $this->pivsign;
206
			for ($j = 0; $j < $this->n; ++$j) {
207
				$d *= $this->LU[$j][$j];
208
			}
209
			return $d;
210
		} else {
211
			throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::MatrixDimensionException);
212
		}
213
	}	//	function det()
214
 
215
 
216
	/**
217
	 *	Solve A*X = B
218
	 *
219
	 *	@param  $B  A Matrix with as many rows as A and any number of columns.
220
	 *	@return  X so that L*U*X = B(piv,:)
221
	 *	@PHPExcel_Calculation_Exception  IllegalArgumentException Matrix row dimensions must agree.
222
	 *	@PHPExcel_Calculation_Exception  RuntimeException  Matrix is singular.
223
	 */
224
	public function solve($B) {
225
		if ($B->getRowDimension() == $this->m) {
226
			if ($this->isNonsingular()) {
227
				// Copy right hand side with pivoting
228
				$nx = $B->getColumnDimension();
229
				$X  = $B->getMatrix($this->piv, 0, $nx-1);
230
				// Solve L*Y = B(piv,:)
231
				for ($k = 0; $k < $this->n; ++$k) {
232
					for ($i = $k+1; $i < $this->n; ++$i) {
233
						for ($j = 0; $j < $nx; ++$j) {
234
							$X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k];
235
						}
236
					}
237
				}
238
				// Solve U*X = Y;
239
				for ($k = $this->n-1; $k >= 0; --$k) {
240
					for ($j = 0; $j < $nx; ++$j) {
241
						$X->A[$k][$j] /= $this->LU[$k][$k];
242
					}
243
					for ($i = 0; $i < $k; ++$i) {
244
						for ($j = 0; $j < $nx; ++$j) {
245
							$X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k];
246
						}
247
					}
248
				}
249
				return $X;
250
			} else {
251
				throw new PHPExcel_Calculation_Exception(self::MatrixSingularException);
252
			}
253
		} else {
254
			throw new PHPExcel_Calculation_Exception(self::MatrixSquareException);
255
		}
256
	}	//	function solve()
257
 
258
}	//	class PHPExcel_Shared_JAMA_LUDecomposition