### LUDecomposition   B last analyzed 2020-05-15 05:48 UTC

#### Complexity

 Total Complexity 43

#### Size/Duplication

 Total Lines 261 Duplicated Lines 0 %

#### Importance

 Changes 0
Metric Value
wmc 43
eloc 93
dl 0
loc 261
rs 8.96
c 0
b 0
f 0

#### 9 Methods

Rating   Name   Duplication   Size   Complexity
A getDoublePivot() 0 3 1
A getPivot() 0 3 1
A getU() 0 14 4
A getSubMatrix() 0 13 3
A det() 0 8 2
A isNonsingular() 0 9 3
A getL() 0 16 5
B solve() 0 36 10
C __construct() 0 62 14

#### Complex Class

Complex classes like LUDecomposition often do a lot of different things. To break such a class down, we need to identify a cohesive component within that class. A common approach to find such a component is to look for fields/methods that share the same prefixes, or suffixes.

Once you have determined the fields that belong together, you can apply the Extract Class refactoring. If the component makes sense as a sub-class, Extract Subclass is also a candidate, and is often faster.

While breaking up the class, it is a good idea to analyze how other classes use LUDecomposition, and based on these observations, apply Extract Interface, too.

 1 `= n, the LU decomposition is an m-by-n` 9 ` * unit lower triangular matrix L, an n-by-n upper triangular matrix U,` 10 ` * and a permutation vector piv of length m so that A(piv,:) = L*U.` 11 ` * If m < n, then L is m-by-m and U is m-by-n.` 12 ` *` 13 ` * The LU decompostion with pivoting always exists, even if the matrix is` 14 ` * singular, so the constructor will never fail. The primary use of the` 15 ` * LU decomposition is in the solution of square systems of simultaneous` 16 ` * linear equations. This will fail if isNonsingular() returns false.` 17 ` *` 18 ` * @author Paul Meagher` 19 ` * @author Bartosz Matosiuk` 20 ` * @author Michael Bommarito` 21 ` *` 22 ` * @version 1.1` 23 ` *` 24 ` * @license PHP v3.0` 25 ` *` 26 ` * Slightly changed to adapt the original code to PHP-ML library` 27 ` * @date 2017/04/24` 28 ` *` 29 ` * @author Mustafa Karabulut` 30 ` */` 31 32 `namespace Phpml\Math\LinearAlgebra;` 33 34 `use Phpml\Exception\MatrixException;` 35 `use Phpml\Math\Matrix;` 36 37 `class LUDecomposition` 38 `{` 39 ` /**` 40 ` * Decomposition storage` 41 ` *` 42 ` * @var array` 43 ` */` 44 ` private \$LU = [];` 45 46 ` /**` 47 ` * Row dimension.` 48 ` *` 49 ` * @var int` 50 ` */` 51 ` private \$m;` 52 53 ` /**` 54 ` * Column dimension.` 55 ` *` 56 ` * @var int` 57 ` */` 58 ` private \$n;` 59 60 ` /**` 61 ` * Pivot sign.` 62 ` *` 63 ` * @var int` 64 ` */` 65 ` private \$pivsign;` 66 67 ` /**` 68 ` * Internal storage of pivot vector.` 69 ` *` 70 ` * @var array` 71 ` */` 72 ` private \$piv = [];` 73 74 ` /**` 75 ` * Constructs Structure to access L, U and piv.` 76 ` *` 77 ` * @param Matrix \$A Rectangular matrix` 78 ` *` 79 ` * @throws MatrixException` 80 ` */` 81 ` public function __construct(Matrix \$A)` 82 ` {` 83 ` if (\$A->getRows() !== \$A->getColumns()) {` 84 ` throw new MatrixException('Matrix is not square matrix');` 85 ` }` 86 87 ` // Use a "left-looking", dot-product, Crout/Doolittle algorithm.` 88 ` \$this->LU = \$A->toArray();` 89 ` \$this->m = \$A->getRows();` 90 ` \$this->n = \$A->getColumns();` 91 ` for (\$i = 0; \$i < \$this->m; ++\$i) {` 92 ` \$this->piv[\$i] = \$i;` 93 ` }` 94 95 ` \$this->pivsign = 1;` 96 ` \$LUcolj = [];` 97 98 ` // Outer loop.` 99 ` for (\$j = 0; \$j < \$this->n; ++\$j) {` 100 ` // Make a copy of the j-th column to localize references.` 101 ` for (\$i = 0; \$i < \$this->m; ++\$i) {` 102 ` \$LUcolj[\$i] = &\$this->LU[\$i][\$j];` 103 ` }` 104 105 ` // Apply previous transformations.` 106 ` for (\$i = 0; \$i < \$this->m; ++\$i) {` 107 ` \$LUrowi = \$this->LU[\$i];` 108 ` // Most of the time is spent in the following dot product.` 109 ` \$kmax = min(\$i, \$j);` 110 ` \$s = 0.0;` 111 ` for (\$k = 0; \$k < \$kmax; ++\$k) {` 112 ` \$s += \$LUrowi[\$k] * \$LUcolj[\$k];` 113 ` }` 114 115 ` \$LUrowi[\$j] = \$LUcolj[\$i] -= \$s;` 116 ` }` 117 118 ` // Find pivot and exchange if necessary.` 119 ` \$p = \$j;` 120 ` for (\$i = \$j + 1; \$i < \$this->m; ++\$i) {` 121 ` if (abs(\$LUcolj[\$i] ?? 0) > abs(\$LUcolj[\$p] ?? 0)) {` 122 ` \$p = \$i;` 123 ` }` 124 ` }` 125 126 ` if (\$p != \$j) {` 127 ` for (\$k = 0; \$k < \$this->n; ++\$k) {` 128 ` \$t = \$this->LU[\$p][\$k];` 129 ` \$this->LU[\$p][\$k] = \$this->LU[\$j][\$k];` 130 ` \$this->LU[\$j][\$k] = \$t;` 131 ` }` 132 133 ` \$k = \$this->piv[\$p];` 134 ` \$this->piv[\$p] = \$this->piv[\$j];` 135 ` \$this->piv[\$j] = \$k;` 136 ` \$this->pivsign *= -1;` 137 ` }` 138 139 ` // Compute multipliers.` 140 ` if ((\$j < \$this->m) && (\$this->LU[\$j][\$j] != 0.0)) {` 141 ` for (\$i = \$j + 1; \$i < \$this->m; ++\$i) {` 142 ` \$this->LU[\$i][\$j] /= \$this->LU[\$j][\$j];` 143 ` }` 144 ` }` 145 ` }` 146 ` }` 147 148 ` /**` 149 ` * Get lower triangular factor.` 150 ` *` 151 ` * @return Matrix Lower triangular factor` 152 ` */` 153 ` public function getL(): Matrix` 154 ` {` 155 ` \$L = [];` 156 ` for (\$i = 0; \$i < \$this->m; ++\$i) {` 157 ` for (\$j = 0; \$j < \$this->n; ++\$j) {` 158 ` if (\$i > \$j) {` 159 ` \$L[\$i][\$j] = \$this->LU[\$i][\$j];` 160 ` } elseif (\$i == \$j) {` 161 ` \$L[\$i][\$j] = 1.0;` 162 ` } else {` 163 ` \$L[\$i][\$j] = 0.0;` 164 ` }` 165 ` }` 166 ` }` 167 168 ` return new Matrix(\$L);` 169 ` }` 170 171 ` /**` 172 ` * Get upper triangular factor.` 173 ` *` 174 ` * @return Matrix Upper triangular factor` 175 ` */` 176 ` public function getU(): Matrix` 177 ` {` 178 ` \$U = [];` 179 ` for (\$i = 0; \$i < \$this->n; ++\$i) {` 180 ` for (\$j = 0; \$j < \$this->n; ++\$j) {` 181 ` if (\$i <= \$j) {` 182 ` \$U[\$i][\$j] = \$this->LU[\$i][\$j];` 183 ` } else {` 184 ` \$U[\$i][\$j] = 0.0;` 185 ` }` 186 ` }` 187 ` }` 188 189 ` return new Matrix(\$U);` 190 ` }` 191 192 ` /**` 193 ` * Return pivot permutation vector.` 194 ` *` 195 ` * @return array Pivot vector` 196 ` */` 197 ` public function getPivot(): array` 198 ` {` 199 ` return \$this->piv;` 200 ` }` 201 202 ` /**` 203 ` * Alias for getPivot` 204 ` *` 205 ` * @see getPivot` 206 ` */` 207 ` public function getDoublePivot(): array` 208 ` {` 209 ` return \$this->getPivot();` 210 ` }` 211 212 ` /**` 213 ` * Is the matrix nonsingular?` 214 ` *` 215 ` * @return bool true if U, and hence A, is nonsingular.` 216 ` */` 217 ` public function isNonsingular(): bool` 218 ` {` 219 ` for (\$j = 0; \$j < \$this->n; ++\$j) {` 220 ` if (\$this->LU[\$j][\$j] == 0) {` 221 ` return false;` 222 ` }` 223 ` }` 224 225 ` return true;` 226 ` }` 227 228 ` public function det(): float` 229 ` {` 230 ` \$d = \$this->pivsign;` 231 ` for (\$j = 0; \$j < \$this->n; ++\$j) {` 232 ` \$d *= \$this->LU[\$j][\$j];` 233 ` }` 234 235 ` return (float) \$d;` 236 ` }` 237 238 ` /**` 239 ` * Solve A*X = B` 240 ` *` 241 ` * @param Matrix \$B A Matrix with as many rows as A and any number of columns.` 242 ` *` 243 ` * @return array X so that L*U*X = B(piv,:)` 244 ` *` 245 ` * @throws MatrixException` 246 ` */` 247 ` public function solve(Matrix \$B): array` 248 ` {` 249 ` if (\$B->getRows() != \$this->m) {` 250 ` throw new MatrixException('Matrix is not square matrix');` 251 ` }` 252 253 ` if (!\$this->isNonsingular()) {` 254 ` throw new MatrixException('Matrix is singular');` 255 ` }` 256 257 ` // Copy right hand side with pivoting` 258 ` \$nx = \$B->getColumns();` 259 ` \$X = \$this->getSubMatrix(\$B->toArray(), \$this->piv, 0, \$nx - 1);` 260 ` // Solve L*Y = B(piv,:)` 261 ` for (\$k = 0; \$k < \$this->n; ++\$k) {` 262 ` for (\$i = \$k + 1; \$i < \$this->n; ++\$i) {` 263 ` for (\$j = 0; \$j < \$nx; ++\$j) {` 264 ` \$X[\$i][\$j] -= \$X[\$k][\$j] * \$this->LU[\$i][\$k];` 265 ` }` 266 ` }` 267 ` }` 268 269 ` // Solve U*X = Y;` 270 ` for (\$k = \$this->n - 1; \$k >= 0; --\$k) {` 271 ` for (\$j = 0; \$j < \$nx; ++\$j) {` 272 ` \$X[\$k][\$j] /= \$this->LU[\$k][\$k];` 273 ` }` 274 275 ` for (\$i = 0; \$i < \$k; ++\$i) {` 276 ` for (\$j = 0; \$j < \$nx; ++\$j) {` 277 ` \$X[\$i][\$j] -= \$X[\$k][\$j] * \$this->LU[\$i][\$k];` 278 ` }` 279 ` }` 280 ` }` 281 282 ` return \$X;` 283 ` }` 284 285 ` protected function getSubMatrix(array \$matrix, array \$RL, int \$j0, int \$jF): array` 286 ` {` 287 ` \$m = count(\$RL);` 288 ` \$n = \$jF - \$j0;` 289 ` \$R = array_fill(0, \$m, array_fill(0, \$n + 1, 0.0));` 290 291 ` for (\$i = 0; \$i < \$m; ++\$i) {` 292 ` for (\$j = \$j0; \$j <= \$jF; ++\$j) {` 293 ` \$R[\$i][\$j - \$j0] = \$matrix[\$RL[\$i]][\$j];` 294 ` }` 295 ` }` 296 297 ` return \$R;` 298 ` }` 299 `}` 300