diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math/linear/LUDecompositionImpl.java')
| -rw-r--r-- | src/main/java/org/apache/commons/math/linear/LUDecompositionImpl.java | 423 |
1 files changed, 423 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math/linear/LUDecompositionImpl.java b/src/main/java/org/apache/commons/math/linear/LUDecompositionImpl.java new file mode 100644 index 0000000..e5c9bbf --- /dev/null +++ b/src/main/java/org/apache/commons/math/linear/LUDecompositionImpl.java @@ -0,0 +1,423 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.math.linear; + +import org.apache.commons.math.MathRuntimeException; +import org.apache.commons.math.exception.util.LocalizedFormats; +import org.apache.commons.math.util.FastMath; + +/** + * Calculates the LUP-decomposition of a square matrix. + * <p>The LUP-decomposition of a matrix A consists of three matrices + * L, U and P that satisfy: PA = LU, L is lower triangular, and U is + * upper triangular and P is a permutation matrix. All matrices are + * m×m.</p> + * <p>As shown by the presence of the P matrix, this decomposition is + * implemented using partial pivoting.</p> + * + * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 août 2010) $ + * @since 2.0 + */ +public class LUDecompositionImpl implements LUDecomposition { + + /** Default bound to determine effective singularity in LU decomposition */ + private static final double DEFAULT_TOO_SMALL = 10E-12; + + /** Entries of LU decomposition. */ + private double lu[][]; + + /** Pivot permutation associated with LU decomposition */ + private int[] pivot; + + /** Parity of the permutation associated with the LU decomposition */ + private boolean even; + + /** Singularity indicator. */ + private boolean singular; + + /** Cached value of L. */ + private RealMatrix cachedL; + + /** Cached value of U. */ + private RealMatrix cachedU; + + /** Cached value of P. */ + private RealMatrix cachedP; + + /** + * Calculates the LU-decomposition of the given matrix. + * @param matrix The matrix to decompose. + * @exception InvalidMatrixException if matrix is not square + */ + public LUDecompositionImpl(RealMatrix matrix) + throws InvalidMatrixException { + this(matrix, DEFAULT_TOO_SMALL); + } + + /** + * Calculates the LU-decomposition of the given matrix. + * @param matrix The matrix to decompose. + * @param singularityThreshold threshold (based on partial row norm) + * under which a matrix is considered singular + * @exception NonSquareMatrixException if matrix is not square + */ + public LUDecompositionImpl(RealMatrix matrix, double singularityThreshold) + throws NonSquareMatrixException { + + if (!matrix.isSquare()) { + throw new NonSquareMatrixException(matrix.getRowDimension(), matrix.getColumnDimension()); + } + + final int m = matrix.getColumnDimension(); + lu = matrix.getData(); + pivot = new int[m]; + cachedL = null; + cachedU = null; + cachedP = null; + + // Initialize permutation array and parity + for (int row = 0; row < m; row++) { + pivot[row] = row; + } + even = true; + singular = false; + + // Loop over columns + for (int col = 0; col < m; col++) { + + double sum = 0; + + // upper + for (int row = 0; row < col; row++) { + final double[] luRow = lu[row]; + sum = luRow[col]; + for (int i = 0; i < row; i++) { + sum -= luRow[i] * lu[i][col]; + } + luRow[col] = sum; + } + + // lower + int max = col; // permutation row + double largest = Double.NEGATIVE_INFINITY; + for (int row = col; row < m; row++) { + final double[] luRow = lu[row]; + sum = luRow[col]; + for (int i = 0; i < col; i++) { + sum -= luRow[i] * lu[i][col]; + } + luRow[col] = sum; + + // maintain best permutation choice + if (FastMath.abs(sum) > largest) { + largest = FastMath.abs(sum); + max = row; + } + } + + // Singularity check + if (FastMath.abs(lu[max][col]) < singularityThreshold) { + singular = true; + return; + } + + // Pivot if necessary + if (max != col) { + double tmp = 0; + final double[] luMax = lu[max]; + final double[] luCol = lu[col]; + for (int i = 0; i < m; i++) { + tmp = luMax[i]; + luMax[i] = luCol[i]; + luCol[i] = tmp; + } + int temp = pivot[max]; + pivot[max] = pivot[col]; + pivot[col] = temp; + even = !even; + } + + // Divide the lower elements by the "winning" diagonal elt. + final double luDiag = lu[col][col]; + for (int row = col + 1; row < m; row++) { + lu[row][col] /= luDiag; + } + } + + } + + /** {@inheritDoc} */ + public RealMatrix getL() { + if ((cachedL == null) && !singular) { + final int m = pivot.length; + cachedL = MatrixUtils.createRealMatrix(m, m); + for (int i = 0; i < m; ++i) { + final double[] luI = lu[i]; + for (int j = 0; j < i; ++j) { + cachedL.setEntry(i, j, luI[j]); + } + cachedL.setEntry(i, i, 1.0); + } + } + return cachedL; + } + + /** {@inheritDoc} */ + public RealMatrix getU() { + if ((cachedU == null) && !singular) { + final int m = pivot.length; + cachedU = MatrixUtils.createRealMatrix(m, m); + for (int i = 0; i < m; ++i) { + final double[] luI = lu[i]; + for (int j = i; j < m; ++j) { + cachedU.setEntry(i, j, luI[j]); + } + } + } + return cachedU; + } + + /** {@inheritDoc} */ + public RealMatrix getP() { + if ((cachedP == null) && !singular) { + final int m = pivot.length; + cachedP = MatrixUtils.createRealMatrix(m, m); + for (int i = 0; i < m; ++i) { + cachedP.setEntry(i, pivot[i], 1.0); + } + } + return cachedP; + } + + /** {@inheritDoc} */ + public int[] getPivot() { + return pivot.clone(); + } + + /** {@inheritDoc} */ + public double getDeterminant() { + if (singular) { + return 0; + } else { + final int m = pivot.length; + double determinant = even ? 1 : -1; + for (int i = 0; i < m; i++) { + determinant *= lu[i][i]; + } + return determinant; + } + } + + /** {@inheritDoc} */ + public DecompositionSolver getSolver() { + return new Solver(lu, pivot, singular); + } + + /** Specialized solver. */ + private static class Solver implements DecompositionSolver { + + /** Entries of LU decomposition. */ + private final double lu[][]; + + /** Pivot permutation associated with LU decomposition. */ + private final int[] pivot; + + /** Singularity indicator. */ + private final boolean singular; + + /** + * Build a solver from decomposed matrix. + * @param lu entries of LU decomposition + * @param pivot pivot permutation associated with LU decomposition + * @param singular singularity indicator + */ + private Solver(final double[][] lu, final int[] pivot, final boolean singular) { + this.lu = lu; + this.pivot = pivot; + this.singular = singular; + } + + /** {@inheritDoc} */ + public boolean isNonSingular() { + return !singular; + } + + /** {@inheritDoc} */ + public double[] solve(double[] b) + throws IllegalArgumentException, InvalidMatrixException { + + final int m = pivot.length; + if (b.length != m) { + throw MathRuntimeException.createIllegalArgumentException( + LocalizedFormats.VECTOR_LENGTH_MISMATCH, b.length, m); + } + if (singular) { + throw new SingularMatrixException(); + } + + final double[] bp = new double[m]; + + // Apply permutations to b + for (int row = 0; row < m; row++) { + bp[row] = b[pivot[row]]; + } + + // Solve LY = b + for (int col = 0; col < m; col++) { + final double bpCol = bp[col]; + for (int i = col + 1; i < m; i++) { + bp[i] -= bpCol * lu[i][col]; + } + } + + // Solve UX = Y + for (int col = m - 1; col >= 0; col--) { + bp[col] /= lu[col][col]; + final double bpCol = bp[col]; + for (int i = 0; i < col; i++) { + bp[i] -= bpCol * lu[i][col]; + } + } + + return bp; + + } + + /** {@inheritDoc} */ + public RealVector solve(RealVector b) + throws IllegalArgumentException, InvalidMatrixException { + try { + return solve((ArrayRealVector) b); + } catch (ClassCastException cce) { + + final int m = pivot.length; + if (b.getDimension() != m) { + throw MathRuntimeException.createIllegalArgumentException( + LocalizedFormats.VECTOR_LENGTH_MISMATCH, b.getDimension(), m); + } + if (singular) { + throw new SingularMatrixException(); + } + + final double[] bp = new double[m]; + + // Apply permutations to b + for (int row = 0; row < m; row++) { + bp[row] = b.getEntry(pivot[row]); + } + + // Solve LY = b + for (int col = 0; col < m; col++) { + final double bpCol = bp[col]; + for (int i = col + 1; i < m; i++) { + bp[i] -= bpCol * lu[i][col]; + } + } + + // Solve UX = Y + for (int col = m - 1; col >= 0; col--) { + bp[col] /= lu[col][col]; + final double bpCol = bp[col]; + for (int i = 0; i < col; i++) { + bp[i] -= bpCol * lu[i][col]; + } + } + + return new ArrayRealVector(bp, false); + + } + } + + /** Solve the linear equation A × X = B. + * <p>The A matrix is implicit here. It is </p> + * @param b right-hand side of the equation A × X = B + * @return a vector X such that A × X = B + * @exception IllegalArgumentException if matrices dimensions don't match + * @exception InvalidMatrixException if decomposed matrix is singular + */ + public ArrayRealVector solve(ArrayRealVector b) + throws IllegalArgumentException, InvalidMatrixException { + return new ArrayRealVector(solve(b.getDataRef()), false); + } + + /** {@inheritDoc} */ + public RealMatrix solve(RealMatrix b) + throws IllegalArgumentException, InvalidMatrixException { + + final int m = pivot.length; + if (b.getRowDimension() != m) { + throw MathRuntimeException.createIllegalArgumentException( + LocalizedFormats.DIMENSIONS_MISMATCH_2x2, + b.getRowDimension(), b.getColumnDimension(), m, "n"); + } + if (singular) { + throw new SingularMatrixException(); + } + + final int nColB = b.getColumnDimension(); + + // Apply permutations to b + final double[][] bp = new double[m][nColB]; + for (int row = 0; row < m; row++) { + final double[] bpRow = bp[row]; + final int pRow = pivot[row]; + for (int col = 0; col < nColB; col++) { + bpRow[col] = b.getEntry(pRow, col); + } + } + + // Solve LY = b + for (int col = 0; col < m; col++) { + final double[] bpCol = bp[col]; + for (int i = col + 1; i < m; i++) { + final double[] bpI = bp[i]; + final double luICol = lu[i][col]; + for (int j = 0; j < nColB; j++) { + bpI[j] -= bpCol[j] * luICol; + } + } + } + + // Solve UX = Y + for (int col = m - 1; col >= 0; col--) { + final double[] bpCol = bp[col]; + final double luDiag = lu[col][col]; + for (int j = 0; j < nColB; j++) { + bpCol[j] /= luDiag; + } + for (int i = 0; i < col; i++) { + final double[] bpI = bp[i]; + final double luICol = lu[i][col]; + for (int j = 0; j < nColB; j++) { + bpI[j] -= bpCol[j] * luICol; + } + } + } + + return new Array2DRowRealMatrix(bp, false); + + } + + /** {@inheritDoc} */ + public RealMatrix getInverse() throws InvalidMatrixException { + return solve(MatrixUtils.createRealIdentityMatrix(pivot.length)); + } + + } + +} |