Skip to content

Commit fd8271c

Browse files
ChrisRackauckas-ClaudeChrisRackauckasclaude
authored
Bound the UMFPACK-backed sparse cachevals to the eltypes UMFPACK has (#1122)
`init(LinearProblem(A, b), nothing)` throws for a Float32 sparse matrix on main: defaultalg(...) = DefaultLinearSolver(DefaultAlgorithmChoice.KLUFactorization, ...) init(...) = MethodError: no method matching SparseArrays.UMFPACK.UmfpackLU(::SparseMatrixCSC{Float32, Int64}) Note the algorithm that throws is not the one selected. `_init_default_cacheval` (src/default.jl, @generated) emits an unconditional `init_cacheval` for every slot - only the three Krylov slots are guarded - so a too-permissive `init_cacheval` on an *unselected* algorithm still breaks `init`. KLU is chosen; the `LUFactorization` slot is built first and throws. `LUFactorization`'s and `UMFPACKFactorization`'s sparse cachevals are bounded by `BLASELTYPES = Union{Float32, Float64, ComplexF32, ComplexF64}`, but UMFPACK's `UmfpackLU` exists only for Float64/ComplexF64 (SuiteSparse's `UMFVTypes`), so the bound is too wide by exactly Float32/ComplexF32. KLU already does this correctly with `T <: KLU.KLUTypes`. Bound them to a local `UMFPACKELTYPES` instead. The additional `BLASELTYPES` method is not redundant: the old methods' `is_cusparse(A)` branch is the only `LUFactorization`+CuSparse `init_cacheval` in the tree, so narrowing alone would have silently dropped Float32 cuDSS support. `UMFPACKELTYPES` is defined locally rather than reusing `SparseArrays.UMFPACK.UMFVTypes`, which is a non-public internal. Float32/ComplexF32 sparse now resolve to the `KLUFactorization` slot, which `algchoice_to_alg` maps to `PureKLUFactorization()` - pure Julia, and the only thing in the tree that supports these eltypes. It is a hard dependency, so no extra `using` is needed. Two separate breakages, both bisected: 0ce03c2 2025-05-22 widened the cachevals to `T <: BLASELTYPES` while adding ComplexF64. Broke Float32 sparse for Sparspak users; everyone else still hit the "using Sparspak required" error in `defaultalg` first, so it stayed hidden. 52ebf25 (#1037) 2026-06-13 removed that early error, exposing the breakage to everyone. Measured before/after, all 10 eltype/index combinations: before: Float32 and ComplexF32 (Int64 and Int32) all `INIT THREW` after: all 10 solve, residual 0, eltype preserved (no Float64 promotion) Explicitly-requested `LUFactorization`/`UMFPACKFactorization` on Float32 sparse remain unsupported, as before this change - they now fail the way SuiteSparse KLU already does for an unsupported eltype rather than with a `MethodError`. Making them work would mean silent Float32->Float64 promotion, which is a behavior addition, not a bug fix. `Test.detect_ambiguities` on the extension: 12 before, 12 after. GROUP=Core passes on the fixed tree and on unmodified main. Claude-Session: https://claude.ai/code/session_017rr42T1vFRRT8D7DMAmJzf Co-authored-by: ChrisRackauckas-Claude <accounts@chrisrackauckas.com> Co-authored-by: Claude Fable 5 <noreply@anthropic.com>
1 parent 2fcfb43 commit fd8271c

2 files changed

Lines changed: 43 additions & 4 deletions

File tree

ext/LinearSolveSparseArraysExt.jl

Lines changed: 24 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -128,6 +128,11 @@ end
128128
)
129129
end # @static if Base.USE_GPL_LIBS
130130

131+
# UMFPACK's `UmfpackLU` only exists for `Float64`/`ComplexF64` (SuiteSparse's
132+
# `UMFVTypes`). The other BLAS eltypes are not a subset of it, so cachevals for
133+
# UMFPACK-backed algorithms must be bounded by this, not by `BLASELTYPES`.
134+
const UMFPACKELTYPES = Union{Float64, ComplexF64}
135+
131136
function LinearSolve.init_cacheval(
132137
alg::LUFactorization, A::AbstractSparseArray{<:Number, <:Integer}, b, u,
133138
Pl, Pr,
@@ -169,7 +174,7 @@ end
169174
Pl, Pr,
170175
maxiters::Int, abstol, reltol,
171176
verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions
172-
) where {T <: BLASELTYPES}
177+
) where {T <: UMFPACKELTYPES}
173178
if LinearSolve.is_cusparse(A)
174179
LinearSolve.cudss_loaded(A) ? ArrayInterface.lu_instance(A) : nothing
175180
else
@@ -185,7 +190,7 @@ end
185190
Pl, Pr,
186191
maxiters::Int, abstol, reltol,
187192
verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions
188-
) where {T <: BLASELTYPES}
193+
) where {T <: UMFPACKELTYPES}
189194
if LinearSolve.is_cusparse(A)
190195
LinearSolve.cudss_loaded(A) ? ArrayInterface.lu_instance(A) : nothing
191196
else
@@ -198,6 +203,21 @@ end
198203
end
199204
end # @static if Base.USE_GPL_LIBS
200205

