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×A = L×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 [EMAIL
PROTECTED]
+ * #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 [EMAIL PROTECTED]
#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×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();
+
+ }
+
+ /** [EMAIL PROTECTED] */
+ 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;
+ }
+
+ /** [EMAIL PROTECTED] */
+ 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;
+ }
+
+ /** [EMAIL PROTECTED] */
+ 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;
+ }
+
+ /** [EMAIL PROTECTED] */
+ public int[] getPivot() {
+ return pivot.clone();
+ }
+
+ /** [EMAIL PROTECTED] */
+ public boolean isNonSingular() {
+ return !singular;
+ }
+
+ /** [EMAIL PROTECTED] */
+ public double getDeterminant() {
+ if (singular) {
+ return 0;
+ } else {
+ double determinant = parity;
+ for (int i = 0; i < m; i++) {
+ determinant *= lu[i][i];
+ }
+ return determinant;
+ }
+ }
+
+ /** [EMAIL PROTECTED] */
+ 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;
+
+ }
+
+ /** [EMAIL PROTECTED] */
+ 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 × 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
+ * @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);
+ }
+
+ /** [EMAIL PROTECTED] */
+ 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
