Added solve methods using double[][] to linear algebra decomposition solvers.

JIRA: MATH-564

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1097730 13f79535-47bb-0310-9956-ffa450edef68

src/main/java/org/apache/commons/math/linear/CholeskyDecompositionImpl.java

@@ -273,15 +273,34 @@ public class CholeskyDecompositionImpl implements CholeskyDecomposition {
             return new ArrayRealVector(solve(b.getDataRef()), false);
         }
 
-        /** {@inheritDoc} */
+        /** Solve the linear equation A × X = B for matrices A.
-        public RealMatrix solve(RealMatrix b) {
+         * <p>The A matrix is implicit, it is provided by the underlying
+         * decomposition algorithm.</p>
+         * @param b right-hand side of the equation A &times; X = B
+         * @param reUseB if true, the b array will be reused and returned,
+         * instead of being copied
+         * @return a matrix X that minimizes the two norm of A &times; X - B
+         * @throws org.apache.commons.math.exception.DimensionMismatchException
+         * if the matrices dimensions do not match.
+         * @throws SingularMatrixException
+         * if the decomposed matrix is singular.
+         */
+        private double[][] solve(double[][] b, boolean reUseB) {
             final int m = lTData.length;
-            if (b.getRowDimension() != m) {
+            if (b.length != m) {
-                throw new DimensionMismatchException(b.getRowDimension(), m);
+                throw new DimensionMismatchException(b.length, m);
             }
 
-            final int nColB = b.getColumnDimension();
+            final int nColB = b[0].length;
-            double[][] x = b.getData();
+            final double[][] x;
+            if (reUseB) {
+                x = b;
+            } else {
+                x = new double[b.length][nColB];
+                for (int i = 0; i < b.length; ++i) {
+                    System.arraycopy(b[i], 0, x[i], 0, nColB);
+                }
+            }
 
             // Solve LY = b
             for (int j = 0; j < m; j++) {
@@ -316,10 +335,20 @@ public class CholeskyDecompositionImpl implements CholeskyDecomposition {
                 }
             }
 
-            return new Array2DRowRealMatrix(x, false);
+            return x;
 
         }
 
+        /** {@inheritDoc} */
+        public double[][] solve(double[][] b) {
+            return solve(b, false);
+        }
+
+        /** {@inheritDoc} */
+        public RealMatrix solve(RealMatrix b) {
+            return new Array2DRowRealMatrix(solve(b.getData(), true), false);
+        }
+
         /** {@inheritDoc} */
         public RealMatrix getInverse() {
             return solve(MatrixUtils.createRealIdentityMatrix(lTData.length));

src/main/java/org/apache/commons/math/linear/DecompositionSolver.java

@@ -59,6 +59,18 @@ public interface DecompositionSolver {
      */
     RealVector solve(final RealVector b);
 
+    /** Solve the linear equation A × X = B for matrices A.
+     * <p>The A matrix is implicit, it is provided by the underlying
+     * decomposition algorithm.</p>
+     * @param b right-hand side of the equation A &times; X = B
+     * @return a matrix X that minimizes the two norm of A &times; X - B
+     * @throws org.apache.commons.math.exception.DimensionMismatchException
+     * if the matrices dimensions do not match.
+     * @throws SingularMatrixException
+     * if the decomposed matrix is singular.
+     */
+    double[][] solve(final double[][] b);
+
     /** Solve the linear equation A &times; X = B for matrices A.
      * <p>The A matrix is implicit, it is provided by the underlying
      * decomposition algorithm.</p>

src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java

@@ -338,37 +338,48 @@ public class EigenDecompositionImpl implements EigenDecomposition {
             return new ArrayRealVector(bp, false);
         }
 
-        /**
+        /** Solve the linear equation A × X = B for matrices A.
-         * Solve the linear equation A × X = B for symmetric matrices A.
+         * <p>The A matrix is implicit, it is provided by the underlying
-         * <p>
+         * decomposition algorithm.</p>
-         * This method only find exact linear solutions, i.e. solutions for
+         * @param b right-hand side of the equation A &times; X = B
-         * which ||A &times; X - B|| is exactly 0.
+         * @param reUseB if true, the b array will be reused and returned,
-         * </p>
+         * instead of being copied
-         * @param b Right-hand side of the equation A &times; X = B
+         * @return a matrix X that minimizes the two norm of A &times; X - B
-         * @return a Matrix X that minimizes the two norm of A &times; X - B
+         * @throws org.apache.commons.math.exception.DimensionMismatchException
-         * @throws DimensionMismatchException if the matrices dimensions do not match.
+         * if the matrices dimensions do not match.
-         * @throws SingularMatrixException if the decomposed matrix is singular.
+         * @throws SingularMatrixException
+         * if the decomposed matrix is singular.
          */
-        public RealMatrix solve(final RealMatrix b) {
+        private double[][] solve(double[][] b, boolean reUseB) {
 
             if (!isNonSingular()) {
                 throw new SingularMatrixException();
             }
 
             final int m = realEigenvalues.length;
-            if (b.getRowDimension() != m) {
+            if (b.length != m) {
-                throw new DimensionMismatchException(b.getRowDimension(), m);
+                throw new DimensionMismatchException(b.length, m);
             }
 
-            final int nColB = b.getColumnDimension();
+            final int nColB = b[0].length;
-            final double[][] bp = new double[m][nColB];
+            final double[][] bp;
+            if (reUseB) {
+                bp = b;
+            } else {
+                bp = new double[m][nColB];
+            }
+            final double[] tmpCol = new double[m];
             for (int k = 0; k < nColB; ++k) {
+                for (int i = 0; i < m; ++i) {
+                    tmpCol[i] = b[i][k];
+                    bp[i][k]  = 0;
+                }
                 for (int i = 0; i < m; ++i) {
                     final ArrayRealVector v = eigenvectors[i];
                     final double[] vData = v.getDataRef();
                     double s = 0;
                     for (int j = 0; j < m; ++j) {
-                        s += v.getEntry(j) * b.getEntry(j, k);
+                        s += v.getEntry(j) * tmpCol[j];
                     }
                     s /= realEigenvalues[i];
                     for (int j = 0; j < m; ++j) {
@@ -377,10 +388,20 @@ public class EigenDecompositionImpl implements EigenDecomposition {
                 }
             }
 
-            return MatrixUtils.createRealMatrix(bp);
+            return bp;
 
         }
 
+        /** {@inheritDoc} */
+        public double[][] solve(double[][] b) {
+            return solve(b, false);
+        }
+
+        /** {@inheritDoc} */
+        public RealMatrix solve(RealMatrix b) {
+            return new Array2DRowRealMatrix(solve(b.getData(), true), false);
+        }
+
         /**
          * Check if the decomposed matrix is non-singular.
          * @return true if the decomposed matrix is non-singular

src/main/java/org/apache/commons/math/linear/LUDecompositionImpl.java

@@ -337,17 +337,17 @@ public class LUDecompositionImpl implements LUDecomposition {
         }
 
         /** {@inheritDoc} */
-        public RealMatrix solve(RealMatrix b) {
+        public double[][] solve(double[][] b) {
 
             final int m = pivot.length;
-            if (b.getRowDimension() != m) {
+            if (b.length != m) {
-                throw new DimensionMismatchException(b.getRowDimension(), m);
+                throw new DimensionMismatchException(b.length, m);
             }
             if (singular) {
                 throw new SingularMatrixException();
             }
 
-            final int nColB = b.getColumnDimension();
+            final int nColB = b[0].length;
 
             // Apply permutations to b
             final double[][] bp = new double[m][nColB];
@@ -355,7 +355,7 @@ public class LUDecompositionImpl implements LUDecomposition {
                 final double[] bpRow = bp[row];
                 final int pRow = pivot[row];
                 for (int col = 0; col < nColB; col++) {
-                    bpRow[col] = b.getEntry(pRow, col);
+                    bpRow[col] = b[pRow][col];
                 }
             }
 
@@ -387,7 +387,13 @@ public class LUDecompositionImpl implements LUDecomposition {
                 }
             }
 
-            return new Array2DRowRealMatrix(bp, false);
+            return bp;
+
+        }
+
+        /** {@inheritDoc} */
+        public RealMatrix solve(RealMatrix b) {
+            return new Array2DRowRealMatrix(solve(b.getData()), false);
         }
 
         /** {@inheritDoc} */

src/main/java/org/apache/commons/math/linear/QRDecompositionImpl.java

@@ -337,6 +337,11 @@ public class QRDecompositionImpl implements QRDecomposition {
             return new ArrayRealVector(solve(b.getDataRef()), false);
         }
 
+        /** {@inheritDoc} */
+        public double[][] solve(double[][] b) {
+            return solve(new BlockRealMatrix(b)).getData();
+        }
+
         /** {@inheritDoc} */
         public RealMatrix solve(RealMatrix b) {
             final int n = qrt.length;

src/main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java

@@ -301,6 +301,22 @@ public class SingularValueDecompositionImpl implements SingularValueDecompositio
             return pseudoInverse.operate(b);
         }
 
+        /**
+         * Solve the linear equation A × X = B in least square sense.
+         * <p>
+         * The m&times;n matrix A may not be square, the solution X is such that
+         * ||A &times; X - B|| is minimal.
+         * </p>
+         *
+         * @param b Right-hand side of the equation A &times; X = B
+         * @return a matrix X that minimizes the two norm of A &times; X - B
+         * @throws org.apache.commons.math.exception.DimensionMismatchException
+         * if the matrices dimensions do not match.
+         */
+        public double[][] solve(final double[][] b) {
+            return pseudoInverse.multiply(MatrixUtils.createRealMatrix(b)).getData();
+        }
+
         /**
          * Solve the linear equation A &times; X = B in least square sense.
          * <p>

src/site/xdoc/changes.xml

@@ -52,6 +52,9 @@ The <action> type attribute can be add,update,fix,remove.
     If the output is not quite correct, check for invisible trailing spaces!
      -->
     <release version="3.0" date="TBD" description="TBD">
+      <action dev="luc" type="add" issue="MATH-564">
+        Added solve methods using double[][] to linear algebra decomposition solvers.
+      </action>
       <action dev="erans" type="add" issue="MATH-562">
         Added an interpolator adapter for data with known period.
       </action>

src/test/java/org/apache/commons/math/linear/CholeskySolverTest.java

@@ -81,6 +81,9 @@ public class CholeskySolverTest {
         // using RealMatrix
         Assert.assertEquals(0, solver.solve(b).subtract(xRef).getNorm(), 1.0e-13);
 
+        // using double[][]
+        Assert.assertEquals(0, MatrixUtils.createRealMatrix(solver.solve(b.getData())).subtract(xRef).getNorm(), 1.0e-13);
+
         // using double[]
         for (int i = 0; i < b.getColumnDimension(); ++i) {
             Assert.assertEquals(0,

src/test/java/org/apache/commons/math/linear/EigenSolverTest.java

@@ -120,6 +120,10 @@ public class EigenSolverTest {
         RealMatrix solution=es.solve(b);
         Assert.assertEquals(0, solution.subtract(xRef).getNorm(), 2.5e-12);
 
+        // using double[][]
+        solution = MatrixUtils.createRealMatrix(es.solve(b.getData()));
+        Assert.assertEquals(0, solution.subtract(xRef).getNorm(), 2.5e-12);
+
         // using double[]
         for (int i = 0; i < b.getColumnDimension(); ++i) {
             Assert.assertEquals(0,

src/test/java/org/apache/commons/math/linear/LUSolverTest.java

@@ -143,6 +143,9 @@ public class LUSolverTest {
         // using RealMatrix
         Assert.assertEquals(0, solver.solve(b).subtract(xRef).getNorm(), 1.0e-13);
 
+        // using double[][]
+        Assert.assertEquals(0, MatrixUtils.createRealMatrix(solver.solve(b.getData())).subtract(xRef).getNorm(), 1.0e-13);
+
         // using double[]
         for (int i = 0; i < b.getColumnDimension(); ++i) {
             Assert.assertEquals(0,

src/test/java/org/apache/commons/math/linear/QRSolverTest.java

@@ -137,6 +137,9 @@ public class QRSolverTest {
         // using RealMatrix
         Assert.assertEquals(0, solver.solve(b).subtract(xRef).getNorm(), 2.0e-16 * xRef.getNorm());
 
+        // using double[][]
+        Assert.assertEquals(0, MatrixUtils.createRealMatrix(solver.solve(b.getData())).subtract(xRef).getNorm(), 2.0e-16 * xRef.getNorm());
+
         // using double[]
         for (int i = 0; i < b.getColumnDimension(); ++i) {
             final double[] x = solver.solve(b.getColumn(i));

src/test/java/org/apache/commons/math/linear/SingularValueSolverTest.java

@@ -101,6 +101,9 @@ public class SingularValueSolverTest {
         // using RealMatrix
         Assert.assertEquals(0, solver.solve(b).subtract(xRef).getNorm(), normTolerance);
 
+        // using double[][]
+        Assert.assertEquals(0, MatrixUtils.createRealMatrix(solver.solve(b.getData())).subtract(xRef).getNorm(), normTolerance);
+
         // using double[]
         for (int i = 0; i < b.getColumnDimension(); ++i) {
             Assert.assertEquals(0,