206+
# Single-precision sparse: no UMFPACK cacheval exists, but a cuDSS-backed
207+
# CuSparse matrix still needs its `lu_instance`.
208+
function LinearSolve.init_cacheval(
209+
alg::LUFactorization, A::AbstractSparseArray{T, <:Union{Int32, Int64}}, b, u,
210+
Pl, Pr,
211+
maxiters::Int, abstol, reltol,
212+
verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions
213+
) where {T <: BLASELTYPES}
214+
if LinearSolve.is_cusparse(A)
215+
return LinearSolve.cudss_loaded(A) ? ArrayInterface.lu_instance(A) : nothing
216+
else
217+
return nothing
218+
end
219+
end
220+
201221
function LinearSolve.init_cacheval(
202222
alg::LUFactorization, A::LinearSolve.GPUArraysCore.AnyGPUArray, b, u,
203223
Pl, Pr,
@@ -231,7 +251,7 @@ end
231251
Pl, Pr,
232252
maxiters::Int, abstol, reltol,
233253
verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions
234-
) where {T <: BLASELTYPES}
254+
) where {T <: UMFPACKELTYPES}
235255
SparseArrays.UMFPACK.UmfpackLU(
236256
SparseMatrixCSC{T, Int64}(
237257
zero(Int64), zero(Int64), [Int64(1)], Int64[], T[]
@@ -244,7 +264,7 @@ end
244264
Pl, Pr,
245265
maxiters::Int, abstol, reltol,
246266
verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions
247-
) where {T <: BLASELTYPES}
267+
) where {T <: UMFPACKELTYPES}
248268
SparseArrays.UMFPACK.UmfpackLU(
249269
SparseMatrixCSC{T, Int32}(
250270
zero(Int32), zero(Int32), [Int32(1)], Int32[], T[]

test/Core/default_algs.jl

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -115,6 +115,25 @@ let As_bf = sparse(BigFloat.([1 2 3; 2 4 6; 1 1 1])), bs_bf = BigFloat.([1, 2, 3
115115
@test all(isfinite, sol.u)
116116
end
117117

118+
# Single-precision sparse has no UMFPACK/KLU (SuiteSparse) support, so the
119+
# default polyalgorithm — which eagerly builds a cacheval for *every* slot, not
120+
# just the selected one — must not try to allocate an UMFPACK cacheval for it.
121+
@testset "Sparse $T with $Ti indices" for T in (Float32, ComplexF32),
122+
Ti in (Int64, Int32)
123+
124+
n = 20
125+
A32 = sparse(Ti.(1:n), Ti.(1:n), fill(T(n), n), n, n) +
126+
sparse(Ti.(1:(n - 1)), Ti.(2:n), fill(T(1), n - 1), n, n)
127+
b32 = ones(T, n)
128+
@test A32 isa SparseMatrixCSC{T, Ti}
129+
@test LinearSolve.defaultalg(A32, b32, LinearSolve.OperatorAssumptions(true)).alg ===
130+
LinearSolve.DefaultAlgorithmChoice.KLUFactorization
131+
sol32 = solve(LinearProblem(A32, b32))
132+
@test SciMLBase.successful_retcode(sol32.retcode)
133+
@test eltype(sol32.u) === T
134+
@test norm(A32 * sol32.u - b32) / norm(b32) < 1.0f-5
135+
end
136+
118137
@test LinearSolve.defaultalg(sprand(10^4, 10^4, 1.0e-5) + I, zeros(1000)).alg ===
119138
LinearSolve.DefaultAlgorithmChoice.KLUFactorization
120139
prob = LinearProblem(sprand(1000, 1000, 0.5), zeros(1000))

0 commit comments

Comments
 (0)