Incorrect derivative of getindex() with repeating indices on CuArrays

[dfdx]

using Zygote
using CUDA

f(A, I) = sum(A[I])
A = rand(4)
I = [1, 3, 1]

# CPU - everything is OK
Zygote.gradient(f, A, I)
# ==> ([2.0, 0.0, 1.0, 0.0], nothing)

# GPU - dA[1] is incorrect
Zygote.gradient(f, cu(A), cu(I))
# => (Float32[1.0, 0.0, 1.0, 0.0], nothing)

I believe CPU version comes from ChainRules.jl which correctly adds several derivatives to dA[1], but I'm not sure what code is used for CUDA version.

Here is how I came to this issue and how I try to resolve it in Yota.

[mcabbott]

Similar to #600, see also JuliaGPU/CUDA.jl#89 and JuliaLang/julia#31407 .

[marius311]

My brief reading of the linked CUDA issue suggests this can't be fixed? Could we add an error here? Just spent a bunch of time discovering subtly wrong results in my code / reducing to a MWE / eventually finding this already-filed issue. 😢

[ToucheSir]

I thought JuliaDiff/ChainRules.jl#655 was the last word on this. If that works then we "just" need to do #1328.