diff --git a/docs/src/solvers/solvers.md b/docs/src/solvers/solvers.md index 5e790a323..cb2dde20d 100644 --- a/docs/src/solvers/solvers.md +++ b/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 diff --git a/ext/LinearSolveSparseArraysExt.jl b/ext/LinearSolveSparseArraysExt.jl index 10cbfdd4f..2264fcfa1 100644 --- a/ext/LinearSolveSparseArraysExt.jl +++ b/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( diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl index 1275e40ae..71a60bd2c 100644 --- a/src/LinearSolve.jl +++ b/src/LinearSolve.jl @@ -310,7 +310,7 @@ EnumX.@enumx DefaultAlgorithmChoice begin DirectLdiv! SparspakFactorization KLUFactorization - UMFPACKFactorization + SupernodalLUFactorization KrylovJL_GMRES GenericLUFactorization RFLUFactorization diff --git a/src/default.jl b/src/default.jl index 1ecbcc0d5..4c55b9ba9 100644 --- a/src/default.jl +++ b/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, diff --git a/test/Core/default_algs.jl b/test/Core/default_algs.jl index 5c741ffd6..73c2085b4 100644 --- a/test/Core/default_algs.jl +++ b/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.