Route the structured-sparse default LU to SupernodalLUFactorization (#1109)

Renames the DefaultAlgorithmChoice enum entry UMFPACKFactorization ->
SupernodalLUFactorization (and the matching DefaultLinearSolverInit
cacheval field), so the structured-sparse slot of the default
polyalgorithm is the pure-Julia supernodal left-right-looking LU.
Explicit UMFPACKFactorization() requests are unchanged; the enum entry
is not public API (users select algorithms by type, not by enum), so
this is not a breaking change.

Both sparse LU default slots now resolve to pure-Julia solvers, so the
LU side of the sparse default no longer depends on Base.USE_GPL_LIBS;
previously the non-GPL build sent structured systems to the scalar KLU
path (131-195s on poisson3d_48 vs 3.8s here). The QR side keeps its
SPQR gating.

Benchmark basis (SciMLBenchmarks SparsePDE fdmatrix sweep 1D/2D/3D to
n~164k, plus 18 SuiteSparse-collection matrices, accuracy-aware with
res <= 1e-6 required to count as a win): UMFPACK won 0 of 45 PDE-sweep
rows and 1 of 18 collection matrices (bayer01: fully unsymmetric
pattern, numerically tame - the one regime where COLAMD-style
unsymmetric ordering beats pattern symmetrization; SupernodalLU still
won its refactorization). KLU keeps banded/1D, circuits, and small n,
which is the klulike side of the split and is unaffected.

Known boundary case: large circuit matrices (scircuit, n=171k) fall on
the structured side of use_klulike_sparse_structure but KLU is ~1.9x
faster there. Density alone cannot separate circuits from 2D meshes
(torso2 is sparser than scircuit yet SupernodalLU wins it 3x); a
predicted-fill probe would fix the routing and is left as a follow-up.

GROUP=Core passes.

Co-authored-by: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-authored-by: Claude Fable 5 <noreply@anthropic.com>

Author: 3 people

docs/src/solvers/solvers.md

@@ -63,8 +63,12 @@ when reduced precision is acceptable for the factorization step.
 
 For sparse LU-factorizations, `PureKLUFactorization` (a pure-Julia KLU with no
 SuiteSparse dependency, the default) if there is less structure to the sparsity
-pattern and `UMFPACKFactorization` if there is more structure. The SuiteSparse-backed
-`KLUFactorization` remains available as an explicit alternative.
+pattern and `SupernodalLUFactorization` (a pure-Julia supernodal
+left–right-looking LU implementing the Schenk & Gärtner method, vendored in
+`src/SupernodalLU`, no binary dependency) if there is more structure — so the
+default sparse LU stack is entirely pure Julia and no longer depends on
+`Base.USE_GPL_LIBS`. The SuiteSparse-backed `KLUFactorization` and
+`UMFPACKFactorization` remain available as explicit alternatives.
 For sparse QR-factorizations (used for non-square, rank-deficient, or least-squares
 systems, and as the fallback when the sparse LU hits a (near-)singular matrix), the
 same less-structure/more-structure split as the LU case applies:
@@ -248,8 +252,9 @@ SparseColumnPivotedQRFactorization
     left–right-looking sparse LU method of Schenk & Gärtner (vendored
     self-contained in `src/SupernodalLU`, no binary dependency). It is the
     strongest choice for "more structured" (PDE-mesh-like) sparse systems,
-    where it outperforms both `UMFPACKFactorization` and `KLUFactorization`.
-    It is not yet wired into the default polyalgorithm; request it explicitly.
+    where it outperforms both `UMFPACKFactorization` and `KLUFactorization`,
+    and it is the default for such systems in the polyalgorithm (the
+    `UMFPACKFactorization` slot resolves to it).
 
 !!! note
 

ext/LinearSolveSparseArraysExt.jl

@@ -1102,42 +1102,35 @@ function LinearSolve.use_klulike_sparse_structure(A::AbstractSparseMatrixCSC, b)
         (size(b, 1) <= 10_000 && length(nonzeros(A)) / length(A) < 2.0e-4)
 end
 
-@static if Base.USE_GPL_LIBS
-    function LinearSolve.defaultalg(
-            A::AbstractSparseMatrixCSC{<:Union{Float64, ComplexF64}, Ti}, b,
-            assump::OperatorAssumptions{Bool}
-        ) where {Ti}
-        klulike = LinearSolve.use_klulike_sparse_structure(A, b)
-        if assump.issq
-            # Less structure → PureKLU (pure-Julia); more structure → UMFPACK.
-            if klulike
-                LinearSolve.DefaultLinearSolver(LinearSolve.DefaultAlgorithmChoice.KLUFactorization)
-            else
-                LinearSolve.DefaultLinearSolver(LinearSolve.DefaultAlgorithmChoice.UMFPACKFactorization)
-            end
+function LinearSolve.defaultalg(
+        A::AbstractSparseMatrixCSC{<:Union{Float64, ComplexF64}, Ti}, b,
+        assump::OperatorAssumptions{Bool}
+    ) where {Ti}
+    klulike = LinearSolve.use_klulike_sparse_structure(A, b)
+    return if assump.issq
+        # Less structure → PureKLU; more structure → the supernodal LU.  Both
+        # are pure Julia, so the sparse LU default does not depend on
+        # Base.USE_GPL_LIBS.
+        if klulike
+            LinearSolve.DefaultLinearSolver(LinearSolve.DefaultAlgorithmChoice.KLUFactorization)
         else
-            # Same split for sparse QR: less structure → SparseColumnPivotedQR
-            # (pure-Julia); more structure → SuiteSparse SPQR (`QRFactorization`).
+            LinearSolve.DefaultLinearSolver(LinearSolve.DefaultAlgorithmChoice.SupernodalLUFactorization)
+        end
+    else
+        @static if Base.USE_GPL_LIBS
+            # Sparse QR: less structure → SparseColumnPivotedQR (pure-Julia);
+            # more structure → SuiteSparse SPQR (`QRFactorization`).
             if klulike
                 LinearSolve.DefaultLinearSolver(LinearSolve.DefaultAlgorithmChoice.SparseColumnPivotedQRFactorization)
             else
                 LinearSolve.DefaultLinearSolver(LinearSolve.DefaultAlgorithmChoice.QRFactorization)
             end
-        end
-    end
-else
-    function LinearSolve.defaultalg(
-            A::AbstractSparseMatrixCSC{<:Union{Float64, ComplexF64}, Ti}, b,
-            assump::OperatorAssumptions{Bool}
-        ) where {Ti}
-        # No SuiteSparse (UMFPACK/SPQR) available: always use the pure-Julia solvers.
-        if assump.issq
-            LinearSolve.DefaultLinearSolver(LinearSolve.DefaultAlgorithmChoice.KLUFactorization)
         else
+            # No SuiteSparse SPQR available: pure-Julia sparse QR throughout.
             LinearSolve.DefaultLinearSolver(LinearSolve.DefaultAlgorithmChoice.SparseColumnPivotedQRFactorization)
         end
     end
-end # @static if Base.USE_GPL_LIBS
+end
 
 # SPQR Handling
 function LinearSolve.init_cacheval(

src/LinearSolve.jl

@@ -310,7 +310,7 @@ EnumX.@enumx DefaultAlgorithmChoice begin
     DirectLdiv!
     SparspakFactorization
     KLUFactorization
-    UMFPACKFactorization
+    SupernodalLUFactorization
     KrylovJL_GMRES
     GenericLUFactorization
     RFLUFactorization

src/default.jl

@@ -10,7 +10,7 @@ mutable struct DefaultLinearSolverInit{
     DirectLdiv!::T4
     SparspakFactorization::T5
     KLUFactorization::T6
-    UMFPACKFactorization::T7
+    SupernodalLUFactorization::T7
     KrylovJL_GMRES::T8
     GenericLUFactorization::T9
     RFLUFactorization::T10
@@ -35,7 +35,7 @@ mutable struct DefaultLinearSolverInit{
     a_backup_allocated::Bool  # true once A_backup has been replaced with a private buffer
     fell_back_to_qr::Bool  # true after QR fallback; reuse QR until matrix is refreshed
     # Persistent-nonstructural-zero reduction state, shared across the sparse
-    # sub-algorithm slots (KLU/UMFPACK + the QR fallback all factor the one reduced
+    # sub-algorithm slots (the two sparse LU slots + the QR fallback all factor
     # matrix). `nothing` when inactive / non-sparse. See `init_sparse_reduction`.
     sparse_reduction::TR
 end
@@ -458,8 +458,8 @@ function algchoice_to_alg(alg::Symbol)
         # PureKLU. The SuiteSparse `KLUFactorization` is unchanged and remains
         # available when requested explicitly.
         PureKLUFactorization()
-    elseif alg === :UMFPACKFactorization
-        UMFPACKFactorization()
+    elseif alg === :SupernodalLUFactorization
+        SupernodalLUFactorization()
     elseif alg === :KrylovJL_GMRES
         KrylovJL_GMRES()
     elseif alg === :GenericLUFactorization
@@ -602,7 +602,7 @@ defaultalg_symbol(::Type{<:QRFactorization{ColumnNorm}}) = :QRFactorizationPivot
 const _SPARSE_ONLY_ALGORITHMS = Symbol.(
     (
         DefaultAlgorithmChoice.KLUFactorization,
-        DefaultAlgorithmChoice.UMFPACKFactorization,
+        DefaultAlgorithmChoice.SupernodalLUFactorization,
         DefaultAlgorithmChoice.SparspakFactorization,
         DefaultAlgorithmChoice.CHOLMODFactorization,
         DefaultAlgorithmChoice.SparseColumnPivotedQRFactorization,
@@ -616,7 +616,7 @@ const _SPARSE_ONLY_ALGORITHMS = Symbol.(
 const _SPARSE_LU_FALLBACK_ALGORITHMS = Symbol.(
     (
         DefaultAlgorithmChoice.KLUFactorization,
-        DefaultAlgorithmChoice.UMFPACKFactorization,
+        DefaultAlgorithmChoice.SupernodalLUFactorization,
         DefaultAlgorithmChoice.SparspakFactorization,
     )
 )
@@ -717,7 +717,7 @@ end
     _do_sparse_qr_fallback(cache::LinearCache, alg, sol, reason::Symbol)
 
 Perform column-pivoted sparse QR (`SparseColumnPivotedQRFactorization`) fallback
-after a sparse LU (`KLUFactorization`, `UMFPACKFactorization`,
+after a sparse LU (`KLUFactorization`, `SupernodalLUFactorization`,
 `SparspakFactorization`) solve failed or produced non-finite output.
 
 Sparse LU does not modify `cache.A` in place — UMFPACK and KLU wrap a
@@ -1125,7 +1125,7 @@ end
                     DefaultAlgorithmChoice.LUFactorization,
                     DefaultAlgorithmChoice.QRFactorization,
                     DefaultAlgorithmChoice.KLUFactorization,
-                    DefaultAlgorithmChoice.UMFPACKFactorization,
+                    DefaultAlgorithmChoice.SupernodalLUFactorization,
                     DefaultAlgorithmChoice.LDLtFactorization,
                     DefaultAlgorithmChoice.SparspakFactorization,
                     DefaultAlgorithmChoice.BunchKaufmanFactorization,

test/Core/default_algs.jl

@@ -121,7 +121,7 @@ prob = LinearProblem(sprand(1000, 1000, 0.5), zeros(1000))
 solve(prob)
 
 @test LinearSolve.defaultalg(sprand(11000, 11000, 0.001), zeros(11000)).alg ===
-    LinearSolve.DefaultAlgorithmChoice.UMFPACKFactorization
+    LinearSolve.DefaultAlgorithmChoice.SupernodalLUFactorization
 prob = LinearProblem(sprand(11000, 11000, 0.5), zeros(11000))
 solve(prob)
 
@@ -491,7 +491,7 @@ let
     # Just past the fast-path boundary on a dense matrix → UMFPACK
     A_past = sprand(1_001, 1_001, 0.5) + I
     @test LinearSolve.defaultalg(A_past, rand(1_001), LinearSolve.OperatorAssumptions(true)).alg ===
-        LinearSolve.DefaultAlgorithmChoice.UMFPACKFactorization
+        LinearSolve.DefaultAlgorithmChoice.SupernodalLUFactorization
 
     # Medium-size, very sparse (density < 2e-4) → density branch picks KLU
     n = 9_000
@@ -503,7 +503,7 @@ let
     # Medium-size, dense sparse → UMFPACK
     A_med_dense = sprand(5_000, 5_000, 0.5) + I
     @test LinearSolve.defaultalg(A_med_dense, rand(5_000), LinearSolve.OperatorAssumptions(true)).alg ===
-        LinearSolve.DefaultAlgorithmChoice.UMFPACKFactorization
+        LinearSolve.DefaultAlgorithmChoice.SupernodalLUFactorization
 end
 
 # === Sparse LU → SPQR fallback ===
@@ -559,7 +559,7 @@ let
     A_u[1, :] .= 0
     A_u = sparse(A_u)
     @test LinearSolve.defaultalg(A_u, ones(n_u), LinearSolve.OperatorAssumptions(true)).alg ===
-        LinearSolve.DefaultAlgorithmChoice.UMFPACKFactorization
+        LinearSolve.DefaultAlgorithmChoice.SupernodalLUFactorization
     sol_u = solve(LinearProblem(A_u, ones(n_u)))
     @test sol_u.retcode === ReturnCode.Success
     # Sanity: a non-singular matrix should NOT fall back.