Use of linear operators on a GPU

Specific Domains GPU
1

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?

1 Like
2

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.

3

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

4

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.

2 Likes