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 java.lang.reflect.Array;
021
022 import org.apache.commons.math.Field;
023 import org.apache.commons.math.FieldElement;
024 import org.apache.commons.math.MathRuntimeException;
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 element.</p>
035 *
036 * @param <T> the type of the field elements
037 * @version $Revision: 903046 $ $Date: 2010-01-25 21:07:26 -0500 (Mon, 25 Jan 2010) $
038 * @since 2.0
039 */
040 public class FieldLUDecompositionImpl<T extends FieldElement<T>> implements FieldLUDecomposition<T> {
041
042 /** Field to which the elements belong. */
043 private final Field<T> field;
044
045 /** Entries of LU decomposition. */
046 private T lu[][];
047
048 /** Pivot permutation associated with LU decomposition */
049 private int[] pivot;
050
051 /** Parity of the permutation associated with the LU decomposition */
052 private boolean even;
053
054 /** Singularity indicator. */
055 private boolean singular;
056
057 /** Cached value of L. */
058 private FieldMatrix<T> cachedL;
059
060 /** Cached value of U. */
061 private FieldMatrix<T> cachedU;
062
063 /** Cached value of P. */
064 private FieldMatrix<T> cachedP;
065
066 /**
067 * Calculates the LU-decomposition of the given matrix.
068 * @param matrix The matrix to decompose.
069 * @exception NonSquareMatrixException if matrix is not square
070 */
071 public FieldLUDecompositionImpl(FieldMatrix<T> matrix)
072 throws NonSquareMatrixException {
073
074 if (!matrix.isSquare()) {
075 throw new NonSquareMatrixException(matrix.getRowDimension(), matrix.getColumnDimension());
076 }
077
078 final int m = matrix.getColumnDimension();
079 field = matrix.getField();
080 lu = matrix.getData();
081 pivot = new int[m];
082 cachedL = null;
083 cachedU = null;
084 cachedP = null;
085
086 // Initialize permutation array and parity
087 for (int row = 0; row < m; row++) {
088 pivot[row] = row;
089 }
090 even = true;
091 singular = false;
092
093 // Loop over columns
094 for (int col = 0; col < m; col++) {
095
096 T sum = field.getZero();
097
098 // upper
099 for (int row = 0; row < col; row++) {
100 final T[] luRow = lu[row];
101 sum = luRow[col];
102 for (int i = 0; i < row; i++) {
103 sum = sum.subtract(luRow[i].multiply(lu[i][col]));
104 }
105 luRow[col] = sum;
106 }
107
108 // lower
109 int nonZero = col; // permutation row
110 for (int row = col; row < m; row++) {
111 final T[] luRow = lu[row];
112 sum = luRow[col];
113 for (int i = 0; i < col; i++) {
114 sum = sum.subtract(luRow[i].multiply(lu[i][col]));
115 }
116 luRow[col] = sum;
117
118 if (lu[nonZero][col].equals(field.getZero())) {
119 // try to select a better permutation choice
120 ++nonZero;
121 }
122 }
123
124 // Singularity check
125 if (nonZero >= m) {
126 singular = true;
127 return;
128 }
129
130 // Pivot if necessary
131 if (nonZero != col) {
132 T tmp = field.getZero();
133 for (int i = 0; i < m; i++) {
134 tmp = lu[nonZero][i];
135 lu[nonZero][i] = lu[col][i];
136 lu[col][i] = tmp;
137 }
138 int temp = pivot[nonZero];
139 pivot[nonZero] = pivot[col];
140 pivot[col] = temp;
141 even = !even;
142 }
143
144 // Divide the lower elements by the "winning" diagonal elt.
145 final T luDiag = lu[col][col];
146 for (int row = col + 1; row < m; row++) {
147 final T[] luRow = lu[row];
148 luRow[col] = luRow[col].divide(luDiag);
149 }
150 }
151
152 }
153
154 /** {@inheritDoc} */
155 public FieldMatrix<T> getL() {
156 if ((cachedL == null) && !singular) {
157 final int m = pivot.length;
158 cachedL = new Array2DRowFieldMatrix<T>(field, m, m);
159 for (int i = 0; i < m; ++i) {
160 final T[] luI = lu[i];
161 for (int j = 0; j < i; ++j) {
162 cachedL.setEntry(i, j, luI[j]);
163 }
164 cachedL.setEntry(i, i, field.getOne());
165 }
166 }
167 return cachedL;
168 }
169
170 /** {@inheritDoc} */
171 public FieldMatrix<T> getU() {
172 if ((cachedU == null) && !singular) {
173 final int m = pivot.length;
174 cachedU = new Array2DRowFieldMatrix<T>(field, m, m);
175 for (int i = 0; i < m; ++i) {
176 final T[] luI = lu[i];
177 for (int j = i; j < m; ++j) {
178 cachedU.setEntry(i, j, luI[j]);
179 }
180 }
181 }
182 return cachedU;
183 }
184
185 /** {@inheritDoc} */
186 public FieldMatrix<T> getP() {
187 if ((cachedP == null) && !singular) {
188 final int m = pivot.length;
189 cachedP = new Array2DRowFieldMatrix<T>(field, m, m);
190 for (int i = 0; i < m; ++i) {
191 cachedP.setEntry(i, pivot[i], field.getOne());
192 }
193 }
194 return cachedP;
195 }
196
197 /** {@inheritDoc} */
198 public int[] getPivot() {
199 return pivot.clone();
200 }
201
202 /** {@inheritDoc} */
203 public T getDeterminant() {
204 if (singular) {
205 return field.getZero();
206 } else {
207 final int m = pivot.length;
208 T determinant = even ? field.getOne() : field.getZero().subtract(field.getOne());
209 for (int i = 0; i < m; i++) {
210 determinant = determinant.multiply(lu[i][i]);
211 }
212 return determinant;
213 }
214 }
215
216 /** {@inheritDoc} */
217 public FieldDecompositionSolver<T> getSolver() {
218 return new Solver<T>(field, lu, pivot, singular);
219 }
220
221 /** Specialized solver. */
222 private static class Solver<T extends FieldElement<T>> implements FieldDecompositionSolver<T> {
223
224 /** Serializable version identifier. */
225 private static final long serialVersionUID = -6353105415121373022L;
226
227 /** Field to which the elements belong. */
228 private final Field<T> field;
229
230 /** Entries of LU decomposition. */
231 private final T lu[][];
232
233 /** Pivot permutation associated with LU decomposition. */
234 private final int[] pivot;
235
236 /** Singularity indicator. */
237 private final boolean singular;
238
239 /**
240 * Build a solver from decomposed matrix.
241 * @param field field to which the matrix elements belong
242 * @param lu entries of LU decomposition
243 * @param pivot pivot permutation associated with LU decomposition
244 * @param singular singularity indicator
245 */
246 private Solver(final Field<T> field, final T[][] lu,
247 final int[] pivot, final boolean singular) {
248 this.field = field;
249 this.lu = lu;
250 this.pivot = pivot;
251 this.singular = singular;
252 }
253
254 /** {@inheritDoc} */
255 public boolean isNonSingular() {
256 return !singular;
257 }
258
259 /** {@inheritDoc} */
260 public T[] solve(T[] b)
261 throws IllegalArgumentException, InvalidMatrixException {
262
263 final int m = pivot.length;
264 if (b.length != m) {
265 throw MathRuntimeException.createIllegalArgumentException(
266 "vector length mismatch: got {0} but expected {1}",
267 b.length, m);
268 }
269 if (singular) {
270 throw new SingularMatrixException();
271 }
272
273 @SuppressWarnings("unchecked") // field is of type T
274 final T[] bp = (T[]) Array.newInstance(field.getZero().getClass(), m);
275
276 // Apply permutations to b
277 for (int row = 0; row < m; row++) {
278 bp[row] = b[pivot[row]];
279 }
280
281 // Solve LY = b
282 for (int col = 0; col < m; col++) {
283 final T bpCol = bp[col];
284 for (int i = col + 1; i < m; i++) {
285 bp[i] = bp[i].subtract(bpCol.multiply(lu[i][col]));
286 }
287 }
288
289 // Solve UX = Y
290 for (int col = m - 1; col >= 0; col--) {
291 bp[col] = bp[col].divide(lu[col][col]);
292 final T bpCol = bp[col];
293 for (int i = 0; i < col; i++) {
294 bp[i] = bp[i].subtract(bpCol.multiply(lu[i][col]));
295 }
296 }
297
298 return bp;
299
300 }
301
302 /** {@inheritDoc} */
303 public FieldVector<T> solve(FieldVector<T> b)
304 throws IllegalArgumentException, InvalidMatrixException {
305 try {
306 return solve((ArrayFieldVector<T>) b);
307 } catch (ClassCastException cce) {
308
309 final int m = pivot.length;
310 if (b.getDimension() != m) {
311 throw MathRuntimeException.createIllegalArgumentException(
312 "vector length mismatch: got {0} but expected {1}",
313 b.getDimension(), m);
314 }
315 if (singular) {
316 throw new SingularMatrixException();
317 }
318
319 @SuppressWarnings("unchecked") // field is of type T
320 final T[] bp = (T[]) Array.newInstance(field.getZero().getClass(), m);
321
322 // Apply permutations to b
323 for (int row = 0; row < m; row++) {
324 bp[row] = b.getEntry(pivot[row]);
325 }
326
327 // Solve LY = b
328 for (int col = 0; col < m; col++) {
329 final T bpCol = bp[col];
330 for (int i = col + 1; i < m; i++) {
331 bp[i] = bp[i].subtract(bpCol.multiply(lu[i][col]));
332 }
333 }
334
335 // Solve UX = Y
336 for (int col = m - 1; col >= 0; col--) {
337 bp[col] = bp[col].divide(lu[col][col]);
338 final T bpCol = bp[col];
339 for (int i = 0; i < col; i++) {
340 bp[i] = bp[i].subtract(bpCol.multiply(lu[i][col]));
341 }
342 }
343
344 return new ArrayFieldVector<T>(bp, false);
345
346 }
347 }
348
349 /** Solve the linear equation A × X = B.
350 * <p>The A matrix is implicit here. It is </p>
351 * @param b right-hand side of the equation A × X = B
352 * @return a vector X such that A × X = B
353 * @exception IllegalArgumentException if matrices dimensions don't match
354 * @exception InvalidMatrixException if decomposed matrix is singular
355 */
356 public ArrayFieldVector<T> solve(ArrayFieldVector<T> b)
357 throws IllegalArgumentException, InvalidMatrixException {
358 return new ArrayFieldVector<T>(solve(b.getDataRef()), false);
359 }
360
361 /** {@inheritDoc} */
362 public FieldMatrix<T> solve(FieldMatrix<T> b)
363 throws IllegalArgumentException, InvalidMatrixException {
364
365 final int m = pivot.length;
366 if (b.getRowDimension() != m) {
367 throw MathRuntimeException.createIllegalArgumentException(
368 "dimensions mismatch: got {0}x{1} but expected {2}x{3}",
369 b.getRowDimension(), b.getColumnDimension(), m, "n");
370 }
371 if (singular) {
372 throw new SingularMatrixException();
373 }
374
375 final int nColB = b.getColumnDimension();
376
377 // Apply permutations to b
378 @SuppressWarnings("unchecked") // field is of type T
379 final T[][] bp = (T[][]) Array.newInstance(field.getZero().getClass(), new int[] { m, nColB });
380 for (int row = 0; row < m; row++) {
381 final T[] bpRow = bp[row];
382 final int pRow = pivot[row];
383 for (int col = 0; col < nColB; col++) {
384 bpRow[col] = b.getEntry(pRow, col);
385 }
386 }
387
388 // Solve LY = b
389 for (int col = 0; col < m; col++) {
390 final T[] bpCol = bp[col];
391 for (int i = col + 1; i < m; i++) {
392 final T[] bpI = bp[i];
393 final T luICol = lu[i][col];
394 for (int j = 0; j < nColB; j++) {
395 bpI[j] = bpI[j].subtract(bpCol[j].multiply(luICol));
396 }
397 }
398 }
399
400 // Solve UX = Y
401 for (int col = m - 1; col >= 0; col--) {
402 final T[] bpCol = bp[col];
403 final T luDiag = lu[col][col];
404 for (int j = 0; j < nColB; j++) {
405 bpCol[j] = bpCol[j].divide(luDiag);
406 }
407 for (int i = 0; i < col; i++) {
408 final T[] bpI = bp[i];
409 final T luICol = lu[i][col];
410 for (int j = 0; j < nColB; j++) {
411 bpI[j] = bpI[j].subtract(bpCol[j].multiply(luICol));
412 }
413 }
414 }
415
416 return new Array2DRowFieldMatrix<T>(bp, false);
417
418 }
419
420 /** {@inheritDoc} */
421 public FieldMatrix<T> getInverse() throws InvalidMatrixException {
422 final int m = pivot.length;
423 final T one = field.getOne();
424 FieldMatrix<T> identity = new Array2DRowFieldMatrix<T>(field, m, m);
425 for (int i = 0; i < m; ++i) {
426 identity.setEntry(i, i, one);
427 }
428 return solve(identity);
429 }
430
431 }
432
433 }