# Cannot solve Ax=B on GPU with A CuSparseMatrixCSC

**URL:** <https://discourse.julialang.org/t/cannot-solve-ax-b-on-gpu-with-a-cusparsematrixcsc/60966>\
**Category:** New to Julia\
**Tags:** cuda, sparse\
**Created:** [May 11, 2021, 4:03pm UTC](https://discourse.julialang.org/t/cannot-solve-ax-b-on-gpu-with-a-cusparsematrixcsc/60966 "2021-05-11T16:03:06Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ludovic\_Dumoulin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ludovic_dumoulin/32/22414_2.png) [@Ludovic\_Dumoulin](https://discourse.julialang.org/u/Ludovic_Dumoulin)\
**Post date:** [May 11, 2021, 4:03pm UTC](https://discourse.julialang.org/t/cannot-solve-ax-b-on-gpu-with-a-cusparsematrixcsc/60966/1 "2021-05-11T16:03:06Z")

</div>

Hello,

I am trying to solve a `A*x=B` matrix equation with `A` being the discretization of a differential operator.  
`A` is a `CuSparseMatrixCSC{Float64}`, `B` and `x` are `CuArray{Float64,1}`.  
I write the following line

> x .= A \ B

and get the following warning:

> Warning: Performing scalar operations on GPU arrays: This is very slow, consider disallowing these operations with `allowscalar(false)`

, and error:

> MethodError: no method matching lu!(::CuSparseMatrixCSC{Float64}, ::Val{true}; check=true)  
> Closest candidates are:  
> lu!(!Matched::StridedArray{T, 2}, ::Union{Val{true}, Val{false}}; check) where T\<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64} at D:\buildbot\worker\package\_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\lu.jl:79

> Stacktrace:  
> [1] lu(::CuSparseMatrixCSC{Float64}, ::Val{true}; check::Bool) at D:\buildbot\worker\package\_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\lu.jl:273  
> [2] lu(::CuSparseMatrixCSC{Float64}, ::Val{true}) at D:\buildbot\worker\package\_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\lu.jl:272 (repeats 2 times)  
> [3] (::CuSparseMatrixCSC{Float64}, ::CuArray{Float64,1}) at D:\buildbot\worker\package\_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\generic.jl:1116

For this test `A` is a small matrix, but once it works `A` will be a `5324800×5324800 CuSparseMatrixCSC{Float64} with 47806720 stored entries:`. Also, I have to do this operation at each time step (500’000) because `B=B(t)`.

Thank you in advance for your help,

Ludo

---

<div class="post-metadata">

**Author:** ![Ludovic\_Dumoulin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ludovic_dumoulin/32/22414_2.png) [@Ludovic\_Dumoulin](https://discourse.julialang.org/u/Ludovic_Dumoulin)\
**Post date:** [May 11, 2021, 4:04pm UTC](https://discourse.julialang.org/t/cannot-solve-ax-b-on-gpu-with-a-cusparsematrixcsc/60966/2 "2021-05-11T16:04:38Z")

</div>

```julia
Julia Version 1.5.2
Commit 539f3ce943 (2020-09-23 23:17 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-7700 CPU @ 3.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-9.0.1 (ORCJIT, skylake)

```

---

<div class="post-metadata">

**Author:** ![roble](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/roble/32/12303_2.png) [@roble](https://discourse.julialang.org/u/roble)\
**Post date:** [May 11, 2021, 10:32pm UTC](https://discourse.julialang.org/t/cannot-solve-ax-b-on-gpu-with-a-cusparsematrixcsc/60966/3 "2021-05-11T22:32:32Z")

</div>

At first glance, I would say the `.=` is wrong. The dot operator is used for broadcasting which is not necessary to save the resulting vector in `x`.

Unfortunately, I don’t know much about Cuda. So I can’t help you there.

---

<div class="post-metadata">

**Author:** ![zhiwei](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zhiwei/32/12385_2.png) [@zhiwei](https://discourse.julialang.org/u/zhiwei)\
**Post date:** [May 11, 2021, 11:06pm UTC](https://discourse.julialang.org/t/cannot-solve-ax-b-on-gpu-with-a-cusparsematrixcsc/60966/4 "2021-05-11T23:06:23Z")

</div>

Try using CSR matrix and csrlsvqr!. Here’s an example:

```julia
using SparseArrays
using LinearAlgebra
using CUDA 

n = 100
a = sprand(n,n,0.2) + sparse(I, n,n)
A = CUSPARSE.CuSparseMatrixCSR(a)
b = CUDA.rand(Float64, n)
x = CuVector{Float64}(undef,n)

tol = 1e-8
CUSOLVER.csrlsvqr!(A, b, x, tol, one(Cint),'O')

```

See [https://github.com/JuliaGPU/CUDA.jl/blob/master/test/cusolver/sparse.jl](https://github.com/JuliaGPU/CUDA.jl/blob/master/test/cusolver/sparse.jl)

---

<div class="post-metadata">

**Author:** ![Ludovic\_Dumoulin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ludovic_dumoulin/32/22414_2.png) [@Ludovic\_Dumoulin](https://discourse.julialang.org/u/Ludovic_Dumoulin)\
**Post date:** [May 12, 2021, 9:09am UTC](https://discourse.julialang.org/t/cannot-solve-ax-b-on-gpu-with-a-cusparsematrixcsc/60966/5 "2021-05-12T09:09:35Z")

</div>

Thank you,

I tried to apply your advice

> [@zhiwei](#):
>
> `CUSOLVER.csrlsvqr!(A, b, x, tol, one(Cint),'O')`

and I got

> LinearAlgebra.SingularException(29)

I just realized that my matrix `A` is not invertible 😢

---

<div class="post-metadata">

**Author:** ![Ludovic\_Dumoulin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ludovic_dumoulin/32/22414_2.png) [@Ludovic\_Dumoulin](https://discourse.julialang.org/u/Ludovic_Dumoulin)\
**Post date:** [May 14, 2021, 12:56pm UTC](https://discourse.julialang.org/t/cannot-solve-ax-b-on-gpu-with-a-cusparsematrixcsc/60966/6 "2021-05-14T12:56:26Z")

</div>

After some changes in `A` it works well.  
However it is very slow compare to a naive iterative scheme.

As I have to use `CUSOLVER.csrlsvqr!(A, b, x, tol, one(Cint),'O')` many times (500’000) I would like to know if I can first compute the QR (or LU) decomposition and then use the solver ?  
Or maybe there is a way to compute only once the decomposition of `A` and to keep it in memory?

---

<div class="post-metadata">

**Author:** ![zhiwei](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zhiwei/32/12385_2.png) [@zhiwei](https://discourse.julialang.org/u/zhiwei)\
**Post date:** [May 19, 2021, 1:15am UTC](https://discourse.julialang.org/t/cannot-solve-ax-b-on-gpu-with-a-cusparsematrixcsc/60966/7 "2021-05-19T01:15:29Z")

</div>

I’m not sure how to do sparse factorization on the GPU directly, maybe someone else could help?

If you’d like to use iterative solvers, there are a few packages that implement Krylov-type methods. Also, maybe you could do the factorization on the CPU and copy arrays to GPU, I haven’t tried this so might not work.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [May 19, 2021, 8:27am UTC](https://discourse.julialang.org/t/cannot-solve-ax-b-on-gpu-with-a-cusparsematrixcsc/60966/8 "2021-05-19T08:27:58Z")

</div>

This is what I do. You have an example here:  
[https://rveltz.github.io/BifurcationKit.jl/dev/tutorialsCGL/#Continuation-of-periodic-orbits-on-the-GPU-(Advanced)-1](https://rveltz.github.io/BifurcationKit.jl/dev/tutorialsCGL/#Continuation-of-periodic-orbits-on-the-GPU-(Advanced)-1)

---

<div class="post-metadata">

**Author:** ![Ludovic\_Dumoulin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ludovic_dumoulin/32/22414_2.png) [@Ludovic\_Dumoulin](https://discourse.julialang.org/u/Ludovic_Dumoulin)\
**Post date:** [May 20, 2021, 8:22am UTC](https://discourse.julialang.org/t/cannot-solve-ax-b-on-gpu-with-a-cusparsematrixcsc/60966/9 "2021-05-20T08:22:22Z")

</div>

Thank you,  
So I tried this method, it is working well !

However, I think my matrix is to large to be solved by direct method. There are 80 times more elements in my `L` matrix than in `A`, and approximately the same ratio in `U`. (I’m using an incomplete LU decomposition, with 1e-7 tolerance. Then I have the same precision than a naive Jacobi iterative scheme. As my system is chaotic, I was looking for a better precision.

---

<div class="post-metadata">

**Author:** ![Ludovic\_Dumoulin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ludovic_dumoulin/32/22414_2.png) [@Ludovic\_Dumoulin](https://discourse.julialang.org/u/Ludovic_Dumoulin)\
**Post date:** [May 20, 2021, 8:22am UTC](https://discourse.julialang.org/t/cannot-solve-ax-b-on-gpu-with-a-cusparsematrixcsc/60966/10 "2021-05-20T08:22:47Z")

</div>

Thank you, I will have a look.