Sparse GPU linear solve from documentation fails

Specific Domains Numerics
1

To solve the Poisson equation on the GPU, I would like to use the interface from LinearSolve.jl to do a sparse linear solve. However, the following minimal reproducible example (also provided in the documentation) gives an error:

The versions I use:

Julia: 1.11.5
LinearSolve: 3.17.0
CUDA: 5.8.2
CUDSS: 0.4.4

The code:

using LinearAlgebra, SparseArrays, LinearSolve
using CUDA, CUDA.CUSPARSE, CUDSS

T = Float32
n = 100
A_cpu = sprand(T, n, n, 0.05) + I
x_cpu = zeros(T, n)
b_cpu = rand(T, n)

A_gpu_csr = CuSparseMatrixCSR(A_cpu)
b_gpu = CuVector(b_cpu)

prob = LinearProblem(A_gpu_csr, b_gpu)
sol = solve(prob, LUFactorization())

The error I got:

MethodError: no method matching SparseArrays.UMFPACK.UmfpackLU(::SparseMatrixCSC{Float32, Int32})
The type `SparseArrays.UMFPACK.UmfpackLU` exists, but no method is defined for this combination of argument types when trying to construct it.

Closest candidates are:
  SparseArrays.UMFPACK.UmfpackLU(::SparseArrays.UMFPACK.Symbolic{Tv, Ti}, ::SparseArrays.UMFPACK.Numeric{Tv, Ti}, ::Int64, ::Int64, ::Vector{Ti}, ::Vector{Ti}, ::Vector{Tv}, ::Int64, ::SparseArrays.UMFPACK.UmfpackWS{Ti}, ::Vector{Float64}, ::Vector{Float64}, ::ReentrantLock) where {Tv<:Union{Float64, ComplexF64}, Ti<:Union{Int32, Int64}}
   @ SparseArrays C:\Users\joerg\AppData\Local\Programs\Julia-1.11.5\share\julia\stdlib\v1.11\SparseArrays\src\solvers\umfpack.jl:227
  SparseArrays.UMFPACK.UmfpackLU(::SparseArrays.AbstractSparseMatrixCSC{Tv, Ti}; control) where {Tv<:Union{Float64, ComplexF64}, Ti<:Union{Int32, Int64}}
   @ SparseArrays C:\Users\joerg\AppData\Local\Programs\Julia-1.11.5\share\julia\stdlib\v1.11\SparseArrays\src\solvers\umfpack.jl:241

Stacktrace:
 [1] init_cacheval(alg::LUFactorization{RowMaximum}, A::CuSparseMatrixCSR{Float32, Int32}, b::CuArray{Float32, 1, CUDA.DeviceMemory}, u::CuArray{Float32, 1, CUDA.DeviceMemory}, Pl::IdentityOperator, Pr::IdentityOperator, maxiters::Int64, abstol::Float32, reltol::Float32, verbose::Bool, assumptions::OperatorAssumptions{Bool})
   @ LinearSolveSparseArraysExt C:\Users\joerg\.julia\packages\LinearSolve\gQrff\ext\LinearSolveSparseArraysExt.jl:111
 [2] init(::LinearProblem{Nothing, true, CuSparseMatrixCSR{Float32, Int32}, CuArray{Float32, 1, CUDA.DeviceMemory}, SciMLBase.NullParameters, Nothing, @Kwargs{}}, ::LUFactorization{RowMaximum}; alias::LinearAliasSpecifier, abstol::Float32, reltol::Float32, maxiters::Int64, verbose::Bool, Pl::Nothing, Pr::Nothing, assumptions::OperatorAssumptions{Bool}, sensealg::LinearSolveAdjoint{Missing}, kwargs::@Kwargs{})
   @ LinearSolve C:\Users\joerg\.julia\packages\LinearSolve\gQrff\src\common.jl:239
 [3] init(::LinearProblem{Nothing, true, CuSparseMatrixCSR{Float32, Int32}, CuArray{Float32, 1, CUDA.DeviceMemory}, SciMLBase.NullParameters, Nothing, @Kwargs{}}, ::LUFactorization{RowMaximum})
   @ LinearSolve C:\Users\joerg\.julia\packages\LinearSolve\gQrff\src\common.jl:140
 [4] solve(::LinearProblem{Nothing, true, CuSparseMatrixCSR{Float32, Int32}, CuArray{Float32, 1, CUDA.DeviceMemory}, SciMLBase.NullParameters, Nothing, @Kwargs{}}, ::LUFactorization{RowMaximum}; kwargs::@Kwargs{})
   @ LinearSolve C:\Users\joerg\.julia\packages\LinearSolve\gQrff\src\common.jl:295
 [5] solve(::LinearProblem{Nothing, true, CuSparseMatrixCSR{Float32, Int32}, CuArray{Float32, 1, CUDA.DeviceMemory}, SciMLBase.NullParameters, Nothing, @Kwargs{}}, ::LUFactorization{RowMaximum})
   @ LinearSolve C:\Users\joerg\.julia\packages\LinearSolve\gQrff\src\common.jl:293
 [6] top-level scope
   @ REPL[11]:1

However, using cuDSS directly works normally:

x_gpu = CuVector(x_cpu)
solver = CudssSolver(A_gpu_csr, "G", 'F')

cudss("analysis", solver, x_gpu, b_gpu)
cudss("factorization", solver, x_gpu, b_gpu)
cudss("solve", solver, x_gpu, b_gpu)

Am I missing an API change or is this an (documentation) issue?

I would be happy to receive an answer, thank you in advance!

1 Like