ILIAS  release_5-4 Revision v5.4.26-12-gabc799a52e6
LU.php
Go to the documentation of this file.
1<?php
2
4
7
8class LU
9{
10 private $luMatrix;
11 private $rows;
12 private $columns;
13
14 private $pivot = [];
15
16 public function __construct(Matrix $matrix)
17 {
18 $this->luMatrix = $matrix->toArray();
19 $this->rows = $matrix->rows;
20 $this->columns = $matrix->columns;
21
22 $this->buildPivot();
23 }
24
30 public function getL(): Matrix
31 {
32 $lower = [];
33
34 $columns = min($this->rows, $this->columns);
35 for ($row = 0; $row < $this->rows; ++$row) {
36 for ($column = 0; $column < $columns; ++$column) {
37 if ($row > $column) {
38 $lower[$row][$column] = $this->luMatrix[$row][$column];
39 } elseif ($row === $column) {
40 $lower[$row][$column] = 1.0;
41 } else {
42 $lower[$row][$column] = 0.0;
43 }
44 }
45 }
46
47 return new Matrix($lower);
48 }
49
55 public function getU(): Matrix
56 {
57 $upper = [];
58
59 $rows = min($this->rows, $this->columns);
60 for ($row = 0; $row < $rows; ++$row) {
61 for ($column = 0; $column < $this->columns; ++$column) {
62 if ($row <= $column) {
63 $upper[$row][$column] = $this->luMatrix[$row][$column];
64 } else {
65 $upper[$row][$column] = 0.0;
66 }
67 }
68 }
69
70 return new Matrix($upper);
71 }
72
78 public function getP(): Matrix
79 {
80 $pMatrix = [];
81
82 $pivots = $this->pivot;
83 $pivotCount = count($pivots);
84 foreach ($pivots as $row => $pivot) {
85 $pMatrix[$row] = array_fill(0, $pivotCount, 0);
86 $pMatrix[$row][$pivot] = 1;
87 }
88
89 return new Matrix($pMatrix);
90 }
91
97 public function getPivot(): array
98 {
99 return $this->pivot;
100 }
101
107 public function isNonsingular(): bool
108 {
109 for ($diagonal = 0; $diagonal < $this->columns; ++$diagonal) {
110 if ($this->luMatrix[$diagonal][$diagonal] === 0.0) {
111 return false;
112 }
113 }
114
115 return true;
116 }
117
118 private function buildPivot(): void
119 {
120 for ($row = 0; $row < $this->rows; ++$row) {
121 $this->pivot[$row] = $row;
122 }
123
124 for ($column = 0; $column < $this->columns; ++$column) {
125 $luColumn = $this->localisedReferenceColumn($column);
126
127 $this->applyTransformations($column, $luColumn);
128
129 $pivot = $this->findPivot($column, $luColumn);
130 if ($pivot !== $column) {
131 $this->pivotExchange($pivot, $column);
132 }
133
134 $this->computeMultipliers($column);
135
136 unset($luColumn);
137 }
138 }
139
140 private function localisedReferenceColumn($column): array
141 {
142 $luColumn = [];
143
144 for ($row = 0; $row < $this->rows; ++$row) {
145 $luColumn[$row] = &$this->luMatrix[$row][$column];
146 }
147
148 return $luColumn;
149 }
150
151 private function applyTransformations($column, array $luColumn): void
152 {
153 for ($row = 0; $row < $this->rows; ++$row) {
154 $luRow = $this->luMatrix[$row];
155 // Most of the time is spent in the following dot product.
156 $kmax = min($row, $column);
157 $sValue = 0.0;
158 for ($kValue = 0; $kValue < $kmax; ++$kValue) {
159 $sValue += $luRow[$kValue] * $luColumn[$kValue];
160 }
161 $luRow[$column] = $luColumn[$row] -= $sValue;
162 }
163 }
164
165 private function findPivot($column, array $luColumn): int
166 {
167 $pivot = $column;
168 for ($row = $column + 1; $row < $this->rows; ++$row) {
169 if (abs($luColumn[$row]) > abs($luColumn[$pivot])) {
170 $pivot = $row;
171 }
172 }
173
174 return $pivot;
175 }
176
177 private function pivotExchange($pivot, $column): void
178 {
179 for ($kValue = 0; $kValue < $this->columns; ++$kValue) {
180 $tValue = $this->luMatrix[$pivot][$kValue];
181 $this->luMatrix[$pivot][$kValue] = $this->luMatrix[$column][$kValue];
182 $this->luMatrix[$column][$kValue] = $tValue;
183 }
184
185 $lValue = $this->pivot[$pivot];
186 $this->pivot[$pivot] = $this->pivot[$column];
187 $this->pivot[$column] = $lValue;
188 }
189
190 private function computeMultipliers($diagonal): void
191 {
192 if (($diagonal < $this->rows) && ($this->luMatrix[$diagonal][$diagonal] != 0.0)) {
193 for ($row = $diagonal + 1; $row < $this->rows; ++$row) {
194 $this->luMatrix[$row][$diagonal] /= $this->luMatrix[$diagonal][$diagonal];
195 }
196 }
197 }
198
199 private function pivotB(Matrix $B): array
200 {
201 $X = [];
202 foreach ($this->pivot as $rowId) {
203 $row = $B->getRows($rowId + 1)->toArray();
204 $X[] = array_pop($row);
205 }
206
207 return $X;
208 }
209
219 public function solve(Matrix $B): Matrix
220 {
221 if ($B->rows !== $this->rows) {
222 throw new Exception('Matrix row dimensions are not equal');
223 }
224
225 if ($this->rows !== $this->columns) {
226 throw new Exception('LU solve() only works on square matrices');
227 }
228
229 if (!$this->isNonsingular()) {
230 throw new Exception('Can only perform operation on singular matrix');
231 }
232
233 // Copy right hand side with pivoting
234 $nx = $B->columns;
235 $X = $this->pivotB($B);
236
237 // Solve L*Y = B(piv,:)
238 for ($k = 0; $k < $this->columns; ++$k) {
239 for ($i = $k + 1; $i < $this->columns; ++$i) {
240 for ($j = 0; $j < $nx; ++$j) {
241 $X[$i][$j] -= $X[$k][$j] * $this->luMatrix[$i][$k];
242 }
243 }
244 }
245
246 // Solve U*X = Y;
247 for ($k = $this->columns - 1; $k >= 0; --$k) {
248 for ($j = 0; $j < $nx; ++$j) {
249 $X[$k][$j] /= $this->luMatrix[$k][$k];
250 }
251 for ($i = 0; $i < $k; ++$i) {
252 for ($j = 0; $j < $nx; ++$j) {
253 $X[$i][$j] -= $X[$k][$j] * $this->luMatrix[$i][$k];
254 }
255 }
256 }
257
258 return new Matrix($X);
259 }
260}
An exception for terminatinating execution or to throw for unit testing.
pivotB(Matrix $B)
Definition: LU.php:199
getU()
Get upper triangular factor.
Definition: LU.php:55
applyTransformations($column, array $luColumn)
Definition: LU.php:151
localisedReferenceColumn($column)
Definition: LU.php:140
getP()
Return pivot permutation vector.
Definition: LU.php:78
computeMultipliers($diagonal)
Definition: LU.php:190
__construct(Matrix $matrix)
Definition: LU.php:16
isNonsingular()
Is the matrix nonsingular?
Definition: LU.php:107
getL()
Get lower triangular factor.
Definition: LU.php:30
pivotExchange($pivot, $column)
Definition: LU.php:177
getPivot()
Return pivot permutation vector.
Definition: LU.php:97
solve(Matrix $B)
Solve A*X = B.
Definition: LU.php:219
findPivot($column, array $luColumn)
Definition: LU.php:165
getRows(int $row, int $rowCount=1)
Return a new matrix as a subset of rows from this matrix, starting at row number $row,...
Definition: Matrix.php:165
rows()
Returns a Generator that will yield each row of the matrix in turn as a vector matrix or the value of...
Definition: Matrix.php:282
columns()
Returns a Generator that will yield each column of the matrix in turn as a vector matrix or the value...
Definition: Matrix.php:297
$i
Definition: disco.tpl.php:19
$matrix
Definition: test.php:18
$X
Definition: test.php:23
$row
Class for the creating "special" Matrices.
Definition: Builder.php:11