Source code

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.legacy.field.linalg;
019
020import org.apache.commons.numbers.field.Field;
021import org.apache.commons.math4.legacy.linear.SingularMatrixException;
022
023/**
024 * Calculates the LUP-decomposition of a square matrix.
025 *
026 * <p>The LUP-decomposition of a matrix A consists of three matrices
027 * L, U and P that satisfy: PA = LU, L is lower triangular, and U is
028 * upper triangular and P is a permutation matrix. All matrices are
029 * m&times;m.</p>
030 *
031 * <p>Since {@link Field field} elements do not provide an ordering
032 * operator, the permutation matrix is computed here only in order to
033 * avoid a zero pivot element, no attempt is done to get the largest
034 * pivot element.</p>
035 *
036 * @param <T> Type of the field elements.
037 *
038 * @see <a href="http://mathworld.wolfram.com/LUDecomposition.html">MathWorld</a>
039 * @see <a href="http://en.wikipedia.org/wiki/LU_decomposition">Wikipedia</a>
040 *
041 * @since 4.0
042 */
043public final class FieldLUDecomposition<T> {
044    /** Field to which the elements belong. */
045    private final Field<T> field;
046    /** Entries of LU decomposition. */
047    private final FieldDenseMatrix<T> mLU;
048    /** Pivot permutation associated with LU decomposition. */
049    private final int[] pivot;
050    /** Singularity indicator. */
051    private final boolean isSingular;
052    /** Parity of the permutation associated with the LU decomposition. */
053    private final boolean isEven;
054
055    /**
056     * Calculates the LU-decomposition of the given {@code matrix}.
057     *
058     * @param matrix Matrix to decompose.
059     * @throws IllegalArgumentException if the matrix is not square.
060     */
061    private FieldLUDecomposition(FieldDenseMatrix<T> matrix) {
062        matrix.checkMultiply(matrix);
063
064        field = matrix.getField();
065        final int m = matrix.getRowDimension();
066        pivot = new int[m];
067
068        // Initialize permutation array and parity.
069        for (int row = 0; row < m; row++) {
070            pivot[row] = row;
071        }
072        mLU = matrix.copy();
073
074        boolean even = true;
075        boolean singular = false;
076        // Loop over columns.
077        for (int col = 0; col < m; col++) {
078            T sum = field.zero();
079
080            // Upper.
081            for (int row = 0; row < col; row++) {
082                sum = mLU.get(row, col);
083                for (int i = 0; i < row; i++) {
084                    sum = field.subtract(sum,
085                                         field.multiply(mLU.get(row, i),
086                                                        mLU.get(i, col)));
087                }
088                mLU.set(row, col, sum);
089            }
090
091            // Lower.
092            int nonZero = col; // Permutation row.
093            for (int row = col; row < m; row++) {
094                sum = mLU.get(row, col);
095                for (int i = 0; i < col; i++) {
096                    sum = field.subtract(sum,
097                                         field.multiply(mLU.get(row, i),
098                                                        mLU.get(i, col)));
099                }
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}