The proper way to do automatic sparse differentiation on GPU

I was trying to use DI+SMC+SCT on a small example on CUDA, not sure I am doing this right, really appreciate some suggestions.

using DifferentiationInterface
using SparseConnectivityTracer, SparseMatrixColorings
using SparseArrays, StableRNGs
using ForwardDiff
using ADTypes: ADTypes
using CUDA
const DI = DifferentiationInterface

n = 1000
J = spdiagm(0 => ones(n))

function MyAutoSparse(backend)
    sparsity_detector = ADTypes.KnownJacobianSparsityDetector(J)
    coloring_algorithm = GreedyColoringAlgorithm(RandomOrder(StableRNG(0)))
    return AutoSparse(backend; sparsity_detector, coloring_algorithm)
end

diffmode = MyAutoSparse(AutoForwardDiff())

x, y = CuArray(ones(n)), CuArray(zeros(n))

diffcache = DI.prepare_jacobian(copyto!, y, diffmode, x)
J = DI.jacobian(copyto!, y, diffcache, diffmode, x)

ERROR: LoadError: Scalar indexing is disallowed.
Invocation of getindex resulted in scalar indexing of a GPU array.
This is typically caused by calling an iterating implementation of a method.
Such implementations *do not* execute on the GPU, but very slowly on the CPU,
and therefore should be avoided.

If you want to allow scalar iteration, use `allowscalar` or `@allowscalar`
to enable scalar iteration globally or for the operations in question.
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:44
  [2] errorscalar(op::String)
    @ GPUArraysCore ~/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:151
  [3] _assertscalar(op::String, behavior::GPUArraysCore.ScalarIndexing)
    @ GPUArraysCore ~/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:124
  [4] assertscalar(op::String)
    @ GPUArraysCore ~/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:112
  [5] getindex
    @ ~/.julia/packages/GPUArrays/3a5jB/src/host/indexing.jl:50 [inlined]
  [6] decompress!(A::SparseMatrixCSC{…}, B::CuArray{…}, result::SparseMatrixColorings.ColumnColoringResult{…})
    @ SparseMatrixColorings ~/.julia/packages/SparseMatrixColorings/cNs01/src/decompression.jl:379
  [7] _sparse_jacobian_aux!(::Tuple{…}, ::SparseMatrixCSC{…}, ::DifferentiationInterfaceSparseMatrixColoringsExt.SMCPushforwardSparseJacobianPrep{…}, ::AutoSparse{…}, ::CuArray{…})
    @ DifferentiationInterfaceSparseMatrixColoringsExt ~/.julia/packages/DifferentiationInterface/afUhd/ext/DifferentiationInterfaceSparseMatrixColoringsExt/jacobian.jl:326
  [8] jacobian!
    @ ~/.julia/packages/DifferentiationInterface/afUhd/ext/DifferentiationInterfaceSparseMatrixColoringsExt/jacobian.jl:237 [inlined]
  [9] jacobian(::typeof(copyto!), ::CuArray{…}, ::DifferentiationInterfaceSparseMatrixColoringsExt.SMCPushforwardSparseJacobianPrep{…}, ::AutoSparse{…}, ::CuArray{…})
    @ DifferentiationInterfaceSparseMatrixColoringsExt ~/.julia/packages/DifferentiationInterface/afUhd/ext/DifferentiationInterfaceSparseMatrixColoringsExt/jacobian.jl:250
 [10] top-level scope
    @ ~/work/demo.jl:35
 [11] include(mapexpr::Function, mod::Module, _path::String)
    @ Base ./Base.jl:307
 [12] top-level scope
    @ REPL[4]:1
in expression starting at /home/featurize/work/demo.jl:35
Some type information was truncated. Use `show(err)` to see complete types.

There are some scalar indexing errors with SparseMatrixColorings.jl when using the greedy coloring algorithm. I wonder if there should be some documentatios for this feature?

Reproduce the CUDA example first, then inspect DifferentiationInterfaceSparseMatrixColoringsExt/jacobian.jl and the SparseMatrixColorings decompression.jl call shown in the stack trace. Determine the supported path for automatic sparse Jacobians on CuArray inputs and document the required usage or limitation, with the example no longer producing an unexplained scalar-indexing error.