@@ -47,29 +47,31 @@ private function decompose()
4747 {
4848 $ scaling = [];
4949 $ this ->parity = 1 ; // start parity at +1 (parity is "even" for zero row interchanges)
50- $ n = $ this ->getRowCount ();
51- $ p =& $ this ->permutations ;
50+ $ rowCount = $ this ->getRowCount ();
5251
5352 // We want to find the largest element in each row for scaling.
54- for ($ i = 0 ; $ i < $ n ; ++$ i ) {
53+ for ($ i = 0 ; $ i < $ rowCount ; ++$ i ) {
5554 $ biggest = 0 ;
56- for ($ j = 0 ; $ j < $ n ; ++$ j ) {
55+
56+ for ($ j = 0 ; $ j < $ rowCount ; ++$ j ) {
5757 $ temp = abs ($ this ->internal [$ i ][$ j ]);
5858 $ biggest = max ($ temp , $ biggest );
5959 }
60+
6061 if ($ biggest == 0 ) {
6162 throw new MatrixException ("Matrix is singular. " );
6263 }
64+
6365 $ scaling [$ i ] = 1 / $ biggest ;
64- $ p [$ i ] = $ i ; // Initialize permutations vector
66+ $ this -> permutations [$ i ] = $ i ; // Initialize permutations vector
6567 }
6668
6769 // Now we find the LU decomposition. This is the outer loop over diagonal elements.
68- for ($ k = 0 ; $ k < $ n ; ++$ k ) {
70+ for ($ k = 0 ; $ k < $ rowCount ; ++$ k ) {
6971 // Search for the best (biggest) pivot element
7072 $ biggest = 0 ;
7173 $ maxRowIndex = $ k ;
72- for ($ i = $ k ; $ i < $ n ; ++$ i ) {
74+ for ($ i = $ k ; $ i < $ rowCount ; ++$ i ) {
7375 $ temp = $ scaling [$ i ] * abs ($ this ->internal [$ i ][$ k ]);
7476 if ($ temp > $ biggest ) {
7577 $ biggest = $ temp ;
@@ -80,9 +82,9 @@ private function decompose()
8082 // Perform the row pivot and store in the permutations vector
8183 if ($ k != $ maxRowIndex ) {
8284 $ this ->rowPivot ($ k , $ maxRowIndex );
83- $ temp = $ p [$ k ];
84- $ p [$ k ] = $ p [$ maxRowIndex ];
85- $ p [$ maxRowIndex ] = $ temp ;
85+ $ temp = $ this -> permutations [$ k ];
86+ $ this -> permutations [$ k ] = $ this -> permutations [$ maxRowIndex ];
87+ $ this -> permutations [$ maxRowIndex ] = $ temp ;
8688 $ this ->parity = -$ this ->parity ; // flip parity
8789 }
8890
@@ -91,13 +93,13 @@ private function decompose()
9193 }
9294
9395 // Crout's algorithm
94- for ($ i = $ k + 1 ; $ i < $ n ; ++$ i ) {
96+ for ($ i = $ k + 1 ; $ i < $ rowCount ; ++$ i ) {
9597
9698 // Divide by the pivot element
9799 $ this ->internal [$ i ][$ k ] = $ this ->internal [$ i ][$ k ] / $ this ->internal [$ k ][$ k ];
98100
99101 // Subtract from each element in the sub-matrix
100- for ($ j = $ k + 1 ; $ j < $ n ; ++$ j ) {
102+ for ($ j = $ k + 1 ; $ j < $ rowCount ; ++$ j ) {
101103 $ this ->internal [$ i ][$ j ] = $ this ->internal [$ i ][$ j ] - $ this ->internal [$ i ][$ k ] * $ this ->internal [$ k ][$ j ];
102104 }
103105 }
@@ -158,29 +160,28 @@ public function solve(array $known)
158160 $ skip = true ;
159161
160162 // Solve L * y = b for y (forward substitution)
161- for ($ i = 0 ; $ i < $ rowCount ; ++ $ i ) {
163+ for ($ i = 0 ; $ i < $ rowCount ; $ i ++ ) {
162164 $ thisB = $ known [$ this ->permutations [$ i ]]; // Unscramble the permutations
163- if ($ skip && $ thisB == 0 ) { // Leading zeroes in b give zeroes in y.
165+
166+ if ($ skip && $ thisB == 0 ) {
167+ // Leading zeroes in b give zeroes in y.
164168 $ y [$ i ] = 0 ;
165169 } else {
166- if ($ skip ) {
167- // We found a non-zero element, so don't skip any more.
168- $ skip = false ;
169- }
170-
170+ // We found a non-zero element, so don't skip any more.
171+ $ skip = false ;
171172 $ y [$ i ] = $ thisB ;
172173
173- for ($ j = 0 ; $ j < $ i ; ++ $ j ) {
174+ for ($ j = 0 ; $ j < $ i ; $ j ++ ) {
174175 $ y [$ i ] = $ y [$ i ] - $ this ->get ($ i , $ j ) * $ y [$ j ];
175176 }
176177 }
177178 }
178179
179180 // Solve U * unknown = y for unknown (backward substitution)
180- for ($ i = $ rowCount - 1 ; $ i >= 0 ; -- $ i ) {
181+ for ($ i = $ rowCount - 1 ; $ i >= 0 ; $ i -- ) {
181182 $ unknown [$ i ] = $ y [$ i ];
182183
183- for ($ j = $ i + 1 ; $ j < $ rowCount ; ++ $ j ) {
184+ for ($ j = $ i + 1 ; $ j < $ rowCount ; $ j ++ ) {
184185 $ unknown [$ i ] = $ unknown [$ i ] - $ this ->get ($ i , $ j ) * $ unknown [$ j ];
185186 }
186187
@@ -201,22 +202,21 @@ public function solve(array $known)
201202 public function inverse ()
202203 {
203204 $ inverse = [];
204-
205- // Get size of matrix
206- $ n = $ this ->getRowCount ();
207- $ b = array_fill (0 , $ n , 0 ); // initialize empty vector
205+ $ rowCount = $ this ->getRowCount ();
206+ $ currentRow = array_fill (0 , $ rowCount , 0 );
208207
209208 // For each j from 0 to n-1
210- for ($ j = 0 ; $ j < $ n ; ++ $ j ) {
209+ for ($ j = 0 ; $ j < $ rowCount ; $ j ++ ) {
211210 // this is the jth column of the identity matrix
212- $ b [$ j ] = 1 ;
211+ $ currentRow [$ j ] = 1 ;
212+
213+ $ solution = $ this ->solve ($ currentRow );
213214
214- // get the solution vector and copy to the jth column of the inverse matrix
215- $ x = $ this ->solve ($ b );
216- for ($ i = 0 ; $ i < $ n ; ++$ i ) {
217- $ inverse [$ i ][$ j ] = $ x [$ i ];
215+ for ($ i = 0 ; $ i < $ rowCount ; ++$ i ) {
216+ $ inverse [$ i ][$ j ] = $ solution [$ i ];
218217 }
219- $ b [$ j ] = 0 ; // Get the vector ready for the next column.
218+
219+ $ currentRow [$ j ] = 0 ; // Get the vector ready for the next column.
220220 }
221221
222222 return new Matrix ($ inverse );
0 commit comments