sparse-rs Internals

sparse-rs is a pure Rust sparse direct solver library providing two independent backends for solving where is a sparse matrix.

Source: sim/sparse-rs/src/

Two backends

KLU

Port of SuiteSparse KLU by Timothy A. Davis. A general-purpose sparse direct solver using:

KLU uses a three-phase pipeline: symbolic analysis (once per matrix pattern), numeric factorization (once per value change), and solve (once per RHS). Refactorization reuses the symbolic structure for 2-5x speedup.

Markowitz

Port of Sparse 1.3 (Kundert, 1988). A circuit-simulation-optimized solver using:

Markowitz is ngspice's default solver. spice-rs uses this backend for the MNA system.

Shared interface

Both backends accept input in CSC (Compressed Sparse Column) format via CscMatrix:

pub struct CscMatrix {
    pub n: usize,
    pub col_ptr: Vec<usize>,
    pub row_idx: Vec<usize>,
    pub values: Vec<f64>,
}

However, in the spice-rs MNA integration, the Markowitz backend is used directly through MarkowitzMatrix rather than through CscMatrix.

KLU Deep Dive

Pipeline overview

Input matrix A (CSC)
        |
   +----v----+
   |   BTF    |  Hopcroft-Karp matching + SCC decomposition
   +----+----+  -> block structure, permutation
        |
   +----v----+
   |   AMD    |  per-block fill-reducing column ordering
   +----+----+  -> column permutation within each block
        |
   +----v----------------+
   |  Gilbert-Peierls     |  left-looking LU with DFS-based sparse triangular solve
   +----+----------------+  -> L, U factors, pivot permutation
        |
   +----v----+
   |  Solve  |  sparse forward/back substitution
   +---------+  -> solution vector x

Phase 1: Symbolic analysis

BTF decomposition

BTF (Block Triangular Form) permutes the matrix into upper block-triangular form. Each diagonal block can be factored independently.

The algorithm:

  1. Hopcroft-Karp maximum matching -- find a permutation that puts nonzeros on the diagonal (maximum transversal).
  2. Strongly connected components (SCC) -- find the SCC decomposition. Each SCC becomes one diagonal block.

Singletons (1x1 blocks) are particularly valuable: they require no LU factorization, just a division. Circuit matrices often have many singletons.

AMD ordering

Within each non-singleton BTF block, AMD computes a column permutation that minimizes fill-in during LU factorization. The algorithm maintains a quotient graph representation and greedily selects the column with minimum degree at each step.

Phase 2: Numeric factorization (Gilbert-Peierls LU)

For each BTF block, the Gilbert-Peierls algorithm computes sparse L and U factors using left-looking factorization:

For each column k:

  1. Sparse triangular solve -- solve where only the nonzero entries of the solution are computed. The nonzero pattern is found by a DFS on the graph of L.
  2. Partial pivoting -- among the entries below the diagonal, select the largest as the pivot.
  3. Split -- entries above the pivot become column k of U; entries at and below become column k of L.

The DFS-based sparse triangular solve is what makes Gilbert-Peierls efficient: it only visits rows that will have nonzero values.

Row scaling

Before factorization, each row is scaled by the reciprocal of its maximum absolute value. This improves numerical stability.

Phase 3: Solve

  1. Apply row scaling
  2. Apply row permutation
  3. Forward substitution with L (unit lower triangular)
  4. Back substitution with U (upper triangular)
  5. Apply column permutation

For multi-block BTF: solve each block's subsystem, then use off-diagonal entries for inter-block back-substitution.

Refactorization

When the sparsity pattern hasn't changed, refactorization reuses the symbolic analysis and pivot permutation. This skips BTF, AMD, DFS-based symbolic structure discovery, and pivot selection. For circuit matrices, refactorization is 2-5x faster than full numeric factorization and 10-100x faster than the full symbolic+numeric pipeline.

Markowitz Deep Dive

Arena-based linked-list matrix

