Scalar indexing error when computing Jacobian with a `CUDA.CuArray`

Consider the following minimal reproducible example:

using CUDA, ForwardDiff
x = cu(ones(3))
f(x) = x.^2
ForwardDiff.jacobian(f, x)

Result on ForwardDiff 1.0.0:

julia> ForwardDiff.jacobian(f, x)
3×3 CuArray{Float32, 2, CUDA.DeviceMemory}:
 2.0  0.0  0.0
 0.0  2.0  0.0
 0.0  0.0  2.0

which is obviously correct.

On ForwardDiff 1.0.1:

ERROR: 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.

Issue introduced in: https://github.com/JuliaDiff/ForwardDiff.jl/pull/739

Run the minimal CUDA and ForwardDiff example from the issue, comparing ForwardDiff 1.0.0 with 1.0.1. Read the changes in PR #739 to identify the regression, then verify that ForwardDiff.jacobian(f, x) works on a CUDA CuArray without scalar-indexing errors.