feat: add solveOnly() for chord Newton — reuse cached LU factorization

New method that skips equilibration, scatter, and factorLU. Just does:
permute RHS → solveLU (reuse cached factors) → iterative refinement.

For chord Newton in NR loops: the caller recomputes only the residual f
with updated voltages, then calls solveOnly to solve J*delta = -f using
the Jacobian factorization from the first NR iteration. Saves ~35ms per
chord iteration on c6288 (skip 40ms model eval + 8ms factorLU).

The csr_data parameter is still needed for the f64 SpMV in iterative
refinement (residual computed with original matrix, not the factored one).

Tested: 3x3 matrix, solveOnly matches fresh solve within 2.2e-16.
C API: sprux_ffi_solve_only() added.

Co-developed-by: Claude Code v2.1.81 (claude-opus-4-6)

Author: robtaylor

sprux/sprux/SpruxFFISolver.cpp

@@ -504,6 +504,132 @@ int SpruxFFISolver::solve(const double* csr_data, const double* rhs, double* x_o
   return itersUsedCpu;
 }
 
+// ---------------------------------------------------------------------------
+// solveOnly(): reuse cached factorization for chord Newton.
+// No equilibrate, no scatter, no factorLU — just permute RHS + solveLU + refine.
+// Uses equilibration scales and factored data from the most recent solve().
+// ---------------------------------------------------------------------------
+
+int SpruxFFISolver::solveOnly(const double* csr_data, const double* rhs, double* x_out) {
+  auto& d = *impl_;
+  const int64_t n = d.n;
+  const auto& perm = *d.perm;
+
+  // Use equilibration scales from the most recent solve() (stored in d.rowScale/colScale
+  // for CPU path, or in the last-used slot for Metal path).
+
+#ifdef SPRUX_USE_METAL
+  if (d.useMetal) {
+    auto& symCtx = d.solver->internalSymbolicContext();
+    auto& metalCtx = MetalContext::instance();
+
+    // The last completed slot has the factored data (solve() calls begin+end, swaps curSlot).
+    int factSlot = 1 - d.curSlot;
+    auto& slot = d.slots[factSlot];
+
+    // Permute RHS using the cached equilibration scales from the factored solve
+    for (int64_t j = 0; j < n; j++) {
+      slot.xGpu.ptr()[perm[j]] = float(slot.rowScale[j] * rhs[d.preproc.rowPerm[j]]);
+    }
+
+    // GPU solve only (no factor) — reuse factored dataGpu + devPivots
+    void* cmdBuf = metalCtx.createCommandBuffer();
+    void* encoder = metalCtx.createComputeEncoder(cmdBuf);
+    symCtx.setExternalEncoder(cmdBuf, encoder);
+
+    d.solver->solveLU(slot.dataGpu.ptr(), slot.devPivots.ptr(), slot.xGpu.ptr(), n, 1,
+                      *d.solveCtx, PivotLocation::Device);
+
+    // Iterative refinement with the NEW csr_data (for accurate f64 SpMV)
+    std::fill(slot.xAccum.begin(), slot.xAccum.end(), 0.0);
+    double bNormSq = 0.0;
+    for (int64_t j = 0; j < n; j++) bNormSq += rhs[j] * rhs[j];
+
+    int itersUsed = 0;
+    for (int iter = 0; iter < d.maxRefine; iter++) {
+      symCtx.clearExternalEncoder();
+
+      for (int64_t j = 0; j < n; j++) {
+        slot.xAccum[j] += slot.colScale[j] * double(slot.xGpu.ptr()[perm[j]]);
+      }
+
+      double resNormSq = 0.0;
+      for (int64_t j = 0; j < n; j++) {
+        int64_t srcRow = d.preproc.rowPerm[j];
+        double sum = 0.0;
+        for (int64_t k = d.csrIndptr[srcRow]; k < d.csrIndptr[srcRow + 1]; k++) {
+          sum += csr_data[k] * slot.xAccum[d.csrIndices[k]];
+        }
+        double residual = rhs[srcRow] - sum;
+        resNormSq += residual * residual;
+        slot.xGpu.ptr()[perm[j]] = float(slot.rowScale[j] * residual);
+      }
+      itersUsed = iter + 1;
+
+      if (d.refineTol > 0 && resNormSq <= d.refineTol * d.refineTol * std::max(bNormSq, 1e-300)) {
+        break;
+      }
+
+      void* newCmdBuf = metalCtx.createCommandBuffer();
+      void* newEncoder = metalCtx.createComputeEncoder(newCmdBuf);
+      symCtx.setExternalEncoder(newCmdBuf, newEncoder);
+
+      d.solver->solveLU(slot.dataGpu.ptr(), slot.devPivots.ptr(), slot.xGpu.ptr(), n, 1,
+                        *d.solveCtx, PivotLocation::Device);
+    }
+
+    symCtx.clearExternalEncoder();
+
+    for (int64_t j = 0; j < n; j++) {
+      x_out[j] = slot.xAccum[j] + slot.colScale[j] * double(slot.xGpu.ptr()[perm[j]]);
+    }
+    return itersUsed;
+  }
+#endif
+
+  // CPU fallback — reuse cached dataCpu + pivotsCpu (from last solve())
+  auto& bp = d.bpCpu;
+  for (int64_t j = 0; j < n; j++) {
+    bp[perm[j]] = float(d.rowScale[j] * rhs[d.preproc.rowPerm[j]]);
+  }
+  d.solver->solveLU(d.dataCpu.data(), d.pivotsCpu.data(), bp.data(), n, 1);
+
+  for (int64_t j = 0; j < n; j++) {
+    x_out[j] = d.colScale[j] * double(bp[perm[j]]);
+  }
+
+  double bNormSq = 0.0;
+  for (int64_t j = 0; j < n; j++) bNormSq += rhs[j] * rhs[j];
+  int itersUsed = 0;
+
+  for (int iter = 0; iter < d.maxRefine; iter++) {
+    std::fill(bp.begin(), bp.end(), 0.0f);
+    double resNormSq = 0.0;
+    for (int64_t j = 0; j < n; j++) {
+      int64_t srcRow = d.preproc.rowPerm[j];
+      double sum = 0.0;
+      for (int64_t k = d.csrIndptr[srcRow]; k < d.csrIndptr[srcRow + 1]; k++) {
+        sum += csr_data[k] * x_out[d.csrIndices[k]];
+      }
+      double residual = rhs[srcRow] - sum;
+      resNormSq += residual * residual;
+      bp[perm[j]] = float(d.rowScale[j] * residual);
+    }
+    itersUsed = iter + 1;
+
+    if (d.refineTol > 0 && resNormSq <= d.refineTol * d.refineTol * std::max(bNormSq, 1e-300)) {
+      break;
+    }
+
+    d.solver->solveLU(d.dataCpu.data(), d.pivotsCpu.data(), bp.data(), n, 1);
+
+    for (int64_t j = 0; j < n; j++) {
+      x_out[j] += d.colScale[j] * double(bp[perm[j]]);
+    }
+  }
+  return itersUsed;
+}
+
 // ---------------------------------------------------------------------------
 // dot(): sparse matrix-vector multiply (CPU, f64, no permutation)
 // ---------------------------------------------------------------------------

