Author: psteitz Date: Fri Aug 12 05:39:55 2011 New Revision: 1156968 URL: http://svn.apache.org/viewvc?rev=1156968&view=rev Log: Added methods to solve upper and lower triangular systems. JIRA: MATH-624. Contributed by Greg Sterijevski.
Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/MatrixUtils.java commons/proper/math/trunk/src/site/xdoc/changes.xml commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/MatrixUtilsTest.java Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/MatrixUtils.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/MatrixUtils.java?rev=1156968&r1=1156967&r2=1156968&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/MatrixUtils.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/MatrixUtils.java Fri Aug 12 05:39:55 2011 @@ -25,6 +25,8 @@ import java.util.Arrays; import org.apache.commons.math.Field; import org.apache.commons.math.FieldElement; +import org.apache.commons.math.exception.MathArithmeticException; +import org.apache.commons.math.exception.MathIllegalArgumentException; import org.apache.commons.math.exception.OutOfRangeException; import org.apache.commons.math.exception.NoDataException; import org.apache.commons.math.exception.NumberIsTooSmallException; @@ -34,6 +36,8 @@ import org.apache.commons.math.exception import org.apache.commons.math.exception.util.LocalizedFormats; import org.apache.commons.math.fraction.BigFraction; import org.apache.commons.math.fraction.Fraction; +import org.apache.commons.math.util.FastMath; +import org.apache.commons.math.util.MathUtils; /** * A collection of static methods that operate on or return matrices. @@ -803,4 +807,84 @@ public class MatrixUtils { throw ioe; } } + + /**Solve a system of composed of a Lower Triangular Matrix + * {@link RealMatrix}. + * <p> + * This method is called to solve systems of equations which are + * of the lower triangular form. The matrix {@link RealMatrix} + * is assumed, though not checked, to be in lower triangular form. + * The vector {@link RealVector} is overwritten with the solution. + * The matrix is checked that it is square and its dimensions match + * the length of the vector. + * </p> + * @param rm RealMatrix which is lower triangular + * @param b RealVector this is overwritten + * @exception IllegalArgumentException if the matrix and vector are not conformable + * @exception ArithmeticException there is a zero or near zero on the diagonal of rm + */ + public static void solveLowerTriangularSystem( RealMatrix rm, RealVector b){ + if ((rm == null) || (b == null) || ( rm.getRowDimension() != b.getDimension())) { + throw new MathIllegalArgumentException(LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, + (rm == null) ? 0 : rm.getRowDimension(), + (b == null) ? 0 : b.getDimension()); + } + if( rm.getColumnDimension() != rm.getRowDimension() ){ + throw new MathIllegalArgumentException(LocalizedFormats.DIMENSIONS_MISMATCH_2x2, + rm.getRowDimension(),rm.getRowDimension(), + rm.getRowDimension(),rm.getColumnDimension()); + } + int rows = rm.getRowDimension(); + for( int i = 0 ; i < rows ; i++ ){ + double diag = rm.getEntry(i, i); + if( FastMath.abs(diag) < MathUtils.SAFE_MIN ){ + throw new MathArithmeticException(LocalizedFormats.ZERO_DENOMINATOR); + } + double bi = b.getEntry(i)/diag; + b.setEntry(i, bi ); + for( int j = i+1; j< rows; j++ ){ + b.setEntry(j, b.getEntry(j)-bi*rm.getEntry(j,i) ); + } + } + } + + /** Solver a system composed of an Upper Triangular Matrix + * {@link RealMatrix}. + * <p> + * This method is called to solve systems of equations which are + * of the lower triangular form. The matrix {@link RealMatrix} + * is assumed, though not checked, to be in upper triangular form. + * The vector {@link RealVector} is overwritten with the solution. + * The matrix is checked that it is square and its dimensions match + * the length of the vector. + * </p> + * @param rm RealMatrix which is upper triangular + * @param b RealVector this is overwritten + * @exception IllegalArgumentException if the matrix and vector are not conformable + * @exception ArithmeticException there is a zero or near zero on the diagonal of rm + */ + public static void solveUpperTriangularSystem( RealMatrix rm, RealVector b){ + if ((rm == null) || (b == null) || ( rm.getRowDimension() != b.getDimension())) { + throw new MathIllegalArgumentException(LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, + (rm == null) ? 0 : rm.getRowDimension(), + (b == null) ? 0 : b.getDimension()); + } + if( rm.getColumnDimension() != rm.getRowDimension() ){ + throw new MathIllegalArgumentException(LocalizedFormats.DIMENSIONS_MISMATCH_2x2, + rm.getRowDimension(),rm.getRowDimension(), + rm.getRowDimension(),rm.getColumnDimension()); + } + int rows = rm.getRowDimension(); + for( int i = rows-1 ; i >-1 ; i-- ){ + double diag = rm.getEntry(i, i); + if( FastMath.abs(diag) < MathUtils.SAFE_MIN ){ + throw new MathArithmeticException(LocalizedFormats.ZERO_DENOMINATOR); + } + double bi = b.getEntry(i)/diag; + b.setEntry(i, bi ); + for( int j = i-1; j>-1; j-- ){ + b.setEntry(j, b.getEntry(j)-bi*rm.getEntry(j,i) ); + } + } + } } Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=1156968&r1=1156967&r2=1156968&view=diff ============================================================================== --- commons/proper/math/trunk/src/site/xdoc/changes.xml (original) +++ commons/proper/math/trunk/src/site/xdoc/changes.xml Fri Aug 12 05:39:55 2011 @@ -52,6 +52,9 @@ The <action> type attribute can be add,u If the output is not quite correct, check for invisible trailing spaces! --> <release version="3.0" date="TBD" description="TBD"> + <action dev="psteitz" type="update" issue="MATH-624" due-to="Greg Sterijevski"> + Added methods to solve upper and lower triangular systems to MatrixUtils. + </action> <action dev="psteitz" type="update" issue="MATH-642"> Improved performance of nextInt(int) in BitsStreamGenerator. </action> Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/MatrixUtilsTest.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/MatrixUtilsTest.java?rev=1156968&r1=1156967&r2=1156968&view=diff ============================================================================== --- commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/MatrixUtilsTest.java (original) +++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/MatrixUtilsTest.java Fri Aug 12 05:39:55 2011 @@ -17,8 +17,7 @@ package org.apache.commons.math.linear; import java.math.BigDecimal; - - +import org.apache.commons.math.TestUtils; import org.apache.commons.math.fraction.BigFraction; import org.apache.commons.math.fraction.Fraction; import org.apache.commons.math.fraction.FractionConversionException; @@ -302,5 +301,29 @@ public final class MatrixUtilsTest { } return d; } + + @Test + public void testSolveLowerTriangularSystem(){ + RealMatrix rm = new Array2DRowRealMatrix( + new double[][] { {2,0,0,0 }, { 1,1,0,0 }, { 3,3,3,0 }, { 3,3,3,4 } }, + false); + RealVector b = new ArrayRealVector(new double[] { 2,3,4,8 }, false); + MatrixUtils.solveLowerTriangularSystem(rm, b); + TestUtils.assertEquals( new double[]{1,2,-1.66666666666667, 1.0} , b.getData() , 1.0e-12); + } + + + /* + * Taken from R manual http://stat.ethz.ch/R-manual/R-patched/library/base/html/backsolve.html + */ + @Test + public void testSolveUpperTriangularSystem(){ + RealMatrix rm = new Array2DRowRealMatrix( + new double[][] { {1,2,3 }, { 0,1,1 }, { 0,0,2 } }, + false); + RealVector b = new ArrayRealVector(new double[] { 8,4,2 }, false); + MatrixUtils.solveUpperTriangularSystem(rm, b); + TestUtils.assertEquals( new double[]{-1,3,1} , b.getData() , 1.0e-12); + } }