How to set the diagonal part of a GPU Array


#1

Let’s see the code

julia> using CuArrays;

julia> s
4-element CuArray{Float32,1}:
 3.2211537
 2.1344962
 0.7451281
 0.3308171

julia> v
4×4 CuArray{Float32,2}:
  0.267688   0.355418  -0.311468  -0.839648
 -0.156514   0.844306  -0.296999   0.417663
  0.658483  -0.260938  -0.623603   0.330803
 -0.685745  -0.304527  -0.65261   -0.105441

julia> v + Diagonal(s)   # I disabled getindex of cuarrays.
ERROR: scalar getindex is disallowed
...

How to make it correct?


#2

This seems a broadcast issue, please file an issue at CuArrays.jl.

julia> v + Diagonal(s)
ERROR: scalar getindex is disallowed
Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] assertscalar at /home/tbesard/Julia/GPUArrays/src/indexing.jl:8 [inlined]
 [3] getindex(::CuArray{Float32,2}, ::Int64) at /home/tbesard/Julia/GPUArrays/src/indexing.jl:44
 [4] _getindex at ./abstractarray.jl:928 [inlined]
 [5] getindex at ./abstractarray.jl:905 [inlined]
 [6] _broadcast_getindex at ./broadcast.jl:540 [inlined]
 [7] _getindex at ./broadcast.jl:570 [inlined]
 [8] _broadcast_getindex at ./broadcast.jl:546 [inlined]
 [9] getindex at ./broadcast.jl:507 [inlined]
 [10] macro expansion at ./broadcast.jl:838 [inlined]
 [11] macro expansion at ./simdloop.jl:73 [inlined]
 [12] copyto! at ./broadcast.jl:837 [inlined]
 [13] copyto!(::Array{Float32,2}, ::Base.Broadcast.Broadcasted{Base.Broadcast.ArrayConflict,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}},typeof(+),Tuple{CuArray{Float32,2},Diagonal{Float32,CuArray{Float32,1}}}}) at ./broadcast.jl:792
 [14] copy at ./broadcast.jl:768 [inlined]
 [15] materialize at ./broadcast.jl:748 [inlined]
 [16] broadcast at ./broadcast.jl:702 [inlined]
 [17] +(::CuArray{Float32,2}, ::Diagonal{Float32,CuArray{Float32,1}}) at ./arraymath.jl:39
 [18] top-level scope at none:0

The Broadcast.ArrayConflict looks iffy here, but I’m no expert of the broadcast internals.


#3

Yes, that’s exactly the issue — both structured matrices and CuArrays define custom broadcasting implementations but haven’t decided which one “wins,” so Broadcast falls back to the totally generic builtin version.


#4

Issue has been filed here. Am I supposed to write CUDAnative to set diagonal elements now?


#5

You can just collect the diagonal and use that, as a temporary workaround:

julia> v + CuArray(Diagonal(s))
4×4 CuArray{Float32,2}:
 1.66536   0.470219  0.71401   0.319409
 0.44635   1.43419   0.358502  0.810913
 0.667319  0.776175  0.720439  0.154112
 0.109931  0.294766  0.292015  1.41413 

julia> CuArray(Diagonal(s))
4×4 CuArray{Float32,2}:
 0.740219  0.0       0.0        0.0     
 0.0       0.920994  0.0        0.0     
 0.0       0.0       0.0390205  0.0     
 0.0       0.0       0.0        0.968963

#6

Or you can assign directly into the diagonal using this linear indexing trick:

v[1:size(v, 1)+1:end] .+= s

#7

v + Diagonal(s)
now this works