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