Subversion Repositories eFlore/Applications.cel

Rev

Rev 2388 | Details | Compare with Previous | 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