# (fast) sparse solve on GPU
**URL:** https://discourse.julialang.org/t/fast-sparse-solve-on-gpu/32297
**Category:** General Usage
**Tags:** question, gpu
**Created:** [December 15, 2019, 9:28am UTC](https://discourse.julialang.org/t/fast-sparse-solve-on-gpu/32297 "2019-12-15T09:28:12Z")
**Posts on this page:** 8
**Page:** 1
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [December 15, 2019, 9:28am UTC](https://discourse.julialang.org/t/fast-sparse-solve-on-gpu/32297/1 "2019-12-15T09:28:12Z")
Hi,
I would like to solve sparse linear systems on the GPU. But when I try the tests [example](https://github.com/JuliaGPU/CuArrays.jl/blob/master/test/sparse_solver.jl) of CuArrays, the GPU version seems way slower than the CPU one. Is it expected?
Additionnally, say I have a LU factorisation of a cpu matrix. Can I transfer this to the GPU so that `\` is overloaded. I may be asking for too much 🤔 but maybe somebody has already done it!
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [December 15, 2019, 9:37am UTC](https://discourse.julialang.org/t/fast-sparse-solve-on-gpu/32297/2 "2019-12-15T09:37:18Z")
> [@rveltz](#):
>
> I would like to solve sparse linear systems on the GPU. But when I try the tests [example](https://github.com/JuliaGPU/CuArrays.jl/blob/master/test/sparse_solver.jl) of CuArrays, the GPU version seems way slower than the CPU one. Is it expected?
GPU sparse solves silently use the CPU for part of it IIRC, kind of like SVD.
> [@rveltz](#):
>
> Can I transfer this to the GPU so that `\` is overloaded. I may be asking for too much 🤔 but maybe somebody has already done it!
You’ll want to upstream our fix in DiffEqBase: [https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/master/src/init.jl#L146-L150](https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/master/src/init.jl#L146-L150)
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [December 15, 2019, 3:03pm UTC](https://discourse.julialang.org/t/fast-sparse-solve-on-gpu/32297/3 "2019-12-15T15:03:20Z")
> [@ChrisRackauckas](#):
>
> You’ll want to upstream our fix in DiffEqBase:
Interesting, I will try that!
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [December 27, 2019, 5:24pm UTC](https://discourse.julialang.org/t/fast-sparse-solve-on-gpu/32297/4 "2019-12-27T17:24:00Z")
Hi,
I tried the above solution but I have a bug I cannot find. It is about using `\` for transpose sparse CuArrays. I perform an iLU decomposition and then try to use the decomposition to solve linear systems on the GPU:
```julia
using SparseArrays, LinearAlgebra, IncompleteLU
n = 1000
A = I + sprand(n,n,0.01)
Precilu = ilu(A, Ï„ = 0.1)
```
Now I check that my formula for the inverse works:
```julia
rhs = rand(n)
sol_0 = Precilu \ rhs
sol_1 = Precilu.U' \ ((I+Precilu.L) \ (rhs))
norm(sol_1-sol_0, Inf64)
```
It returns 0!! Hence, I have perform these solves in the GPU. The first one works well:
```julia
using CuArrays
CuArrays.allowscalar(false)
sol_0 = (I+Precilu.L) \ rhs
sol_1 = LowerTriangular(CuArrays.CUSPARSE.CuSparseMatrixCSR(I+Precilu.L)) \ CuArray(rhs)
norm(sol_0-Array(sol_1), Inf64)
```
and this returns the `4.492221705731936e-9`. However, the following fails and I don’t understand why 🙄
```julia
sol_0 = (Precilu.U)' \ rhs
sol_1 = LowerTriangular(CuArrays.CUSPARSE.CuSparseMatrixCSR(sparse(Precilu.U'))) \ CuArray(rhs)
norm(sol_0-Array(sol_1), Inf64)
```
and this returns `1137.9857325312835`. Maybe the fact that
```julia
rhs = rand(n)
sol_0 = Precilu \ rhs
sol_1 = LowerTriangular(Precilu.U') \ (LowerTriangular(I+Precilu.L) \ (rhs))
norm(sol_1-sol_0, Inf64)
```
does not return zero whereas this one does is key 🤔
```julia
rhs = rand(n)
sol_0 = Precilu \ rhs
sol_1 =(Precilu.U') \ (LowerTriangular(I+Precilu.L) \ (rhs))
norm(sol_1-sol_0, Inf64)
```
---
### Author:  [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)
#### Post date: [December 27, 2019, 6:12pm UTC](https://discourse.julialang.org/t/fast-sparse-solve-on-gpu/32297/5 "2019-12-27T18:12:56Z")
It seems that `Precilu.U` is not `UpperTriangular` but lower. So `LowerTriangular(Precilu.U')` in your code should be `UpperTriangular(Precilu.U')` which makes sense because the second linear system solve should be using the upper triangular matrix from the factorization. I don’t know why @stabbles went with that convention in the package but fixing the above fixes your problem.
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [December 28, 2019, 7:52am UTC](https://discourse.julialang.org/t/fast-sparse-solve-on-gpu/32297/6 "2019-12-28T07:52:30Z")
Wow, this was nicely spotted, thank you!
---
### Author:  [@stabbles](https://discourse.julialang.org/u/stabbles)
#### Post date: [December 28, 2019, 9:18pm UTC](https://discourse.julialang.org/t/fast-sparse-solve-on-gpu/32297/7 "2019-12-28T21:18:15Z")
It’s a technical issue; you get the U factor row by row in crout ILU, but SparseMatrixCSC is created column by column.
---
### Author:  [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)
#### Post date: [December 28, 2019, 9:25pm UTC](https://discourse.julialang.org/t/fast-sparse-solve-on-gpu/32297/8 "2019-12-28T21:25:30Z")
Might be worth returning a `Transpose` then, it’s a bit counter-intuitive that `U` is lower triangular 🙂