Issue with using IterativeSolvers, LinearMaps and CuArrays

General Usage
1

Dear all,

I am trying to solve big linear systems with iterative methods on the GPU using IterativeSolvers.jl and CuArrays.jl. But I can’t get it to work. My problem is very related to this one and I put my error on this link. My code is a bit involved, so I tried to find a simpler example of bug. This one already produces a bug although it is not completely related to link.

using CuArrays,IterativeSolvers, LinearMaps
A = LinearMap{Float32}(cumsum, 100; ismutating=false)
b = rand(100) |> cu
tol = 1e-5
x = gmres(A, b; tol=tol, maxiter=2000)

returns

ERROR: conversion to pointer not defined for CuArray{Float32,1}
Stacktrace:
 [1] copy!(::SubArray{Float32,1,Array{Float32,2},Tuple{Base.Slice{Base.OneTo{Int64}},Int64},true}, ::CuArray{Float32,1}) at /user/rveltz/home/.julia/v0.6/GPUArrays/src/abstractarray.jl:120
 [2] #init!#55(::Bool, ::Function, ::IterativeSolvers.ArnoldiDecomp{Float32,LinearMaps.FunctionMap{Float32,Base.#cumsum,Void}}, ::Array{Float32,1}, ::CuArray{Float32,1}, ::IterativeSolvers.Identity, ::Array{Float32,1}) at /user/rveltz/home/.julia/v0.6/IterativeSolvers/src/gmres.jl:219
 [3] (::IterativeSolvers.#kw##init!)(::Array{Any,1}, ::IterativeSolvers.#init!, ::IterativeSolvers.ArnoldiDecomp{Float32,LinearMaps.FunctionMap{Float32,Base.#cumsum,Void}}, ::Array{Float32,1}, ::CuArray{Float32,1}, ::IterativeSolvers.Identity, ::Array{Float32,1}) at ./<missing>:0
 [4] #gmres_iterable!#52(::IterativeSolvers.Identity, ::IterativeSolvers.Identity, ::Float64, ::Int64, ::Int64, ::Bool, ::Function, ::Array{Float32,1}, ::LinearMaps.FunctionMap{Float32,Base.#cumsum,Void}, ::CuArray{Float32,1}) at /user/rveltz/home/.julia/v0.6/IterativeSolvers/src/gmres.jl:118
 [5] (::IterativeSolvers.#kw##gmres_iterable!)(::Array{Any,1}, ::IterativeSolvers.#gmres_iterable!, ::Array{Float32,1}, ::LinearMaps.FunctionMap{Float32,Base.#cumsum,Void}, ::CuArray{Float32,1}) at ./<missing>:0
 [6] #gmres!#54(::IterativeSolvers.Identity, ::IterativeSolvers.Identity, ::Float64, ::Int64, ::Int64, ::Bool, ::Bool, ::Bool, ::IterativeSolvers.#gmres!, ::Array{Float32,1}, ::LinearMaps.FunctionMap{Float32,Base.#cumsum,Void}, ::CuArray{Float32,1}) at /user/rveltz/home/.julia/v0.6/IterativeSolvers/src/gmres.jl:185
 [7] (::IterativeSolvers.#kw##gmres!)(::Array{Any,1}, ::IterativeSolvers.#gmres!, ::Array{Float32,1}, ::LinearMaps.FunctionMap{Float32,Base.#cumsum,Void}, ::CuArray{Float32,1}) at ./<missing>:0
 [8] #gmres#53(::Array{Any,1}, ::Function, ::LinearMaps.FunctionMap{Float32,Base.#cumsum,Void}, ::CuArray{Float32,1}) at /user/rveltz/home/.julia/v0.6/IterativeSolvers/src/gmres.jl:134
 [9] (::IterativeSolvers.#kw##gmres)(::Array{Any,1}, ::IterativeSolvers.#gmres, ::LinearMaps.FunctionMap{Float32,Base.#cumsum,Void}, ::CuArray{Float32,1}) at ./<missing>:100:

I hope that this is the reason of the bug.

Thank you for your help

2

IterativeSolvers.jl doesn’t support arbitrary array types (yet) AFAIK.

3

but link seems to imply that it used to work. @samuelpowell confirmed this to me.

4

I saw IterativeSolvers.jl/gmres.jl at 30062518106d28f7a52e35637e68f273e41fcee5 · JuliaLinearAlgebra/IterativeSolvers.jl · GitHub and assumed it doesn’t. Use of Matrix and zeros(n) is also common all over IterativeSolvers.jl. I would be happy to see it working but it may need a bit of work. @stabbles may have a comment on this.

5

The PR to which you have linked is for the CUDAdrv.jl implementation of CuArray which I believe to be different to that provided by CuArrays.jl - see comments by @maleadt in said PR.

Further, I can only confirm that this worked for my particular use case. As I recall this involved a linear operator exposed through the AbstractArray interface which undertook its own memory transfers and launched kernels on the GPU. The only missing piece for me was the copy! function for the views provided by the solver.

6

Oh that’s my mistake, I did not pay attention to the fact that it used the CUDAdrv.jl implementation of CuArray.

Then I am stuck…

7

It seems unfortunate that IterativeSolvers.jl does not support GPU computations. I have no idea how much work it is to allow it. What is the list of priorities for IterativeSolvers.jl?

Does anyone know of an alternative?