sprux/sprux/SpruxFFISolver.h

@@ -91,6 +91,27 @@ class SpruxFFISolver {
   void beginSolve(const double* csr_data, const double* rhs);
   int endSolve(double* x_out);
 
+  /**
+   * Solve with the previously factored matrix (chord Newton).
+   *
+   * Reuses the factored LU data from the most recent solve() or beginSolve()
+   * call. No equilibration, no scatter, no refactorization — just:
+   *   1. Permute RHS by BTF + AMD ordering + cached equilibration scales
+   *   2. GPU solveLU (forward/backward substitution)
+   *   3. CPU f64 iterative refinement
+   *   4. Unpermute result
+   *
+   * For chord Newton: the caller recomputes the residual f with updated
+   * voltages but reuses the Jacobian factorization from the first NR iteration.
+   * The csr_data is still needed for the f64 SpMV in iterative refinement.
+   *
+   * @param csr_data  CSR non-zero values [nnz], f64 (for refinement SpMV)
+   * @param rhs       Right-hand side vector [n], f64
+   * @param x_out     Solution vector [n], f64 (output)
+   * @return Number of refinement iterations performed.
+   */
+  int solveOnly(const double* csr_data, const double* rhs, double* x_out);
+
   /**
    * Sparse matrix-vector multiply: b_out = A @ x (CPU, f64).
    *

sprux/sprux/sprux_c_api.cpp

@@ -197,6 +197,15 @@ int sprux_ffi_dot(sprux_ffi_solver_t h, const double* csr_data, const double* x,
   }
 }
 
+int sprux_ffi_solve_only(sprux_ffi_solver_t h, const double* csr_data, const double* rhs,
+                         double* x_out) {
+  try {
+    return h->solver->solveOnly(csr_data, rhs, x_out);
+  } catch (...) {
+    return -1;
+  }
+}
+
 int sprux_ffi_begin_solve(sprux_ffi_solver_t h, const double* csr_data, const double* rhs) {
   try {
     h->solver->beginSolve(csr_data, rhs);

sprux/sprux/sprux_c_api.h

@@ -197,6 +197,25 @@ int sprux_ffi_solve(sprux_ffi_solver_t h, const double* csr_data, const double*
  */
 int sprux_ffi_dot(sprux_ffi_solver_t h, const double* csr_data, const double* x, double* b_out);
 
+/**
+ * Solve using the previously factored matrix (chord Newton).
+ *
+ * Reuses the LU factorization from the most recent sprux_ffi_solve() call.
+ * No equilibration, no scatter, no refactorization — just permute RHS,
+ * forward/backward substitution, and iterative refinement.
+ *
+ * For chord Newton: the caller recomputes the residual f with updated
+ * voltages but reuses the Jacobian factorization.
+ *
+ * @param h         Solver handle
+ * @param csr_data  CSR non-zero values [nnz], f64 (needed for refinement SpMV)
+ * @param rhs       Right-hand side vector [n], f64
+ * @param x_out     Solution vector [n], f64 (output)
+ * @return Number of refinement iterations, or -1 on error
+ */
+int sprux_ffi_solve_only(sprux_ffi_solver_t h, const double* csr_data, const double* rhs,
+                         double* x_out);
+
 /**
  * Split-phase solve: submit GPU factor+solve asynchronously.
  *

sprux/tests/SpruxFFISolverTest.cpp

@@ -223,6 +223,50 @@ TEST(SpruxFFISolver, SmallMatrix) {
   EXPECT_LT(relRes, 1e-5) << "Small matrix solve failed";
 }
 
+// Test: solveOnly() reuses factored data for chord Newton
+TEST(SpruxFFISolver, SolveOnly) {
+  int32_t n = 3, nnz = 7;
+  vector<int32_t> indptr = {0, 2, 5, 7};
+  vector<int32_t> indices = {0, 1, 0, 1, 2, 1, 2};
+  vector<double> data = {4.0, 1.0, 1.0, 3.0, 1.0, 1.0, 4.0};
+
+  SpruxFFISolver solver(n, nnz, indptr.data(), indices.data(), data.data(), 10);
+
+  // First: full solve to establish factorization
+  vector<double> rhs1 = {5.0, 5.0, 5.0};
+  vector<double> x1(n);
+  solver.solve(data.data(), rhs1.data(), x1.data());
+
+  // Now: solveOnly with a different RHS but same matrix (chord Newton scenario)
+  vector<double> rhs2 = {1.0, 2.0, 3.0};
+  vector<double> x2(n);
+  solver.solveOnly(data.data(), rhs2.data(), x2.data());
+
+  // Verify x2 = A^-1 * rhs2
+  double relRes = 0, bNorm = 0;
+  for (int i = 0; i < n; i++) {
+    double Ax_i = 0;
+    for (int k = indptr[i]; k < indptr[i + 1]; k++) {
+      Ax_i += data[k] * x2[indices[k]];
+    }
+    double r = Ax_i - rhs2[i];
+    relRes += r * r;
+    bNorm += rhs2[i] * rhs2[i];
+  }
+  relRes = sqrt(relRes / bNorm);
+  cout << "  solveOnly: x = [" << x2[0] << ", " << x2[1] << ", " << x2[2] << "]" << endl;
+  cout << "  solveOnly: relative residual = " << scientific << setprecision(3) << relRes << endl;
+  EXPECT_LT(relRes, 1e-5) << "solveOnly failed";
+
+  // Also test solveOnly vs fresh solve — should give same result
+  vector<double> x2_ref(n);
+  solver.solve(data.data(), rhs2.data(), x2_ref.data());
+  double maxDiff = 0;
+  for (int i = 0; i < n; i++) maxDiff = max(maxDiff, abs(x2[i] - x2_ref[i]));
+  cout << "  solveOnly vs solve diff: " << scientific << setprecision(3) << maxDiff << endl;
+  EXPECT_LT(maxDiff, 1e-10) << "solveOnly doesn't match solve";
+}
+
 // Test: dot() computes correct SpMV
 TEST(SpruxFFISolver, DotProduct) {
   string dir = findTestDataDir("c6288_sequence");