commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r687519 - in /commons/proper/math/branches/MATH_2_0/src: java/org/apache/commons/math/linear/LUDecomposition.java java/org/apache/commons/math/linear/LUDecompositionImpl.java test/org/apache/commons/math/linear/LUDecompositionImplTest.java
Date Thu, 21 Aug 2008 00:19:48 GMT
Author: luc
Date: Wed Aug 20 17:19:48 2008
New Revision: 687519

URL: http://svn.apache.org/viewvc?rev=687519&view=rev
Log:
added JAMA-like LU-decomposition

Added:
    commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/LUDecomposition.java
  (with props)
    commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/LUDecompositionImpl.java
  (with props)
    commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/linear/LUDecompositionImplTest.java
  (with props)

Added: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/LUDecomposition.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/LUDecomposition.java?rev=687519&view=auto
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/LUDecomposition.java
(added)
+++ commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/LUDecomposition.java
Wed Aug 20 17:19:48 2008
@@ -0,0 +1,85 @@
+/*
+ * 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;
+
+/**
+ * An interface to classes that implement a algorithm to calculate the 
+ * LU-decomposition of a real matrix.
+ * <p>The LU-decomposition of matrix A is a set of three matrices: P, L and U
+ * such that P&times;A = L&times;U. P is a rows permutation matrix that is used
+ * to rearrange the rows of A before so that it can be decomposed. L is a lower
+ * triangular matrix with unit diagonal terms and U is an upper triangular matrix.</p>
+ * <p>This interface is similar to the class with similar name from the now defunct
+ * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library, with the
+ * exception of the <code>det</code> method which has been renamed as {@link
+ * #getDeterminant() getDeterminant}.</p>
+ *   
+ * @see <a href="http://mathworld.wolfram.com/LUDecomposition.html">MathWorld</a>
+ * @see <a href="http://en.wikipedia.org/wiki/LU_decomposition">Wikipedia</a>
+ * @version $Revision$ $Date$
+ * @since 2.0
+ */
+public interface LUDecomposition extends DecompositionSolver {
+
+    /**
+     * Returns the matrix L of the decomposition. 
+     * <p>L is an lower-triangular matrix</p>
+     * @return the L matrix (or null if decomposed matrix is singular)
+     */
+    RealMatrix getL();
+
+    /**
+     * Returns the matrix U of the decomposition. 
+     * <p>U is an upper-triangular matrix</p>
+     * @return the U matrix (or null if decomposed matrix is singular)
+     */
+    RealMatrix getU();
+
+    /**
+     * Returns the P rows permutation matrix.
+     * <p>P is a sparse matrix with exactly one element set to 1.0 in
+     * each row and each column, all other elements being set to 0.0.</p>
+     * <p>The positions of the 1 elements are given by the {@link #getPivot()
+     * pivot permutation vector}.</p>
+     * @return the P rows permutation matrix (or null if decomposed matrix is singular)
+     * @see #getPivot()
+     */
+    RealMatrix getP();
+
+    /**
+     * Returns the pivot permutation vector.
+     * @return the pivot permutation vector
+     * @see #getPermutation()
+     */
+    int[] getPivot();
+
+    /**
+     * Check if the decomposed matrix is non-singular.
+     * @return true if the decomposed matrix is non-singular
+     * @see #getDeterminant()
+     */
+    boolean isNonSingular();
+
+    /**
+     * Return the determinant of the matrix
+     * @return determinant of the matrix
+     * @see #isNonSingular()
+     */
+    double getDeterminant();
+
+}

Propchange: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/LUDecomposition.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/LUDecomposition.java
------------------------------------------------------------------------------
    svn:keywords = Author Date Id Revision

