Add support of `CuSparseMatrix`-type matrices for `QRFactorization()`

[ytdHuang]

I tried to solve some sparse matrix with the type of CuSparseMatrixCSC. When I switch to the algorithm QRFactorization(), it seems to automatically convert the sparse matrix to dense matrix CuArray.

I'm using LinearSolve v2.12.1

[ChrisRackauckas]

Is this the behavior of qr(A) in CUDA.jl itself?

[ytdHuang]

I'm using CUDA v5.0.0

Therefore, if I want to solve large sparse problems using GPU (CUDA), the Krylov-based algorithms are the only choices so far ?

[ytdHuang]

By the way, I also want to report that there might be a small issue in the Krylov algorithms. But I'm not sure whether the issue is from LinearSolve or it's also the behavior of Krylov.jl itself.

When the sparse matrix is given as the type of CuSparseMatrixCSC{ComplexF32, Int32}, we will get an TypeError:

using LinearSolve
using CUDA
using CUDA.CUSPARSE
A::CuSparseMatrixCSC{ComplexF32, Int32}
b::CuArray{ComplexF32}
cache = init(LinearProblem(A, b, solver=KrylovJL_GMRES())
sol = solve!(cache)

TypeError: in keyword argument atol, expected Float32, got a value of type ComplexF32

Of course we can provide the atol and rtol, but we have to convert it to Float32 by ourselves:

LinearProblem(A, b, 
    solver=KrylovJL_GMRES(rtol=1f-6, atol=1f-8)
)

Otherwise we will still get an error

LinearProblem(A, b, 
    solver=KrylovJL_GMRES(rtol=1e-6, atol=1e-8)
)

TypeError: in keyword argument atol, expected Float32, got a value of type Float64

[ChrisRackauckas]

I'm asking about something requires no LinearSolve at all. If you run qr(A) on a CuSparseMatrixCSC, do it give you a dense factorization?

[ytdHuang]

Yeah, I think so.

using SparseArrays
using LinearAlgebra
using CUDA
using CUDA.CUSPARSE

A = sprand(100, 100, 0.1)
colptr = CuArray{Int32}(A.colptr)
rowval = CuArray{Int32}(A.rowval)
nzval  = CuArray{ComplexF32}(A.nzval)
A_gpu  = CuSparseMatrixCSC{ComplexF32, Int32}(colptr, rowval, nzval, size(A))

qr(A_gpu)

output:

LinearAlgebra.QRCompactWY{ComplexF32, Matrix{ComplexF32}, Matrix{ComplexF32}}
Q factor:
bla bla bla ...

the output is something in the type of Matrix{ComplexF32}

[ChrisRackauckas]

Okay yeah so it's an upstream issue with CUDA.jl. File an issue about getting that handled and it should automatically fix the LinearSolve interface here.

[ytdHuang]

I think this also relates to this issue https://github.com/JuliaGPU/CUDA.jl/issues/1396

[ytdHuang]

Maybe this issue is solved in here https://github.com/JuliaGPU/CUDA.jl/pull/2121

[amontoison]

For the issue with Krylov.jl, it was solved with https://github.com/SciML/LinearSolve.jl/pull/397.

[ChrisRackauckas]

Yeah that piece was fixed earlier last week so I'm not worried there. The CUDA part is the bigger issue though.