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.math3.linear;
019
020 import java.lang.reflect.Array;
021
022 import org.apache.commons.math3.Field;
023 import org.apache.commons.math3.FieldElement;
024 import org.apache.commons.math3.exception.DimensionMismatchException;
025
026 /**
027 * Calculates the LUP-decomposition of a square matrix.
028 * <p>The LUP-decomposition of a matrix A consists of three matrices
029 * L, U and P that satisfy: PA = LU, L is lower triangular, and U is
030 * upper triangular and P is a permutation matrix. All matrices are
031 * m×m.</p>
032 * <p>Since {@link FieldElement field elements} do not provide an ordering
033 * operator, the permutation matrix is computed here only in order to avoid
034 * a zero pivot element, no attempt is done to get the largest pivot
035 * element.</p>
036 * <p>This class is based on the class with similar name from the
037 * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library.</p>
038 * <ul>
039 * <li>a {@link #getP() getP} method has been added,</li>
040 * <li>the {@code det} method has been renamed as {@link #getDeterminant()
041 * getDeterminant},</li>
042 * <li>the {@code getDoublePivot} method has been removed (but the int based
043 * {@link #getPivot() getPivot} method has been kept),</li>
044 * <li>the {@code solve} and {@code isNonSingular} methods have been replaced
045 * by a {@link #getSolver() getSolver} method and the equivalent methods
046 * provided by the returned {@link DecompositionSolver}.</li>
047 * </ul>
048 *
049 * @param <T> the type of the field elements
050 * @see <a href="http://mathworld.wolfram.com/LUDecomposition.html">MathWorld</a>
051 * @see <a href="http://en.wikipedia.org/wiki/LU_decomposition">Wikipedia</a>
052 * @version $Id: FieldLUDecomposition.java 1416643 2012-12-03 19:37:14Z tn $
053 * @since 2.0 (changed to concrete class in 3.0)
054 */
055 public class FieldLUDecomposition<T extends FieldElement<T>> {
056
057 /** Field to which the elements belong. */
058 private final Field<T> field;
059
060 /** Entries of LU decomposition. */
061 private T[][] lu;
062
063 /** Pivot permutation associated with LU decomposition. */
064 private int[] pivot;
065
066 /** Parity of the permutation associated with the LU decomposition. */
067 private boolean even;
068
069 /** Singularity indicator. */
070 private boolean singular;
071
072 /** Cached value of L. */
073 private FieldMatrix<T> cachedL;
074
075 /** Cached value of U. */
076 private FieldMatrix<T> cachedU;
077
078 /** Cached value of P. */
079 private FieldMatrix<T> cachedP;
080
081 /**
082 * Calculates the LU-decomposition of the given matrix.
083 * @param matrix The matrix to decompose.
084 * @throws NonSquareMatrixException if matrix is not square
085 */
086 public FieldLUDecomposition(FieldMatrix<T> matrix) {
087 if (!matrix.isSquare()) {
088 throw new NonSquareMatrixException(matrix.getRowDimension(),
089 matrix.getColumnDimension());
090 }
091
092 final int m = matrix.getColumnDimension();
093 field = matrix.getField();
094 lu = matrix.getData();
095 pivot = new int[m];
096 cachedL = null;
097 cachedU = null;
098 cachedP = null;
099
100 // Initialize permutation array and parity
101 for (int row = 0; row < m; row++) {
102 pivot[row] = row;
103 }
104 even = true;
105 singular = false;
106
107 // Loop over columns
108 for (int col = 0; col < m; col++) {
109
110 T sum = field.getZero();
111
112 // upper
113 for (int row = 0; row < col; row++) {
114 final T[] luRow = lu[row];
115 sum = luRow[col];
116 for (int i = 0; i < row; i++) {
117 sum = sum.subtract(luRow[i].multiply(lu[i][col]));
118 }
119 luRow[col] = sum;
120 }
121
122 // lower
123 int nonZero = col; // permutation row
124 for (int row = col; row < m; row++) {
125 final T[] luRow = lu[row];
126 sum = luRow[col];
127 for (int i = 0; i < col; i++) {
128 sum = sum.subtract(luRow[i].multiply(lu[i][col]));
129 }
130 luRow[col] = sum;
131
132 if (lu[nonZero][col].equals(field.getZero())) {
133 // try to select a better permutation choice
134 ++nonZero;
135 }
136 }
137
138 // Singularity check
139 if (nonZero >= m) {
140 singular = true;
141 return;
142 }
143
144 // Pivot if necessary
145 if (nonZero != col) {
146 T tmp = field.getZero();
147 for (int i = 0; i < m; i++) {
148 tmp = lu[nonZero][i];
149 lu[nonZero][i] = lu[col][i];
150 lu[col][i] = tmp;
151 }
152 int temp = pivot[nonZero];
153 pivot[nonZero] = pivot[col];
154 pivot[col] = temp;
155 even = !even;
156 }
157
158 // Divide the lower elements by the "winning" diagonal elt.
159 final T luDiag = lu[col][col];
160 for (int row = col + 1; row < m; row++) {
161 final T[] luRow = lu[row];
162 luRow[col] = luRow[col].divide(luDiag);
163 }
164 }
165
166 }
167
168 /**
169 * Returns the matrix L of the decomposition.
170 * <p>L is a lower-triangular matrix</p>
171 * @return the L matrix (or null if decomposed matrix is singular)
172 */
173 public FieldMatrix<T> getL() {
174 if ((cachedL == null) && !singular) {
175 final int m = pivot.length;
176 cachedL = new Array2DRowFieldMatrix<T>(field, m, m);
177 for (int i = 0; i < m; ++i) {
178 final T[] luI = lu[i];
179 for (int j = 0; j < i; ++j) {
180 cachedL.setEntry(i, j, luI[j]);
181 }
182 cachedL.setEntry(i, i, field.getOne());
183 }
184 }
185 return cachedL;
186 }
187
188 /**
189 * Returns the matrix U of the decomposition.
190 * <p>U is an upper-triangular matrix</p>
191 * @return the U matrix (or null if decomposed matrix is singular)
192 */
193 public FieldMatrix<T> getU() {
194 if ((cachedU == null) && !singular) {
195 final int m = pivot.length;
196 cachedU = new Array2DRowFieldMatrix<T>(field, m, m);
197 for (int i = 0; i < m; ++i) {
198 final T[] luI = lu[i];
199 for (int j = i; j < m; ++j) {
200 cachedU.setEntry(i, j, luI[j]);
201 }
202 }
203 }
204 return cachedU;
205 }
206
207 /**
208 * Returns the P rows permutation matrix.
209 * <p>P is a sparse matrix with exactly one element set to 1.0 in
210 * each row and each column, all other elements being set to 0.0.</p>
211 * <p>The positions of the 1 elements are given by the {@link #getPivot()
212 * pivot permutation vector}.</p>
213 * @return the P rows permutation matrix (or null if decomposed matrix is singular)
214 * @see #getPivot()
215 */
216 public FieldMatrix<T> getP() {
217 if ((cachedP == null) && !singular) {
218 final int m = pivot.length;
219 cachedP = new Array2DRowFieldMatrix<T>(field, m, m);
220 for (int i = 0; i < m; ++i) {
221 cachedP.setEntry(i, pivot[i], field.getOne());
222 }
223 }
224 return cachedP;
225 }
226
227 /**
228 * Returns the pivot permutation vector.
229 * @return the pivot permutation vector
230 * @see #getP()
231 */
232 public int[] getPivot() {
233 return pivot.clone();
234 }
235
236 /**
237 * Return the determinant of the matrix.
238 * @return determinant of the matrix
239 */
240 public T getDeterminant() {
241 if (singular) {
242 return field.getZero();
243 } else {
244 final int m = pivot.length;
245 T determinant = even ? field.getOne() : field.getZero().subtract(field.getOne());
246 for (int i = 0; i < m; i++) {
247 determinant = determinant.multiply(lu[i][i]);
248 }
249 return determinant;
250 }
251 }
252
253 /**
254 * Get a solver for finding the A × X = B solution in exact linear sense.
255 * @return a solver
256 */
257 public FieldDecompositionSolver<T> getSolver() {
258 return new Solver<T>(field, lu, pivot, singular);
259 }
260
261 /** Specialized solver. */
262 private static class Solver<T extends FieldElement<T>> implements FieldDecompositionSolver<T> {
263
264 /** Field to which the elements belong. */
265 private final Field<T> field;
266
267 /** Entries of LU decomposition. */
268 private final T[][] lu;
269
270 /** Pivot permutation associated with LU decomposition. */
271 private final int[] pivot;
272
273 /** Singularity indicator. */
274 private final boolean singular;
275
276 /**
277 * Build a solver from decomposed matrix.
278 * @param field field to which the matrix elements belong
279 * @param lu entries of LU decomposition
280 * @param pivot pivot permutation associated with LU decomposition
281 * @param singular singularity indicator
282 */
283 private Solver(final Field<T> field, final T[][] lu,
284 final int[] pivot, final boolean singular) {
285 this.field = field;
286 this.lu = lu;
287 this.pivot = pivot;
288 this.singular = singular;
289 }
290
291 /** {@inheritDoc} */
292 public boolean isNonSingular() {
293 return !singular;
294 }
295
296 /** {@inheritDoc} */
297 public FieldVector<T> solve(FieldVector<T> b) {
298 try {
299 return solve((ArrayFieldVector<T>) b);
300 } catch (ClassCastException cce) {
301
302 final int m = pivot.length;
303 if (b.getDimension() != m) {
304 throw new DimensionMismatchException(b.getDimension(), m);
305 }
306 if (singular) {
307 throw new SingularMatrixException();
308 }
309
310 @SuppressWarnings("unchecked") // field is of type T
311 final T[] bp = (T[]) Array.newInstance(field.getRuntimeClass(), m);
312
313 // Apply permutations to b
314 for (int row = 0; row < m; row++) {
315 bp[row] = b.getEntry(pivot[row]);
316 }
317
318 // Solve LY = b
319 for (int col = 0; col < m; col++) {
320 final T bpCol = bp[col];
321 for (int i = col + 1; i < m; i++) {
322 bp[i] = bp[i].subtract(bpCol.multiply(lu[i][col]));
323 }
324 }
325
326 // Solve UX = Y
327 for (int col = m - 1; col >= 0; col--) {
328 bp[col] = bp[col].divide(lu[col][col]);
329 final T bpCol = bp[col];
330 for (int i = 0; i < col; i++) {
331 bp[i] = bp[i].subtract(bpCol.multiply(lu[i][col]));
332 }
333 }
334
335 return new ArrayFieldVector<T>(field, bp, false);
336
337 }
338 }
339
340 /** Solve the linear equation A × X = B.
341 * <p>The A matrix is implicit here. It is </p>
342 * @param b right-hand side of the equation A × X = B
343 * @return a vector X such that A × X = B
344 * @throws DimensionMismatchException if the matrices dimensions do not match.
345 * @throws SingularMatrixException if the decomposed matrix is singular.
346 */
347 public ArrayFieldVector<T> solve(ArrayFieldVector<T> b) {
348 final int m = pivot.length;
349 final int length = b.getDimension();
350 if (length != m) {
351 throw new DimensionMismatchException(length, m);
352 }
353 if (singular) {
354 throw new SingularMatrixException();
355 }
356
357 @SuppressWarnings("unchecked")
358 // field is of type T
359 final T[] bp = (T[]) Array.newInstance(field.getRuntimeClass(),
360 m);
361
362 // Apply permutations to b
363 for (int row = 0; row < m; row++) {
364 bp[row] = b.getEntry(pivot[row]);
365 }
366
367 // Solve LY = b
368 for (int col = 0; col < m; col++) {
369 final T bpCol = bp[col];
370 for (int i = col + 1; i < m; i++) {
371 bp[i] = bp[i].subtract(bpCol.multiply(lu[i][col]));
372 }
373 }
374
375 // Solve UX = Y
376 for (int col = m - 1; col >= 0; col--) {
377 bp[col] = bp[col].divide(lu[col][col]);
378 final T bpCol = bp[col];
379 for (int i = 0; i < col; i++) {
380 bp[i] = bp[i].subtract(bpCol.multiply(lu[i][col]));
381 }
382 }
383
384 return new ArrayFieldVector<T>(bp, false);
385 }
386
387 /** {@inheritDoc} */
388 public FieldMatrix<T> solve(FieldMatrix<T> b) {
389 final int m = pivot.length;
390 if (b.getRowDimension() != m) {
391 throw new DimensionMismatchException(b.getRowDimension(), m);
392 }
393 if (singular) {
394 throw new SingularMatrixException();
395 }
396
397 final int nColB = b.getColumnDimension();
398
399 // Apply permutations to b
400 @SuppressWarnings("unchecked") // field is of type T
401 final T[][] bp = (T[][]) Array.newInstance(field.getRuntimeClass(), new int[] { m, nColB });
402 for (int row = 0; row < m; row++) {
403 final T[] bpRow = bp[row];
404 final int pRow = pivot[row];
405 for (int col = 0; col < nColB; col++) {
406 bpRow[col] = b.getEntry(pRow, col);
407 }
408 }
409
410 // Solve LY = b
411 for (int col = 0; col < m; col++) {
412 final T[] bpCol = bp[col];
413 for (int i = col + 1; i < m; i++) {
414 final T[] bpI = bp[i];
415 final T luICol = lu[i][col];
416 for (int j = 0; j < nColB; j++) {
417 bpI[j] = bpI[j].subtract(bpCol[j].multiply(luICol));
418 }
419 }
420 }
421
422 // Solve UX = Y
423 for (int col = m - 1; col >= 0; col--) {
424 final T[] bpCol = bp[col];
425 final T luDiag = lu[col][col];
426 for (int j = 0; j < nColB; j++) {
427 bpCol[j] = bpCol[j].divide(luDiag);
428 }
429 for (int i = 0; i < col; i++) {
430 final T[] bpI = bp[i];
431 final T luICol = lu[i][col];
432 for (int j = 0; j < nColB; j++) {
433 bpI[j] = bpI[j].subtract(bpCol[j].multiply(luICol));
434 }
435 }
436 }
437
438 return new Array2DRowFieldMatrix<T>(field, bp, false);
439
440 }
441
442 /** {@inheritDoc} */
443 public FieldMatrix<T> getInverse() {
444 final int m = pivot.length;
445 final T one = field.getOne();
446 FieldMatrix<T> identity = new Array2DRowFieldMatrix<T>(field, m, m);
447 for (int i = 0; i < m; ++i) {
448 identity.setEntry(i, i, one);
449 }
450 return solve(identity);
451 }
452 }
453 }