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    }