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
18 package org.apache.commons.math4.legacy.field.linalg;
19
20 import org.apache.commons.numbers.field.Field;
21 import org.apache.commons.math4.legacy.linear.SingularMatrixException;
22
23 /**
24 * Calculates the LUP-decomposition of a square matrix.
25 *
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×m.</p>
30 *
31 * <p>Since {@link Field field} elements do not provide an ordering
32 * operator, the permutation matrix is computed here only in order to
33 * avoid a zero pivot element, no attempt is done to get the largest
34 * pivot element.</p>
35 *
36 * @param <T> Type of the field elements.
37 *
38 * @see <a href="http://mathworld.wolfram.com/LUDecomposition.html">MathWorld</a>
39 * @see <a href="http://en.wikipedia.org/wiki/LU_decomposition">Wikipedia</a>
40 *
41 * @since 4.0
42 */
43 public final class FieldLUDecomposition<T> {
44 /** Field to which the elements belong. */
45 private final Field<T> field;
46 /** Entries of LU decomposition. */
47 private final FieldDenseMatrix<T> mLU;
48 /** Pivot permutation associated with LU decomposition. */
49 private final int[] pivot;
50 /** Singularity indicator. */
51 private final boolean isSingular;
52 /** Parity of the permutation associated with the LU decomposition. */
53 private final boolean isEven;
54
55 /**
56 * Calculates the LU-decomposition of the given {@code matrix}.
57 *
58 * @param matrix Matrix to decompose.
59 * @throws IllegalArgumentException if the matrix is not square.
60 */
61 private FieldLUDecomposition(FieldDenseMatrix<T> matrix) {
62 matrix.checkMultiply(matrix);
63
64 field = matrix.getField();
65 final int m = matrix.getRowDimension();
66 pivot = new int[m];
67
68 // Initialize permutation array and parity.
69 for (int row = 0; row < m; row++) {
70 pivot[row] = row;
71 }
72 mLU = matrix.copy();
73
74 boolean even = true;
75 boolean singular = false;
76 // Loop over columns.
77 for (int col = 0; col < m; col++) {
78 T sum = field.zero();
79
80 // Upper.
81 for (int row = 0; row < col; row++) {
82 sum = mLU.get(row, col);
83 for (int i = 0; i < row; i++) {
84 sum = field.subtract(sum,
85 field.multiply(mLU.get(row, i),
86 mLU.get(i, col)));
87 }
88 mLU.set(row, col, sum);
89 }
90
91 // Lower.
92 int nonZero = col; // Permutation row.
93 for (int row = col; row < m; row++) {
94 sum = mLU.get(row, col);
95 for (int i = 0; i < col; i++) {
96 sum = field.subtract(sum,
97 field.multiply(mLU.get(row, i),
98 mLU.get(i, col)));
99 }
100 mLU.set(row, col, sum);
101
102 if (mLU.get(nonZero, col).equals(field.zero())) {
103 // try to select a better permutation choice
104 ++nonZero;
105 }
106 }
107
108 // Singularity check.
109 if (nonZero >= m) {
110 singular = true;
111 } else {
112 // Pivot if necessary.
113 if (nonZero != col) {
114 T tmp = field.zero();
115 for (int i = 0; i < m; i++) {
116 tmp = mLU.get(nonZero, i);
117 mLU.set(nonZero, i, mLU.get(col, i));
118 mLU.set(col, i, tmp);
119 }
120 int temp = pivot[nonZero];
121 pivot[nonZero] = pivot[col];
122 pivot[col] = temp;
123 even = !even;
124 }
125
126 // Divide the lower elements by the "winning" diagonal element.
127 final T luDiag = mLU.get(col, col);
128 for (int row = col + 1; row < m; row++) {
129 mLU.set(row, col, field.divide(mLU.get(row, col),
130 luDiag));
131 }
132 }
133 }
134
135 isSingular = singular;
136 isEven = even;
137 }
138
139 /**
140 * Factory method.
141 *
142 * @param <T> Type of the field elements.
143 * @param m Matrix to decompose.
144 * @return a new instance.
145 */
146 public static <T> FieldLUDecomposition<T> of(FieldDenseMatrix<T> m) {
147 return new FieldLUDecomposition<>(m);
148 }
149
150 /**
151 * @return {@code true} if the matrix is singular.
152 */
153 public boolean isSingular() {
154 return isSingular;
155 }
156
157 /**
158 * Builds the "L" matrix of the decomposition.
159 *
160 * @return the lower triangular matrix.
161 * @throws SingularMatrixException if the matrix is singular.
162 */
163 public FieldDenseMatrix<T> getL() {
164 if (isSingular) {
165 throw new SingularMatrixException();
166 }
167
168 final int m = pivot.length;
169 final FieldDenseMatrix<T> mL = FieldDenseMatrix.zero(field, m, m);
170 for (int i = 0; i < m; i++) {
171 for (int j = 0; j < i; j++) {
172 mL.set(i, j, mLU.get(i, j));
173 }
174 mL.set(i, i, field.one());
175 }
176
177 return mL;
178 }
179
180 /**
181 * Builds the "U" matrix of the decomposition.
182 *
183 * @return the upper triangular matrix.
184 * @throws SingularMatrixException if the matrix is singular.
185 */
186 public FieldDenseMatrix<T> getU() {
187 if (isSingular) {
188 throw new SingularMatrixException();
189 }
190
191 final int m = pivot.length;
192 final FieldDenseMatrix<T> mU = FieldDenseMatrix.zero(field, m, m);
193 for (int i = 0; i < m; i++) {
194 for (int j = i; j < m; j++) {
195 mU.set(i, j, mLU.get(i, j));
196 }
197 }
198
199 return mU;
200 }
201
202 /**
203 * Builds the "P" matrix.
204 *
205 * <p>P is a matrix with exactly one element set to {@link Field#one() one} in
206 * each row and each column, all other elements being set to {@link Field#zero() zero}.
207 * The positions of the "one" elements are given by the {@link #getPivot()
208 * pivot permutation vector}.</p>
209 * @return the "P" rows permutation matrix.
210 * @throws SingularMatrixException if the matrix is singular.
211 *
212 * @see #getPivot()
213 */
214 public FieldDenseMatrix<T> getP() {
215 if (isSingular) {
216 throw new SingularMatrixException();
217 }
218
219 final int m = pivot.length;
220 final FieldDenseMatrix<T> mP = FieldDenseMatrix.zero(field, m, m);
221
222 for (int i = 0; i < m; i++) {
223 mP.set(i, pivot[i], field.one());
224 }
225
226 return mP;
227 }
228
229 /**
230 * Gets the pivot permutation vector.
231 *
232 * @return the pivot permutation vector.
233 *
234 * @see #getP()
235 */
236 public int[] getPivot() {
237 return pivot.clone();
238 }
239
240 /**
241 * Return the determinant of the matrix.
242 * @return determinant of the matrix
243 */
244 public T getDeterminant() {
245 if (isSingular) {
246 return field.zero();
247 } else {
248 final int m = pivot.length;
249 T determinant = isEven ?
250 field.one() :
251 field.negate(field.one());
252
253 for (int i = 0; i < m; i++) {
254 determinant = field.multiply(determinant,
255 mLU.get(i, i));
256 }
257
258 return determinant;
259 }
260 }
261
262 /**
263 * Creates a solver for finding the solution {@code X} of the linear
264 * system of equations {@code A X = B}.
265 *
266 * @return a solver.
267 * @throws SingularMatrixException if the matrix is singular.
268 */
269 public FieldDecompositionSolver<T> getSolver() {
270 if (isSingular) {
271 throw new SingularMatrixException();
272 }
273
274 return new Solver<>(mLU, pivot);
275 }
276
277 /**
278 * Specialized solver.
279 *
280 * @param <T> Type of the field elements.
281 */
282 private static final class Solver<T> implements FieldDecompositionSolver<T> {
283 /** Field to which the elements belong. */
284 private final Field<T> field;
285 /** LU decomposition. */
286 private final FieldDenseMatrix<T> mLU;
287 /** Pivot permutation associated with LU decomposition. */
288 private final int[] pivot;
289
290 /**
291 * Builds a solver from a LU-decomposed matrix.
292 *
293 * @param mLU LU matrix.
294 * @param pivot Pivot permutation associated with the decomposition.
295 */
296 private Solver(final FieldDenseMatrix<T> mLU,
297 final int[] pivot) {
298 field = mLU.getField();
299 this.mLU = mLU.copy();
300 this.pivot = pivot.clone();
301 }
302
303 /** {@inheritDoc} */
304 @Override
305 public FieldDenseMatrix<T> solve(final FieldDenseMatrix<T> b) {
306 mLU.checkMultiply(b);
307
308 final FieldDenseMatrix<T> bp = b.copy();
309 final int nColB = b.getColumnDimension();
310 final int m = pivot.length;
311
312 // Apply permutations.
313 for (int row = 0; row < m; row++) {
314 final int pRow = pivot[row];
315 for (int col = 0; col < nColB; col++) {
316 bp.set(row, col,
317 b.get(row, col));
318 }
319 }
320
321 // Solve LY = b
322 for (int col = 0; col < m; col++) {
323 for (int i = col + 1; i < m; i++) {
324 for (int j = 0; j < nColB; j++) {
325 bp.set(i, j,
326 field.subtract(bp.get(i, j),
327 field.multiply(bp.get(col, j),
328 mLU.get(i, col))));
329 }
330 }
331 }
332
333 // Solve UX = Y
334 for (int col = m - 1; col >= 0; col--) {
335 for (int j = 0; j < nColB; j++) {
336 bp.set(col, j,
337 field.divide(bp.get(col, j),
338 mLU.get(col, col)));
339 }
340 for (int i = 0; i < col; i++) {
341 for (int j = 0; j < nColB; j++) {
342 bp.set(i, j,
343 field.subtract(bp.get(i, j),
344 field.multiply(bp.get(col, j),
345 mLU.get(i, col))));
346 }
347 }
348 }
349
350 return bp;
351 }
352
353 /** {@inheritDoc} */
354 @Override
355 public FieldDenseMatrix<T> getInverse() {
356 return solve(FieldDenseMatrix.identity(field, pivot.length));
357 }
358 }
359 }