Optimize Solver

src/arcade/core/util/MatrixArray.java

@@ -0,0 +1,166 @@
+package arcade.core.util;
+
+import java.util.ArrayList;
+import java.util.Collections;
+import arcade.core.util.Matrix.Value;
+
+/**
+ * Container class for array-based sparse matrix representation.
+ *
+ * <p>Class provides a subset of matrix operations needed for solving a system of linear equations
+ * using the successive over-relaxation method in {@link arcade.core.util.Solver}.
+ */
+public class MatrixArray {
+    /** Number of rows in the matrix. */
+    int nRows;
+
+    /** Number of columns in the matrix. */
+    int nColumns;
+
+    /** List of non-zero values in the matrix. */
+    Value[] values;
+
+    /**
+     * Constructs a sparse matrix from a dense two-dimensional array.
+     *
+     * <p>All non-zero elements of the input matrix are stored in the internal sparse
+     * representation. The input matrix must be rectangular; that is, every row must have the same
+     * number of columns.
+     *
+     * @param a the dense matrix to convert
+     * @throws IllegalArgumentException if the matrix is empty, contains empty rows, or is not
+     *     rectangular
+     */
+    public MatrixArray(double[][] a) {
+        nRows = a.length;
+        if (nRows == 0) {
+            throw new IllegalArgumentException("MatrixArray ctor: input is empty");
+        }
+        nColumns = a[0].length;
+
+        if (nColumns == 0) {
+            throw new IllegalArgumentException("Matrix array ctor: empty columns");
+        }
+
+        ArrayList<Value> alValues = new ArrayList<>();
+        int nNonZero = 0;
+        for (int row = 0; row < nRows; row++) {
+            if (a[row].length != nColumns) {
+                throw new IllegalArgumentException(
+                        "Matrix array ctor: not all columns are the same length");
+            }
+
+            for (int column = 0; column < nColumns; column++) {
+                if (a[row][column] != 0.) {
+                    alValues.add(new Value(row, column, a[row][column]));
+                    nNonZero++;
+                }
+            }
+        }
+
+        //
+        // No need to sort, since we built them in order.
+        //
+
+        values = new Value[nNonZero];
+        int i = 0;
+        for (Value v : alValues) {
+            values[i] = v;
+            i++;
+        }
+    } // ctor
+
+    /**
+     * Constructs a sparse matrix from a collection of matrix entries.
+     *
+     * <p>The entries are sorted by row and column index and copied into an internal array to
+     * support efficient sparse matrix operations such as matrix-vector multiplication.
+     *
+     * @param alValues the non-zero matrix entries
+     * @param rows the number of rows in the matrix
+     * @param columns the number of columns in the matrix
+     */
+    public MatrixArray(ArrayList<Value> alValues, int rows, int columns) {
+        nRows = rows;
+        nColumns = columns;
+        values = new Value[alValues.size()];
+
+        Collections.sort(
+                alValues,
+                (v1, v2) -> (v1.i == v2.i ? Integer.compare(v1.j, v2.j) : (v1.i > v2.i ? 1 : -1)));
+
+        int i = 0;
+        for (Value v : alValues) {
+            values[i] = v;
+            i++;
+        }
+    }
+
+    /**
+     * Solves the equation {@code Lx = b} using forward substitution for an array-based sparse
+     * matrix.
+     *
+     * @param b the right-hand side vector
+     * @return the left-hand side vector
+     */
+    public double[] forwardSubstitution(double[] b) {
+
+        int n = b.length;
+        double[] subbed = new double[n];
+        double[] diag = new double[n];
+
+        // Group lower diagonal by row.
+        ArrayList<ArrayList<Value>> rowsL = new ArrayList<ArrayList<Value>>();
+        for (int r = 0; r < n; r++) {
+            rowsL.add(new ArrayList<>());
+        }
+        for (Value v : values) {
+            rowsL.get(v.i).add(v);
+        }
+
+        // Get values along diagonal.
+        for (Value v : values) {
+            if (v.i == v.j) {
+                diag[v.i] = v.v;
+            }
+        }
+
+        // Iterate only through non-zero entries in the lower diagonal matrix.
+        for (int i = 0; i < n; i++) {
+            ArrayList<Value> rowL = rowsL.get(i);
+            double val = 0;
+            for (Value v : rowL) {
+                val += subbed[v.j] * v.v;
+            }
+            val = b[i] - val;
+            subbed[i] = val / diag[i];
+        }
+
+        return subbed;
+    }
+
+    /**
+     * Multiplies this matrix by a vector.
+     *
+     * <p>The matrix is stored in sparse form, so only non-zero entries are processed during the
+     * multiplication.
+     *
+     * @param b the vector to multiply by
+     * @return the resulting vector {@code A * b}
+     * @throws IllegalArgumentException if the vector length is not compatible with the matrix
+     *     dimensions
+     */
+    public double[] multiply(double[] b) {
+        if (b.length != nRows) {
+            throw new IllegalArgumentException("MatrixArray.multiply (by a vector): conformation");
+        }
+        double[] multiplied = new double[nRows];
+
+        // Iterate through all entries and multiply.
+        for (Value a : values) {
+            multiplied[a.i] += a.v * b[a.j];
+        }
+
+        return multiplied;
+    }
+}

src/arcade/core/util/Solver.java

@@ -1,7 +1,6 @@
 package arcade.core.util;
 
 import java.util.ArrayList;
-import java.util.logging.Logger;
 import arcade.core.util.Matrix.Value;
 import static arcade.core.util.Matrix.*;
 
