ILIAS  release_5-2 Revision v5.2.25-18-g3f80b828510
LUDecomposition.php
Go to the documentation of this file.
1<?php
22
23 const MatrixSingularException = "Can only perform operation on singular matrix.";
24 const MatrixSquareException = "Mismatched Row dimension";
25
30 private $LU = array();
31
36 private $m;
37
42 private $n;
43
48 private $pivsign;
49
54 private $piv = array();
55
56
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 {
119 }
120 } // function __construct()
121
122
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
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
168 public function getPivot() {
169 return $this->piv;
170 } // function getPivot()
171
172
178 public function getDoublePivot() {
179 return $this->getPivot();
180 } // function getDoublePivot()
181
182
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
203 public function det() {
204 if ($this->m == $this->n) {
206 for ($j = 0; $j < $this->n; ++$j) {
207 $d *= $this->LU[$j][$j];
208 }
209 return $d;
210 } else {
212 }
213 } // function det()
214
215
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 {
252 }
253 } else {
254 throw new PHPExcel_Calculation_Exception(self::MatrixSquareException);
255 }
256 } // function solve()
257
258} // class PHPExcel_Shared_JAMA_LUDecomposition
for($col=0; $col< 50; $col++) $d
const MatrixSingularException
Definition: Error.php:54
An exception for terminatinating execution or to throw for unit testing.
__construct($A)
LU Decomposition constructor.
getPivot()
Return pivot permutation vector.
getL()
Get lower triangular factor.
getU()
Get upper triangular factor.
isNonsingular()
Is the matrix nonsingular?
PHPExcel root directory.
Definition: Matrix.php:27