The matrix is stored as doubly-linked lists (by row and by column) with all elements in a flat arena:

pub struct MarkowitzMatrix {
    size: usize,
    elements: Vec<Element>,       // arena, index 0 = dummy sentinel
    first_in_row: Vec<u32>,       // 1-indexed head pointers
    first_in_col: Vec<u32>,
    diag: Vec<u32>,               // direct pointers to diagonal elements
    markowitz_row: Vec<i32>,      // nonzero count per row
    markowitz_col: Vec<i32>,      // nonzero count per column
    markowitz_prod: Vec<i64>,     // row_count * col_count per row/col
    int_to_ext_row: Vec<usize>,   // permutation: internal -> external
    int_to_ext_col: Vec<usize>,
    ext_to_int_row: Vec<usize>,   // permutation: external -> internal
    ext_to_int_col: Vec<usize>,
}

pub struct Element {
    pub real: f64,
    pub imag: f64,                // for AC analysis
    pub row: u32,
    pub col: u32,
    pub next_in_row: u32,         // arena index (NONE = u32::MAX for NULL)
    pub next_in_col: u32,
}

All arrays are 1-indexed to match the C code. NONE = u32::MAX replaces NULL pointers.

Markowitz pivot criterion

The Markowitz criterion selects the pivot that minimizes , where counts are in the unreduced submatrix. This product approximates the fill-in the pivot will cause. A product of 0 means a singleton -- zero fill-in.

The four-level pivot cascade

Level 1: Search for singletons. Row singletons and column singletons. Zero fill-in guaranteed.

Level 2: Quick diagonal search. Scan diagonal elements for smallest Markowitz product with quick magnitude threshold check. Strong diagonal preference.

Level 3: Careful diagonal search. Same as Level 2 but with full column-maximum threshold check.

Level 4: Full matrix search. Scan every element. Last resort, guarantees finding a pivot if matrix is non-singular.

Factorization

order_and_factor() (first time): For each elimination step: select pivot, swap rows/columns, update Markowitz counts, eliminate column and row, create fill-in as needed.

L stores 1/pivot on the diagonal (not the pivot itself). U has implicit unit diagonal.

factor() (refactorization): Reuses the pivot ordering from order_and_factor(). Only recomputes numeric values. This is the normal path after the first NR iteration.

Solve

Forward substitution (L) then backward substitution (U), operating on the RHS vector in-place. The permutation arrays map between external and internal ordering.

Complex mode (AC analysis)

Elements have both real and imag fields. The factorization and solve operate on complex values for AC analysis. In DC/transient mode, imag is always zero.

Benchmarks

The sparse-eval crate benchmarks sparse-rs against SuiteSparse C implementations.

sparse-rs KLU vs SuiteSparse C KLU

Correctness is validated: solutions match to machine precision on all test matrices.

Performance: sparse-rs is typically within 2x of SuiteSparse C. The gap comes from hand-optimized C with careful cache layout, bounds checking in safe Rust, and decades of micro-optimization in SuiteSparse.

For circuit simulation, the solver is not the bottleneck -- device model evaluation dominates. The 2x overhead on the solver translates to a much smaller overhead on total simulation time.

Markowitz vs KLU

Markowitz KLU
Diagonal preference Strong None (relies on BTF)
Optimized for Circuit matrices General sparse matrices
Refactorization Reuse pivot ordering Reuse symbolic structure
Larger matrices Slower Faster (BTF + cache-friendly)

For typical SPICE circuits (tens to hundreds of nodes), both converge to the same solution.

Refactorization performance

Both backends support refactorization -- typically 2-5x faster than full numeric factorization. This is the critical optimization for SPICE, where the NR loop calls the solver many times with the same matrix structure.

Companion benchmarks

cargo run --release --bin amd-compare    # AMD ordering quality comparison
cargo run --release --bin solver-compare # Full solve pipeline comparison
cargo run --release --bin sparse-eval    # Main benchmark suite