source: src/main/java/agents/org/apache/commons/math/linear/LUDecompositionImpl.java

Last change on this file was 1, checked in by Wouter Pasman, 7 years ago

Initial import : Genius 9.0.0
File size: 13.7 KB
Line 
1/*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18package agents.org.apache.commons.math.linear;
19
20import agents.org.apache.commons.math.MathRuntimeException;
21import agents.org.apache.commons.math.exception.util.LocalizedFormats;
22import agents.org.apache.commons.math.util.FastMath;
23
24/**
25 * Calculates the LUP-decomposition of a square matrix.
26 * <p>The LUP-decomposition of a matrix A consists of three matrices
27 * L, U and P that satisfy: PA = LU, L is lower triangular, and U is
28 * upper triangular and P is a permutation matrix. All matrices are
29 * m&times;m.</p>
30 * <p>As shown by the presence of the P matrix, this decomposition is
31 * implemented using partial pivoting.</p>
32 *
33 * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 août 2010) $
34 * @since 2.0
35 */
36public class LUDecompositionImpl implements LUDecomposition {
37
38 /** Default bound to determine effective singularity in LU decomposition */
39 private static final double DEFAULT_TOO_SMALL = 10E-12;
40
41 /** Entries of LU decomposition. */
42 private double lu[][];
43
44 /** Pivot permutation associated with LU decomposition */
45 private int[] pivot;
46
47 /** Parity of the permutation associated with the LU decomposition */
48 private boolean even;
49
50 /** Singularity indicator. */
51 private boolean singular;
52
53 /** Cached value of L. */
54 private RealMatrix cachedL;
55
56 /** Cached value of U. */
57 private RealMatrix cachedU;
58
59 /** Cached value of P. */
60 private RealMatrix cachedP;
61
62 /**
63 * Calculates the LU-decomposition of the given matrix.
64 * @param matrix The matrix to decompose.
65 * @exception InvalidMatrixException if matrix is not square
66 */
67 public LUDecompositionImpl(RealMatrix matrix)
68 throws InvalidMatrixException {
69 this(matrix, DEFAULT_TOO_SMALL);
70 }
71
72 /**
73 * Calculates the LU-decomposition of the given matrix.
74 * @param matrix The matrix to decompose.
75 * @param singularityThreshold threshold (based on partial row norm)
76 * under which a matrix is considered singular
77 * @exception NonSquareMatrixException if matrix is not square
78 */
79 public LUDecompositionImpl(RealMatrix matrix, double singularityThreshold)
80 throws NonSquareMatrixException {
81
82 if (!matrix.isSquare()) {
83 throw new NonSquareMatrixException(matrix.getRowDimension(), matrix.getColumnDimension());
84 }
85
86 final int m = matrix.getColumnDimension();
87 lu = matrix.getData();
88 pivot = new int[m];
89 cachedL = null;
90 cachedU = null;
91 cachedP = null;
92
93 // Initialize permutation array and parity
94 for (int row = 0; row < m; row++) {
95 pivot[row] = row;
96 }
97 even = true;
98 singular = false;
99
100 // Loop over columns
101 for (int col = 0; col < m; col++) {
102
103 double sum = 0;
104
105 // upper
106 for (int row = 0; row < col; row++) {
107 final double[] luRow = lu[row];
108 sum = luRow[col];
109 for (int i = 0; i < row; i++) {
110 sum -= luRow[i] * lu[i][col];
111 }
112 luRow[col] = sum;
113 }
114
115 // lower
116 int max = col; // permutation row
117 double largest = Double.NEGATIVE_INFINITY;
118 for (int row = col; row < m; row++) {
119 final double[] luRow = lu[row];
120 sum = luRow[col];
121 for (int i = 0; i < col; i++) {
122 sum -= luRow[i] * lu[i][col];
123 }
124 luRow[col] = sum;
125
126 // maintain best permutation choice
127 if (FastMath.abs(sum) > largest) {
128 largest = FastMath.abs(sum);
129 max = row;
130 }
131 }
132
133 // Singularity check
134 if (FastMath.abs(lu[max][col]) < singularityThreshold) {
135 singular = true;
136 return;
137 }
138
139 // Pivot if necessary
140 if (max != col) {
141 double tmp = 0;
142 final double[] luMax = lu[max];
143 final double[] luCol = lu[col];
144 for (int i = 0; i < m; i++) {
145 tmp = luMax[i];
146 luMax[i] = luCol[i];
147 luCol[i] = tmp;
148 }
149 int temp = pivot[max];
150 pivot[max] = pivot[col];
151 pivot[col] = temp;
152 even = !even;
153 }
154
155 // Divide the lower elements by the "winning" diagonal elt.
156 final double luDiag = lu[col][col];
157 for (int row = col + 1; row < m; row++) {
158 lu[row][col] /= luDiag;
159 }
160 }
161
162 }
163
164 /** {@inheritDoc} */
165 public RealMatrix getL() {
166 if ((cachedL == null) && !singular) {
167 final int m = pivot.length;
168 cachedL = MatrixUtils.createRealMatrix(m, m);
169 for (int i = 0; i < m; ++i) {
170 final double[] luI = lu[i];
171 for (int j = 0; j < i; ++j) {
172 cachedL.setEntry(i, j, luI[j]);
173 }
174 cachedL.setEntry(i, i, 1.0);
175 }
176 }
177 return cachedL;
178 }
179
180 /** {@inheritDoc} */
181 public RealMatrix getU() {
182 if ((cachedU == null) && !singular) {
183 final int m = pivot.length;
184 cachedU = MatrixUtils.createRealMatrix(m, m);
185 for (int i = 0; i < m; ++i) {
186 final double[] luI = lu[i];
187 for (int j = i; j < m; ++j) {
188 cachedU.setEntry(i, j, luI[j]);
189 }
190 }
191 }
192 return cachedU;
193 }
194
195 /** {@inheritDoc} */
196 public RealMatrix getP() {
197 if ((cachedP == null) && !singular) {
198 final int m = pivot.length;
199 cachedP = MatrixUtils.createRealMatrix(m, m);
200 for (int i = 0; i < m; ++i) {
201 cachedP.setEntry(i, pivot[i], 1.0);
202 }
203 }
204 return cachedP;
205 }
206
207 /** {@inheritDoc} */
208 public int[] getPivot() {
209 return pivot.clone();
210 }
211
212 /** {@inheritDoc} */
213 public double getDeterminant() {
214 if (singular) {
215 return 0;
216 } else {
217 final int m = pivot.length;
218 double determinant = even ? 1 : -1;
219 for (int i = 0; i < m; i++) {
220 determinant *= lu[i][i];
221 }
222 return determinant;
223 }
224 }
225
226 /** {@inheritDoc} */
227 public DecompositionSolver getSolver() {
228 return new Solver(lu, pivot, singular);
229 }
230
231 /** Specialized solver. */
232 private static class Solver implements DecompositionSolver {
233
234 /** Entries of LU decomposition. */
235 private final double lu[][];
236
237 /** Pivot permutation associated with LU decomposition. */
238 private final int[] pivot;
239
240 /** Singularity indicator. */
241 private final boolean singular;
242
243 /**
244 * Build a solver from decomposed matrix.
245 * @param lu entries of LU decomposition
246 * @param pivot pivot permutation associated with LU decomposition
247 * @param singular singularity indicator
248 */
249 private Solver(final double[][] lu, final int[] pivot, final boolean singular) {
250 this.lu = lu;
251 this.pivot = pivot;
252 this.singular = singular;
253 }
254
255 /** {@inheritDoc} */
256 public boolean isNonSingular() {
257 return !singular;
258 }
259
260 /** {@inheritDoc} */
261 public double[] solve(double[] b)
262 throws IllegalArgumentException, InvalidMatrixException {
263
264 final int m = pivot.length;
265 if (b.length != m) {
266 throw MathRuntimeException.createIllegalArgumentException(
267 LocalizedFormats.VECTOR_LENGTH_MISMATCH, b.length, m);
268 }
269 if (singular) {
270 throw new SingularMatrixException();
271 }
272
273 final double[] bp = new double[m];
274
275 // Apply permutations to b
276 for (int row = 0; row < m; row++) {
277 bp[row] = b[pivot[row]];
278 }
279
280 // Solve LY = b
281 for (int col = 0; col < m; col++) {
282 final double bpCol = bp[col];
283 for (int i = col + 1; i < m; i++) {
284 bp[i] -= bpCol * lu[i][col];
285 }
286 }
287
288 // Solve UX = Y
289 for (int col = m - 1; col >= 0; col--) {
290 bp[col] /= lu[col][col];
291 final double bpCol = bp[col];
292 for (int i = 0; i < col; i++) {
293 bp[i] -= bpCol * lu[i][col];
294 }
295 }
296
297 return bp;
298
299 }
300
301 /** {@inheritDoc} */
302 public RealVector solve(RealVector b)
303 throws IllegalArgumentException, InvalidMatrixException {
304 try {
305 return solve((ArrayRealVector) b);
306 } catch (ClassCastException cce) {
307
308 final int m = pivot.length;
309 if (b.getDimension() != m) {
310 throw MathRuntimeException.createIllegalArgumentException(
311 LocalizedFormats.VECTOR_LENGTH_MISMATCH, b.getDimension(), m);
312 }
313 if (singular) {
314 throw new SingularMatrixException();
315 }
316
317 final double[] bp = new double[m];
318
319 // Apply permutations to b
320 for (int row = 0; row < m; row++) {
321 bp[row] = b.getEntry(pivot[row]);
322 }
323
324 // Solve LY = b
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 // Solve UX = Y
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 /** Solve the linear equation A &times; X = B.
347 * <p>The A matrix is implicit here. It is </p>
348 * @param b right-hand side of the equation A &times; X = B
349 * @return a vector X such that A &times; X = B
350 * @exception IllegalArgumentException if matrices dimensions don't match
351 * @exception InvalidMatrixException if decomposed matrix is singular
352 */
353 public ArrayRealVector solve(ArrayRealVector b)
354 throws IllegalArgumentException, InvalidMatrixException {
355 return new ArrayRealVector(solve(b.getDataRef()), false);
356 }
357
358 /** {@inheritDoc} */
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 LocalizedFormats.DIMENSIONS_MISMATCH_2x2,
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 // Apply permutations to b
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 // Solve LY = b
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 // Solve UX = Y
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 /** {@inheritDoc} */
417 public RealMatrix getInverse() throws InvalidMatrixException {
418 return solve(MatrixUtils.createRealIdentityMatrix(pivot.length));
419 }
420
421 }
422
423}
Note: See TracBrowser for help on using the repository browser.