Backends and Factorization Reuse

The Y-bus admittance matrix and Newton-Raphson Jacobian are stored and factored as sparse matrices, not dense ones. This page covers why that matters, how factorization work is reused across solves, and what the five interchangeable linear-solver backends are.

Why sparse

At 2,605 nodes the original dense solver was over 100,000x slower than power-grid-model; the sparse rewrite closed that to roughly an order of magnitude on a cold, single-shot solve. The underlying grid is sparse — each bus only connects to a handful of neighbors — so a dense representation was doing asymptotically unnecessary work regardless of how fast the constant-factor arithmetic was.

Two things make this work, rather than just "swap in a sparse matrix type":

Sparse-aware assembly. Both the Jacobian build (solver::build_jacobian_triplets) and the linear initial-guess warm start (network::linear_initial_guess) walk each bus's actual admittance neighbors (via network::YBusSparse::row) instead of looping over every possible bus pair. The O(n²)/O(m²) assembly cost has to go too, or a sparse solve alone doesn't fix the bottleneck.
Symbolic factorization reuse. A Newton-Raphson Jacobian has the same sparsity pattern every iteration — same bus topology, only numeric values change — so solver::newton_raphson computes the symbolic factorization (fill-reducing ordering) once and reuses it for a cheap numeric-only refactorization on each iteration (sparse::RealSparseSystem), mirroring what PGM's own solver does internally.

Inside KLU walks through exactly what "symbolic factorization" and "numeric-only refactorization" mean, step by step, on a real Jacobian.

Reusing factorization across repeated solves

A single newton_raphson/newton_raphson_with_backend call reuses its symbolic factorization across its own NR iterations, but starts cold on every call — re-deriving the fill-reducing ordering from scratch. That is fine for a genuinely one-off solve and wasteful for anything that solves the same topology repeatedly: a time series, a batch of scenarios, contingency analysis. In those, only bus values (p_spec, q_spec, voltage guess) change between calls, not the topology.

solver::PersistentSolver extends the reuse across calls. Construct one per topology, then solve as many times as needed; only the first call pays for symbolic factorization.

#![allow(unused)]
fn main() {
use gridoxide::solver::{JacobianBackend, PersistentSolver};

let mut solver = PersistentSolver::new(JacobianBackend::Klu);
for scenario in scenarios {
    apply_scenario(&mut buses, scenario); // changes p_spec/q_spec only
    solver.solve(&mut buses, &ybus, 1e-6, 20);
}
}

Call .reset() (or construct a new solver) if the topology itself changes between solves. This is a meaningful win on real-world grids, since a cold solve otherwise redoes COLAMD/AMD/BTF ordering from scratch every call. examples/bench_network.rs exposes it as an optional warm mode; its default cold mode still measures "N independent flat-start solves with no shared state," a different and also legitimate number. See Benchmarking for the measured warm-vs-cold figures.

The backend interface

src/sparse.rs is the thin backend wrapper around faer, and intentionally the only file that imports faer types directly, so a different sparse-solver backend can be swapped in behind the same interface without touching the rest of the codebase.

solver::newton_raphson always uses the default Scalar (faer-backed) path. newton_raphson_with_backend additionally accepts four alternatives via solver::JacobianBackend:

`Block`

src/block_sparse.rs, no extra build requirements. Groups each bus's own (angle, voltage-magnitude) unknowns into one dense 2×2 block, mirroring power-grid-model's block-per-bus matrix structure, with a hand-written Gilbert-Peierls sparse block LU (block_sparse::BlockLu). Symmetric power flow only. Consistently faster than Scalar.

`Klu`

src/sparse_klu.rs, needs cargo build --features klu. The same scalar Jacobian as Scalar, solved by SuiteSparse's KLU instead of faer, vendored and compiled from source (vendor/suitesparse/) rather than depending on a third-party Rust wrapper crate. Needs a C compiler and libclang (for bindgen) at build time.

KLU and BTF (one of KLU's own dependencies) are LGPL-2.1-or-later. See Provenance and Licensing.

`KluNative`

src/klu_native/, no extra build requirements. A from-scratch Rust translation of the same KLU algorithm Klu links over FFI: BTF block-triangular preprocessing, per-block AMD ordering, a partial-pivoting Gilbert-Peierls LU kernel with Eisenstat-Liu pruning, and cheap numeric-only refactorization — all faithfully ported, not a simplified reimplementation. No C compiler or libclang needed, so unlike Klu it is always built.

Validated end-to-end against real KLU on all 13 real MATPOWER benchmark cases (14 to 9,241 buses): identical iteration counts and identical converged voltages on every case. One known, documented gap — row scaling (klu_native::scale) is ported and independently tested but not yet wired into the factor/refactor path (see src/klu_native/mod.rs's module doc comment). That is a numerical-stability preconditioning step, not a correctness one, so it doesn't affect the results above.

Inside KLU is a walkthrough of this port specifically.

`Pardiso`

src/sparse_pardiso.rs, needs cargo build --features pardiso and MKLROOT set at build time. The same scalar Jacobian as Scalar, solved by Intel oneMKL's PARDISO direct solver. Unlike Klu, nothing is vendored — MKL is proprietary (Intel Simplified Software License, not OSS), so this only dynamically links a locally-installed oneMKL (libmkl_rt.so, discovered via MKLROOT, e.g. source /opt/intel/oneapi/setvars.sh) and generates FFI bindings via bindgen against that install's own mkl_pardiso.h. No MKL header or source is copied into this repo.

PARDISO's C API is one function called repeatedly with different phase values against a persistent opaque handle, rather than KLU's separate analyze/factor/refactor/solve functions. mtype = 11 (real, nonsymmetric) and iparm[34] = 1 (0-based indexing) are the two settings that matter for matching gridoxide's CSR/CSC conventions.

Not built or tested in CI — no CI runner has MKL installed — so this is a local/manual-verification-only backend.

How they compare

All four alternatives are strictly parallel to Scalar, not replacements. A bug in any of them can't affect newton_raphson's default behavior, and every existing test keeps using Scalar unless it explicitly opts into a different backend.

All five backends produce identical converged voltages at every scale — these are purely performance comparisons, not correctness trade-offs. In rough terms:

Block, Klu, and KluNative are all meaningfully faster than Scalar.
Klu and KluNative land close to each other, slightly ahead of Block.
Pardiso carries a largely size-independent fixed setup cost from its default matching/scaling preprocessing, making it the slowest backend at small problem sizes — even behind Scalar — though it scales better than Scalar as node count grows.

PGM is clearly faster than any gridoxide backend on synthetic radial-distribution/LV topology, a real, standing gap this project hasn't closed. That gap doesn't hold universally, though: on real-world transmission-topology grids, gridoxide's Klu backend is frequently faster than lightsim2grid's own KLU-backed C++ solver. The comparison depends on topology, not just implementation language.

Benchmarking points at the full measured numbers, exact ratios, and how to reproduce them.

Power Grid Book