Seed dense arrays of isbits duals without scalar indexing (fixes GPU jacobians)

[ChrisRackauckas-Claude]

What changed and why

ForwardDiff GPU Jacobians currently fail because the seeding helpers scalar-index GPU arrays. An independent clean-checkout bisect found two regression points:

• fce7f76cf8885ded3e150f11645c772961ce8dc5 (released in v1.0.1) replaced broadcast seeding with scalar seed! loops.
• a337ee6658ed2a26bfa7251f1580da59c0b36625 (merged in PR 821 and released in v1.4.4) added the chunk-tail seed_zero_partials! call, which introduced a second scalar-indexing path.

The older seed! failure normally masks the newer tail-clear failure. This PR therefore fixes both paths: dense arrays with isbits values use broadcast for full-buffer writes and map! over contiguous views for chunk writes. A concrete isbits callable carries the seed tuple and offset into map!, avoiding captured closures. Structural wrappers, offset axes, and non-isbits values retain the existing structural fallback.

The runtime implementation uses only Base array interfaces and has no GPUArraysCore dependency or AbstractGPUArray dispatch. JLArrays is added only as a test dependency; it is MIT-licensed, matching ForwardDiff's license.

Downstream fallout

DeepEquilibriumNetworks PR 243 works around the v1.4.4 failure by selecting autodiff=false for GPU steady-state adjoints. Do not revert that workaround before DeepEquilibriumNetworks requires a ForwardDiff release containing this PR: reverting earlier restores the CuArray failure.

After such a ForwardDiff release is required, the regression-specific need for PR 243 is gone. Whether to revert it is then a policy choice:

• Revert it if DeepEquilibriumNetworks intends to use ForwardDiff Jacobians on GPU, after an actual CUDA steady-state adjoint run passes; the direct GPUArraysCore dependency added by PR 243 can then move back to a test extra or be removed if unused.
• Keep it if DeepEquilibriumNetworks intends to continue matching SciMLSensitivity's own GPU-specific autodiff=false default.

In short: treat PR 243 as release-gated cleanup, not an immediate revert.

Regression evidence

The final GPU-like test file was applied to clean master and run with scalar indexing disabled:

zero chunk tail: Error During Test
  Scalar indexing is disallowed.
  [6] _seed_zero_partials!(...)
      @ ForwardDiff src/apiutils.jl:93

Test Summary:                               | Error  Total
ForwardDiff seeding on GPU arrays           |    16     16
ERROR: Some tests did not pass: 0 passed, 0 failed, 16 errored, 0 broken.

The same test on this PR passes:

Test Summary:                     | Pass  Total   Time
ForwardDiff seeding on GPU arrays |   17     17  44.4s

The focused command was:

TMPDIR="$HOME/tmp" ~/.juliaup/bin/julia +release --startup-file=no -e \
  'using Pkg; Pkg.activate(; temp=true); Pkg.develop(path=pwd()); Pkg.add(name="JLArrays", version="0.2"); include("test/GPUArraysTest.jl")'

Full verification

TMPDIR="$HOME/tmp" ~/.juliaup/bin/julia +lts \
  --startup-file=no --project=. -e 'using Pkg; Pkg.test()'

Test Summary:  | Pass  Total     Time
ForwardDiff.jl | 9309   9309  8m15.3s
Testing ForwardDiff tests passed

TMPDIR="$HOME/tmp" ~/.juliaup/bin/julia +release \
  --startup-file=no --project=. -e 'using Pkg; Pkg.test()'

Test Summary:  | Pass  Total      Time
ForwardDiff.jl | 9327   9327  12m35.4s
Testing ForwardDiff tests passed

Both full runs include allocations and JET QA. These checks also exited successfully with no output:

~/.juliaup/bin/julia +release -m Runic --check test/GPUArraysTest.jl test/runtests.jl
typos Project.toml src/apiutils.jl test/GPUArraysTest.jl test/runtests.jl
git diff --check

Fresh Julia 1.12.6 process-local medians against current master after PR 821:

benchmark	master	this PR
gradient!, n=1000	277.2 µs	222.5 µs
jacobian!, n=100	39.3 µs	27.2 µs
jacobian!, n=1000	4.13 ms	3.30 ms

