1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math.linear;
19
20 import org.apache.commons.math.MathRuntimeException;
21
22
23
24
25
26
27
28
29
30
31
32
33
34 public class LUDecompositionImpl implements LUDecomposition {
35
36
37 private double lu[][];
38
39
40 private int[] pivot;
41
42
43 private boolean even;
44
45
46 private boolean singular;
47
48
49 private RealMatrix cachedL;
50
51
52 private RealMatrix cachedU;
53
54
55 private RealMatrix cachedP;
56
57
58 private static final double DEFAULT_TOO_SMALL = 10E-12;
59
60
61
62
63
64
65 public LUDecompositionImpl(RealMatrix matrix)
66 throws InvalidMatrixException {
67 this(matrix, DEFAULT_TOO_SMALL);
68 }
69
70
71
72
73
74
75
76
77 public LUDecompositionImpl(RealMatrix matrix, double singularityThreshold)
78 throws NonSquareMatrixException {
79
80 if (!matrix.isSquare()) {
81 throw new NonSquareMatrixException(matrix.getRowDimension(), matrix.getColumnDimension());
82 }
83
84 final int m = matrix.getColumnDimension();
85 lu = matrix.getData();
86 pivot = new int[m];
87 cachedL = null;
88 cachedU = null;
89 cachedP = null;
90
91
92 for (int row = 0; row < m; row++) {
93 pivot[row] = row;
94 }
95 even = true;
96 singular = false;
97
98
99 for (int col = 0; col < m; col++) {
100
101 double sum = 0;
102
103
104 for (int row = 0; row < col; row++) {
105 final double[] luRow = lu[row];
106 sum = luRow[col];
107 for (int i = 0; i < row; i++) {
108 sum -= luRow[i] * lu[i][col];
109 }
110 luRow[col] = sum;
111 }
112
113
114 int max = col;
115 double largest = Double.NEGATIVE_INFINITY;
116 for (int row = col; row < m; row++) {
117 final double[] luRow = lu[row];
118 sum = luRow[col];
119 for (int i = 0; i < col; i++) {
120 sum -= luRow[i] * lu[i][col];
121 }
122 luRow[col] = sum;
123
124
125 if (Math.abs(sum) > largest) {
126 largest = Math.abs(sum);
127 max = row;
128 }
129 }
130
131
132 if (Math.abs(lu[max][col]) < singularityThreshold) {
133 singular = true;
134 return;
135 }
136
137
138 if (max != col) {
139 double tmp = 0;
140 final double[] luMax = lu[max];
141 final double[] luCol = lu[col];
142 for (int i = 0; i < m; i++) {
143 tmp = luMax[i];
144 luMax[i] = luCol[i];
145 luCol[i] = tmp;
146 }
147 int temp = pivot[max];
148 pivot[max] = pivot[col];
149 pivot[col] = temp;
150 even = !even;
151 }
152
153
154 final double luDiag = lu[col][col];
155 for (int row = col + 1; row < m; row++) {
156 lu[row][col] /= luDiag;
157 }
158 }
159
160 }
161
162
163 public RealMatrix getL() {
164 if ((cachedL == null) && !singular) {
165 final int m = pivot.length;
166 cachedL = MatrixUtils.createRealMatrix(m, m);
167 for (int i = 0; i < m; ++i) {
168 final double[] luI = lu[i];
169 for (int j = 0; j < i; ++j) {
170 cachedL.setEntry(i, j, luI[j]);
171 }
172 cachedL.setEntry(i, i, 1.0);
173 }
174 }
175 return cachedL;
176 }
177
178
179 public RealMatrix getU() {
180 if ((cachedU == null) && !singular) {
181 final int m = pivot.length;
182 cachedU = MatrixUtils.createRealMatrix(m, m);
183 for (int i = 0; i < m; ++i) {
184 final double[] luI = lu[i];
185 for (int j = i; j < m; ++j) {
186 cachedU.setEntry(i, j, luI[j]);
187 }
188 }
189 }
190 return cachedU;
191 }
192
193
194 public RealMatrix getP() {
195 if ((cachedP == null) && !singular) {
196 final int m = pivot.length;
197 cachedP = MatrixUtils.createRealMatrix(m, m);
198 for (int i = 0; i < m; ++i) {
199 cachedP.setEntry(i, pivot[i], 1.0);
200 }
201 }
202 return cachedP;
203 }
204
205
206 public int[] getPivot() {
207 return pivot.clone();
208 }
209
210
211 public double getDeterminant() {
212 if (singular) {
213 return 0;
214 } else {
215 final int m = pivot.length;
216 double determinant = even ? 1 : -1;
217 for (int i = 0; i < m; i++) {
218 determinant *= lu[i][i];
219 }
220 return determinant;
221 }
222 }
223
224
225 public DecompositionSolver getSolver() {
226 return new Solver(lu, pivot, singular);
227 }
228
229
230 private static class Solver implements DecompositionSolver {
231
232
233 private final double lu[][];
234
235
236 private final int[] pivot;
237
238
239 private final boolean singular;
240
241
242
243
244
245
246
247 private Solver(final double[][] lu, final int[] pivot, final boolean singular) {
248 this.lu = lu;
249 this.pivot = pivot;
250 this.singular = singular;
251 }
252
253
254 public boolean isNonSingular() {
255 return !singular;
256 }
257
258
259 public double[] solve(double[] b)
260 throws IllegalArgumentException, InvalidMatrixException {
261
262 final int m = pivot.length;
263 if (b.length != m) {
264 throw MathRuntimeException.createIllegalArgumentException(
265 "vector length mismatch: got {0} but expected {1}",
266 b.length, m);
267 }
268 if (singular) {
269 throw new SingularMatrixException();
270 }
271
272 final double[] bp = new double[m];
273
274
275 for (int row = 0; row < m; row++) {
276 bp[row] = b[pivot[row]];
277 }
278
279
280 for (int col = 0; col < m; col++) {
281 final double bpCol = bp[col];
282 for (int i = col + 1; i < m; i++) {
283 bp[i] -= bpCol * lu[i][col];
284 }
285 }
286
287
288 for (int col = m - 1; col >= 0; col--) {
289 bp[col] /= lu[col][col];
290 final double bpCol = bp[col];
291 for (int i = 0; i < col; i++) {
292 bp[i] -= bpCol * lu[i][col];
293 }
294 }
295
296 return bp;
297
298 }
299
300
301 public RealVector solve(RealVector b)
302 throws IllegalArgumentException, InvalidMatrixException {
303 try {
304 return solve((ArrayRealVector) b);
305 } catch (ClassCastException cce) {
306
307 final int m = pivot.length;
308 if (b.getDimension() != m) {
309 throw MathRuntimeException.createIllegalArgumentException(
310 "vector length mismatch: got {0} but expected {1}",
311 b.getDimension(), m);
312 }
313 if (singular) {
314 throw new SingularMatrixException();
315 }
316
317 final double[] bp = new double[m];
318
319
320 for (int row = 0; row < m; row++) {
321 bp[row] = b.getEntry(pivot[row]);
322 }
323
324
325 for (int col = 0; col < m; col++) {
326 final double bpCol = bp[col];
327 for (int i = col + 1; i < m; i++) {
328 bp[i] -= bpCol * lu[i][col];
329 }
330 }
331
332
333 for (int col = m - 1; col >= 0; col--) {
334 bp[col] /= lu[col][col];
335 final double bpCol = bp[col];
336 for (int i = 0; i < col; i++) {
337 bp[i] -= bpCol * lu[i][col];
338 }
339 }
340
341 return new ArrayRealVector(bp, false);
342
343 }
344 }
345
346
347
348
349
350
351
352
353 public ArrayRealVector solve(ArrayRealVector b)
354 throws IllegalArgumentException, InvalidMatrixException {
355 return new ArrayRealVector(solve(b.getDataRef()), false);
356 }
357
358
359 public RealMatrix solve(RealMatrix b)
360 throws IllegalArgumentException, InvalidMatrixException {
361
362 final int m = pivot.length;
363 if (b.getRowDimension() != m) {
364 throw MathRuntimeException.createIllegalArgumentException(
365 "dimensions mismatch: got {0}x{1} but expected {2}x{3}",
366 b.getRowDimension(), b.getColumnDimension(), m, "n");
367 }
368 if (singular) {
369 throw new SingularMatrixException();
370 }
371
372 final int nColB = b.getColumnDimension();
373
374
375 final double[][] bp = new double[m][nColB];
376 for (int row = 0; row < m; row++) {
377 final double[] bpRow = bp[row];
378 final int pRow = pivot[row];
379 for (int col = 0; col < nColB; col++) {
380 bpRow[col] = b.getEntry(pRow, col);
381 }
382 }
383
384
385 for (int col = 0; col < m; col++) {
386 final double[] bpCol = bp[col];
387 for (int i = col + 1; i < m; i++) {
388 final double[] bpI = bp[i];
389 final double luICol = lu[i][col];
390 for (int j = 0; j < nColB; j++) {
391 bpI[j] -= bpCol[j] * luICol;
392 }
393 }
394 }
395
396
397 for (int col = m - 1; col >= 0; col--) {
398 final double[] bpCol = bp[col];
399 final double luDiag = lu[col][col];
400 for (int j = 0; j < nColB; j++) {
401 bpCol[j] /= luDiag;
402 }
403 for (int i = 0; i < col; i++) {
404 final double[] bpI = bp[i];
405 final double luICol = lu[i][col];
406 for (int j = 0; j < nColB; j++) {
407 bpI[j] -= bpCol[j] * luICol;
408 }
409 }
410 }
411
412 return new Array2DRowRealMatrix(bp, false);
413
414 }
415
416
417 public RealMatrix getInverse() throws InvalidMatrixException {
418 return solve(MatrixUtils.createRealIdentityMatrix(pivot.length));
419 }
420
421 }
422
423 }