ILIAS  release_5-4 Revision v5.4.26-12-gabc799a52e6
LUDecomposition.php
Go to the documentation of this file.
1<?php
2
4
5use PhpOffice\PhpSpreadsheet\Calculation\Exception as CalculationException;
6
25{
26 const MATRIX_SINGULAR_EXCEPTION = 'Can only perform operation on singular matrix.';
27 const MATRIX_SQUARE_EXCEPTION = 'Mismatched Row dimension';
28
34 private $LU = [];
35
41 private $m;
42
48 private $n;
49
55 private $pivsign;
56
62 private $piv = [];
63
69 public function __construct($A)
70 {
71 if ($A instanceof Matrix) {
72 // Use a "left-looking", dot-product, Crout/Doolittle algorithm.
73 $this->LU = $A->getArray();
74 $this->m = $A->getRowDimension();
75 $this->n = $A->getColumnDimension();
76 for ($i = 0; $i < $this->m; ++$i) {
77 $this->piv[$i] = $i;
78 }
79 $this->pivsign = 1;
80 $LUrowi = $LUcolj = [];
81
82 // Outer loop.
83 for ($j = 0; $j < $this->n; ++$j) {
84 // Make a copy of the j-th column to localize references.
85 for ($i = 0; $i < $this->m; ++$i) {
86 $LUcolj[$i] = &$this->LU[$i][$j];
87 }
88 // Apply previous transformations.
89 for ($i = 0; $i < $this->m; ++$i) {
90 $LUrowi = $this->LU[$i];
91 // Most of the time is spent in the following dot product.
92 $kmax = min($i, $j);
93 $s = 0.0;
94 for ($k = 0; $k < $kmax; ++$k) {
95 $s += $LUrowi[$k] * $LUcolj[$k];
96 }
97 $LUrowi[$j] = $LUcolj[$i] -= $s;
98 }
99 // Find pivot and exchange if necessary.
100 $p = $j;
101 for ($i = $j + 1; $i < $this->m; ++$i) {
102 if (abs($LUcolj[$i]) > abs($LUcolj[$p])) {
103 $p = $i;
104 }
105 }
106 if ($p != $j) {
107 for ($k = 0; $k < $this->n; ++$k) {
108 $t = $this->LU[$p][$k];
109 $this->LU[$p][$k] = $this->LU[$j][$k];
110 $this->LU[$j][$k] = $t;
111 }
112 $k = $this->piv[$p];
113 $this->piv[$p] = $this->piv[$j];
114 $this->piv[$j] = $k;
115 $this->pivsign = $this->pivsign * -1;
116 }
117 // Compute multipliers.
118 if (($j < $this->m) && ($this->LU[$j][$j] != 0.0)) {
119 for ($i = $j + 1; $i < $this->m; ++$i) {
120 $this->LU[$i][$j] /= $this->LU[$j][$j];
121 }
122 }
123 }
124 } else {
125 throw new CalculationException(Matrix::ARGUMENT_TYPE_EXCEPTION);
126 }
127 }
128
129 // function __construct()
130
136 public function getL()
137 {
138 $L = [];
139 for ($i = 0; $i < $this->m; ++$i) {
140 for ($j = 0; $j < $this->n; ++$j) {
141 if ($i > $j) {
142 $L[$i][$j] = $this->LU[$i][$j];
143 } elseif ($i == $j) {
144 $L[$i][$j] = 1.0;
145 } else {
146 $L[$i][$j] = 0.0;
147 }
148 }
149 }
150
151 return new Matrix($L);
152 }
153
154 // function getL()
155
161 public function getU()
162 {
163 $U = [];
164 for ($i = 0; $i < $this->n; ++$i) {
165 for ($j = 0; $j < $this->n; ++$j) {
166 if ($i <= $j) {
167 $U[$i][$j] = $this->LU[$i][$j];
168 } else {
169 $U[$i][$j] = 0.0;
170 }
171 }
172 }
173
174 return new Matrix($U);
175 }
176
177 // function getU()
178
184 public function getPivot()
185 {
186 return $this->piv;
187 }
188
189 // function getPivot()
190
196 public function getDoublePivot()
197 {
198 return $this->getPivot();
199 }
200
201 // function getDoublePivot()
202
208 public function isNonsingular()
209 {
210 for ($j = 0; $j < $this->n; ++$j) {
211 if ($this->LU[$j][$j] == 0) {
212 return false;
213 }
214 }
215
216 return true;
217 }
218
219 // function isNonsingular()
220
226 public function det()
227 {
228 if ($this->m == $this->n) {
230 for ($j = 0; $j < $this->n; ++$j) {
231 $d *= $this->LU[$j][$j];
232 }
233
234 return $d;
235 }
236
237 throw new CalculationException(Matrix::MATRIX_DIMENSION_EXCEPTION);
238 }
239
240 // function det()
241
249 public function solve(Matrix $B)
250 {
251 if ($B->getRowDimension() == $this->m) {
252 if ($this->isNonsingular()) {
253 // Copy right hand side with pivoting
254 $nx = $B->getColumnDimension();
255 $X = $B->getMatrix($this->piv, 0, $nx - 1);
256 // Solve L*Y = B(piv,:)
257 for ($k = 0; $k < $this->n; ++$k) {
258 for ($i = $k + 1; $i < $this->n; ++$i) {
259 for ($j = 0; $j < $nx; ++$j) {
260 $X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k];
261 }
262 }
263 }
264 // Solve U*X = Y;
265 for ($k = $this->n - 1; $k >= 0; --$k) {
266 for ($j = 0; $j < $nx; ++$j) {
267 $X->A[$k][$j] /= $this->LU[$k][$k];
268 }
269 for ($i = 0; $i < $k; ++$i) {
270 for ($j = 0; $j < $nx; ++$j) {
271 $X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k];
272 }
273 }
274 }
275
276 return $X;
277 }
278
279 throw new CalculationException(self::MATRIX_SINGULAR_EXCEPTION);
280 }
281
282 throw new CalculationException(self::MATRIX_SQUARE_EXCEPTION);
283 }
284}
An exception for terminatinating execution or to throw for unit testing.
For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n unit lower triangular matrix L,...
getColumnDimension()
getColumnDimension.
Definition: Matrix.php:137
for( $i=6;$i< 13;$i++) for($i=1; $i< 13; $i++) $d
Definition: date.php:296
$i
Definition: disco.tpl.php:19
$X
Definition: test.php:23
Class for the creating "special" Matrices.
Definition: Builder.php:11
$s
Definition: pwgen.php:45