Differences Between: [Versions 310 and 311] [Versions 310 and 400] [Versions 310 and 401]
1 <?php 2 3 namespace PhpOffice\PhpSpreadsheet\Shared\JAMA; 4 5 use PhpOffice\PhpSpreadsheet\Calculation\Exception as CalculationException; 6 7 /** 8 * For an m-by-n matrix A with m >= 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 class LUDecomposition 25 { 26 const MATRIX_SINGULAR_EXCEPTION = 'Can only perform operation on singular matrix.'; 27 const MATRIX_SQUARE_EXCEPTION = 'Mismatched Row dimension'; 28 29 /** 30 * Decomposition storage. 31 * 32 * @var array 33 */ 34 private $LU = []; 35 36 /** 37 * Row dimension. 38 * 39 * @var int 40 */ 41 private $m; 42 43 /** 44 * Column dimension. 45 * 46 * @var int 47 */ 48 private $n; 49 50 /** 51 * Pivot sign. 52 * 53 * @var int 54 */ 55 private $pivsign; 56 57 /** 58 * Internal storage of pivot vector. 59 * 60 * @var array 61 */ 62 private $piv = []; 63 64 /** 65 * LU Decomposition constructor. 66 * 67 * @param Matrix $A Rectangular matrix 68 */ 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 131 /** 132 * Get lower triangular factor. 133 * 134 * @return Matrix Lower triangular factor 135 */ 136 public function getL() 137 { 138 for ($i = 0; $i < $this->m; ++$i) { 139 for ($j = 0; $j < $this->n; ++$j) { 140 if ($i > $j) { 141 $L[$i][$j] = $this->LU[$i][$j]; 142 } elseif ($i == $j) { 143 $L[$i][$j] = 1.0; 144 } else { 145 $L[$i][$j] = 0.0; 146 } 147 } 148 } 149 150 return new Matrix($L); 151 } 152 153 // function getL() 154 155 /** 156 * Get upper triangular factor. 157 * 158 * @return Matrix Upper triangular factor 159 */ 160 public function getU() 161 { 162 for ($i = 0; $i < $this->n; ++$i) { 163 for ($j = 0; $j < $this->n; ++$j) { 164 if ($i <= $j) { 165 $U[$i][$j] = $this->LU[$i][$j]; 166 } else { 167 $U[$i][$j] = 0.0; 168 } 169 } 170 } 171 172 return new Matrix($U); 173 } 174 175 // function getU() 176 177 /** 178 * Return pivot permutation vector. 179 * 180 * @return array Pivot vector 181 */ 182 public function getPivot() 183 { 184 return $this->piv; 185 } 186 187 // function getPivot() 188 189 /** 190 * Alias for getPivot. 191 * 192 * @see getPivot 193 */ 194 public function getDoublePivot() 195 { 196 return $this->getPivot(); 197 } 198 199 // function getDoublePivot() 200 201 /** 202 * Is the matrix nonsingular? 203 * 204 * @return bool true if U, and hence A, is nonsingular 205 */ 206 public function isNonsingular() 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 // function isNonsingular() 218 219 /** 220 * Count determinants. 221 * 222 * @return array d matrix deterninat 223 */ 224 public function det() 225 { 226 if ($this->m == $this->n) { 227 $d = $this->pivsign; 228 for ($j = 0; $j < $this->n; ++$j) { 229 $d *= $this->LU[$j][$j]; 230 } 231 232 return $d; 233 } 234 235 throw new CalculationException(Matrix::MATRIX_DIMENSION_EXCEPTION); 236 } 237 238 // function det() 239 240 /** 241 * Solve A*X = B. 242 * 243 * @param mixed $B a Matrix with as many rows as A and any number of columns 244 * 245 * @throws CalculationException illegalArgumentException Matrix row dimensions must agree 246 * @throws CalculationException runtimeException Matrix is singular 247 * 248 * @return Matrix X so that L*U*X = B(piv,:) 249 */ 250 public function solve($B) 251 { 252 if ($B->getRowDimension() == $this->m) { 253 if ($this->isNonsingular()) { 254 // Copy right hand side with pivoting 255 $nx = $B->getColumnDimension(); 256 $X = $B->getMatrix($this->piv, 0, $nx - 1); 257 // Solve L*Y = B(piv,:) 258 for ($k = 0; $k < $this->n; ++$k) { 259 for ($i = $k + 1; $i < $this->n; ++$i) { 260 for ($j = 0; $j < $nx; ++$j) { 261 $X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k]; 262 } 263 } 264 } 265 // Solve U*X = Y; 266 for ($k = $this->n - 1; $k >= 0; --$k) { 267 for ($j = 0; $j < $nx; ++$j) { 268 $X->A[$k][$j] /= $this->LU[$k][$k]; 269 } 270 for ($i = 0; $i < $k; ++$i) { 271 for ($j = 0; $j < $nx; ++$j) { 272 $X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k]; 273 } 274 } 275 } 276 277 return $X; 278 } 279 280 throw new CalculationException(self::MATRIX_SINGULAR_EXCEPTION); 281 } 282 283 throw new CalculationException(self::MATRIX_SQUARE_EXCEPTION); 284 } 285 }
title
Description
Body
title
Description
Body
title
Description
Body
title
Body