@@ -19,8 +18,6 @@
  * </ul>
  */
 public class Solver {
-    /** Logger for {@code Solver}. */
-    private static final Logger LOGGER = Logger.getLogger(Solver.class.getName());
 
     /** Error tolerance for Cash-Karp. */
     private static final double ERROR = 1E-5;
@@ -38,13 +35,13 @@ public class Solver {
     private static final double OMEGA = 1.4;
 
     /** Maximum number of iterations. */
-    private static final int MAX_ITERS = 10000;
+    private static final int MAX_ITERS = 20000;
 
     /** Error tolerance for SOR. */
-    private static final double TOLERANCE = 1E-8;
+    private static final double TOLERANCE = 1E-5;
 
     /** Convergence delta for bisection method. */
-    private static final double DELTA = 1E-5;
+    private static final double DELTA = 1E-6;
 
     /** Matrix size threshold for dense representation. */
     private static final int MATRIX_THRESHOLD = 100;
@@ -370,7 +367,6 @@ private static double[] denseSOR(
             double[] r = subtract(vec, multiply(mat, xCurr));
             error = normalize(r);
         }
-
         return xCurr;
     }
 
@@ -396,6 +392,7 @@ private static double[] sparseSOR(
         double[] c = forwardSubstitution(sparseA, vec);
         ArrayList<Value> t = forwardSubstitution(sparseA);
         t = scale(t, -1);
+        MatrixArray tArray = new MatrixArray(t, mat.length, mat[0].length);
 
         // Set initial guess.
         double[] xCurr = x0;
@@ -404,7 +401,7 @@ private static double[] sparseSOR(
         // Iterate until convergence.
         while (i < maxIters && error > tolerance) {
             // Calculate new guess for x.
-            xCurr = add(scale(add(multiply(t, xPrev), c), OMEGA), scale(xPrev, 1 - OMEGA));
+            xCurr = add(scale(add(tArray.multiply(xPrev), c), OMEGA), scale(xPrev, 1 - OMEGA));
 
             // Set previous to copy of current and increment iteration count.
             xPrev = xCurr;
@@ -414,7 +411,6 @@ private static double[] sparseSOR(
             double[] r = subtract(vec, multiply(sparseA, xCurr));
             error = normalize(r);
         }
-
         return xCurr;
     }
 
@@ -428,9 +424,11 @@ private static double[] sparseSOR(
      * @param a the lower bound on the interval
      * @param b the upper bound on the interval
      * @param maxIters the maximum number of iterations
+     * @param tolerance the error
      * @return the root of the function
      */
-    public static double bisection(Function func, double a, double b, int maxIters) {
+    public static double bisection(
+            Function func, double a, double b, int maxIters, double tolerance) {
         double c;
         double fc;
         int i = 0;
@@ -452,7 +450,7 @@ public static double bisection(Function func, double a, double b, int maxIters)
             fc = func.f(c);
 
             // Check for exit conditions.
-            if (fc == 0 || (b - a) / 2 < DELTA) {
+            if (fc == 0 || (b - a) / 2 < tolerance) {
                 return c;
             } else {
                 if (Math.signum(fc) == Math.signum(func.f(a))) {
@@ -469,7 +467,7 @@ public static double bisection(Function func, double a, double b, int maxIters)
     }
 
     /**
-     * Finds root using bisection method with default maximum iterations.
+     * Finds root using bisection method with default maximum iterations and tolerance.
      *
      * <p>Root is found by repeatedly bisecting the interval and selecting the interval in which the
      * function changes sign. If no root is found, the simulation will throw an ArithmeticException.
@@ -480,6 +478,38 @@ public static double bisection(Function func, double a, double b, int maxIters)
      * @return the root of the function
      */
     public static double bisection(Function func, double a, double b) {
-        return bisection(func, a, b, MAX_ITERS);
+        return bisection(func, a, b, MAX_ITERS, DELTA);
+    }
+
+    /**
+     * Finds root using bisection method with default maximum iterations.
+     *
+     * <p>Root is found by repeatedly bisecting the interval and selecting the interval in which the
+     * function changes sign. If no root is found, the simulation will throw an ArithmeticException.
+     *
+     * @param func the function
+     * @param a the lower bound on the interval
+     * @param b the upper bound on the interval
+     * @param delta the error tolerance
+     * @return the root of the function
+     */
+    public static double bisection(Function func, double a, double b, double delta) {
+        return bisection(func, a, b, MAX_ITERS, delta);
+    }
+
+    /**
+     * Finds root using bisection method with default tolerance.
+     *
+     * <p>Root is found by repeatedly bisecting the interval and selecting the interval in which the
+     * function changes sign. If no root is found, the simulation will throw an ArithmeticException.
+     *
+     * @param func the function
+     * @param a the lower bound on the interval
+     * @param b the upper bound on the interval
+     * @param maxIters the maximum number of iterations
+     * @return the root of the function
+     */
+    public static double bisection(Function func, double a, double b, int maxIters) {
+        return bisection(func, a, b, maxIters, DELTA);
     }
 }

test/arcade/core/util/SolverTest.java

@@ -236,4 +236,21 @@ public void testBisection_quadraticFunctionAndSwappedInputs_returnsAnswer() {
 
         assertEquals(1.41421, result, 0.0001);
     }
+
+    @Test
+    public void testBisection_exceedsMaxIterations_returnsNaN() {
+        Function f = (x) -> x * x - 2;
+        double result = Solver.bisection(f, 2, 0, 2);
+
+        assertEquals(Double.NaN, result, 0.001);
+    }
+
+    @Test
+    public void testBisection_givenPrecision_returnsAnswersWithinPrecision() {
+        Function f = (x) -> x * x - 2;
+        double result = Solver.bisection(f, 2, 0, 1E-3);
+
+        assertEquals(1.41421, result, 1E-3);
+        assertNotEquals(1.41421, result, 1E-4);
+    }
 }