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
259: