Add SupernodalLUFactorization (pure-Julia Schenk–Gärtner supernodal LU via PurePardiso.jl)

[ChrisRackauckas-Claude]

⚠️ Draft — please ignore until reviewed by @ChrisRackauckas.

Prerequisite: the backing package is not yet registered (and pending a possible rename away from PurePardiso — the algorithm struct name here is already trademark-free). CI will fail on dependency resolution until registration; everything below was verified locally with the package dev'd.

What

New sparse LU algorithm SupernodalLUFactorization, backed by PurePardiso.jl as a hard dependency (same pattern as PureKLU): the supernodal left–right-looking LU method of O. Schenk & K. Gärtner (FGCS 20(3) 2004; ETNA 23 2006) — supernodal BLAS-3 LU on the symmetrized pattern, restricted-block static pivoting with iterative-refinement recovery, MC64-style matching + scaling. Pure Julia, no binary dependencies, generic over element types (BigFloat verified to 1e-78 residual through the LinearSolve interface).

Why

It is the strongest solver for "more structured" (PDE-mesh-like) sparse systems: fastest single-threaded full factorization on every mesh-like problem in PurePardiso's benchmark suite — ahead of UMFPACK (0.5–0.8×) and KLU (up to 50× on 3D Poisson), and ahead of MKL PARDISO at 1 thread. The cache-reuse path here maps to PurePardiso's allocation-free numeric refactorization (measured 176 bytes per refactor+solve through the full LinearSolve interface).

Semantics

• reuse_symbolic/check_pattern mirror the KLU/UMFPACK algs; pattern change re-analyzes.
• Static pivoting never aborts, so numerically singular systems would otherwise return finite garbage with Success — a residual check on the (rare) perturbed-pivot path surfaces them as ReturnCode.Infeasible.
• using RecursiveFactorization (e.g. for RFLUFactorization) automatically activates PurePardiso's speed extension (RecursiveFactorization + TriangularSolve panel kernels, −8–18 % refactorization) — no wiring needed here, the extension triggers on package presence.
• threaded=true opt-in supernodal-etree parallel factorization (1.5–3.2× at 8 threads).

Not in this PR

Default-polyalgorithm routing. The natural follow-up: route the "more structure" branch of the sparse default (currently UMFPACKFactorization, gated on Base.USE_GPL_LIBS) to this algorithm — which would (a) speed up the GPL build on mesh-like systems and (b) rescue the non-GPL build from PureKLU-on-3D collapse, making the entire default sparse LU stack pure Julia. That is a defaults decision deserving its own PR and review.

Local verification

• End-to-end through LinearProblem/init/solve!: basic solve, cache-reuse refactorization (allocation-free), pattern-change re-analysis, weak-diagonal matching path, BigFloat, singular→Infeasible, ordering=:nd, threaded=true — all pass.
• GROUP=Core Pkg.test() run locally with the new testset (result pending at PR time; will update).

🤖 Generated with Claude Code

https://claude.ai/code/session_017rr42T1vFRRT8D7DMAmJzf

[ChrisRackauckas-Claude]

Local test status: the initial GROUP=Core run surfaced 4 errors, all in test/Core/resolve.jl — its re-solve sweep iterates every AbstractSparseFactorization subtype and only sparsifies the test matrix for algorithms in its explicit lists, so the new algorithm received a dense Matrix. Fixed by adding SupernodalLUFactorization to both lists (pushed). resolve.jl re-run locally with the full dep set (including RecursiveFactorization loaded, which also exercises PurePardiso's speed extension through the LinearSolve path): full algorithm sweep completes with zero failures. No other Core-group failures involved the new algorithm.

[ChrisRackauckas-Claude]

Updated per direction: the algorithm is now also the default on the structured sparse branch, and the dense-kernel composition with LinearSolve's default dense stack is settled.

Defaults change (new commit): the default polyalgorithm's :UMFPACKFactorization slot now resolves to SupernodalLUFactorization() — the exact pattern already used for :KLUFactorization → PureKLUFactorization(). Explicit UMFPACKFactorization() requests are untouched. Because the slot symbol is unchanged, the established sparse-LU → column-pivoted-QR fallback chain applies automatically (verified: consistent singular system falls back and solves to residual 0.0; inconsistent returns the chain's usual least-squares behavior; explicit alg requests still return Infeasible). Side effect worth calling out: the sparse LU default no longer depends on Base.USE_GPL_LIBS — the non-GPL build previously sent structured systems to scalar-kernel KLU (131–195 s on poisson3d_48); now both builds get the BLAS-3 solver (3.8 s). Honest trade-off: circuit-like large-n systems are ~1.6× slower than UMFPACK on that row of the benchmark suite.

"Use the LinearSolve default dense solver" for the panel kernels — measured result: a plain LAPACK-getrf fast path for the base kernel was benchmarked and not adopted (±5 %, the mandatory block backup cancels getrf's gains at supernode sizes; recorded in the PurePardiso commit history as a negative result). RecursiveFactorization is a genuine win (−8–18 % refactorization) and its kernels activate automatically whenever RF is loaded — verified that using RecursiveFactorization alone triggers the extension (TriangularSolve arrives transitively). That is exactly the condition under which LinearSolve's own dense default prefers RFLU, so the sparse solver's dense panels and the dense default stay on the same engine with zero additional wiring.

Local verification on top of the previous rounds: test/Core/default_algs.jl passes unchanged (slot enum untouched), default-path solve/refactor on structured systems routes through PardisoLU cacheval at 1e-15 residuals, and resolve.jl still passes.

[ChrisRackauckas-Claude]

Regime-mapping data relevant to the defaults change in this PR — two benchmark campaigns using SciMLBenchmarks problem setups (SparsePDE fdmatrix scaling sweep 1D/2D/3D to n≈164k, plus 18 curated SuiteSparse-collection matrices via MatrixDepot), accuracy-aware (res ≤ 1e-6 to count as a win). Full analysis: results/REGIMES.md.

Where KLU is better (validates the existing use_klulike_sparse_structure split): banded/1D at every size; circuits/BTF-reducible at every size tested (to n=171k); anything below n≈1k (3D-like) to n≈5–10k (2D-like); and refactorization in those same regimes.

Where UMFPACK is better: one regime — fully-unsymmetric-pattern but numerically tame matrices (bayer01-class, symfrac≈0), where COLAMD avoids the A+Aᵀ symmetrization penalty (1.5×; SupernodalLU still wins the refactorization there). It won 0 of 45 PDE-sweep rows and 1 of 18 collection matrices. The replacement this PR makes is safe.

Boundary caveat for the current heuristic: large circuits (scircuit n=171k) fall on the "structured" side of use_klulike_sparse_structure but KLU is ~1.9× faster there; density alone can't separate circuits from 2D meshes (torso2 is sparser than scircuit yet SupernodalLU wins it 3×). A predicted-fill probe (KLU-like iff est. nnz(L)/nnz(A) ≲ 4 from the AMD analysis) would fix the routing — suggested follow-up, not in this PR.

Also fixed in PurePardiso while benchmarking: ASIC_100ks exposed a mass-perturbation policy bug (96 % pivots perturbed → res 0.76); matching=:auto now self-heals by retrying with matching (res 4e-12, faster than KLU/UMFPACK/PureKLU on that matrix even counting the discarded attempt). And a robustness note: on goodwin, SupernodalLU is the only accurate solver of the five (KLU res 6e2, MKL PARDISO 0.31, UMFPACK 6e-4).