add selectable sparse matrix solver

src/com/lushprojects/circuitjs1/client/EditOptions.java

@@ -96,12 +96,18 @@ public EditInfo getEditInfo(int n) {
 		    ei.checkbox = new Checkbox("Auto-Adjust Timestep", sim.adjustTimeStep);
 		    return ei;
 		}
-		if (n == 16 && sim.adjustTimeStep)
+		if (n == 16) {
+		    EditInfo ei = new EditInfo("Matrix Solver", 0, -1, -1);
+		    ei.choice = new Choice();
+		    ei.choice.add(Locale.LS("Auto"));
+		    ei.choice.add(Locale.LS("Dense (Crout's LU)"));
+		    ei.choice.add(Locale.LS("Sparse (CSC LU)"));
+		    ei.choice.select(sim.solverType);
+		    return ei;
+		}
+		if (n == 17 && sim.adjustTimeStep)
 		    return new EditInfo("Minimum time step size (s)", sim.minTimeStep, 0, 0).setPositive();
 
-		// don't add new options here.  they are only visible if sim.adjustTimeStemp is set, and it isn't by default.
-		// add them before the "Auto-Adjust Timestep" checkbox.
-
 		return null;
 	}
 
@@ -197,7 +203,19 @@ public void setEditValue(int n, EditInfo ei) {
 		    sim.adjustTimeStep = ei.checkbox.getState();
 		    ei.newDialog = true;
 		}
-		if (n == 16)
+		if (n == 16) {
+		    int newType = ei.choice.getSelectedIndex();
+		    if (newType != sim.solverType) {
+			sim.solverType = newType;
+			// Save as default for new circuits; existing circuit will persist its own
+			// solverType via the "st" XML attr when saved.
+			Storage stor = Storage.getLocalStorageIfSupported();
+			if (stor != null)
+			    stor.setItem("solverType", Integer.toString(sim.solverType));
+			app.needAnalyze();
+		    }
+		}
+		if (n == 17 && ei.value > 0)
 		    sim.minTimeStep = ei.value;
 	}
 

src/com/lushprojects/circuitjs1/client/SimulationManager.java

@@ -7,11 +7,14 @@
 import java.util.Map;
 import java.util.Vector;
 
+import com.google.gwt.storage.client.Storage;
 import com.google.gwt.user.client.Window;
 import com.google.gwt.xml.client.Document;
 import com.google.gwt.xml.client.Element;
 import com.google.gwt.xml.client.Node;
 import com.google.gwt.xml.client.XMLParser;