Small n=10 gradient timings varied between process runs, so no directional claim is made for that case.

The review-requested replacement of captured closures with the concrete callable was also compared in fresh Julia processes. Against the preceding PR head, medians were 223.9 µs vs 222.5 µs for gradient! at n=1000, 38.6 µs vs 27.2 µs for jacobian! at n=100, and 3.86 ms vs 3.30 ms for jacobian! at n=1000. Allocation tests in the full suite remain green.

CI

All 27 checks on current head 8697dff pass, including the full Linux/macOS/Windows matrix across minimum-patch, LTS, Julia 1, and prerelease versions with NaN-safe mode enabled and disabled: https://github.com/JuliaDiff/ForwardDiff.jl/actions/runs/33137219763. Codecov reports 90.94% project coverage (+0.25% from master) and covers every modified line: https://app.codecov.io/gh/JuliaDiff/ForwardDiff.jl/pull/816.

The preceding tree-identical CI attempt had one randomized gamma_inc tolerance failure on Windows Julia 1.12.7. The exact values reproduce on clean master with that job's seed and do not involve this PR's code; the failing job and clean-master reproduction are documented at #816 (comment). No assertion was loosened or silenced.

Not verified locally

• No physical CUDA GPU was available. JLArrays exercises GPU array dispatch and the scalar-indexing ban on CPU; CI still needs to cover actual device compilation.
• NaN-safe-enabled, prerelease Julia, Windows, and macOS jobs were not run locally.
• Documentation was not built because this changes no public API, docstring, or documentation source.

Note

This remains a draft opened by an agent on behalf of @ChrisRackauckas. Please ignore it until reviewed by @ChrisRackauckas.

Links

• ForwardDiff PR: Seed dense arrays of isbits duals without scalar indexing (fixes GPU jacobians) #816
• Chunk-tail change: Unseed only the chunk, and split unseed! out of seed! #821
• DeepEquilibriumNetworks workaround: Avoid ForwardDiff Jacobians for GPU steady-state adjoints SciML/DeepEquilibriumNetworks.jl#243
• First scalar-seeding regression: fce7f76
• Chunk-tail scalar-zeroing regression: a337ee6

🤖 Generated with Codex CLI 0.150.1 (model: unknown; session: 01a04603-0ac8-7570-9713-851acf5b8f5d)

[ChrisRackauckas]

Because the broadcasts were purposefully removed due to performance things, we should at least recover the GPU functionality by handling AbstractGPUArray.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 90.94%. Comparing base (b742809) to head (8697dff).

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #816      +/-   ##
==========================================
+ Coverage   90.68%   90.94%   +0.25%     
==========================================
  Files          11       11              
  Lines        1052     1082      +30     
==========================================
+ Hits          954      984      +30     
  Misses         98       98              

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[devmotion]

GPU arrays are always
dense, one-based, and carry isbits element types, so the
structural_eachindex / unset-element handling of the generic methods is
unnecessary on this path and plain broadcast suffices.

It's difficult to verify this claim given that the AbstractGPUArray interface is not documented: https://juliagpu.github.io/GPUArrays.jl/dev/functionality/host/

As an alternative to an extension, I wonder whether special-casing DenseArrays with isbits types to use e.g. map! (broadcasting caused performance regressions IIRC but maybe that would perform better? should be compared with eachindex though) could fix GPU support as well.

[ChrisRackauckas-Claude]

Good suggestion — I benchmarked it, and the DenseArray special-case is both viable and faster than the current scalar loops on the CPU, so I've reworked the PR to that approach and dropped the extension entirely.

