1<?php 2/** 3 * @package JAMA 4 * 5 * For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n 6 * unit lower triangular matrix L, an n-by-n upper triangular matrix U, 7 * and a permutation vector piv of length m so that A(piv,:) = L*U. 8 * If m < n, then L is m-by-m and U is m-by-n. 9 * 10 * The LU decompostion with pivoting always exists, even if the matrix is 11 * singular, so the constructor will never fail. The primary use of the 12 * LU decomposition is in the solution of square systems of simultaneous 13 * linear equations. This will fail if isNonsingular() returns false. 14 * 15 * @author Paul Meagher 16 * @author Bartosz Matosiuk 17 * @author Michael Bommarito 18 * @version 1.1 19 * @license PHP v3.0 20 */ 21class PHPExcel_Shared_JAMA_LUDecomposition { 22 23 const MatrixSingularException = "Can only perform operation on singular matrix."; 24 const MatrixSquareException = "Mismatched Row dimension"; 25 26 /** 27 * Decomposition storage 28 * @var array 29 */ 30 private $LU = array(); 31 32 /** 33 * Row dimension. 34 * @var int 35 */ 36 private $m; 37 38 /** 39 * Column dimension. 40 * @var int 41 */ 42 private $n; 43 44 /** 45 * Pivot sign. 46 * @var int 47 */ 48 private $pivsign; 49 50 /** 51 * Internal storage of pivot vector. 52 * @var array 53 */ 54 private $piv = array(); 55 56 57 /** 58 * LU Decomposition constructor. 59 * 60 * @param $A Rectangular matrix 61 * @return Structure to access L, U and piv. 62 */ 63 public function __construct($A) { 64 if ($A instanceof PHPExcel_Shared_JAMA_Matrix) { 65 // Use a "left-looking", dot-product, Crout/Doolittle algorithm. 66 $this->LU = $A->getArray(); 67 $this->m = $A->getRowDimension(); 68 $this->n = $A->getColumnDimension(); 69 for ($i = 0; $i < $this->m; ++$i) { 70 $this->piv[$i] = $i; 71 } 72 $this->pivsign = 1; 73 $LUrowi = $LUcolj = array(); 74 75 // Outer loop. 76 for ($j = 0; $j < $this->n; ++$j) { 77 // Make a copy of the j-th column to localize references. 78 for ($i = 0; $i < $this->m; ++$i) { 79 $LUcolj[$i] = &$this->LU[$i][$j]; 80 } 81 // Apply previous transformations. 82 for ($i = 0; $i < $this->m; ++$i) { 83 $LUrowi = $this->LU[$i]; 84 // Most of the time is spent in the following dot product. 85 $kmax = min($i,$j); 86 $s = 0.0; 87 for ($k = 0; $k < $kmax; ++$k) { 88 $s += $LUrowi[$k] * $LUcolj[$k]; 89 } 90 $LUrowi[$j] = $LUcolj[$i] -= $s; 91 } 92 // Find pivot and exchange if necessary. 93 $p = $j; 94 for ($i = $j+1; $i < $this->m; ++$i) { 95 if (abs($LUcolj[$i]) > abs($LUcolj[$p])) { 96 $p = $i; 97 } 98 } 99 if ($p != $j) { 100 for ($k = 0; $k < $this->n; ++$k) { 101 $t = $this->LU[$p][$k]; 102 $this->LU[$p][$k] = $this->LU[$j][$k]; 103 $this->LU[$j][$k] = $t; 104 } 105 $k = $this->piv[$p]; 106 $this->piv[$p] = $this->piv[$j]; 107 $this->piv[$j] = $k; 108 $this->pivsign = $this->pivsign * -1; 109 } 110 // Compute multipliers. 111 if (($j < $this->m) && ($this->LU[$j][$j] != 0.0)) { 112 for ($i = $j+1; $i < $this->m; ++$i) { 113 $this->LU[$i][$j] /= $this->LU[$j][$j]; 114 } 115 } 116 } 117 } else { 118 throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::ArgumentTypeException); 119 } 120 } // function __construct() 121 122 123 /** 124 * Get lower triangular factor. 125 * 126 * @return array Lower triangular factor 127 */ 128 public function getL() { 129 for ($i = 0; $i < $this->m; ++$i) { 130 for ($j = 0; $j < $this->n; ++$j) { 131 if ($i > $j) { 132 $L[$i][$j] = $this->LU[$i][$j]; 133 } elseif ($i == $j) { 134 $L[$i][$j] = 1.0; 135 } else { 136 $L[$i][$j] = 0.0; 137 } 138 } 139 } 140 return new PHPExcel_Shared_JAMA_Matrix($L); 141 } // function getL() 142 143 144 /** 145 * Get upper triangular factor. 146 * 147 * @return array Upper triangular factor 148 */ 149 public function getU() { 150 for ($i = 0; $i < $this->n; ++$i) { 151 for ($j = 0; $j < $this->n; ++$j) { 152 if ($i <= $j) { 153 $U[$i][$j] = $this->LU[$i][$j]; 154 } else { 155 $U[$i][$j] = 0.0; 156 } 157 } 158 } 159 return new PHPExcel_Shared_JAMA_Matrix($U); 160 } // function getU() 161 162 163 /** 164 * Return pivot permutation vector. 165 * 166 * @return array Pivot vector 167 */ 168 public function getPivot() { 169 return $this->piv; 170 } // function getPivot() 171 172 173 /** 174 * Alias for getPivot 175 * 176 * @see getPivot 177 */ 178 public function getDoublePivot() { 179 return $this->getPivot(); 180 } // function getDoublePivot() 181 182 183 /** 184 * Is the matrix nonsingular? 185 * 186 * @return true if U, and hence A, is nonsingular. 187 */ 188 public function isNonsingular() { 189 for ($j = 0; $j < $this->n; ++$j) { 190 if ($this->LU[$j][$j] == 0) { 191 return false; 192 } 193 } 194 return true; 195 } // function isNonsingular() 196 197 198 /** 199 * Count determinants 200 * 201 * @return array d matrix deterninat 202 */ 203 public function det() { 204 if ($this->m == $this->n) { 205 $d = $this->pivsign; 206 for ($j = 0; $j < $this->n; ++$j) { 207 $d *= $this->LU[$j][$j]; 208 } 209 return $d; 210 } else { 211 throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::MatrixDimensionException); 212 } 213 } // function det() 214 215 216 /** 217 * Solve A*X = B 218 * 219 * @param $B A Matrix with as many rows as A and any number of columns. 220 * @return X so that L*U*X = B(piv,:) 221 * @PHPExcel_Calculation_Exception IllegalArgumentException Matrix row dimensions must agree. 222 * @PHPExcel_Calculation_Exception RuntimeException Matrix is singular. 223 */ 224 public function solve($B) { 225 if ($B->getRowDimension() == $this->m) { 226 if ($this->isNonsingular()) { 227 // Copy right hand side with pivoting 228 $nx = $B->getColumnDimension(); 229 $X = $B->getMatrix($this->piv, 0, $nx-1); 230 // Solve L*Y = B(piv,:) 231 for ($k = 0; $k < $this->n; ++$k) { 232 for ($i = $k+1; $i < $this->n; ++$i) { 233 for ($j = 0; $j < $nx; ++$j) { 234 $X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k]; 235 } 236 } 237 } 238 // Solve U*X = Y; 239 for ($k = $this->n-1; $k >= 0; --$k) { 240 for ($j = 0; $j < $nx; ++$j) { 241 $X->A[$k][$j] /= $this->LU[$k][$k]; 242 } 243 for ($i = 0; $i < $k; ++$i) { 244 for ($j = 0; $j < $nx; ++$j) { 245 $X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k]; 246 } 247 } 248 } 249 return $X; 250 } else { 251 throw new PHPExcel_Calculation_Exception(self::MatrixSingularException); 252 } 253 } else { 254 throw new PHPExcel_Calculation_Exception(self::MatrixSquareException); 255 } 256 } // function solve() 257 258} // class PHPExcel_Shared_JAMA_LUDecomposition 259