Added: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/LUDecompositionImpl.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/LUDecompositionImpl.java?rev=687519&view=auto
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/LUDecompositionImpl.java
(added)
+++ commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/LUDecompositionImpl.java
Wed Aug 20 17:19:48 2008
@@ -0,0 +1,394 @@
+/*
+ * 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;
+
+/**
+ * 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: A = LUP, L is lower triangular, and U is
+ * upper triangular and P is a permutation matrix. All matrices are
+ * m&times;m.</p>
+ * <p>As shown by the presence of the P matrix, this decomposition is
+ * implemented using partial pivoting.</p>
+ *
+ * @version $Revision$ $Date$
+ * @since 2.0
+ */
+public class LUDecompositionImpl implements LUDecomposition {
+
+    /** Serializable version identifier. */
+    private static final long serialVersionUID = -1606789599960880183L;
+
+    /** Bound to determine effective singularity in LU decomposition */
+    private final double singularityThreshold;
+
+    /** Size of the matrix. */
+    private final int m;
+
+    /** Entries of LU decomposition. */
+    private final double lu[][];
+
+    /** Pivot permutation associated with LU decomposition */
+    private final int[] pivot;
+
+    /** Parity of the permutation associated with the LU decomposition */
+    private int parity;
+
+    /** 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;
+
+    /** Default bound to determine effective singularity in LU decomposition */
+    private static final double DEFAULT_TOO_SMALL = 10E-12;
+
+    /**
+     * 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 InvalidMatrixException if matrix is not square
+     */
+    public LUDecompositionImpl(RealMatrix matrix, double singularityThreshold)
+        throws InvalidMatrixException {
+        if (!matrix.isSquare()) {
+            throw new InvalidMatrixException("LU decomposition requires that the matrix be
square.");
+        }
+        this.singularityThreshold = singularityThreshold;
+        m = matrix.getColumnDimension();
+        lu = matrix.getData();
+        pivot = new int[m];
+        cachedL = null;
+        cachedU = null;
+        cachedP = null;
+
+        // perform decomposition
+        luDecompose();
+
+    }
+
+    /** {@inheritDoc} */
+    public RealMatrix getL() {
+        if ((cachedL == null) && !singular) {
+            final double[][] lData = new double[m][m];
+            for (int i = 0; i < m; ++i) {
+                System.arraycopy(lu[i], 0, lData[i], 0, i);
+                lData[i][i] = 1.0;
+            }
+            cachedL = new RealMatrixImpl(lData, false);
+        }
+        return cachedL;
+    }
+
+    /** {@inheritDoc} */
+    public RealMatrix getU() {
+        if ((cachedU == null) && !singular) {
+            final double[][] uData = new double[m][m];
+            for (int i = 0; i < m; ++i) {
+                System.arraycopy(lu[i], i, uData[i], i, m - i);
+            }
+            cachedU = new RealMatrixImpl(uData, false);
+        }
+        return cachedU;
+    }
+
+    /** {@inheritDoc} */
+    public RealMatrix getP() {
+        if ((cachedP == null) && !singular) {
+            final double[][] pData = new double[m][m];
+            for (int i = 0; i < m; ++i) {
+                pData[i][pivot[i]] = 1.0;
+            }
+            cachedP = new RealMatrixImpl(pData, false);
+        }
+        return cachedP;
+    }
+
+    /** {@inheritDoc} */
+    public int[] getPivot() {
+        return pivot.clone();
+    }
+
+    /** {@inheritDoc} */
+    public boolean isNonSingular() {
+        return !singular;
+    }
+
+    /** {@inheritDoc} */
+    public double getDeterminant() {
+        if (singular) {
+            return 0;
+        } else {
+            double determinant = parity;
+            for (int i = 0; i < m; i++) {
+                determinant *= lu[i][i];
+            }
+            return determinant;
+        }
+    }
+
+    /** {@inheritDoc} */
+    public double[] solve(double[] b)
+        throws IllegalArgumentException, InvalidMatrixException {
+
+        if (b.length != m) {
+            throw new IllegalArgumentException("constant vector has wrong length");
+        }
+        if (singular) {
+            throw new InvalidMatrixException("Matrix is singular.");
+        }
+
+        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++) {
+            for (int i = col + 1; i < m; i++) {
+                bp[i] -= bp[col] * lu[i][col];
+            }
+        }
+
+        // Solve UX = Y
+        for (int col = m - 1; col >= 0; col--) {
+            bp[col] /= lu[col][col];
+            for (int i = 0; i < col; i++) {
+                bp[i] -= bp[col] * lu[i][col];
+            }
+        }
+
+        return bp;
+
+    }
+
+    /** {@inheritDoc} */
+    public RealVector solve(RealVector b)
+        throws IllegalArgumentException, InvalidMatrixException {
+        try {
+            return solve((RealVectorImpl) b);
+        } catch (ClassCastException cce) {
+
+            if (b.getDimension() != m) {
+                throw new IllegalArgumentException("constant vector has wrong length");
+            }
+            if (singular) {
+                throw new InvalidMatrixException("Matrix is singular.");
+            }
+
+            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++) {
+                for (int i = col + 1; i < m; i++) {
+                    bp[i] -= bp[col] * lu[i][col];
+                }
+            }
+
+            // Solve UX = Y
+            for (int col = m - 1; col >= 0; col--) {
+                bp[col] /= lu[col][col];
+                for (int i = 0; i < col; i++) {
+                    bp[i] -= bp[col] * lu[i][col];
+                }
+            }
+
+            return new RealVectorImpl(bp, false);
+
+        }
+    }
+
+    /** Solve the linear equation A &times; X = B.
+     * <p>The A matrix is implicit here. It is </p>
+     * @param b right-hand side of the equation A &times; X = B
+     * @return a vector X such that A &times; X = B
+     * @throws IllegalArgumentException if matrices dimensions don't match
+     * @throws InvalidMatrixException if decomposed matrix is singular
+     */
+    public RealVectorImpl solve(RealVectorImpl b)
+        throws IllegalArgumentException, InvalidMatrixException {
+        return new RealVectorImpl(solve(b.getDataRef()), false);
+    }
+
+    /** {@inheritDoc} */
+    public RealMatrix solve(RealMatrix b)
+        throws IllegalArgumentException, InvalidMatrixException {
+        if (b.getRowDimension() != m) {
+            throw new IllegalArgumentException("Incorrect row dimension");
+        }
+        if (singular) {
+            throw new InvalidMatrixException("Matrix is singular.");
+        }
+
+        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 RealMatrixImpl(bp, false);
+
+    }
+
+    /**
+     * Computes a new
+     * <a href="http://www.math.gatech.edu/~bourbaki/math2601/Web-notes/2num.pdf">
+     * LU decomposition</a> for this matrix, storing the result for use by other methods.
+     * <p>
+     * <strong>Implementation Note</strong>:<br>
+     * Uses <a href="http://www.damtp.cam.ac.uk/user/fdl/people/sd/lectures/nummeth98/linear.htm">
+     * Crout's algorithm</a>, with partial pivoting.</p>
+     * <p>
+     * <strong>Usage Note</strong>:<br>
+     * This method should rarely be invoked directly. Its only use is
+     * to force recomputation of the LU decomposition when changes have been
+     * made to the underlying data using direct array references. Changes
+     * made using setXxx methods will trigger recomputation when needed
+     * automatically.</p>
+     */
+    private void luDecompose() {
+
+        // Initialize permutation array and parity
+        for (int row = 0; row < m; row++) {
+            pivot[row] = row;
+        }
+        parity = 1;
+        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 (Math.abs(sum) > largest) {
+                    largest = Math.abs(sum);
+                    max = row;
+                }
+            }
+
+            // Singularity check
+            if (Math.abs(lu[max][col]) < singularityThreshold) {
+                singular = true;
+                return;
+            }
+
+            // Pivot if necessary
+            if (max != col) {
+                double tmp = 0;
+                for (int i = 0; i < m; i++) {
+                    tmp = lu[max][i];
+                    lu[max][i] = lu[col][i];
+                    lu[col][i] = tmp;
+                }
+                int temp = pivot[max];
+                pivot[max] = pivot[col];
+                pivot[col] = temp;
+                parity = -parity;
+            }
+
+            // 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;
+            }
+        }
+    }
+
+}

