# Cuda cg + ilu?

**URL:** <https://discourse.julialang.org/t/cuda-cg-ilu/53069>\
**Category:** GPU\
**Created:** [January 8, 2021, 10:41pm UTC](https://discourse.julialang.org/t/cuda-cg-ilu/53069 "2021-01-08T22:41:24Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![Omega-xyZac](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/omega-xyzac/32/47908_2.png) [@Omega-xyZac](https://discourse.julialang.org/u/Omega-xyZac)\
**Post date:** [January 8, 2021, 10:41pm UTC](https://discourse.julialang.org/t/cuda-cg-ilu/53069/1 "2021-01-08T22:41:24Z")

</div>

Hi all,  
I’m looking to use a GPU CG method with an ILU preconditioner. I’m currently using a CPU CG method through IterativeSolvers and ILU through **[IncompleteLU.jl](https://github.com/haampie/IncompleteLU.jl)**. The CG method from IterativeSolvers is easily extended to GPUs however, some work is required on IncompleteLU.jl to implement it on a GPU. I think only `forward_substitution` and `backward_substitution` methods are required to extend this. I’ve tried using `sv2!` and `ilu02` with little success. Does anyone have any ideas?

---

<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:** [January 8, 2021, 11:46pm UTC](https://discourse.julialang.org/t/cuda-cg-ilu/53069/2 "2021-01-08T23:46:53Z")

</div>

do you need the ILU on GPU? You could just pass the decomposition to the GPU, see [https://rveltz.github.io/BifurcationKit.jl/dev/tutorialsCGL/#Complex-Ginzburg-Landau-2d-1](https://rveltz.github.io/BifurcationKit.jl/dev/tutorialsCGL/#Complex-Ginzburg-Landau-2d-1)

---

<div class="post-metadata">

**Author:** ![Omega-xyZac](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/omega-xyzac/32/47908_2.png) [@Omega-xyZac](https://discourse.julialang.org/u/Omega-xyZac)\
**Post date:** [January 8, 2021, 11:57pm UTC](https://discourse.julialang.org/t/cuda-cg-ilu/53069/3 "2021-01-08T23:57:11Z")

</div>

The issue is that the preconditioner requires the ‘ldiv!’ method as per [https://julialinearalgebra.github.io/IterativeSolvers.jl/dev/preconditioning/](https://julialinearalgebra.github.io/IterativeSolvers.jl/dev/preconditioning/). Unless I am mistaken?

---

<div class="post-metadata">

**Author:** ![Omega-xyZac](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/omega-xyzac/32/47908_2.png) [@Omega-xyZac](https://discourse.julialang.org/u/Omega-xyZac)\
**Post date:** [January 8, 2021, 11:58pm UTC](https://discourse.julialang.org/t/cuda-cg-ilu/53069/4 "2021-01-08T23:58:53Z")

</div>

The link you sent does look helpful though. Thank you for that

---

<div class="post-metadata">

**Author:** ![Omega-xyZac](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/omega-xyzac/32/47908_2.png) [@Omega-xyZac](https://discourse.julialang.org/u/Omega-xyZac)\
**Post date:** [January 9, 2021, 12:38am UTC](https://discourse.julialang.org/t/cuda-cg-ilu/53069/5 "2021-01-09T00:38:24Z")

</div>

Here is some code:

```julia
using CUDA, CUDA.CUSPARSE, CUDA.CUSOLVER
using LinearAlgebra
using SparseArrays
using IncompleteLU
using IterativeSolvers
CUDA.allowscalar(false)

val = sprand(200,200,0.05);
A_cpu = val*val'
b_cpu = rand(200)

A_gpu = CuSparseMatrixCSC(A_cpu)
b_gpu = CuArray(b_cpu)

Precilu = ilu(A_cpu, τ = 3.0)

import Base: ldiv!
function LinearAlgebra.ldiv!(y::CuArray, P::LUgpu, x::CuArray)
  copyto!(y, x)                       
  sv2!('N', 'L', 1.0, P, y, 'O') 
  sv2!('N', 'U', 1.0, P, y, 'O')  
  return y
end

function LinearAlgebra.ldiv!(P::LUgpu, x::CuArray)                        
  sv2!('N', 'L', 1.0, P, x, 'O') 
  sv2!('N', 'U', 1.0, P, x, 'O')  
  return x
end

# function LinearAlgebra.ldiv!(_lu::LUperso, rhs::CUDA.CuArray)
# _x = UpperTriangular(_lu.Ut) \ (LowerTriangular(_lu.L) \ rhs)
# rhs .= vec(_x)
# # CUDA.unsafe_free!(_x)
# rhs
# end
#
# function LinearAlgebra.ldiv!(yrhs::CuArray,_lu::LUperso, rhs::CuArray)
# copyto!(yrhs,rhs)
# _x = UpperTriangular(_lu.Ut) \ (LowerTriangular(_lu.L) \ rhs)
# rhs .= vec(_x)
# # CUDA.unsafe_free!(_x)
# rhs
# end

struct LUgpu
	L
	Ut	# transpose of U in LU decomposition
end

P = LUperso(LowerTriangular(CuSparseMatrixCSR(I+Precilu.L)), UpperTriangular(CuSparseMatrixCSR(sparse(Precilu.U'))));
val = cg(A_gpu,b_gpu,verbose=true,Pl=P,tol=10^-7,maxiter=1000)

P_cpu = ilu(A_cpu,τ=3.0)
val = cg(A_cpu,b_cpu;verbose=true,Pl=P_cpu,tol=10^-7,maxiter=1000)
A_cpu\b_cpu

```