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 */
21class 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
259