+import com.lushprojects.circuitjs1.client.matrix.DMatrixSparseCSC;
+import com.lushprojects.circuitjs1.client.matrix.SparseLU;
 import com.lushprojects.circuitjs1.client.util.Locale;
 
 public class SimulationManager {
@@ -35,6 +38,15 @@ public class SimulationManager {
     CircuitElm elmArr[];
     double t;
 
+    // Solver type: 0=Auto, 1=Dense, 2=Sparse
+    static final int SOLVER_AUTO = 0;
+    static final int SOLVER_DENSE = 1;
+    static final int SOLVER_SPARSE = 2;
+    int solverType = SOLVER_AUTO;
+    static final int SPARSE_THRESHOLD = 150;
+    boolean usingSparse;
+    SparseLU sparseLU;
+
     // current timestep (time between iterations)
     double timeStep;
 
@@ -56,7 +68,17 @@ public class SimulationManager {
     public static native void console(String text) /*-{ console.log(text) }-*/;
     public static native void debugger() /*-{ debugger; }-*/;
 
-    SimulationManager(CirSim app_) { theSim = this; app = app_; }
+    SimulationManager(CirSim app_) {
+	theSim = this; app = app_;
+	// load solver preference from localStorage
+	try {
+	    Storage stor = Storage.getLocalStorageIfSupported();
+	    if (stor != null) {
+		String s = stor.getItem("solverType");
+		if (s != null) solverType = Integer.parseInt(s);
+	    }
+	} catch (Exception e) {}
+    }
 
     void resetTime() {
     	t = timeStepAccum = 0;
@@ -1003,16 +1025,26 @@ void stampCircuit() {
 	    return;
 
 	// save original matrices for restoring during nonlinear iterations
+	int maxMatrixSize = 0;
 	for (int mi = 0; mi != matrices.length; mi++) {
 	    CircuitMatrix m = matrices[mi];
 	    int sz = m.size;
+	    if (sz > maxMatrixSize) maxMatrixSize = sz;
 	    for (i = 0; i != sz; i++)
 		m.origRightSide[i] = m.rightSide[i];
 	    for (i = 0; i != sz; i++)
 		for (int j = 0; j != sz; j++)
 		    m.origMatrix[i][j] = m.matrix[i][j];
 	    //CirSim.console("matrix " + mi + " size: " + sz);
 	}
+	// determine which solver to use (AUTO uses the largest matrix size to decide)
+	if (solverType == SOLVER_SPARSE)
+	    usingSparse = true;
+	else if (solverType == SOLVER_DENSE)
+	    usingSparse = false;
+	else // SOLVER_AUTO
+	    usingSparse = (maxMatrixSize >= SPARSE_THRESHOLD);
+	sparseLU = null; // reset on re-stamp
 
 	// if a matrix is linear, we can do the lu_factor here instead of
 	// needing to do it every frame
@@ -1755,31 +1787,42 @@ public CustomCompositeModel getCircuitAsComposite() {
 
 	static void invertMatrix(double a[][], int n) {
 	    int ipvt[] = new int[n];
-	    lu_factor(a, n, ipvt);
+	    lu_factor_dense(a, n, ipvt);
 	    int i, j;
 	    double b[] = new double[n];
 	    double inva[][] = new double[n][n];
-	    
+
 	    // solve for each column of identity matrix
 	    for (i = 0; i != n; i++) {
 		for (j = 0; j != n; j++)
 		    b[j] = 0;
 		b[i] = 1;
-		lu_solve(a, n, ipvt, b);
+		lu_solve_dense(a, n, ipvt, b);
 		for (j = 0; j != n; j++)
 		    inva[j][i] = b[j];
 	    }
-	    
+
 	    // return in original matrix
 	    for (i = 0; i != n; i++)
 		for (j = 0; j != n; j++)
 		    a[i][j] = inva[i][j];
 	}
+    // Dispatching lu_factor: uses sparse or dense solver based on solverType setting
+    static boolean lu_factor(double a[][], int n, int ipvt[]) {
+	SimulationManager sm = theSim;
+	if (sm != null && sm.usingSparse) {
+	    DMatrixSparseCSC csc = DMatrixSparseCSC.convert(a, DMatrixSparseCSC.EPS);
+	    if (sm.sparseLU == null) sm.sparseLU = new SparseLU();
+	    return sm.sparseLU.setA(csc);
+	}
+	return lu_factor_dense(a, n, ipvt);
+    }
+
     // factors a matrix into upper and lower triangular matrices by
     // gaussian elimination.  On entry, a[0..n-1][0..n-1] is the
     // matrix to be factored.  ipvt[] returns an integer vector of pivot
     // indices, used in the lu_solve() routine.
-    static boolean lu_factor(double a[][], int n, int ipvt[]) {
+    static boolean lu_factor_dense(double a[][], int n, int ipvt[]) {
 	int i,j,k;
 
 	// check for a possible singular matrix by scanning for rows that
@@ -1859,10 +1902,20 @@ static boolean lu_factor(double a[][], int n, int ipvt[]) {
 	return true;
     }
 
+    // Dispatching lu_solve: uses sparse or dense solver based on solverType setting
+    static void lu_solve(double a[][], int n, int ipvt[], double b[]) {
+	SimulationManager sm = theSim;
+	if (sm != null && sm.usingSparse) {
+	    sm.sparseLU.solve(b, b);
+	    return;
+	}
+	lu_solve_dense(a, n, ipvt, b);
+    }
+
     // Solves the set of n linear equations using a LU factorization
     // previously performed by lu_factor.  On input, b[0..n-1] is the right
     // hand side of the equations, and on output, contains the solution.
-    static void lu_solve(double a[][], int n, int ipvt[], double b[]) {
+    static void lu_solve_dense(double a[][], int n, int ipvt[], double b[]) {
 	int i;
 
 	// find first nonzero b element

src/com/lushprojects/circuitjs1/client/XMLDeserializer.java

@@ -69,6 +69,7 @@ void readCircuit(String text, int readFlags) {
 	    CircuitElm.voltageRange = parseDoubleAttr("vr", CircuitElm.voltageRange);
 	    ui.powerBar.setValue(parseIntAttr("pb", ui.powerBar.getValue()));
 	    sim.minTimeStep = parseDoubleAttr("mts", sim.minTimeStep);
+	    sim.solverType = parseIntAttr("st", sim.solverType);
 	    app.setGrid();
 	}
 

src/com/lushprojects/circuitjs1/client/XMLSerializer.java

@@ -126,6 +126,8 @@ String dumpCircuit() {
 	XMLSerializer.dumpAttr(root, "pb", ui.powerBar.getValue());
 	XMLSerializer.dumpAttr(root, "vr", CircuitElm.voltageRange);
 	XMLSerializer.dumpAttr(root, "mts", sim.minTimeStep);
+	if (sim.solverType != 0)
+	    XMLSerializer.dumpAttr(root, "st", sim.solverType);
 
 	for (CircuitElm ce: app.elmList) {
 	    Element elem = doc.createElement(ce.getXmlDumpType());

src/com/lushprojects/circuitjs1/client/matrix/DGrowArray.java

@@ -0,0 +1,99 @@
+/*
+ * Derived from EJML (Efficient Java Matrix Library)
+ * Copyright (c) 2009-2020, Peter Abeles. All Rights Reserved.
+ *
+ * Licensed 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.
+ *
+ * Adapted for CircuitJS1 by H-Dynamite (sharpie7/circuitjs1 PR #920).
+ */
+
+package com.lushprojects.circuitjs1.client.matrix;
+
+public class DGrowArray {
+    public double[] data;
+    public int length;
+
+    public DGrowArray( int length ) {
+        this.data = new double[length];
+        this.length = length;
+    }
+
+    public DGrowArray() {
+        this(0);
+    }
+
+    public int length() {
+        return length;
+    }
+
+    public void reset() {reshape(0);}
+
+    /**
+     * Changes the array's length and doesn't attempt to preserve previous values if a new array is required
+     *
+     * @param length New array length
+     */
+    public DGrowArray reshape( int length ) {
+        if (data.length < length) {
+            data = new double[length];
+        }
+        this.length = length;
+        return this;
+    }
+
+    /**
+     * Increases the internal array's length by the specified amount. Previous values are preserved.
+     * The length value is not modified since this does not change the 'meaning' of the array, just
+     * increases the amount of data which can be stored in it.
+     *
+     * this.data = new data_type[ data.length + amount ]
+     *
+     * @param amount Number of elements added to the internal array's length
+     */
+    public void growInternal( int amount ) {
+        double[] tmp = new double[data.length + amount];
+
+        System.arraycopy(data, 0, tmp, 0, data.length);
+        this.data = tmp;
+    }
+
+    public void setTo( DGrowArray original ) {
+        reshape(original.length);
+        System.arraycopy(original.data, 0, data, 0, original.length);
+    }
+
+    public void add( double value ) {
+        if (length >= data.length) {
+            growInternal(Math.min(500_000, data.length + 10));
+        }
+
+        data[length++] = value;
+    }
+
+    public double get( int index ) {
+        if (index < 0 || index >= length)
+            throw new IllegalArgumentException("Out of bounds");
+        return data[index];
+    }
+
+    public void set( int index, double value ) {
+        if (index < 0 || index >= length)
+            throw new IllegalArgumentException("Out of bounds");
+        data[index] = value;
+    }
+
+    public void free() {
+        data = new double[0];
+        length = 0;
+    }
+}