minres_qlp performs scalar operations on GPU

[mtanneau]

Something is happening inside minres_qlp that seems to not be GPU-friendly

using Random, LinearAlgebra, CUDA, Krylov

Random.seed!(0)
n = 64
A = Matrix(Symmetric(rand(n, n)))
cA = CuArray(A)

b = ones(n)
cb = CuArray(b)

x, stats = minres_qlp(A, b);
cx, cstats = minres_qlp(cA, cb);

The last call leads to this warning:

┌ Warning: Performing scalar operations on GPU arrays: This is very slow, consider disallowing these operations with `allowscalar(false)`
└ @ GPUArrays ~/.julia/packages/GPUArrays/eVYIC/src/host/indexing.jl:43

Subsequent timings look like

julia> @time x, stats = minres_qlp(A, b);
  0.000091 seconds (28 allocations: 8.391 KiB)

julia> CUDA.@time cx, stats = minres_qlp(cA, cb);
  0.088764 seconds (91.31 k CPU allocations: 3.673 MiB) (6 GPU allocations: 3.000 KiB, 0.06% gc time of which 79.26% spent allocating)

---

julia> versioninfo()
Julia Version 1.5.1
Commit 697e782ab8 (2020-08-25 20:08 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i5-6600 CPU @ 3.30GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-9.0.1 (ORCJIT, skylake)

[mtanneau]

Digging into it, this seems to come from some internal tools.
Setting CUDA.allowscalar(false) disables the warning and throws an error instead.

julia> CUDA.allowscalar(false);

julia> cx, stats = minres_qlp(cA, cb);
ERROR: scalar getindex is disallowed
Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] assertscalar(::String) at /home/mtanneau/.julia/packages/GPUArrays/eVYIC/src/host/indexing.jl:41
 [3] getindex(::CuArray{Float64,1}, ::Int64) at /home/mtanneau/.julia/packages/GPUArrays/eVYIC/src/host/indexing.jl:96
 [4] macro expansion at /home/mtanneau/.julia/packages/Krylov/axdoL/src/krylov_aux.jl:153 [inlined]
 [5] macro expansion at ./simdloop.jl:77 [inlined]
 [6] krylov_ref! at /home/mtanneau/.julia/packages/Krylov/axdoL/src/krylov_aux.jl:152 [inlined]
 [7] minres_qlp(::LinearOperators.PreallocatedLinearOperator{Float64}, ::CuArray{Float64,1}; M::LinearOperators.opEye, atol::Float64, rtol::Float64, λ::Float64, itmax::Int64, verbose::Bool) at /home/mtanneau/.julia/packages/Krylov/axdoL/src/minres_qlp.jl:252
 [8] minres_qlp at /home/mtanneau/.julia/packages/Krylov/axdoL/src/minres_qlp.jl:32 [inlined]
 [9] #minres_qlp#100 at /home/mtanneau/.julia/packages/Krylov/axdoL/src/variants.jl:49 [inlined]
 [10] minres_qlp(::CuArray{Float64,2}, ::CuArray{Float64,1}) at /home/mtanneau/.julia/packages/Krylov/axdoL/src/variants.jl:49
 [11] top-level scope at REPL[11]:1

[amontoison]

The problem comes from the function krylov_ref!, I added it in Julia some months ago and it's now available with julia 1.5 (reflect! and rotate! functions). I just need to add specialized version in CuArrays.jl. With the multiple dispatch the optimized function will be used.