001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017
018 package org.apache.commons.math.linear;
019
020 import org.apache.commons.math.MathRuntimeException;
021
022 /**
023 * Calculates the LUP-decomposition of a square matrix.
024 * <p>The LUP-decomposition of a matrix A consists of three matrices
025 * L, U and P that satisfy: PA = LU, L is lower triangular, and U is
026 * upper triangular and P is a permutation matrix. All matrices are
027 * m×m.</p>
028 * <p>As shown by the presence of the P matrix, this decomposition is
029 * implemented using partial pivoting.</p>
030 *
031 * @version $Revision: 885278 $ $Date: 2009-11-29 16:47:51 -0500 (Sun, 29 Nov 2009) $
032 * @since 2.0
033 */
034 public class LUDecompositionImpl implements LUDecomposition {
035
036 /** Default bound to determine effective singularity in LU decomposition */
037 private static final double DEFAULT_TOO_SMALL = 10E-12;
038
039 /** Message for vector length mismatch. */
040 private static final String VECTOR_LENGTH_MISMATCH_MESSAGE =
041 "vector length mismatch: got {0} but expected {1}";
042
043 /** Entries of LU decomposition. */
044 private double lu[][];
045
046 /** Pivot permutation associated with LU decomposition */
047 private int[] pivot;
048
049 /** Parity of the permutation associated with the LU decomposition */
050 private boolean even;
051
052 /** Singularity indicator. */
053 private boolean singular;
054
055 /** Cached value of L. */
056 private RealMatrix cachedL;
057
058 /** Cached value of U. */
059 private RealMatrix cachedU;
060
061 /** Cached value of P. */
062 private RealMatrix cachedP;
063
064 /**
065 * Calculates the LU-decomposition of the given matrix.
066 * @param matrix The matrix to decompose.
067 * @exception InvalidMatrixException if matrix is not square
068 */
069 public LUDecompositionImpl(RealMatrix matrix)
070 throws InvalidMatrixException {
071 this(matrix, DEFAULT_TOO_SMALL);
072 }
073
074 /**
075 * Calculates the LU-decomposition of the given matrix.
076 * @param matrix The matrix to decompose.
077 * @param singularityThreshold threshold (based on partial row norm)
078 * under which a matrix is considered singular
079 * @exception NonSquareMatrixException if matrix is not square
080 */
081 public LUDecompositionImpl(RealMatrix matrix, double singularityThreshold)
082 throws NonSquareMatrixException {
083
084 if (!matrix.isSquare()) {
085 throw new NonSquareMatrixException(matrix.getRowDimension(), matrix.getColumnDimension());
086 }
087
088 final int m = matrix.getColumnDimension();
089 lu = matrix.getData();
090 pivot = new int[m];
091 cachedL = null;
092 cachedU = null;
093 cachedP = null;
094
095 // Initialize permutation array and parity
096 for (int row = 0; row < m; row++) {
097 pivot[row] = row;
098 }
099 even = true;
100 singular = false;
101
102 // Loop over columns
103 for (int col = 0; col < m; col++) {
104
105 double sum = 0;
106
107 // upper
108 for (int row = 0; row < col; row++) {
109 final double[] luRow = lu[row];
110 sum = luRow[col];
111 for (int i = 0; i < row; i++) {
112 sum -= luRow[i] * lu[i][col];
113 }
114 luRow[col] = sum;
115 }
116
117 // lower
118 int max = col; // permutation row
119 double largest = Double.NEGATIVE_INFINITY;
120 for (int row = col; row < m; row++) {
121 final double[] luRow = lu[row];
122 sum = luRow[col];
123 for (int i = 0; i < col; i++) {
124 sum -= luRow[i] * lu[i][col];
125 }
126 luRow[col] = sum;
127
128 // maintain best permutation choice
129 if (Math.abs(sum) > largest) {
130 largest = Math.abs(sum);
131 max = row;
132 }
133 }
134
135 // Singularity check
136 if (Math.abs(lu[max][col]) < singularityThreshold) {
137 singular = true;
138 return;
139 }
140
141 // Pivot if necessary
142 if (max != col) {
143 double tmp = 0;
144 final double[] luMax = lu[max];
145 final double[] luCol = lu[col];
146 for (int i = 0; i < m; i++) {
147 tmp = luMax[i];
148 luMax[i] = luCol[i];
149 luCol[i] = tmp;
150 }
151 int temp = pivot[max];
152 pivot[max] = pivot[col];
153 pivot[col] = temp;
154 even = !even;
155 }
156
157 // Divide the lower elements by the "winning" diagonal elt.
158 final double luDiag = lu[col][col];
159 for (int row = col + 1; row < m; row++) {
160 lu[row][col] /= luDiag;
161 }
162 }
163
164 }
165
166 /** {@inheritDoc} */
167 public RealMatrix getL() {
168 if ((cachedL == null) && !singular) {
169 final int m = pivot.length;
170 cachedL = MatrixUtils.createRealMatrix(m, m);
171 for (int i = 0; i < m; ++i) {
172 final double[] luI = lu[i];
173 for (int j = 0; j < i; ++j) {
174 cachedL.setEntry(i, j, luI[j]);
175 }
176 cachedL.setEntry(i, i, 1.0);
177 }
178 }
179 return cachedL;
180 }
181
182 /** {@inheritDoc} */
183 public RealMatrix getU() {
184 if ((cachedU == null) && !singular) {
185 final int m = pivot.length;
186 cachedU = MatrixUtils.createRealMatrix(m, m);
187 for (int i = 0; i < m; ++i) {
188 final double[] luI = lu[i];
189 for (int j = i; j < m; ++j) {
190 cachedU.setEntry(i, j, luI[j]);
191 }
192 }
193 }
194 return cachedU;
195 }
196
197 /** {@inheritDoc} */
198 public RealMatrix getP() {
199 if ((cachedP == null) && !singular) {
200 final int m = pivot.length;
201 cachedP = MatrixUtils.createRealMatrix(m, m);
202 for (int i = 0; i < m; ++i) {
203 cachedP.setEntry(i, pivot[i], 1.0);
204 }
205 }
206 return cachedP;
207 }
208
209 /** {@inheritDoc} */
210 public int[] getPivot() {
211 return pivot.clone();
212 }
213
214 /** {@inheritDoc} */
215 public double getDeterminant() {
216 if (singular) {
217 return 0;
218 } else {
219 final int m = pivot.length;
220 double determinant = even ? 1 : -1;
221 for (int i = 0; i < m; i++) {
222 determinant *= lu[i][i];
223 }
224 return determinant;
225 }
226 }
227
228 /** {@inheritDoc} */
229 public DecompositionSolver getSolver() {
230 return new Solver(lu, pivot, singular);
231 }
232
233 /** Specialized solver. */
234 private static class Solver implements DecompositionSolver {
235
236 /** Entries of LU decomposition. */
237 private final double lu[][];
238
239 /** Pivot permutation associated with LU decomposition. */
240 private final int[] pivot;
241
242 /** Singularity indicator. */
243 private final boolean singular;
244
245 /**
246 * Build a solver from decomposed matrix.
247 * @param lu entries of LU decomposition
248 * @param pivot pivot permutation associated with LU decomposition
249 * @param singular singularity indicator
250 */
251 private Solver(final double[][] lu, final int[] pivot, final boolean singular) {
252 this.lu = lu;
253 this.pivot = pivot;
254 this.singular = singular;
255 }
256
257 /** {@inheritDoc} */
258 public boolean isNonSingular() {
259 return !singular;
260 }
261
262 /** {@inheritDoc} */
263 public double[] solve(double[] b)
264 throws IllegalArgumentException, InvalidMatrixException {
265
266 final int m = pivot.length;
267 if (b.length != m) {
268 throw MathRuntimeException.createIllegalArgumentException(
269 VECTOR_LENGTH_MISMATCH_MESSAGE, b.length, m);
270 }
271 if (singular) {
272 throw new SingularMatrixException();
273 }
274
275 final double[] bp = new double[m];
276
277 // Apply permutations to b
278 for (int row = 0; row < m; row++) {
279 bp[row] = b[pivot[row]];
280 }
281
282 // Solve LY = b
283 for (int col = 0; col < m; col++) {
284 final double bpCol = bp[col];
285 for (int i = col + 1; i < m; i++) {
286 bp[i] -= bpCol * lu[i][col];
287 }
288 }
289
290 // Solve UX = Y
291 for (int col = m - 1; col >= 0; col--) {
292 bp[col] /= lu[col][col];
293 final double bpCol = bp[col];
294 for (int i = 0; i < col; i++) {
295 bp[i] -= bpCol * lu[i][col];
296 }
297 }
298
299 return bp;
300
301 }
302
303 /** {@inheritDoc} */
304 public RealVector solve(RealVector b)
305 throws IllegalArgumentException, InvalidMatrixException {
306 try {
307 return solve((ArrayRealVector) b);
308 } catch (ClassCastException cce) {
309
310 final int m = pivot.length;
311 if (b.getDimension() != m) {
312 throw MathRuntimeException.createIllegalArgumentException(
313 VECTOR_LENGTH_MISMATCH_MESSAGE, b.getDimension(), m);
314 }
315 if (singular) {
316 throw new SingularMatrixException();
317 }
318
319 final double[] bp = new double[m];
320
321 // Apply permutations to b
322 for (int row = 0; row < m; row++) {
323 bp[row] = b.getEntry(pivot[row]);
324 }
325
326 // Solve LY = b
327 for (int col = 0; col < m; col++) {
328 final double bpCol = bp[col];
329 for (int i = col + 1; i < m; i++) {
330 bp[i] -= bpCol * lu[i][col];
331 }
332 }
333
334 // Solve UX = Y
335 for (int col = m - 1; col >= 0; col--) {
336 bp[col] /= lu[col][col];
337 final double bpCol = bp[col];
338 for (int i = 0; i < col; i++) {
339 bp[i] -= bpCol * lu[i][col];
340 }
341 }
342
343 return new ArrayRealVector(bp, false);
344
345 }
346 }
347
348 /** Solve the linear equation A × X = B.
349 * <p>The A matrix is implicit here. It is </p>
350 * @param b right-hand side of the equation A × X = B
351 * @return a vector X such that A × X = B
352 * @exception IllegalArgumentException if matrices dimensions don't match
353 * @exception InvalidMatrixException if decomposed matrix is singular
354 */
355 public ArrayRealVector solve(ArrayRealVector b)
356 throws IllegalArgumentException, InvalidMatrixException {
357 return new ArrayRealVector(solve(b.getDataRef()), false);
358 }
359
360 /** {@inheritDoc} */
361 public RealMatrix solve(RealMatrix b)
362 throws IllegalArgumentException, InvalidMatrixException {
363
364 final int m = pivot.length;
365 if (b.getRowDimension() != m) {
366 throw MathRuntimeException.createIllegalArgumentException(
367 "dimensions mismatch: got {0}x{1} but expected {2}x{3}",
368 b.getRowDimension(), b.getColumnDimension(), m, "n");
369 }
370 if (singular) {
371 throw new SingularMatrixException();
372 }
373
374 final int nColB = b.getColumnDimension();
375
376 // Apply permutations to b
377 final double[][] bp = new double[m][nColB];
378 for (int row = 0; row < m; row++) {
379 final double[] bpRow = bp[row];
380 final int pRow = pivot[row];
381 for (int col = 0; col < nColB; col++) {
382 bpRow[col] = b.getEntry(pRow, col);
383 }
384 }
385
386 // Solve LY = b
387 for (int col = 0; col < m; col++) {
388 final double[] bpCol = bp[col];
389 for (int i = col + 1; i < m; i++) {
390 final double[] bpI = bp[i];
391 final double luICol = lu[i][col];
392 for (int j = 0; j < nColB; j++) {
393 bpI[j] -= bpCol[j] * luICol;
394 }
395 }
396 }
397
398 // Solve UX = Y
399 for (int col = m - 1; col >= 0; col--) {
400 final double[] bpCol = bp[col];
401 final double luDiag = lu[col][col];
402 for (int j = 0; j < nColB; j++) {
403 bpCol[j] /= luDiag;
404 }
405 for (int i = 0; i < col; i++) {
406 final double[] bpI = bp[i];
407 final double luICol = lu[i][col];
408 for (int j = 0; j < nColB; j++) {
409 bpI[j] -= bpCol[j] * luICol;
410 }
411 }
412 }
413
414 return new Array2DRowRealMatrix(bp, false);
415
416 }
417
418 /** {@inheritDoc} */
419 public RealMatrix getInverse() throws InvalidMatrixException {
420 return solve(MatrixUtils.createRealIdentityMatrix(pivot.length));
421 }
422
423 }
424
425 }