# Use of linear operators on a GPU

**URL:** <https://discourse.julialang.org/t/use-of-linear-operators-on-a-gpu/92899>\
**Category:** GPU\
**Tags:** question\
**Created:** [January 12, 2023, 9:25pm UTC](https://discourse.julialang.org/t/use-of-linear-operators-on-a-gpu/92899 "2023-01-12T21:25:12Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![shakedregev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shakedregev/32/45821_2.png) [@shakedregev](https://discourse.julialang.org/u/shakedregev)\
**Post date:** [January 12, 2023, 9:25pm UTC](https://discourse.julialang.org/t/use-of-linear-operators-on-a-gpu/92899/1 "2023-01-12T21:25:12Z")

</div>

Say I have a sparse matrices `$A$` and `$B=D*A^T$` and their product is dense. I want to solve systems with them using an iterative method from Krylov.jl. On a CPU, I can do  
`opC =LinearOperator(A)*LinearOperator(B)`  
`cg(opC,r)`  
This works fine, but how do I do this with CUDA?

---

<div class="post-metadata">

**Author:** ![geoffroyleconte](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/geoffroyleconte/32/21661_2.png) [@geoffroyleconte](https://discourse.julialang.org/u/geoffroyleconte)\
**Post date:** [January 15, 2023, 1:08am UTC](https://discourse.julialang.org/t/use-of-linear-operators-on-a-gpu/92899/2 "2023-01-15T01:08:27Z")

</div>

Hi! If A and B are CUDA sparse matrices in Float64 (also works in Float32) you can use  
`opC = LinearOperator(A, S = CuArray{Float64, 1, CUDA.Mem.DeviceBuffer}) * LinearOperator(B, S = CuArray{Float64, 1, CUDA.Mem.DeviceBuffer})`.  
If `r` is a CUDA vector this should work.

---

<div class="post-metadata">

**Author:** ![shakedregev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shakedregev/32/45821_2.png) [@shakedregev](https://discourse.julialang.org/u/shakedregev)\
**Post date:** [January 18, 2023, 12:01am UTC](https://discourse.julialang.org/t/use-of-linear-operators-on-a-gpu/92899/3 "2023-01-18T00:01:12Z")

</div>

Hi and thanks for answering! This code I posted is pretty ridiculous, but it’s just to show the idea. The error follows the code block.

```julia
using Pkg
using NPZ
using SparseArrays
using Krylov
using LinearAlgebra
using LinearOperators
using LDLFactorizations
import LinearAlgebra.ldiv!
using KrylovKit
using CUDA
using CUDA.CUSPARSE
using CUDAKernels
vec::Vector{Float64} = collect(LinRange(1,2,100))
rows::Vector{Int64} = collect(1:length(vec))
cols::Vector{Int64} = collect(1:length(vec))
Dvec = sparse(rows, cols, vec)
opC = LinearOperator(Dvec, S = CuArray{Float64, 1, CUDA.Mem.DeviceBuffer}) * LinearOperator(Dvec, S = CuArray{Float64, 1, CUDA.Mem.DeviceBuffer})
(x,stats) = cg(opC,vec)
println(stats.niter)

```

LoadError: Scalar indexing is disallowed.  
Invocation of getindex resulted in scalar indexing of a GPU array.  
This is typically caused by calling an iterating implementation of a method.  
Such implementations _do not_ execute on the GPU, but very slowly on the CPU

---

<div class="post-metadata">

**Author:** ![skleinbo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/skleinbo/32/36080_2.png) [@skleinbo](https://discourse.julialang.org/u/skleinbo)\
**Post date:** [January 18, 2023, 7:12am UTC](https://discourse.julialang.org/t/use-of-linear-operators-on-a-gpu/92899/4 "2023-01-18T07:12:19Z")

</div>

> [@shakedregev](#):
>
> `Dvec = sparse(rows, cols, vec)`

This needs to be a CUSPARSE array, i.e

```julia
Dvec = CuSparseMatrixCSC(sparse(rows, cols, vec))

```

By the way, `vec` is a symbol in `Base`. Trying to reassign it should throw an error.