On the interface concern: dispatching on DenseArray sidesteps it. GPUArraysCore.AbstractGPUArray <: DenseArray (it's declared that way in GPUArraysCore itself), so a dense fast path in the main package catches all GPU array types without relying on any undocumented AbstractGPUArray property. The gate is

duals isa DenseArray && isbitstype(V) && !Base.has_offset_axes(duals, x)

so nothing needs to be assumed about GPU element types either — non-isbits V (which needs the unset-element handling) and offset axes fall through to the existing structural path.

map! vs broadcast vs the eachindex loop, per seed! method

Isolated seed! benchmarks, Julia 1.12.4 (1.10.11 agrees on all trends), Vector{Float64} input, chunk N = 12:

method	writes	loop (master)	broadcast	map!
seed!(duals, x, seed), n=1000	all n	4.1 μs	2.7 μs	4.1 μs
seed!(duals, x, seeds), n=1000	first N	56 ns	56 ns	52 ns
seed!(duals, x, 500, seed), n=1000	index:end	3.0 μs	1.4 μs	2.0 μs
seed!(duals, x, 500, seeds), n=1000	N at index	458 ns	1.6 μs (7 allocs)	62 ns
seed!(duals, x, 50000, seeds), n=100000	N at index	40 μs	1.6 μs (7 allocs)	68 ns

So the answer to "map! or broadcast?" is: both, per method (bold = what the PR now uses).

• For the full-array write (method 1), broadcast wins; on 1.10 map! is actually ~25% slower than the loop there, so the PR uses broadcast.
• For method 3 (index:end), broadcast is fastest but its dotview allocates 112 bytes under --check-bounds=yes on Julia 1.10 (which Pkg.test uses, and test/AllocationsTest.jl rightly rejects), so the PR uses map! over views — allocation-free on both versions and still 1.5× faster than the loop.
• For the chunk writes (methods 2 and 4), map! wins — but the seeds NTuple can't be a map! source (there's no map!(f, dest, ::Tuple) method), and slicing it at runtime (seeds[1:chunksize], what both 0.10 and the extension in the first iteration of this PR did) is what costs the 1.6 μs + allocations in the broadcast column. map! over the index range with a closure over the tuple avoids both, and GPUArrays' map! handles the range argument fine.
• The last row is the real story for chunk mode: the structural path pays Iterators.drop(structural_eachindex(duals, x), offset) — O(index) per chunk, so O(n²/N) per full chunked jacobian/gradient sweep. The dense path indexes the chunk directly.

