BandedMatrix does not support GPU?

[chaopengshen]

Dear BandedMatrix devs,

I'm trying to put some GPU matrices ((CuArray)) on BandedMatrix. But it seems it has turned the cuArray into CPU arrays. Please see the following.
Do you think you can add GPU support? The GPU computation is so important nowadays. And I was told it should not be that hard to fix this (although I don't know how to).

julia> using BandedMatrices

julia> using CUDA

julia> A0 = cu(rand(10)) #alternatively, do "using Flux" and then "A0=gpu(rand(10))"
10-element CuArray{Float32, 1}:
0.9238173
0.39963627
0.74860376
0.5007854
0.41714704
0.41270784
0.36631182
0.75199664
0.14709264
0.17270018

julia> A = BandedMatrix(0 => A0)
┌ Warning: Performing scalar operations on GPU arrays: This is very slow, consider disallowing these operations with allowscalar(false)
└ @ GPUArrays C:\Users\cxs1024.julia\packages\GPUArrays\Z5nPF\src\host\indexing.jl:64
10×10 BandedMatrix{Float32} with bandwidths (0, 0):
0.923817 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ 0.399636 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ 0.748604 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ 0.500785 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ 0.417147 ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ 0.412708 ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.366312 ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.751997 ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.147093 ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.1727

julia> typeof(A.data)
Matrix{Float32} (alias for Array{Float32, 2})

julia> A.data
1×10 Matrix{Float32}:
0.923817 0.399636 0.748604 0.500785 0.417147 0.412708 0.366312 0.751997 0.147093 0.1727

Also, regarding ”┌ Warning: Performing scalar operations“ above. if I take off allowscalar, it showed where the scalar op is.

julia> CUDA.allowscalar(true)
julia> A = BandedMatrix(0 => A0)
ERROR: scalar getindex is disallowed
Stacktrace:
[1] error(s::String)
@ Base .\error.jl:33
[2] assertscalar(op::String)
@ GPUArrays C:\Users\cxs1024.julia\packages\GPUArrays\Z5nPF\src\host\indexing.jl:62
[3] getindex
@ C:\Users\cxs1024.julia\packages\GPUArrays\Z5nPF\src\host\indexing.jl:104 [inlined]
[4] iterate
@ .\abstractarray.jl:1096 [inlined]
[5] iterate
@ .\abstractarray.jl:1094 [inlined]
[6] _unsafe_setindex!(::IndexLinear, ::Matrix{Float32}, ::CuArray{Float32, 1}, ::Int64, ::UnitRange{Int64})
@ Base .\multidimensional.jl:897
[7] _setindex!
@ .\multidimensional.jl:887 [inlined]
[8] setindex!
@ .\abstractarray.jl:1267 [inlined]
[9] (BandedMatrix{Float32, CONTAINER, RAXIS} where {CONTAINER, RAXIS})(kv::Tuple{Pair{Int64, CuArray{Float32, 1}}}, mn::Tuple{Int64, Int64}, lu::Tuple{Int64, Int64})
@ BandedMatrices C:\Users\cxs1024.julia\packages\BandedMatrices\Gbqeq\src\banded\BandedMatrix.jl:254
[10] BandedMatrix
@ C:\Users\cxs1024.julia\packages\BandedMatrices\Gbqeq\src\banded\BandedMatrix.jl:286 [inlined]
[11] BandedMatrix(kv::Pair{Int64, CuArray{Float32, 1}})
@ BandedMatrices C:\Users\cxs1024.julia\packages\BandedMatrices\Gbqeq\src\banded\BandedMatrix.jl:295
[12] top-level scope
@ REPL[10]:1
[13] top-level scope
@

[dlfivefifty]

You probably want to create the data directly. An (ancient) example is here:

https://github.com/JuliaMatrices/BandedMatrices.jl/blob/master/examples/clarrays.jl

[chaopengshen]

Is it possible to overload lu with CuArray so everything could carry automatically? That will be a lot painless.

[dlfivefifty]

It just calls LAPACK

https://github.com/JuliaMatrices/BandedMatrices.jl/blob/0d506eb5bfeb51eb9477763e6152370517cdf03b/src/banded/BandedLU.jl#L87

so it should work find for CuArray provided the relevant LAPACK call is overloaded