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?

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.

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.

```
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

This needs to be a CUSPARSE array, i.e

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

By the way, `vec`

is a symbol in `Base`

. Trying to reassign it should throw an error.