On the 2018 revert (#354) that removed broadcast for performance: the regression jrevels measured there was pre-Julia-1.0 broadcast plus an allocating x[dual_inds] slice; with views and Ref none of these forms allocate today.

End-to-end (master → this PR, Julia 1.12.4)

• gradient! n=1000 (chunk 12): 464 μs → 344 μs
• gradient! n=100000 (chunk 12, median of 7 runs): 4.80 s → 3.96 s
• jacobian! n=100: 63 μs → 46 μs
• jacobian! n=1000: 5.92 ms → 4.84 ms
• gradient! n=10 (vector mode): 141 ns → 158 ns — the one regression; seed! itself measures identical (46 ns both), so it's an inlining side effect of the added branch.

One more thing that turned up while measuring: the 1.x rewrite also changed the chunk-mode "unseed" call seed!(xdual, x, i) to write from i to the end of the array each chunk — 0.10 wrote exactly the N chunk elements — which is O(n²) dual writes per sweep (~40 GB of memory traffic for the n=100000 gradient above). Restoring the narrow unseed on top of this PR takes the n=100000 gradient! from ~3.9 s to ~2.4 s. I've kept that out of this PR to keep it focused; follow-up PR coming.

JLArray tests (vector/chunk mode, jacobian!, f!, matrix and view inputs) now exercise the main-package code path directly; the extension and weak dep are gone.

[ChrisRackauckas-Claude]

CI status: all 19 Julia 1 / lts / min-patch jobs (3 OSes × NaN-safe on/off) plus Documentation are green. The 6 "Julia pre" jobs fail, which needs a note:

• Master's pre jobs already fail (every run since ≥2025-12): the JET.@test_opt tests in test/QATest.jl report runtime dispatches on Julia 1.13.0-rc1 whose leaf frames are all in Base internals (broadcast preprocess, simd_index/CartesianIndex, SubArray construction, mightalias/_in_tuple). Reproduced locally on a clean master clone with 1.13.0-rc1: same failures, byte-identical reports. JET's own Julia 1.13 support is still in progress (Add Julia 1.13 support while maintaining Julia 1.12 compatibility aviatesk/JET.jl#815, Declare documented API public (for downstream ExplicitImports) #819), so these look like upstream Julia/JET issues, not ForwardDiff regressions.
• This PR adds one more instance of that same class: the gradient QA test (QATest.jl:11) now also reports a Base-broadcast preprocess dispatch on 1.13-rc1, coming from the new duals .= Dual{T,V,N}.(x, Ref(seed)) in seed!. It is clean on 1.10/1.12 (those QA jobs pass).
• I checked the obvious workarounds on 1.13-rc1: both a closure broadcast ((xi -> Dual{T,V,N}(xi, seed)).(x)) and map! are JET-clean — but both measure 35–50% slower than the Type broadcast for the full-array seed on 1.10 and 1.12 (the Type constructor broadcast SIMD-vectorizes better). Trading permanent seeding performance for a prerelease-only JET false positive seemed like the wrong call, so I left the broadcast in. If the reports persist once JET's 1.13 support settles, that decision can be revisited.

The chunk-unseed O(n²) follow-up mentioned above is now open as #821 (stacked on this branch).

[ChrisRackauckas]

@maleadt maybe should comment here on what is expected to be public and used.

[ChrisRackauckas-Claude]

Rebased onto current master after #821 and resolved the conflict by covering both scalar-indexing sites: the older seed! loops and the v1.4.4 chunk-tail seed_zero_partials! path. The final JLArray regression fails 16/16 testsets on clean master and passes 17/17 assertions here; full local suites pass 9309/9309 on Julia 1.10.12 and 9325/9325 on Julia 1.12.6. The PR body now records the independent bisects and treats SciML/DeepEquilibriumNetworks.jl#243 as release-gated cleanup: do not revert it until DeepEquilibriumNetworks requires a ForwardDiff release containing both fixes.

[ChrisRackauckas-Claude]

All existing human review comments are addressed on current head c91c578:

• The broadcast-performance concern from Seed dense arrays of isbits duals without scalar indexing (fixes GPU jacobians) #816 (comment) is covered by retaining broadcast only for the full-buffer write and using map! for bounded chunk writes. Fresh Julia 1.12.6 medians against current master are recorded in the PR body; the measured n=1000 gradient and Jacobian cases are faster, not regressions.
• The undocumented-AbstractGPUArray concern from Seed dense arrays of isbits duals without scalar indexing (fixes GPU jacobians) #816 (comment) and the public-interface question at Seed dense arrays of isbits duals without scalar indexing (fixes GPU jacobians) #816 (comment) are sidestepped: runtime code neither imports GPUArraysCore nor dispatches on AbstractGPUArray. It uses Base's dense-array/isbits/axes interfaces, while JLArrays remains an MIT-licensed test-only dependency that verifies the GPU-array path with scalar indexing disabled.
• The captured-closure suggestion at Seed dense arrays of isbits duals without scalar indexing (fixes GPU jacobians) #816 (comment) is implemented with a concrete isbits callable carrying seeds and offset. The thread has the focused/full test totals and before/after callable benchmarks and is resolved.

The post-review full suites pass 9309/9309 tests on Julia 1.10.12 and 9327/9327 on Julia 1.12.6, including allocation and QA checks. CI for c91c578 is now running.

The downstream conclusion remains release-gated: do not revert SciML/DeepEquilibriumNetworks.jl#243 until DeepEquilibriumNetworks requires a ForwardDiff release containing this fix. After that, reverting is appropriate only if the package intends to restore ForwardDiff GPU Jacobians and an actual CUDA steady-state-adjoint regression passes; otherwise keeping autodiff=false continues to match SciMLSensitivity's GPU policy.

[ChrisRackauckas-Claude]

CI attempt 1 for tree b992631de5526ebc36ca2c76e6e095f74b76347a finished 26/27 green. The sole failure was an unrelated randomized gamma_inc tolerance case in test/DualTest.jl:661 on Windows Julia 1.12.7: https://github.com/JuliaDiff/ForwardDiff.jl/actions/runs/33136312489/job/98736955860. That job's GPUArrays test passed 17/17, and its seeding, allocation, and QA sections all passed.

I reproduced the same four values and failures byte-for-byte on clean master with Julia 1.12.7 and the CI seed 1787884652. I did not loosen or silence the assertion. GitHub rejected gh run rerun --failed because this fork token lacks repository admin rights, so I pushed tree-identical commit 8697dff solely to trigger a fresh matrix and random seed. A separate clean-master bisect/root-cause investigation required by the repository instructions is running in parallel; its result will be reported independently from this PR's GPU fix.

The clean-master investigation is now complete:

• A fresh-environment git bisect run identified 61a7c5b (Define SpecialFunctions.gamma_inc for ForwardDiff.Dual) as the first bad commit. Its parent has neither the method nor this test; the first bad commit introduces both the method and the cancellation-prone Float64 finite-difference oracle.
• For the failing input a=0.009112825048747886, x=10, the ForwardDiff derivative is 4.2469380334661066e-8, matching the closed-form derivative bit-for-bit. The Float64 Calculus oracle is 4.246857314027132e-8, about 1.90e-5 relative error; a BigFloat Calculus oracle has about 3.58e-46 relative error. This is a test-oracle cancellation problem, not a ForwardDiff source defect.
• A separate test-only correction computes that oracle with BigFloat while retaining the existing rtol. The exact CI seed passes all Partials/Dual assertions, and full clean-master Pkg.test() passes 9308/9308 on Julia 1.12.7 in 9m36.7s. Runic on the changed lines, typos, and git diff --check also pass.

That unrelated test correction is intentionally not mixed into this GPU-seeding PR. Current PR head 8697dff remains tree-identical to the locally validated code and has 27/27 green checks: https://github.com/JuliaDiff/ForwardDiff.jl/actions/runs/33137219763.

[ChrisRackauckas-Claude]

I checked whether the dense map! fast path is close enough to direct indexing that the scalar-indexing path can be removed.

On current head 8697dffef2633cc798afd185874fe818b5a5baed, I ran isolated chunk-seeding benchmarks locally with Julia 1.10.11 and 1.12.7. Both implementations were allocation-free:

runtime	dense map!	direct dense indexing
Julia 1.10.11	90–110 ns	80–100 ns
Julia 1.12.7	90–130 ns	80–110 ns

Chunk zeroing was effectively identical. For a 100,000-element full-buffer zero, broadcast, map!, and direct indexing were also effectively tied at about 310–323 μs.

The dense map! path is therefore fast enough, and it avoids the offset-dependent cost of the generic structural iterator. At index 99,989 of a 100,000-element vector, the structural path took 79.45 μs on Julia 1.10.11 and 61.58 μs on Julia 1.12.7, versus 100–120 ns for map!.

However, the structural/indexing fallback cannot be removed:

• Chunk map! over contiguous linear views fails for Diagonal and UpperTriangular: the second linear position is off-structure, and assigning a dual with nonzero partials there throws ArgumentError.
• Full zeroing with a generic map! happens to be valid for those wrappers because the off-structure values are zero, but it is substantially slower. For 100×100 inputs, I measured 83.62 μs versus 1.37 μs for Diagonal (61× slower), and 68.83 μs versus 18.87 μs for UpperTriangular (3.65× slower).
• For non-isbits arrays containing unassigned elements, map! throws UndefRefError; the existing fallback checks isassigned and preserves the holes.

Conclusion: keep the current split. Dense, one-based arrays with isbits values should use the map! fast path, while structured wrappers and non-isbits values still require the structural/indexing fallback. There is no remaining scalar-indexing path for the dense/GPU case.

Commands used:

~/.juliaup/bin/julia +1.10 --project=../bench_env_110 ../bench_env/map_vs_index.jl
~/.juliaup/bin/julia +1.12 --project=../bench_env ../bench_env/map_vs_index.jl

Links

• PR: Seed dense arrays of isbits duals without scalar indexing (fixes GPU jacobians) #816
• Earlier benchmark discussion: Seed dense arrays of isbits duals without scalar indexing (fixes GPU jacobians) #816 (comment)

[ChrisRackauckas]

@devmotion I think this is good for another review now.