How to set the diagonal part of a GPU Array

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?

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.

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.

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

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
2 Likes

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

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

v + Diagonal(s)
now this works

1 Like