These results are based on our legacy PHP analysis, consider migrating to our new PHP analysis engine instead. Learn more
1  <?php 

2  
3  declare(strict_types=1); 

4  
5  /** 

6  * @package JAMA 

7  * 

8  * For an mbyn matrix A with m >= n, the LU decomposition is an mbyn 

9  * unit lower triangular matrix L, an nbyn 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 mbym and U is mbyn. 

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 PHPML 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 MatrixException::notSquareMatrix(); 

85  } 

86  
87  // Use a "leftlooking", dotproduct, 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 jth column to localize references. 

101  View Code Duplication  for ($i = 0; $i < $this>m; ++$i) { 

0 ignored issues
–
show


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]) > abs($LUcolj[$p])) { 

122  $p = $i; 

123  } 

124  } 

125  
126  if ($p != $j) { 

127  View Code Duplication  for ($k = 0; $k < $this>n; ++$k) { 

0 ignored issues
–
show
This code seems to be duplicated across your project.
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository.
Loading history...


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 = $this>pivsign * 1; 

137  } 

138  
139  // Compute multipliers. 

140  if (($j < $this>m) && ($this>LU[$j][$j] != 0.0)) { 

141  View Code Duplication  for ($i = $j + 1; $i < $this>m; ++$i) { 

0 ignored issues
–
show
This code seems to be duplicated across your project.
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository.
Loading history...


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() 

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  View Code Duplication  for ($j = 0; $j < $this>n; ++$j) { 

0 ignored issues
–
show
This code seems to be duplicated across your project.
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository.
Loading history...


220  if ($this>LU[$j][$j] == 0) { 

221  return false; 

222  } 

223  } 

224  
225  return true; 

226  } 

227  
228  /** 

229  * Count determinants 

230  * 

231  * @return floatint d matrix determinant 

232  * 

233  * @throws MatrixException 

234  */ 

235  public function det() 

236  { 

237  if ($this>m !== $this>n) { 

238  throw MatrixException::notSquareMatrix(); 

239  } 

240  
241  $d = $this>pivsign; 

242  View Code Duplication  for ($j = 0; $j < $this>n; ++$j) { 

0 ignored issues
–
show
This code seems to be duplicated across your project.
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository.
Loading history...


243  $d *= $this>LU[$j][$j]; 

244  } 

245  
246  return $d; 

247  } 

248  
249  /** 

250  * Solve A*X = B 

251  * 

252  * @param Matrix $B A Matrix with as many rows as A and any number of columns. 

253  * 

254  * @return array X so that L*U*X = B(piv,:) 

255  * 

256  * @throws MatrixException 

257  */ 

258  public function solve(Matrix $B): array 

259  { 

260  if ($B>getRows() != $this>m) { 

261  throw MatrixException::notSquareMatrix(); 

262  } 

263  
264  if (!$this>isNonsingular()) { 

265  throw MatrixException::singularMatrix(); 

266  } 

267  
268  // Copy right hand side with pivoting 

269  $nx = $B>getColumns(); 

270  $X = $this>getSubMatrix($B>toArray(), $this>piv, 0, $nx  1); 

271  // Solve L*Y = B(piv,:) 

272  for ($k = 0; $k < $this>n; ++$k) { 

273  View Code Duplication  for ($i = $k + 1; $i < $this>n; ++$i) { 

0 ignored issues
–
show
This code seems to be duplicated across your project.
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository.
Loading history...


274  for ($j = 0; $j < $nx; ++$j) { 

275  $X[$i][$j] = $X[$k][$j] * $this>LU[$i][$k]; 

276  } 

277  } 

278  } 

279  
280  // Solve U*X = Y; 

281  for ($k = $this>n  1; $k >= 0; $k) { 

282  for ($j = 0; $j < $nx; ++$j) { 

283  $X[$k][$j] /= $this>LU[$k][$k]; 

284  } 

285  
286  View Code Duplication  for ($i = 0; $i < $k; ++$i) { 

0 ignored issues
–
show
This code seems to be duplicated across your project.
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository.
Loading history...


287  for ($j = 0; $j < $nx; ++$j) { 

288  $X[$i][$j] = $X[$k][$j] * $this>LU[$i][$k]; 

289  } 

290  } 

291  } 

292  
293  return $X; 

294  } 

295  
296  protected function getSubMatrix(array $matrix, array $RL, int $j0, int $jF): array 

297  { 

298  $m = count($RL); 

299  $n = $jF  $j0; 

300  $R = array_fill(0, $m, array_fill(0, $n + 1, 0.0)); 

301  
302  for ($i = 0; $i < $m; ++$i) { 

303  for ($j = $j0; $j <= $jF; ++$j) { 

304  $R[$i][$j  $j0] = $matrix[$RL[$i]][$j]; 

305  } 

306  } 

307  
308  return $R; 

309  } 

310  } 

311 
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation.
You can also find more detailed suggestions in the “Code” section of your repository.