Internals

KLUWrapper

A small, allocation-aware wrapper over libklu (provided by SuiteSparse_jll) designed for the access patterns of PowerNetworkMatrices:

• Cache the symbolic and numeric factorizations of an SPD/asymmetric sparse matrix and reuse them across many solves.
• Refactor (numeric only, or full) without re-allocating.
• Solve dense and sparse right-hand sides without materializing N×N intermediates when the RHS is structurally sparse.

This module is intentionally lighter than KLU.jl: it owns no Julia-side copies of the matrix values, exposes the symbolic/numeric split directly, and binds only the SuiteSparselong (`klul*,kluzl_*`) entry points used by the package.

KLU_POOL_DEBUG :: Bool

Compile-time gate for the KLU wrapper's runtime diagnostics — currently just the precondition snapshot in solve! that fires when klu_l_solve returns FALSE. When false, the snapshot is folded out by @static if and contributes zero runtime cost. Flip to true when reproducing the libklu cross-cache concurrency failure documented on _LIBKLU_LOCK.

This is a const rather than a Preferences.@load_preference flag so that toggling it forces a precompile rebuild and the production binary never carries the diagnostic code.

_LIBKLU_LOCK :: ReentrantLock

Process-wide lock that serializes every libklu ccall. libklu corrupts internal state under concurrent access even when the caller hands out distinct Numeric/Symbolic/Common triples per thread (i.e. the access pattern KLU's user guide implies is supported). The corruption manifests two ways: an intermittent KLU_INVALID return with all input pointers still valid both pre- and post-call, and a SIGSEGV inside klu_l_solve (klu_solve.c:118 in v7.8.3, the row-permutation read in the nrhs == 1 chunk). The pre-call snapshot dump in solve! made both modes reproducible on macOS and confirmed the deterministic Windows-MinGW failure was the same bug. This lock is the only mechanism we have evidence for that prevents both modes.

A cached KLU linear solver designed for repeated solves against the same sparse matrix structure. numeric_refactor! and solve! allocate nothing once the cache is built.

Type parameters:

• Tv ∈ {Float64, ComplexF64} selects the real/complex KLU path (klu_*_factor/klu_z*_factor).
• Ti ∈ {Int32, Int64} selects the index-type entry-point family (klu_* for int/Int32, klu_l_* for SuiteSparse_long/Int64). The cache's colptr/rowval/col_map are stored in this type.

reuse_symbolic controls whether symbolic_refactor! keeps the analysis; check_pattern adds a structural-equality check on refactor calls and is only consulted when reusing.

KLULinSolveCache(A; reuse_symbolic=true, check_pattern=true)

Build a cache for the sparse matrix A. The cache's index type is taken from A: SparseMatrixCSC{Tv, Int32} ⇒ KLULinSolveCache{Tv, Int32}, SparseMatrixCSC{Tv, Int64} ⇒ KLULinSolveCache{Tv, Int64}. Allocates structural arrays and runs the corresponding klu_defaults/klu_l_defaults initializer, but does not factorize. Call full_factor! (or symbolic_factor! followed by numeric_refactor!) before solve!.

A finalizer frees libklu handles on GC; call Base.finalize(cache) to release them eagerly. Releasing the handles leaves Julia-side state intact, so the cache can be re-factorized via symbolic_factor!/numeric_refactor! or full_factor!.

\(cache::KLULinSolveCache, B) -> X

Allocating solve, mirroring LinearAlgebra.Factorization's API.

_ensure_scratch!(cache, block) -> Nothing

Ensure cache.scratch is at least n × block and cache.col_map length block. Grows in place; reuses across solve_sparse! calls.

Release the libklu numeric and symbolic handles held by cache, leaving the Julia-side fields (colptr, rowval, common, scratch, col_map) intact so the cache remains structurally valid and re-factorable. Idempotent: a second call hits the C_NULL guards. Used both by symbolic_factor! mid-life (drop old handles before re-analyzing) and by the GC finalizer.

Drop a corrupted numeric handle and rebuild it from the values cached in cache.nzval. Used by solve_sparse! on the KLU_INVALID retry path, where libklu state has been observed to corrupt without the caller having A in scope. Requires that a numeric factor has been built before (so cache.nzval is populated) and the symbolic factor is still valid.

condest!(cache) -> Float64

Compute the 1-norm condition-number estimate of the cached factorization via libklu's klu_condest. The result lands in cache.common[].condest and is also returned. Cost is roughly two extra solves; use sparingly.

Useful when deciding whether iterative refinement is worth running, or for flagging near-singular Jacobians in Newton-Raphson loops.

Float64 only.

full_factor!(cache, A) -> cache

Run a fresh symbolic analysis followed by a numeric factorization on A. Equivalent to symbolic_factor!(cache, A); numeric_refactor!(cache, A). Use this on a freshly constructed cache, or after _free_klu_handles! has cleared the handles, to bring the cache to a factored state.

full_refactor!(cache, A) -> cache

Refresh both the symbolic and numeric factorizations on A. Defers to symbolic_refactor! (which reuses the existing analysis when cache.reuse_symbolic is set) followed by numeric_refactor!. Use this when the matrix values have changed; if the structure has also changed and the cache was built with reuse_symbolic = false, the symbolic analysis is rerun as well.

is_factored(cache::KLULinSolveCache) -> Bool

Return true when cache holds both a symbolic and a numeric factorization ready for solve! / tsolve! / solve_sparse!. Returns false after construction (before full_factor!) or after the libklu handles have been finalized.

klu_factorize(A; reuse_symbolic=true, check_pattern=true) -> KLULinSolveCache

Build a cache for A and immediately compute the full factorization.

numeric_refactor!(cache, A)

Compute (or refresh) the numeric factorization. The first call after symbolic_factor! invokes klu_*_factor; subsequent calls invoke klu_*_refactor and reuse the existing numeric struct.

rcond!(cache) -> Float64

Compute the cheap reciprocal-condition estimate min(|diag(U)|)/max(|diag(U)|) via libklu's klu_rcond. The result lands in cache.common[].rcond and is also returned. Faster than condest! but less reliable as a conditioning indicator.

Float64 only.

solve!(cache, B) -> B

Solve A · X = B in place. B::StridedVecOrMat{Tv} must have first-dimension size equal to cache.n and unit stride in the first dimension. Multiple columns of B are handled in a single libklu call.

solve_sparse!(cache, B, out; block=64) -> out

Solve A · X = B for a SparseMatrixCSC right-hand side, writing the result into out. Empty columns of B are not solved — out's corresponding columns are filled with zeros. Non-empty columns within each chunk of block consecutive RHS columns are packed into a dense scratch and solved in a single libklu call. out may be any AbstractMatrix{Tv} with shape (cache.n, size(B,2)), including a view into a larger matrix. For an allocating variant, see solve_sparse.

The block chunk size bounds the working set so that processing an n × nrhs sparse RHS requires only O(n · block) extra memory regardless of nrhs. The cache reuses its packing buffer across calls; warm calls allocate nothing in the solver.

Not thread-safe (mutates per-cache scratch). Callers should serialize access through a per-cache lock if invoked from multiple threads.

Allocating wrapper around solve_sparse!.

sort_factors!(cache) -> cache

Sort the columns of the cached L and U factors in place via libklu's klu_sort / klu_l_sort. KLU's numeric phase stores factor columns in arbitrary order; sorting once after the first factor improves cache locality on every subsequent solve! / tsolve!. The cost is O(nnz_factor) and is amortized over many repeated solves — a win whenever the cache is used for ≥ a few solves on the same factorization.

Idempotent and refactor-stable: sorting after the initial factor persists through numeric_refactor! because refactor preserves the column layout. Only call this if the cache will be reused for multiple solves; for a one-shot solve it is pure overhead.

Float64 only.

symbolic_factor!(cache, A)

Free any cached symbolic/numeric factor, replace the structural arrays with A's pattern, and run klu_analyze / klu_l_analyze.

symbolic_refactor!(cache, A)

If cache.reuse_symbolic, optionally verify the structure matches and reuse the existing analysis. Otherwise, rerun symbolic_factor!.

tsolve!(cache, B; conjugate=false) -> B

In-place solve Aᵀ · X = B (or Aᴴ · X = B when conjugate=true on the complex path). Same shape requirements as solve!. The conjugate keyword is ignored on the real path.

@klu_lock expr

Evaluate expr while holding _LIBKLU_LOCK. Wrap every libklu ccall so that no two libklu entries can run concurrently in the process.

AccelerateWrapper

A small, allocation-aware wrapper over Apple's libSparse.dylib (/System/Library/Frameworks/Accelerate.framework/.../libSparse.dylib) designed for the access patterns of PowerNetworkMatrices:

• Cache the symbolic and numeric factorizations of a general (unsymmetric) sparse matrix (LU with threshold partial pivoting + Inf-norm equilibration scaling) and reuse them across many solves.
• Refresh the numeric factor (numeric_refactor!) while keeping the symbolic analysis, without re-allocating the structural arrays.
• Solve dense and sparse right-hand sides in place, with the sparse path packing only non-empty RHS columns into a bounded scratch block.
• Compute A·X and A·x directly via libSparse's SparseMultiply.

This module is intentionally lighter than the upstream AppleAccelerate.jl package: it owns no high-level Julia wrappers over libSparse, exposes the symbolic/numeric split directly, binds only the entry points used by PNM, and is compile-gated to macOS so non-Apple builds never codegen the @ccall sites.