Source code: jmat/data/matrixDecompositions/LUDecomposition.java
1 package jmat.data.matrixDecompositions;
2
3 import jmat.data.Matrix;
4
5
6 /** LU Decomposition.
7 <P>
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 In other words, assuming P the permutation Matrix, P*A = L*U.
12 If m < n, then L is m-by-m and U is m-by-n.
13 <P>
14 The LU decompostion with pivoting always exists, even if the matrix is
15 singular, so the constructor will never fail. The primary use of the
16 LU decomposition is in the solution of square systems of simultaneous
17 linear equations. This will fail if isNonsingular() returns false.
18 */
19 public class LUDecomposition implements java.io.Serializable
20 {
21 //~ Instance fields ////////////////////////////////////////////////////////
22
23 /* ------------------------
24 Class variables
25 * ------------------------ */
26
27 /** Array for internal storage of decomposition.
28 @serial internal array storage.
29 */
30 private double[][] LU;
31
32 /** Internal storage of pivot vector.
33 @serial pivot vector.
34 */
35 private int[] piv;
36
37 /** Row and column dimensions, and pivot sign.
38 @serial column dimension.
39 @serial row dimension.
40 @serial pivot sign.
41 */
42 private int m;
43
44 /** Row and column dimensions, and pivot sign.
45 @serial column dimension.
46 @serial row dimension.
47 @serial pivot sign.
48 */
49 private int n;
50
51 /** Row and column dimensions, and pivot sign.
52 @serial column dimension.
53 @serial row dimension.
54 @serial pivot sign.
55 */
56 private int pivsign;
57
58 //~ Constructors ///////////////////////////////////////////////////////////
59
60 /* ------------------------
61 Constructor
62 * ------------------------ */
63
64 /** LU Decomposition
65 @param A Rectangular matrix
66 @return Structure to access L, U and piv.
67 */
68 public LUDecomposition(Matrix A)
69 {
70 // Use a "left-looking", dot-product, Crout/Doolittle algorithm.
71 LU = A.getArrayCopy();
72 m = A.getRowDimension();
73 n = A.getColumnDimension();
74 piv = new int[m];
75
76 for (int i = 0; i < m; i++)
77 {
78 piv[i] = i;
79 }
80
81 pivsign = 1;
82
83 double[] LUrowi;
84 double[] LUcolj = new double[m];
85
86 // Outer loop.
87 for (int j = 0; j < n; j++)
88 {
89 // Make a copy of the j-th column to localize references.
90 for (int i = 0; i < m; i++)
91 {
92 LUcolj[i] = LU[i][j];
93 }
94
95 // Apply previous transformations.
96 for (int i = 0; i < m; i++)
97 {
98 LUrowi = LU[i];
99
100 // Most of the time is spent in the following dot product.
101 int kmax = Math.min(i, j);
102 double s = 0.0;
103
104 for (int k = 0; k < kmax; k++)
105 {
106 s += (LUrowi[k] * LUcolj[k]);
107 }
108
109 LUrowi[j] = LUcolj[i] -= s;
110 }
111
112 // Find pivot and exchange if necessary.
113 int p = j;
114
115 for (int i = j + 1; i < m; i++)
116 {
117 if (Math.abs(LUcolj[i]) > Math.abs(LUcolj[p]))
118 {
119 p = i;
120 }
121 }
122
123 if (p != j)
124 {
125 for (int k = 0; k < n; k++)
126 {
127 double t = LU[p][k];
128 LU[p][k] = LU[j][k];
129 LU[j][k] = t;
130 }
131
132 int k = piv[p];
133 piv[p] = piv[j];
134 piv[j] = k;
135 pivsign = -pivsign;
136 }
137
138 // Compute multipliers.
139 if ((j < m) & (LU[j][j] != 0.0))
140 {
141 for (int i = j + 1; i < m; i++)
142 {
143 LU[i][j] /= LU[j][j];
144 }
145 }
146 }
147 }
148
149 //~ Methods ////////////////////////////////////////////////////////////////
150
151 /** Return lower triangular factor
152 @return L
153 */
154 public Matrix getL()
155 {
156 Matrix X = new Matrix(m, n);
157 double[][] L = X.getArray();
158
159 for (int i = 0; i < m; i++)
160 {
161 for (int j = 0; j < n; j++)
162 {
163 if (i > j)
164 {
165 L[i][j] = LU[i][j];
166 }
167 else if (i == j)
168 {
169 L[i][j] = 1.0;
170 }
171 else
172 {
173 L[i][j] = 0.0;
174 }
175 }
176 }
177
178 return X;
179 }
180
181 /* ------------------------
182 Temporary, experimental code.
183 ------------------------ *\
184
185 \** LU Decomposition, computed by Gaussian elimination.
186 <P>
187 This constructor computes L and U with the "daxpy"-based elimination
188 algorithm used in LINPACK and MATLAB. In Java, we suspect the dot-product,
189 Crout algorithm will be faster. We have temporarily included this
190 constructor until timing experiments confirm this suspicion.
191 <P>
192 @param A Rectangular matrix
193 @param linpackflag Use Gaussian elimination. Actual value ignored.
194 @return Structure to access L, U and piv.
195 *\
196
197 public LUDecomposition (Matrix A, int linpackflag) {
198 // Initialize.
199 LU = A.getArrayCopy();
200 m = A.getRowDimension();
201 n = A.getColumnDimension();
202 piv = new int[m];
203 for (int i = 0; i < m; i++) {
204 piv[i] = i;
205 }
206 pivsign = 1;
207 // Main loop.
208 for (int k = 0; k < n; k++) {
209 // Find pivot.
210 int p = k;
211 for (int i = k+1; i < m; i++) {
212 if (Math.abs(LU[i][k]) > Math.abs(LU[p][k])) {
213 p = i;
214 }
215 }
216 // Exchange if necessary.
217 if (p != k) {
218 for (int j = 0; j < n; j++) {
219 double t = LU[p][j]; LU[p][j] = LU[k][j]; LU[k][j] = t;
220 }
221 int t = piv[p]; piv[p] = piv[k]; piv[k] = t;
222 pivsign = -pivsign;
223 }
224 // Compute multipliers and eliminate k-th column.
225 if (LU[k][k] != 0.0) {
226 for (int i = k+1; i < m; i++) {
227 LU[i][k] /= LU[k][k];
228 for (int j = k+1; j < n; j++) {
229 LU[i][j] -= LU[i][k]*LU[k][j];
230 }
231 }
232 }
233 }
234 }
235
236 \* ------------------------
237 End of temporary code.
238 * ------------------------ */
239 /* ------------------------
240 Public Methods
241 * ------------------------ */
242
243 /** Is the matrix nonsingular?
244 @return true if U, and hence A, is nonsingular.
245 */
246 public boolean isNonsingular()
247 {
248 for (int j = 0; j < n; j++)
249 {
250 if (LU[j][j] == 0)
251 {
252 return false;
253 }
254 }
255
256 return true;
257 }
258
259 /* ----------
260 JAMA code.
261 ---------- *\
262
263 \** Return pivot permutation vector
264 @return piv
265 *\
266
267 public int[] getPivot () {
268 int[] p = new int[m];
269 for (int i = 0; i < m; i++) {
270 p[i] = piv[i];
271 }
272 return p;
273 }
274
275 \** Return pivot permutation vector as a one-dimensional double array
276 @return (double) piv
277 *\
278
279 public double[] getDoublePivot () {
280 double[] vals = new double[m];
281 for (int i = 0; i < m; i++) {
282 vals[i] = (double) piv[i];
283 }
284 return vals;
285 }
286
287 \* -----------------
288 End of JAMA code.
289 * ----------------- */
290
291 /** Return pivot permutation vector
292 @return piv
293 */
294 public Matrix getP()
295 {
296 Matrix p = new Matrix(m, m);
297
298 for (int i = 0; i < m; i++)
299 {
300 p.set(i, piv[i], 1);
301 }
302
303 return p;
304 }
305
306 /** Return upper triangular factor
307 @return U
308 */
309 public Matrix getU()
310 {
311 Matrix X = new Matrix(n, n);
312 double[][] U = X.getArray();
313
314 for (int i = 0; i < n; i++)
315 {
316 for (int j = 0; j < n; j++)
317 {
318 if (i <= j)
319 {
320 U[i][j] = LU[i][j];
321 }
322 else
323 {
324 U[i][j] = 0.0;
325 }
326 }
327 }
328
329 return X;
330 }
331
332 /** Determinant
333 @return det(A)
334 @exception IllegalArgumentException Matrix must be square
335 */
336 public double det()
337 {
338 if (m != n)
339 {
340 throw new IllegalArgumentException("Matrix must be square.");
341 }
342
343 double d = (double) pivsign;
344
345 for (int j = 0; j < n; j++)
346 {
347 d *= LU[j][j];
348 }
349
350 return d;
351 }
352
353 /** Solve A*X = B
354 @param B A Matrix with as many rows as A and any number of columns.
355 @return X so that L*U*X = B(piv,:)
356 @exception IllegalArgumentException Matrix row dimensions must agree.
357 @exception RuntimeException Matrix is singular.
358 */
359 public Matrix solve(Matrix B)
360 {
361 if (B.getRowDimension() != m)
362 {
363 throw new IllegalArgumentException(
364 "Matrix row dimensions must agree.");
365 }
366
367 if (!this.isNonsingular())
368 {
369 throw new RuntimeException("Matrix is singular.");
370 }
371
372 // Copy right hand side with pivoting
373 int nx = B.getColumnDimension();
374 Matrix Xmat = B.getRows(piv); //B.getMatrix(piv,0,nx-1);
375 double[][] X = Xmat.getArray();
376
377 // Solve L*Y = B(piv,:)
378 for (int k = 0; k < n; k++)
379 {
380 for (int i = k + 1; i < n; i++)
381 {
382 for (int j = 0; j < nx; j++)
383 {
384 X[i][j] -= (X[k][j] * LU[i][k]);
385 }
386 }
387 }
388
389 // Solve U*X = Y;
390 for (int k = n - 1; k >= 0; k--)
391 {
392 for (int j = 0; j < nx; j++)
393 {
394 X[k][j] /= LU[k][k];
395 }
396
397 for (int i = 0; i < k; i++)
398 {
399 for (int j = 0; j < nx; j++)
400 {
401 X[i][j] -= (X[k][j] * LU[i][k]);
402 }
403 }
404 }
405
406 return Xmat;
407 }
408 }
409 ///////////////////////////////////////////////////////////////////////////////
410 // END OF FILE.
411 ///////////////////////////////////////////////////////////////////////////////