How to solve Ax=B using cudss, B is a n*m Array

[crazyfireji]

Like the title, I want to solve Ax=B using the code:

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

T = Float64
n = 1000
A_cpu = sprand(T, n, n, 0.05) + I;
x_cpu = zeros(T, n, 2);
b_cpu = rand(T, n, 2);

A_gpu = CuSparseMatrixCSR(A_cpu);
x_gpu =x_cpu |> cu;
b_gpu = b_cpu |> cu;

solver = CudssSolver(A_gpu, "G", 'F')

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

But it doesnt work!

[amontoison]

@crazyfireji
Please provide some logs instead of just saying it doesn't work.

We can't help you with this kind of comment...

[crazyfireji]

@crazyfireji Please provide some logs instead of just saying it doesn't work.

We can't help you with this kind of comment...

the log as following:

julia> cudss("analysis", solver, x_gpu, b_gpu)
ERROR: MethodError: no method matching cudss(::String, ::CudssSolver{Float64}, ::CuArray{Float32, 2, CUDA.DeviceMemory}, ::CuArray{Float32, 2, CUDA.DeviceMemory})
The function `cudss` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  cudss(::String, ::CudssSolver{T}, ::CuArray{T, 2}, ::CuArray{T, 2}) where T<:Union{Float32, Float64, ComplexF64, ComplexF32}
   @ CUDSS ~/.julia/packages/CUDSS/RoKCe/src/interfaces.jl:219
  cudss(::String, ::CudssSolver{T}, ::CuArray{T, 1}, ::CuArray{T, 1}) where T<:Union{Float32, Float64, ComplexF64, ComplexF32}
   @ CUDSS ~/.julia/packages/CUDSS/RoKCe/src/interfaces.jl:213
  cudss(::String, ::CudssSolver{T}, ::CudssMatrix{T}, ::CudssMatrix{T}) where T<:Union{Float32, Float64, ComplexF64, ComplexF32}
   @ CUDSS ~/.julia/packages/CUDSS/RoKCe/src/interfaces.jl:209

Stacktrace:
 [1] top-level scope
   @ REPL[78]:1

[amontoison]

Ok, I know the issue.
You didn't copy-pasted the example of the README and did some local modifications.
The error comes from theses two lines:

x_gpu =x_cpu |> cu;
b_gpu = b_cpu |> cu;

cu changed the eltype of the data.
x_gpu and b_gpu are CuArray of type Float32 and not Float64.

If you use this code, it should work:

x_gpu =CuMatrix(x_cpu)
b_gpu = CuMatrix(b_cpu)

because it will not change the precision of the data.

[crazyfireji]

Ok, I know the issue. You didn't copy-pasted the example of the README and did some local modifications. The error comes from theses two lines:

x_gpu =x_cpu |> cu;
b_gpu = b_cpu |> cu;

cu changed the eltype of the data. x_gpu and b_gpu are CuArray of type Float32 and not Float64.

If you use this code, it should work:

x_gpu =CuMatrix(x_cpu)
b_gpu = CuMatrix(b_cpu)

because it will not change the precision of the data.

Oh,Oh, Thank you for your help!!!!!!!

[crazyfireji]

but how to solve this equetion if A is unsym matrix Float64, b and x are ComplexF64, like:

T = ComplexF64
R = real(T)
n = 50
p = 5

A_cpu = sprand(R, n, n, 0.01)+ I
X_cpu = zeros(T, n, p)
B_cpu = rand(T, n, p)

A_gpu = CuSparseMatrixCSR(A_cpu)
X_gpu = CuMatrix(X_cpu)
B_gpu = CuMatrix(B_cpu)

solver = CudssSolver(A_gpu, "G", 'F')

cudss("analysis", solver, X_gpu, B_gpu)
cudss("factorization", solver, X_gpu, B_gpu)
cudss("solve", solver, X_gpu, B_gpu)

the log as follows:

julia> cudss("analysis", solver, X_gpu, B_gpu)
ERROR: MethodError: no method matching cudss(::String, ::CudssSolver{…}, ::CuArray{…}, ::CuArray{…})
The function `cudss` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  cudss(::String, ::CudssSolver{T}, ::CuArray{T, 2}, ::CuArray{T, 2}) where T<:Union{Float32, Float64, ComplexF64, ComplexF32}
   @ CUDSS ~/.julia/packages/CUDSS/RoKCe/src/interfaces.jl:219
  cudss(::String, ::CudssSolver{T}, ::CuArray{T, 1}, ::CuArray{T, 1}) where T<:Union{Float32, Float64, ComplexF64, ComplexF32}
   @ CUDSS ~/.julia/packages/CUDSS/RoKCe/src/interfaces.jl:213
  cudss(::String, ::CudssSolver{T}, ::CudssMatrix{T}, ::CudssMatrix{T}) where T<:Union{Float32, Float64, ComplexF64, ComplexF32}
   @ CUDSS ~/.julia/packages/CUDSS/RoKCe/src/interfaces.jl:209

Stacktrace:
 [1] top-level scope
   @ REPL[49]:1
Some type information was truncated. Use `show(err)` to see complete types.

[amontoison]

You just need to extract the real and complex parts of your right-hand side b_gpu, solve two systems, and then combine the solutions.

using CUDA, CUDA.CUSPARSE
using CUDSS

T = ComplexF64
R = real(T)
n = 50
p = 5

A_cpu = sprand(R, n, n, 0.01)+ I
X_cpu = zeros(T, n, p)
B_cpu = rand(T, n, p)

A_gpu = CuSparseMatrixCSR(A_cpu)
X_gpu = CuMatrix(X_cpu)
B_gpu = CuMatrix(B_cpu)

C_gpu = real.(B_gpu)
D_gpu = imag.(B_gpu)

Y_gpu = copy(C_gpu)
Z_gpu = copy(D_gpu)

solver = CudssSolver(A_gpu, "G", 'F')

cudss("analysis", solver, Y_gpu, C_gpu)
cudss("factorization", solver, Y_gpu, C_gpu)
cudss("solve", solver, Y_gpu, C_gpu)
cudss("solve", solver, Z_gpu, D_gpu)
X_gpu .= Y_gpu .+ im .* Z_gpu

# Check that we have the correct solution
norm(B_gpu - ComplexF64.(A_gpu) * X_gpu)