Propchange: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/LUDecompositionImpl.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/linear/LUDecompositionImpl.java
------------------------------------------------------------------------------
    svn:keywords = Author Date Id Revision

Added: commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/linear/LUDecompositionImplTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/linear/LUDecompositionImplTest.java?rev=687519&view=auto
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/linear/LUDecompositionImplTest.java
(added)
+++ commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/linear/LUDecompositionImplTest.java
Wed Aug 20 17:19:48 2008
@@ -0,0 +1,422 @@
+/*
+ * 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 junit.framework.Test;
+import junit.framework.TestCase;
+import junit.framework.TestSuite;
+
+public class LUDecompositionImplTest extends TestCase {
+    private double[][] testData = {
+            { 1.0, 2.0, 3.0},
+            { 2.0, 5.0, 3.0},
+            { 1.0, 0.0, 8.0}
+    };
+    private double[][] testDataMinus = {
+            { -1.0, -2.0, -3.0},
+            { -2.0, -5.0, -3.0},
+            { -1.0,  0.0, -8.0}
+    };
+    private double[][] luData = {
+            { 2.0, 3.0, 3.0 },
+            { 0.0, 5.0, 7.0 },
+            { 6.0, 9.0, 8.0 }
+    };
+    
+    // singular matrices
+    private double[][] singular = {
+            { 2.0, 3.0 },
+            { 2.0, 3.0 }
+    };
+    private double[][] bigSingular = {
+            { 1.0, 2.0,   3.0,    4.0 },
+            { 2.0, 5.0,   3.0,    4.0 },
+            { 7.0, 3.0, 256.0, 1930.0 },
+            { 3.0, 7.0,   6.0,    8.0 }
+    }; // 4th row = 1st + 2nd
+
+    private static final double entryTolerance = 10e-16;
+
+    private static final double normTolerance = 10e-14;
+
+    public LUDecompositionImplTest(String name) {
+        super(name);
+    }
+
+    public static Test suite() {
+        TestSuite suite = new TestSuite(LUDecompositionImplTest.class);
+        suite.setName("LUDecompositionImpl Tests");
+        return suite;
+    }
+
+    /** test dimensions */
+    public void testDimensions() {
+        RealMatrixImpl matrix = new RealMatrixImpl(testData, false);
+        LUDecomposition LU = new LUDecompositionImpl(matrix);
+        assertEquals(testData.length, LU.getL().getRowDimension());
+        assertEquals(testData.length, LU.getL().getColumnDimension());
+        assertEquals(testData.length, LU.getU().getRowDimension());
+        assertEquals(testData.length, LU.getU().getColumnDimension());
+        assertEquals(testData.length, LU.getP().getRowDimension());
+        assertEquals(testData.length, LU.getP().getColumnDimension());
+
+    }
+
+    /** test non-square matrix */
+    public void testNonSquare() {
+        try {
+            new LUDecompositionImpl(new RealMatrixImpl(new double[3][2], false));
+        } catch (InvalidMatrixException ime) {
+            // expected behavior
+        } catch (Exception e) {
+            fail("wrong exception caught");
+        }
+    }
+
+    /** test PA = LU */
+    public void testPAEqualLU() {
+        RealMatrix matrix = new RealMatrixImpl(testData, false);
+        LUDecomposition lu = new LUDecompositionImpl(matrix);
+        RealMatrix l = lu.getL();
+        RealMatrix u = lu.getU();
+        RealMatrix p = lu.getP();
+        double norm = l.multiply(u).subtract(p.multiply(matrix)).getNorm();
+        assertEquals(0, norm, normTolerance);
+
+        matrix = new RealMatrixImpl(testDataMinus, false);
+        lu = new LUDecompositionImpl(matrix);
+        l = lu.getL();
+        u = lu.getU();
+        p = lu.getP();
+        norm = l.multiply(u).subtract(p.multiply(matrix)).getNorm();
+        assertEquals(0, norm, normTolerance);
+
+        matrix = MatrixUtils.createRealIdentityMatrix(17);
+        lu = new LUDecompositionImpl(matrix);
+        l = lu.getL();
+        u = lu.getU();
+        p = lu.getP();
+        norm = l.multiply(u).subtract(p.multiply(matrix)).getNorm();
+        assertEquals(0, norm, normTolerance);
+
+        matrix = new RealMatrixImpl(singular, false);
+        lu = new LUDecompositionImpl(matrix);
+        assertFalse(lu.isNonSingular());
+        assertNull(lu.getL());
+        assertNull(lu.getU());
+        assertNull(lu.getP());
+
+        matrix = new RealMatrixImpl(bigSingular, false);
+        lu = new LUDecompositionImpl(matrix);
+        assertFalse(lu.isNonSingular());
+        assertNull(lu.getL());
+        assertNull(lu.getU());
+        assertNull(lu.getP());
+
+    }
+
+    /** test that L is lower triangular with unit diagonal */
+    public void testLLowerTriangular() {
+        RealMatrixImpl matrix = new RealMatrixImpl(testData, false);
+        RealMatrix l = new LUDecompositionImpl(matrix).getL();
+        for (int i = 0; i < l.getRowDimension(); i++) {
+            assertEquals(l.getEntry(i, i), 1, entryTolerance);
+            for (int j = i + 1; j < l.getColumnDimension(); j++) {
+                assertEquals(l.getEntry(i, j), 0, entryTolerance);
+            }
+        }
+    }
+
+    /** test that U is upper triangular */
+    public void testUUpperTriangular() {
+        RealMatrixImpl matrix = new RealMatrixImpl(testData, false);
+        RealMatrix u = new LUDecompositionImpl(matrix).getU();
+        for (int i = 0; i < u.getRowDimension(); i++) {
+            for (int j = 0; j < i; j++) {
+                assertEquals(u.getEntry(i, j), 0, entryTolerance);
+            }
+        }
+    }
+
+    /** test that P is a permutation matrix */
+    public void testPPermutation() {
+        RealMatrixImpl matrix = new RealMatrixImpl(testData, false);
+        RealMatrix p   = new LUDecompositionImpl(matrix).getP();
+
+        RealMatrix ppT = p.multiply(p.transpose());
+        RealMatrix id  = MatrixUtils.createRealIdentityMatrix(p.getRowDimension());
+        assertEquals(0, ppT.subtract(id).getNorm(), normTolerance);
+
+        for (int i = 0; i < p.getRowDimension(); i++) {
+            int zeroCount  = 0;
+            int oneCount   = 0;
+            int otherCount = 0;
+            for (int j = 0; j < p.getColumnDimension(); j++) {
+                final double e = p.getEntry(i, j);
+                if (e == 0) {
+                    ++zeroCount;
+                } else if (e == 1) {
+                    ++oneCount;
+                } else {
+                    ++otherCount;
+                }
+            }
+            assertEquals(p.getColumnDimension() - 1, zeroCount);
+            assertEquals(1, oneCount);
+            assertEquals(0, otherCount);
+        }
+
+        for (int j = 0; j < p.getColumnDimension(); j++) {
+            int zeroCount  = 0;
+            int oneCount   = 0;
+            int otherCount = 0;
+            for (int i = 0; i < p.getRowDimension(); i++) {
+                final double e = p.getEntry(i, j);
+                if (e == 0) {
+                    ++zeroCount;
+                } else if (e == 1) {
+                    ++oneCount;
+                } else {
+                    ++otherCount;
+                }
+            }
+            assertEquals(p.getRowDimension() - 1, zeroCount);
+            assertEquals(1, oneCount);
+            assertEquals(0, otherCount);
+        }
+
+    }
+
+
+    /** test singular */
+    public void testSingular() {
+        LUDecomposition lu =
+            new LUDecompositionImpl(new RealMatrixImpl(testData, false));
+        assertTrue(lu.isNonSingular());
+        lu = new LUDecompositionImpl(new RealMatrixImpl(singular, false));
+        assertFalse(lu.isNonSingular());
+        lu = new LUDecompositionImpl(new RealMatrixImpl(bigSingular, false));
+        assertFalse(lu.isNonSingular());
+    }
+
+    /** test solve dimension errors */
+    public void testSolveDimensionErrors() {
+        LUDecomposition lu =
+            new LUDecompositionImpl(new RealMatrixImpl(testData, false));
+        RealMatrix b = new RealMatrixImpl(new double[2][2]);
+        try {
+            lu.solve(b);
+            fail("an exception should have been thrown");
+        } catch (IllegalArgumentException iae) {
+            // expected behavior
+        } catch (Exception e) {
+            fail("wrong exception caught");
+        }
+        try {
+            lu.solve(b.getColumn(0));
+            fail("an exception should have been thrown");
+        } catch (IllegalArgumentException iae) {
+            // expected behavior
+        } catch (Exception e) {
+            fail("wrong exception caught");
+        }
+        try {
+            lu.solve(new RealVectorImplTest.RealVectorTestImpl(b.getColumn(0)));
+            fail("an exception should have been thrown");
+        } catch (IllegalArgumentException iae) {
+            // expected behavior
+        } catch (Exception e) {
+            fail("wrong exception caught");
+        }
+    }
+
+    /** test solve singularity errors */
+    public void testSolveSingularityErrors() {
+        LUDecomposition lu =
+            new LUDecompositionImpl(new RealMatrixImpl(singular, false));
+        RealMatrix b = new RealMatrixImpl(new double[2][2]);
+        try {
+            lu.solve(b);
+            fail("an exception should have been thrown");
+        } catch (InvalidMatrixException ime) {
+            // expected behavior
+        } catch (Exception e) {
+            fail("wrong exception caught");
+        }
+        try {
+            lu.solve(b.getColumn(0));
+            fail("an exception should have been thrown");
+        } catch (InvalidMatrixException ime) {
+            // expected behavior
+        } catch (Exception e) {
+            fail("wrong exception caught");
+        }
+        try {
+            lu.solve(b.getColumnVector(0));
+            fail("an exception should have been thrown");
+        } catch (InvalidMatrixException ime) {
+            // expected behavior
+        } catch (Exception e) {
+            fail("wrong exception caught");
+        }
+        try {
+            lu.solve(new RealVectorImplTest.RealVectorTestImpl(b.getColumn(0)));
+            fail("an exception should have been thrown");
+        } catch (InvalidMatrixException ime) {
+            // expected behavior
+        } catch (Exception e) {
+            fail("wrong exception caught");
+        }
+    }
+
+    /** test solve */
+    public void testSolve() {
+        LUDecomposition lu =
+            new LUDecompositionImpl(new RealMatrixImpl(testData, false));
+        RealMatrix b = new RealMatrixImpl(new double[][] {
+                { 1, 0 }, { 2, -5 }, { 3, 1 }
+        });
+        RealMatrix xRef = new RealMatrixImpl(new double[][] {
+                { 19, -71 }, { -6, 22 }, { -2, 9 }
+        });
+
+        // using RealMatrix
+        assertEquals(0, lu.solve(b).subtract(xRef).getNorm(), 1.0e-13);
+
+        // using double[]
+        for (int i = 0; i < b.getColumnDimension(); ++i) {
+            assertEquals(0,
+                         new RealVectorImpl(lu.solve(b.getColumn(i))).subtract(xRef.getColumnVector(i)).getNorm(),
+                         1.0e-13);
+        }
+
+        // using RealVectorImpl
+        for (int i = 0; i < b.getColumnDimension(); ++i) {
+            assertEquals(0,
+                         lu.solve(b.getColumnVector(i)).subtract(xRef.getColumnVector(i)).getNorm(),
+                         1.0e-13);
+        }
+
+        // using RealVector with an alternate implementation
+        for (int i = 0; i < b.getColumnDimension(); ++i) {
+            RealVectorImplTest.RealVectorTestImpl v =
+                new RealVectorImplTest.RealVectorTestImpl(b.getColumn(i));
+            assertEquals(0,
+                         lu.solve(v).subtract(xRef.getColumnVector(i)).getNorm(),
+                         1.0e-13);
+        }
+
+    }
+
+    /** test matrices values */
+    public void testMatricesValues1() {
+       LUDecomposition lu =
+            new LUDecompositionImpl(new RealMatrixImpl(testData, false));
+        RealMatrix lRef = new RealMatrixImpl(new double[][] {
+                { 1.0, 0.0, 0.0 },
+                { 0.5, 1.0, 0.0 },
+                { 0.5, 0.2, 1.0 }
+        });
+        RealMatrix uRef = new RealMatrixImpl(new double[][] {
+                { 2.0,  5.0, 3.0 },
+                { 0.0, -2.5, 6.5 },
+                { 0.0,  0.0, 0.2 }
+        });
+        RealMatrix pRef = new RealMatrixImpl(new double[][] {
+                { 0.0, 1.0, 0.0 },
+                { 0.0, 0.0, 1.0 },
+                { 1.0, 0.0, 0.0 }
+        });
+        int[] pivotRef = { 1, 2, 0 };
+
+        // check values against known references
+        RealMatrix l = lu.getL();
+        assertEquals(0, l.subtract(lRef).getNorm(), 1.0e-13);
+        RealMatrix u = lu.getU();
+        assertEquals(0, u.subtract(uRef).getNorm(), 1.0e-13);
+        RealMatrix p = lu.getP();
+        assertEquals(0, p.subtract(pRef).getNorm(), 1.0e-13);
+        int[] pivot = lu.getPivot();
+        for (int i = 0; i < pivotRef.length; ++i) {
+            assertEquals(pivotRef[i], pivot[i]);
+        }
+
+        // check the same cached instance is returned the second time
+        assertTrue(l == lu.getL());
+        assertTrue(u == lu.getU());
+        assertTrue(p == lu.getP());
+        
+    }
+
+    /** test matrices values */
+    public void testMatricesValues2() {
+       LUDecomposition lu =
+            new LUDecompositionImpl(new RealMatrixImpl(luData, false));
+        RealMatrix lRef = new RealMatrixImpl(new double[][] {
+                {    1.0,    0.0, 0.0 },
+                {    0.0,    1.0, 0.0 },
+                { 1.0 / 3.0, 0.0, 1.0 }
+        });
+        RealMatrix uRef = new RealMatrixImpl(new double[][] {
+                { 6.0, 9.0,    8.0    },
+                { 0.0, 5.0,    7.0    },
+                { 0.0, 0.0, 1.0 / 3.0 }
+        });
+        RealMatrix pRef = new RealMatrixImpl(new double[][] {
+                { 0.0, 0.0, 1.0 },
+                { 0.0, 1.0, 0.0 },
+                { 1.0, 0.0, 0.0 }
+        });
+        int[] pivotRef = { 2, 1, 0 };
+
+        // check values against known references
+        RealMatrix l = lu.getL();
+        assertEquals(0, l.subtract(lRef).getNorm(), 1.0e-13);
+        RealMatrix u = lu.getU();
+        assertEquals(0, u.subtract(uRef).getNorm(), 1.0e-13);
+        RealMatrix p = lu.getP();
+        assertEquals(0, p.subtract(pRef).getNorm(), 1.0e-13);
+        int[] pivot = lu.getPivot();
+        for (int i = 0; i < pivotRef.length; ++i) {
+            assertEquals(pivotRef[i], pivot[i]);
+        }
+
+        // check the same cached instance is returned the second time
+        assertTrue(l == lu.getL());
+        assertTrue(u == lu.getU());
+        assertTrue(p == lu.getP());
+        
+    }
+
+    /** test determinant */
+    public void testDeterminant() {
+        assertEquals(-1,
+                     new LUDecompositionImpl(new RealMatrixImpl(testData, false)).getDeterminant(),
+                     1.0e-15);
+        assertEquals(-10,
+                     new LUDecompositionImpl(new RealMatrixImpl(luData, false)).getDeterminant(),
+                     1.0e-14);
+        assertEquals(0,
+                     new LUDecompositionImpl(new RealMatrixImpl(singular, false)).getDeterminant(),
+                     1.0e-17);
+        assertEquals(0,
+                     new LUDecompositionImpl(new RealMatrixImpl(bigSingular, false)).getDeterminant(),
+                     1.0e-17);
+    }
+
+}

Propchange: commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/linear/LUDecompositionImplTest.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/linear/LUDecompositionImplTest.java
------------------------------------------------------------------------------
    svn:keywords = Author Date Id Revision



Mime
View raw message