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