Code Coverage |
||||||||||
Lines |
Functions and Methods |
Classes and Traits |
||||||||
| Total | |
88.46% |
69 / 78 |
|
85.71% |
6 / 7 |
CRAP | |
0.00% |
0 / 1 |
| LUDecomposition | |
88.46% |
69 / 78 |
|
85.71% |
6 / 7 |
40.22 | |
0.00% |
0 / 1 |
| __construct | |
71.88% |
23 / 32 |
|
0.00% |
0 / 1 |
16.76 | |||
| getL | |
100.00% |
11 / 11 |
|
100.00% |
1 / 1 |
5 | |||
| getU | |
100.00% |
7 / 7 |
|
100.00% |
1 / 1 |
4 | |||
| getPivot | |
100.00% |
1 / 1 |
|
100.00% |
1 / 1 |
1 | |||
| isNonsingular | |
100.00% |
4 / 4 |
|
100.00% |
1 / 1 |
3 | |||
| det | |
100.00% |
4 / 4 |
|
100.00% |
1 / 1 |
2 | |||
| solve | |
100.00% |
19 / 19 |
|
100.00% |
1 / 1 |
10 | |||
| 1 | <?php |
| 2 | /** |
| 3 | * Jingga |
| 4 | * |
| 5 | * PHP Version 8.1 |
| 6 | * |
| 7 | * @package phpOMS\Math\Matrix |
| 8 | * @copyright Dennis Eichhorn |
| 9 | * @copyright JAMA - https://math.nist.gov/javanumerics/jama/ |
| 10 | * @license OMS License 2.0 |
| 11 | * @version 1.0.0 |
| 12 | * @link https://jingga.app |
| 13 | */ |
| 14 | declare(strict_types=1); |
| 15 | |
| 16 | namespace phpOMS\Math\Matrix; |
| 17 | |
| 18 | use phpOMS\Math\Matrix\Exception\InvalidDimensionException; |
| 19 | |
| 20 | /** |
| 21 | * LU decomposition |
| 22 | * |
| 23 | * A(piv,:) = L*U |
| 24 | * |
| 25 | * @package phpOMS\Math\Matrix |
| 26 | * @license OMS License 2.0 |
| 27 | * @link https://jingga.app |
| 28 | * @since 1.0.0 |
| 29 | */ |
| 30 | final class LUDecomposition |
| 31 | { |
| 32 | /** |
| 33 | * LU matrix. |
| 34 | * |
| 35 | * @var array |
| 36 | * @since 1.0.0 |
| 37 | */ |
| 38 | private array $LU = []; |
| 39 | |
| 40 | /** |
| 41 | * Dimension m |
| 42 | * |
| 43 | * @var int |
| 44 | * @since 1.0.0 |
| 45 | */ |
| 46 | private int $m = 0; |
| 47 | |
| 48 | /** |
| 49 | * Dimension n |
| 50 | * |
| 51 | * @var int |
| 52 | * @since 1.0.0 |
| 53 | */ |
| 54 | private int $n = 0; |
| 55 | |
| 56 | /** |
| 57 | * Pivot sign |
| 58 | * |
| 59 | * @var int |
| 60 | * @since 1.0.0 |
| 61 | */ |
| 62 | private int $pivSign = 1; |
| 63 | |
| 64 | /** |
| 65 | * Pivot |
| 66 | * |
| 67 | * @var array |
| 68 | * @since 1.0.0 |
| 69 | */ |
| 70 | private array $piv = []; |
| 71 | |
| 72 | /** |
| 73 | * Constructor. |
| 74 | * |
| 75 | * @param Matrix $M Matrix |
| 76 | * |
| 77 | * @since 1.0.0 |
| 78 | */ |
| 79 | public function __construct(Matrix $M) |
| 80 | { |
| 81 | $this->LU = $M->toArray(); |
| 82 | $this->m = $M->getM(); |
| 83 | $this->n = $M->getN(); |
| 84 | |
| 85 | for ($i = 0; $i < $this->m; ++$i) { |
| 86 | $this->piv[$i] = $i; |
| 87 | } |
| 88 | |
| 89 | $LUrowi = $LUcolj = []; |
| 90 | |
| 91 | for ($j = 0; $j < $this->n; ++$j) { |
| 92 | for ($i = 0; $i < $this->m; ++$i) { |
| 93 | $LUcolj[$i] = &$this->LU[$i][$j]; |
| 94 | } |
| 95 | |
| 96 | for ($i = 0; $i < $this->m; ++$i) { |
| 97 | $LUrowi = $this->LU[$i]; |
| 98 | $kmax = \min($i, $j); |
| 99 | $s = 0.0; |
| 100 | |
| 101 | for ($k = 0; $k < $kmax; ++$k) { |
| 102 | $s += $LUrowi[$k] * $LUcolj[$k]; |
| 103 | } |
| 104 | $LUrowi[$j] = $LUcolj[$i] -= $s; |
| 105 | } |
| 106 | |
| 107 | $p = $j; |
| 108 | for ($i = $j + 1; $i < $this->m; ++$i) { |
| 109 | if (\abs($LUcolj[$i]) > \abs($LUcolj[$p])) { |
| 110 | $p = $i; |
| 111 | } |
| 112 | } |
| 113 | |
| 114 | if ($p !== $j) { |
| 115 | for ($k = 0; $k < $this->n; ++$k) { |
| 116 | $t = $this->LU[$p][$k]; |
| 117 | $this->LU[$p][$k] = $this->LU[$j][$k]; |
| 118 | $this->LU[$j][$k] = $t; |
| 119 | } |
| 120 | |
| 121 | $k = $this->piv[$p]; |
| 122 | $this->piv[$p] = $this->piv[$j]; |
| 123 | $this->piv[$j] = $k; |
| 124 | $this->pivSign *= -1; |
| 125 | } |
| 126 | |
| 127 | if (($j < $this->m) && ($this->LU[$j][$j] != 0.0)) { |
| 128 | for ($i = $j + 1; $i < $this->m; ++$i) { |
| 129 | $this->LU[$i][$j] /= $this->LU[$j][$j]; |
| 130 | } |
| 131 | } |
| 132 | } |
| 133 | } |
| 134 | |
| 135 | /** |
| 136 | * Get L matrix |
| 137 | * |
| 138 | * @return Matrix |
| 139 | * |
| 140 | * @since 1.0.0 |
| 141 | */ |
| 142 | public function getL() : Matrix |
| 143 | { |
| 144 | $L = [[]]; |
| 145 | |
| 146 | for ($i = 0; $i < $this->m; ++$i) { |
| 147 | for ($j = 0; $j < $this->n; ++$j) { |
| 148 | if ($i > $j) { |
| 149 | $L[$i][$j] = $this->LU[$i][$j]; |
| 150 | } elseif ($i === $j) { |
| 151 | $L[$i][$j] = 1.0; |
| 152 | } else { |
| 153 | $L[$i][$j] = 0.0; |
| 154 | } |
| 155 | } |
| 156 | } |
| 157 | |
| 158 | $matrix = new Matrix(); |
| 159 | $matrix->setMatrix($L); |
| 160 | |
| 161 | return $matrix; |
| 162 | } |
| 163 | |
| 164 | /** |
| 165 | * Get U matrix |
| 166 | * |
| 167 | * @return Matrix |
| 168 | * |
| 169 | * @since 1.0.0 |
| 170 | */ |
| 171 | public function getU() : Matrix |
| 172 | { |
| 173 | $U = [[]]; |
| 174 | |
| 175 | for ($i = 0; $i < $this->n; ++$i) { |
| 176 | for ($j = 0; $j < $this->n; ++$j) { |
| 177 | $U[$i][$j] = $i <= $j ? $this->LU[$i][$j] ?? 0 : 0.0; |
| 178 | } |
| 179 | } |
| 180 | |
| 181 | $matrix = new Matrix(); |
| 182 | $matrix->setMatrix($U); |
| 183 | |
| 184 | return $matrix; |
| 185 | } |
| 186 | |
| 187 | /** |
| 188 | * Get pivot |
| 189 | * |
| 190 | * @return array |
| 191 | * |
| 192 | * @since 1.0.0 |
| 193 | */ |
| 194 | public function getPivot() : array |
| 195 | { |
| 196 | return $this->piv; |
| 197 | } |
| 198 | |
| 199 | /** |
| 200 | * Is matrix nonsingular |
| 201 | * |
| 202 | * @return bool |
| 203 | * |
| 204 | * @since 1.0.0 |
| 205 | */ |
| 206 | public function isNonsingular() : bool |
| 207 | { |
| 208 | for ($j = 0; $j < $this->n; ++$j) { |
| 209 | if ($this->LU[$j][$j] == 0) { |
| 210 | return false; |
| 211 | } |
| 212 | } |
| 213 | |
| 214 | return true; |
| 215 | } |
| 216 | |
| 217 | /** |
| 218 | * Get determinant |
| 219 | * |
| 220 | * @return float |
| 221 | * |
| 222 | * @since 1.0.0 |
| 223 | */ |
| 224 | public function det() : float |
| 225 | { |
| 226 | $d = $this->pivSign; |
| 227 | for ($j = 0; $j < $this->n; ++$j) { |
| 228 | $d *= $this->LU[$j][$j]; |
| 229 | } |
| 230 | |
| 231 | return (float) $d; |
| 232 | } |
| 233 | |
| 234 | /** |
| 235 | * Solve Ax = b |
| 236 | * |
| 237 | * @param Matrix $B Matrix |
| 238 | * |
| 239 | * @return Matrix |
| 240 | * |
| 241 | * @throws InvalidDimensionException |
| 242 | * @throws \Exception |
| 243 | * |
| 244 | * @since 1.0.0 |
| 245 | */ |
| 246 | public function solve(Matrix $B) : Matrix |
| 247 | { |
| 248 | if ($B->getM() !== $this->m) { |
| 249 | throw new InvalidDimensionException((string) $B->getM()); |
| 250 | } |
| 251 | |
| 252 | if (!$this->isNonsingular()) { |
| 253 | throw new \Exception(); |
| 254 | } |
| 255 | |
| 256 | $n = $B->getN(); |
| 257 | $X = $B->getSubMatrixByRows($this->piv, 0, $n - 1)->toArray(); |
| 258 | |
| 259 | // Solve L*Y = B(piv,:) |
| 260 | for ($k = 0; $k < $this->n; ++$k) { |
| 261 | for ($i = $k + 1; $i < $this->n; ++$i) { |
| 262 | for ($j = 0; $j < $n; ++$j) { |
| 263 | $X[$i][$j] -= $X[$k][$j] * $this->LU[$i][$k]; |
| 264 | } |
| 265 | } |
| 266 | } |
| 267 | |
| 268 | // Solve U*X = Y; |
| 269 | for ($k = $this->n - 1; $k >= 0; --$k) { |
| 270 | for ($j = 0; $j < $n; ++$j) { |
| 271 | $X[$k][$j] /= $this->LU[$k][$k]; |
| 272 | } |
| 273 | for ($i = 0; $i < $k; ++$i) { |
| 274 | for ($j = 0; $j < $n; ++$j) { |
| 275 | $X[$i][$j] -= $X[$k][$j] * $this->LU[$i][$k]; |
| 276 | } |
| 277 | } |
| 278 | } |
| 279 | |
| 280 | $solution = new Matrix(); |
| 281 | $solution->setMatrix($X); |
| 282 | |
| 283 | return $solution; |
| 284 | } |
| 285 | } |