CUDA solver paths error with ForwardDiff >= 1: seed! scalar-indexes GPU dual arrays

[ChrisRackauckas-Claude]

Summary

With ForwardDiff ≥ 1 in the environment, every solver that computes jacobians through ForwardDiff (e.g. NonlinearReformulation / BokhovenIterativeAlgorithm with the default SimpleNewtonRaphson() nonlinear solver, whose default AD is AutoForwardDiff) errors on CUDA arrays with:

Scalar indexing is disallowed.
Invocation of getindex resulted in scalar indexing of a GPU array.
...
[6] seed!(duals::CuArray{ForwardDiff.Dual{...}}, x::CuArray{Float32, 1, ...}, seeds::Tuple{...})
  @ ForwardDiff .../ForwardDiff/src/apiutils.jl:95
[7] vector_mode_dual_eval!
[8] vector_mode_jacobian
...
[11] compute_jacobian!! @ SimpleNonlinearSolve/src/utils.jl:153

This is what made the CUDA Core CI job fail on #63 (55 errored tests) once the test environment started resolving ForwardDiff 1.4 (previously held back at 0.10.39 by ParametricMCPs).

Root cause (upstream, ForwardDiff)

ForwardDiff 1.x rewrote all four seed! methods in src/apiutils.jl to use scalar setindex! loops over structural_eachindex:

for (i, idx) in zip(1:N, structural_eachindex(duals, x))
    duals[idx] = Dual{T,V,N}(x[idx], seeds[i])
end

ForwardDiff 0.10 used broadcast-based seeding, which worked on GPU arrays:

dual_inds = 1:N
duals[dual_inds] .= Dual{T,V,N}.(view(x, dual_inds), seeds)

So this is an upstream ForwardDiff regression for GPU arrays, not something specific to this package. The last green CUDA Core run on main (2026-04-13) resolved ForwardDiff 0.10.39; with the dependency updates in #63 the environment resolves ForwardDiff 1.4 (required by NonlinearSolve ≥ 4.17), and the GPU jacobian paths break.

Verified workaround (not applied)

Broadcast-based seed! overloads for GPUArraysCore.AbstractGPUArray{<:Dual} restore the 0.10 behavior and make the full CUDA Core suite pass (verified on the GPU CI runners in #63 before being reverted, and locally on JLArrays in both vector and chunk mode). This was prototyped as type piracy in #63 (commit d1d3d70) and removed per review (c3a8fa9, 06b7c1d). The proper home for that fix is ForwardDiff itself (e.g. a GPUArraysCore package extension upstream).

Note ForwardDiff.gradient on GPU arrays has an additional scalar-indexing site in the gradient extraction path (extract_gradient_chunk!), so fixing seed! alone covers jacobians (what this package's solvers use) but not gradients.

Consequences here until fixed upstream

• The CUDA Core CI job fails when the environment resolves ForwardDiff ≥ 1 (which is forced once NonlinearSolve ≥ 4.17 is in the environment).
• GPU usage of ForwardDiff-jacobian-based solvers requires ForwardDiff 0.10, which is incompatible with NonlinearSolve ≥ 4.17 / NNlib ≥ 0.9.32 / CUDA ≥ 6.

🤖 Generated with Claude Code

[ChrisRackauckas-Claude]

Upstream fix opened: JuliaDiff/ForwardDiff.jl#816 adds a GPUArraysCore extension with broadcast-based seed! methods for AbstractGPUArray{<:Dual}, restoring the pre-1.0 GPU-compatible seeding. This covers the jacobian paths (vector and chunk mode) that this package's solvers use. Tested via JLArrays in ForwardDiff's CI.

Note (as flagged in the issue): this fixes seed! only — ForwardDiff.gradient on GPU arrays still has a separate scalar-indexing site in extract_gradient_chunk!, which is out of scope for that PR.

Once it merges and tags, the broadcast seed! overloads prototyped here can be dropped.