rslab (Python bindings)

NumPy/SciPy bindings for RSLAB, a pure-Rust sparse direct solver and preconditioner: complex/real symmetric LDL^T (Bunch-Kaufman), unsymmetric LU, and a KLU-style path for circuit-shaped matrices. A thin wrapper, all numeric work happens in Rust.

Install

pip install rslab

Usage

import numpy as np
import scipy.sparse as sp
import rslab

# Symmetric system (real or complex; the dtype selects the path).
A = sp.random(5000, 5000, density=1e-3, format="csc") + sp.eye(5000) * 10
A = A + A.T
b = np.random.rand(5000)

# One-shot solve.
x = rslab.spsolve(A, b)

# Factor once, solve many right-hand sides.
f = rslab.ldlt(A)
x1 = f.solve(b)
X = f.solve_many(np.random.rand(5000, 8))   # n x nrhs

print(f.n, f.factor_nnz, f.inertia, f.dtype)

Complex-symmetric matrices (EM/FEM, PARDISO mtype 6) work identically:

A = A.astype(np.complex128); A.data += 1j * 0.3 * A.data.real
x = rslab.ldlt(A).solve(np.ones(A.shape[0], dtype=np.complex128))

Unsymmetric matrices use the LU path:

f = rslab.lu(A_general)
x = f.solve(b)

Circuit-shaped matrices (MNA / SPICE-class: very sparse, unsymmetric, near-triangularizable) use the KLU path, bit-deterministic, with a numeric-only refactor for fixed-pattern sweeps:

f = rslab.klu(A_circuit)
x = f.solve(b)
A_circuit.data *= 1.5            # frequency sweep: same pattern, new values
f.refactor(A_circuit.data)       # no symbolic work, no pivot search
x2 = f.solve(b)
y = f.solve_transpose(b)         # A.T @ y = b on the same factors (adjoint)

solve_transpose is the plain transpose; for the conjugate-transpose adjoint use f.solve_transpose(b.conj()).conj().

Preconditioner mode

Never-fail static pivoting plus iterative refinement for hard/indefinite systems:

f = rslab.ldlt(A, preconditioner=1e-4)
x = f.solve(b, refine=20)        # refine against the original A

Configuration

ldlt, lu and spsolve take a Settings object (or the same keywords directly), klu a KluSettings. By default the factor uses RSLAB's deterministic heuristic pick: the adaptive ordering plus an exact nested-dissection bakeoff on large systems, and at most 4 workers (the calibrated count after a one-time rslab.install_diagnose()). Keywords override the pick:

s = rslab.Settings(ordering="metis", threads=2, preconditioner=1e-4)
f = rslab.ldlt(A, settings=s)
f = rslab.lu(A, ordering="amd", pivot_u=0.5)        # the same, as keywords
print(s.to_dict())

KluSettings: pivot_tol (1e-3), row_scaling (on), btf (on), matching (on: MC64 row matching as the BTF transversal), parallel (None = structural auto gate, True / False force), interrupt.

Settings a path ignores are reported under diagnostics()["warnings"]. solve and solve_many on Ldlt and Lu handles are supernodal and tree-parallel (leaf subtrees of the elimination tree in parallel, parallel sections inside the wide top separators), bit-identical for every thread count. Throughput is always reported: diagnostics()["rates"] holds the analysis, factorization and solve rates in million unknowns per second (MDOF/s), the factorization flop rate (GFlop/s) and the factor-entry rate (Mnnz/s); the summary line and the info log carry the factor rate too.

Symbolic reuse

The analysis depends only on the sparsity pattern. Pay it once and factor each value set of a sweep on it:

sym = rslab.analyze(A, path="lu")                # LdltSymbolic / LuSymbolic / KluSymbolic
print(sym.factor_nnz, sym.estimate_memory()["factor_mb"])
for omega in frequencies:
    f = sym.factor((K + 1j * omega * C).data)    # same pattern, new values
    x = f.solve(b)

Krylov solvers

gmres, gmres_block, cocg and cocr run on any matrix, optionally preconditioned by any factor handle (an incomplete factor, or the factor of a nearby matrix):

M = rslab.lu(A, drop_tol=1e-3)                   # ILU-style preconditioner
x, converged, iters, res, stop = rslab.gmres(A, b, M, tol=1e-10)
r = rslab.gmres(A_next, b, M, recycle=M.recycle(8))   # reuse M, deflate across solves

Logging

rslab.set_log_level("info")                      # or the RLA_LOG environment variable
log = logging.getLogger("rslab")
rslab.set_log_sink(lambda level, msg: log.log(logging.getLevelName(level.upper()), msg))

API

The full reference, generated from the docstrings, is in docs/api.md (help(rslab.ldlt) etc. show the same text).

Factor handles share solve(b, refine=0, target=None, measure="normwise"), solve_many(B), gmres, gmres_block, cocg, cocr, recycle(k), diagnostics() and the attributes n, factor_nnz, n_perturbed, dtype; Ldlt adds inertia, Klu adds n_blocks, solve_transpose(b) and the numeric-only refactor(data).

Supported dtypes: float64, float32, complex128, complex64.

License

MIT.