From 05839035378bb595decfcddf2cbecb49007756aa Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 8 Jun 2026 22:13:17 +0000 Subject: [PATCH 1/7] Bump the all-julia-packages group across 1 directory with 5 updates Updates the requirements on [CUDA](https://github.com/JuliaGPU/CUDA.jl), [Statistics](https://github.com/JuliaStats/Statistics.jl), [Optimization](https://github.com/SciML/Optimization.jl), [OrdinaryDiffEq](https://github.com/SciML/OrdinaryDiffEq.jl) and [DiffEqBase](https://github.com/SciML/OrdinaryDiffEq.jl) to permit the latest version. Updates `CUDA` to 6.1.0 - [Release notes](https://github.com/JuliaGPU/CUDA.jl/releases) - [Commits](https://github.com/JuliaGPU/CUDA.jl/commits/v6.1.0) Updates `Statistics` to 1.11.1 - [Release notes](https://github.com/JuliaStats/Statistics.jl/releases) - [Commits](https://github.com/JuliaStats/Statistics.jl/compare/v1.4.0...v1.11.1) Updates `Optimization` to 5.6.1 - [Release notes](https://github.com/SciML/Optimization.jl/releases) - [Changelog](https://github.com/SciML/Optimization.jl/blob/master/NEWS.md) - [Commits](https://github.com/SciML/Optimization.jl/commits/v5.6.1) Updates `OrdinaryDiffEq` to 7.0.0 - [Release notes](https://github.com/SciML/OrdinaryDiffEq.jl/releases) - [Changelog](https://github.com/SciML/OrdinaryDiffEq.jl/blob/master/NEWS.md) - [Commits](https://github.com/SciML/OrdinaryDiffEq.jl/compare/v6.0.0...v7.0.0) Updates `DiffEqBase` to 7.5.5 - [Release notes](https://github.com/SciML/OrdinaryDiffEq.jl/releases) - [Changelog](https://github.com/SciML/OrdinaryDiffEq.jl/blob/master/NEWS.md) - [Commits](https://github.com/SciML/OrdinaryDiffEq.jl/commits) --- updated-dependencies: - dependency-name: CUDA dependency-version: 6.1.0 dependency-type: direct:production dependency-group: all-julia-packages - dependency-name: DiffEqBase dependency-version: 7.5.5 dependency-type: direct:production dependency-group: all-julia-packages - dependency-name: Optimization dependency-version: 5.6.1 dependency-type: direct:production dependency-group: all-julia-packages - dependency-name: OrdinaryDiffEq dependency-version: 7.0.0 dependency-type: direct:production dependency-group: all-julia-packages - dependency-name: Statistics dependency-version: 1.11.1 dependency-type: direct:production dependency-group: all-julia-packages ... Signed-off-by: dependabot[bot] --- test/Project.toml | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/test/Project.toml b/test/Project.toml index 3902bb7..6fe19d0 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -29,23 +29,24 @@ Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [compat] Aqua = "0.8" BenchmarkTools = "1" -CUDA = "4, 5" +CUDA = "4, 5, 6.1" ChainRulesCore = "1" ComponentArrays = "0.15" -DiffEqBase = "6" +DiffEqBase = "6, 7.5" ExplicitImports = "1" FiniteDifferences = "0.12" ForwardDiff = "0.10, 1" NonlinearSolve = "1, 2, 3, 4" Optimisers = "0.3, 0.4" -Optimization = "3, 4" +Optimization = "3, 4, 5.6" OptimizationOptimisers = "0.2, 0.3" -OrdinaryDiffEq = "6" +OrdinaryDiffEq = "6, 7.0" PATHSolver = "1" ParametricMCPs = "0.1" SafeTestsets = "0.1" SciMLSensitivity = "7" SimpleNonlinearSolve = "0.1, 1, 2" StableRNGs = "1" +Statistics = "1.11.1" SteadyStateDiffEq = "2" Zygote = "0.6, 0.7" From 3c1764915f1318b3d86790b1106cae389d43d657 Mon Sep 17 00:00:00 2001 From: "Chris Rackauckas (Claude)" Date: Tue, 9 Jun 2026 21:54:51 -0400 Subject: [PATCH 2/7] Remove ParametricMCPs test dependency to unblock dependency updates ParametricMCPs (including its latest release, 0.1.17) caps ForwardDiff at 0.10 through a weak dependency, which makes it impossible to resolve a test environment containing NonlinearSolve >= 4.17 (requires ForwardDiff 1), NNlib >= 0.9.32, or CUDA >= 6. This is what made the GPU CI jobs on the dependabot update PR fail at version resolution, and it silently held the CPU test jobs back on DiffEqBase 6 / OrdinaryDiffEq 6 / CUDA 5 / ForwardDiff 0.10. Dropping the ParametricMCPs comparison benchmark lets the test environment resolve the full updated stack (CUDA 6.1, DiffEqBase 7.5, OrdinaryDiffEq 7.0, Optimization 5.6, ForwardDiff 1.4, NonlinearSolve 4.19, SciMLSensitivity 7.112, SciMLBase 3.18). The ComplementaritySolve.jl side of the benchmark testset is kept since it exercises forward solves and Zygote adjoints for both in-place and out-of-place problems. Core, Applications, and QA test groups all pass locally against the updated stack with no source changes required. Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Fable 5 --- test/Project.toml | 2 -- test/core/cpu/mcp.jl | 27 ++++----------------------- 2 files changed, 4 insertions(+), 25 deletions(-) diff --git a/test/Project.toml b/test/Project.toml index 6fe19d0..09d39d9 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -15,7 +15,6 @@ Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba" OptimizationOptimisers = "42dfb2eb-d2b4-4451-abcd-913932933ac1" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" PATHSolver = "f5f7c340-0bb3-5c69-969a-41884d311d1b" -ParametricMCPs = "9b992ff8-05bb-4ea1-b9d2-5ef72d82f7ad" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1" SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7" @@ -42,7 +41,6 @@ Optimization = "3, 4, 5.6" OptimizationOptimisers = "0.2, 0.3" OrdinaryDiffEq = "6, 7.0" PATHSolver = "1" -ParametricMCPs = "0.1" SafeTestsets = "0.1" SciMLSensitivity = "7" SimpleNonlinearSolve = "0.1, 1, 2" diff --git a/test/core/cpu/mcp.jl b/test/core/cpu/mcp.jl index 0c05945..56ff819 100644 --- a/test/core/cpu/mcp.jl +++ b/test/core/cpu/mcp.jl @@ -1,6 +1,5 @@ using BenchmarkTools, ComplementaritySolve, ComponentArrays, FiniteDifferences using ForwardDiff, NonlinearSolve, SimpleNonlinearSolve, StableRNGs, Test, Zygote -import ParametricMCPs rng = StableRNG(0) @@ -70,7 +69,10 @@ rng = StableRNG(0) @test ∂θ_zygote ≈ ∂θ_finitediff atol = 1.0e-3 rtol = 1.0e-3 end - @testset "Benchmarking against ParametricMCPs.jl" begin + # ParametricMCPs.jl was previously benchmarked here as well, but its latest + # release caps ForwardDiff at 0.10 (via a weak dependency), which is + # incompatible with NonlinearSolve >= 4.17, NNlib >= 0.9.32, and CUDA >= 6. + @testset "Benchmarking" begin u0 = randn(rng, Float64, 4) lb = Float64[-Inf, -Inf, 0, 0] ub = Float64[Inf, Inf, Inf, Inf] @@ -88,18 +90,6 @@ rng = StableRNG(0) return sum(abs2, sol.u) end - prob_ext = ParametricMCPs.ParametricMCP(f, lb, ub, length(θ)) - function loss_function_parametric_mcp(θ) - sol = ParametricMCPs.solve(prob_ext, θ) - return sum(abs2, sol.z) - end - - function loss_function_parametric_mcp_total(θ) - prob = Zygote.@ignore ParametricMCPs.ParametricMCP(f, lb, ub, length(θ)) - sol = ParametricMCPs.solve(prob, θ) - return sum(abs2, sol.z) - end - loss_function_path_oop = Base.Fix2(loss_function_oop, PATHSolverAlgorithm()) loss_function_nr_oop = Base.Fix2( loss_function_oop, @@ -121,15 +111,6 @@ rng = StableRNG(0) t₂ = @belapsed only(Zygote.gradient($loss_function, $θ)) @info "ComplementaritySolve.jl: $(loss_function)" fwd_time = t₁ with_adjoint_time = t₂ end - - for loss_function in ( - loss_function_parametric_mcp, - loss_function_parametric_mcp_total, - ) - t₁ = @belapsed $loss_function($θ) - t₂ = @belapsed only(Zygote.gradient($loss_function, $θ)) - @info "ParametricMCPs.jl: $(loss_function)" fwd_time = t₁ with_adjoint_time = t₂ - end end end end From d1d3d70680e694f5db2501d13206e734aed2f5c8 Mon Sep 17 00:00:00 2001 From: "Chris Rackauckas (Claude)" Date: Tue, 9 Jun 2026 23:32:31 -0400 Subject: [PATCH 3/7] Fix GPU ForwardDiff seeding and Julia 1.10 Statistics compat ForwardDiff 1.x rewrote seed! to use scalar setindex! loops (via structural_eachindex), which errors on GPU arrays with scalar indexing disallowed; ForwardDiff 0.10 used broadcast and worked. This broke every NonlinearReformulation/BokhovenIterativeAlgorithm CUDA test through SimpleNonlinearSolve's ForwardDiff jacobians. Add broadcast-based seed! overloads for AbstractGPUArray in the existing type-piracy section (should eventually live upstream in ForwardDiff as a GPUArraysCore extension), with a JLArrays-based regression test covering vector and chunk mode. Also relax the Statistics compat to "1.10, 1.11": the stdlib is version 1.10.0 on Julia 1.10, so "1.11.1" was unsatisfiable on the downgrade CI jobs which run on the LTS. Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Fable 5 --- src/ComplementaritySolve.jl | 41 +++++++++++++++++++++++++++++++ test/Project.toml | 4 ++- test/core/cpu/forwarddiff_seed.jl | 22 +++++++++++++++++ test/runtests.jl | 3 +++ 4 files changed, 69 insertions(+), 1 deletion(-) create mode 100644 test/core/cpu/forwarddiff_seed.jl diff --git a/src/ComplementaritySolve.jl b/src/ComplementaritySolve.jl index 3042a3c..c3b740d 100644 --- a/src/ComplementaritySolve.jl +++ b/src/ComplementaritySolve.jl @@ -50,6 +50,47 @@ const DEFAULT_NLSOLVER = SimpleNewtonRaphson() ArrayInterfaceCore.can_setindex(::Type{<:AbstractFill}) = false ArrayInterfaceCore.can_setindex(::Zygote.OneElement) = false +# ForwardDiff 1.x seeds dual arrays with scalar `setindex!` loops (via +# `structural_eachindex`), which errors on GPU arrays with scalar indexing +# disallowed. ForwardDiff 0.10 used broadcast and worked on GPU arrays, so +# restore broadcast-based seeding for them. TODO: upstream to ForwardDiff as a +# GPUArraysCore extension. +function ForwardDiff.seed!( + duals::GPUArraysCore.AbstractGPUArray{ForwardDiff.Dual{T, V, N}}, x, + seed::ForwardDiff.Partials{N, V} = zero(ForwardDiff.Partials{N, V}) + ) where {T, V, N} + duals .= ForwardDiff.Dual{T, V, N}.(x, Ref(seed)) + return duals +end + +function ForwardDiff.seed!( + duals::GPUArraysCore.AbstractGPUArray{ForwardDiff.Dual{T, V, N}}, x, + seeds::NTuple{N, ForwardDiff.Partials{N, V}} + ) where {T, V, N} + dual_inds = 1:N + duals[dual_inds] .= ForwardDiff.Dual{T, V, N}.(view(x, dual_inds), seeds) + return duals +end + +function ForwardDiff.seed!( + duals::GPUArraysCore.AbstractGPUArray{ForwardDiff.Dual{T, V, N}}, x, index, + seed::ForwardDiff.Partials{N, V} = zero(ForwardDiff.Partials{N, V}) + ) where {T, V, N} + dual_inds = index:length(duals) + duals[dual_inds] .= ForwardDiff.Dual{T, V, N}.(view(x, dual_inds), Ref(seed)) + return duals +end + +function ForwardDiff.seed!( + duals::GPUArraysCore.AbstractGPUArray{ForwardDiff.Dual{T, V, N}}, x, index, + seeds::NTuple{N, ForwardDiff.Partials{N, V}}, chunksize = N + ) where {T, V, N} + offset = index - 1 + dual_inds = (1 + offset):(offset + chunksize) + duals[dual_inds] .= ForwardDiff.Dual{T, V, N}.(view(x, dual_inds), seeds[1:chunksize]) + return duals +end + ### ------ Type Piracy Ends ------ ### # NOTE: LinearSolve.defaultalg for AbstractSciMLOperator + AbstractGPUArray was upstreamed diff --git a/test/Project.toml b/test/Project.toml index 09d39d9..7c51ea3 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -8,6 +8,7 @@ DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7" FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +JLArrays = "27aeb0d3-9eb9-45fb-866b-73c2ecf80fcb" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" Optimisers = "3bd65402-5787-11e9-1adc-39752487f4e2" @@ -35,6 +36,7 @@ DiffEqBase = "6, 7.5" ExplicitImports = "1" FiniteDifferences = "0.12" ForwardDiff = "0.10, 1" +JLArrays = "0.1, 0.2, 0.3" NonlinearSolve = "1, 2, 3, 4" Optimisers = "0.3, 0.4" Optimization = "3, 4, 5.6" @@ -45,6 +47,6 @@ SafeTestsets = "0.1" SciMLSensitivity = "7" SimpleNonlinearSolve = "0.1, 1, 2" StableRNGs = "1" -Statistics = "1.11.1" +Statistics = "1.10, 1.11" SteadyStateDiffEq = "2" Zygote = "0.6, 0.7" diff --git a/test/core/cpu/forwarddiff_seed.jl b/test/core/cpu/forwarddiff_seed.jl new file mode 100644 index 0000000..b722cb9 --- /dev/null +++ b/test/core/cpu/forwarddiff_seed.jl @@ -0,0 +1,22 @@ +using ComplementaritySolve, ForwardDiff, JLArrays, Test + +# ForwardDiff 1.x seeds dual arrays with scalar indexing, which errors on GPU +# arrays. ComplementaritySolve pirates broadcast-based `seed!` methods for +# `AbstractGPUArray`; JLArrays emulates GPU array semantics on the CPU so the +# overloads can be exercised without a GPU. +JLArrays.allowscalar(false) + +@testset "ForwardDiff seeding on GPU arrays" begin + f(x) = x .^ 2 .+ 2 .* x + + @testset "vector mode (length $(n))" for n in (4, 8) + x = collect(Float64, 1:n) + @test Array(ForwardDiff.jacobian(f, JLArray(x))) == ForwardDiff.jacobian(f, x) + end + + # lengths above the default chunk size exercise the chunked `seed!` methods + @testset "chunk mode (length $(n))" for n in (16, 20, 27) + x = collect(Float64, 1:n) + @test Array(ForwardDiff.jacobian(f, JLArray(x))) == ForwardDiff.jacobian(f, x) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 383c36d..9b31163 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -25,6 +25,9 @@ end @safetestset "Mixed Complementarity Problems" begin include("core/cpu/mcp.jl") end + @safetestset "ForwardDiff Seeding on GPU Arrays" begin + include("core/cpu/forwarddiff_seed.jl") + end end @testif BACKEND_GROUP "CUDA" begin From 4ba8a52553b40ee2be23c03f4248e7bb7b8639eb Mon Sep 17 00:00:00 2001 From: "Chris Rackauckas (Claude)" Date: Wed, 10 Jun 2026 00:35:41 -0400 Subject: [PATCH 4/7] Retrigger CI: CUDA Core job hit out-of-GPU-memory at context creation (runner infrastructure issue) Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Fable 5 From c3a8fa960ebef5377cc6c534063827a978f48cb1 Mon Sep 17 00:00:00 2001 From: "Chris Rackauckas (Claude)" Date: Wed, 10 Jun 2026 04:56:52 -0400 Subject: [PATCH 5/7] Remove ForwardDiff seed! type piracy, document the GPU limitation instead Per review feedback on #63: drop the pirated broadcast-based seed! overloads for AbstractGPUArray (and the JLArrays regression test that exercised them) and document the limitation instead, in the README CUDA footnote and next to DEFAULT_NLSOLVER. With ForwardDiff >= 1, ForwardDiff-based jacobians seed dual arrays with scalar setindex! loops, so solvers that compute jacobians through ForwardDiff error on CUDA arrays until that is fixed upstream in ForwardDiff. The CUDA CI groups will fail with ForwardDiff 1.x as a consequence. Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Fable 5 --- README.md | 6 +++- src/ComplementaritySolve.jl | 49 +++++-------------------------- test/Project.toml | 2 -- test/core/cpu/forwarddiff_seed.jl | 22 -------------- test/runtests.jl | 3 -- 5 files changed, 13 insertions(+), 69 deletions(-) delete mode 100644 test/core/cpu/forwarddiff_seed.jl diff --git a/README.md b/README.md index 794ea78..094e5ab 100644 --- a/README.md +++ b/README.md @@ -42,7 +42,11 @@ All `LCP`s, `MLCP`s, and `NCP`s can be converted to `MCP`s, and these solvers ca | MixedComplementarityAdjoint | MCP | | :heavy_check_mark: | :heavy_check_mark: | | [^1]: Solvers internally using `NonlinearSolve.jl` need to use a CUDA compatible solver -(like `SimpleNewtonRaphson()`). +(like `SimpleNewtonRaphson()`). Additionally, with ForwardDiff 1.x these solvers +currently error on CUDA arrays: ForwardDiff >= 1 seeds dual arrays with scalar +`setindex!` loops (ForwardDiff 0.10 used broadcast), which is disallowed on GPU +arrays, so ForwardDiff-based jacobians (the default AD) fail until this is fixed +upstream in ForwardDiff. ## Usage diff --git a/src/ComplementaritySolve.jl b/src/ComplementaritySolve.jl index c3b740d..e5b3905 100644 --- a/src/ComplementaritySolve.jl +++ b/src/ComplementaritySolve.jl @@ -44,53 +44,20 @@ const AV = AbstractVector const AM = AbstractMatrix const AA3 = AbstractArray{T, 3} where {T} +# KNOWN LIMITATION: with ForwardDiff >= 1, ForwardDiff-based jacobians (the +# default AD for SimpleNewtonRaphson and friends) seed dual arrays with scalar +# `setindex!` loops (via `structural_eachindex`), which errors on GPU arrays +# with scalar indexing disallowed. ForwardDiff 0.10 used broadcast-based +# seeding and worked on GPU arrays. Until this is fixed upstream in +# ForwardDiff, solvers that compute jacobians via ForwardDiff (e.g. +# NonlinearReformulation with the default nonlinear solver) error on CUDA +# arrays when run with ForwardDiff >= 1. const DEFAULT_NLSOLVER = SimpleNewtonRaphson() ### ----- Type Piracy Starts ----- ### ArrayInterfaceCore.can_setindex(::Type{<:AbstractFill}) = false ArrayInterfaceCore.can_setindex(::Zygote.OneElement) = false -# ForwardDiff 1.x seeds dual arrays with scalar `setindex!` loops (via -# `structural_eachindex`), which errors on GPU arrays with scalar indexing -# disallowed. ForwardDiff 0.10 used broadcast and worked on GPU arrays, so -# restore broadcast-based seeding for them. TODO: upstream to ForwardDiff as a -# GPUArraysCore extension. -function ForwardDiff.seed!( - duals::GPUArraysCore.AbstractGPUArray{ForwardDiff.Dual{T, V, N}}, x, - seed::ForwardDiff.Partials{N, V} = zero(ForwardDiff.Partials{N, V}) - ) where {T, V, N} - duals .= ForwardDiff.Dual{T, V, N}.(x, Ref(seed)) - return duals -end - -function ForwardDiff.seed!( - duals::GPUArraysCore.AbstractGPUArray{ForwardDiff.Dual{T, V, N}}, x, - seeds::NTuple{N, ForwardDiff.Partials{N, V}} - ) where {T, V, N} - dual_inds = 1:N - duals[dual_inds] .= ForwardDiff.Dual{T, V, N}.(view(x, dual_inds), seeds) - return duals -end - -function ForwardDiff.seed!( - duals::GPUArraysCore.AbstractGPUArray{ForwardDiff.Dual{T, V, N}}, x, index, - seed::ForwardDiff.Partials{N, V} = zero(ForwardDiff.Partials{N, V}) - ) where {T, V, N} - dual_inds = index:length(duals) - duals[dual_inds] .= ForwardDiff.Dual{T, V, N}.(view(x, dual_inds), Ref(seed)) - return duals -end - -function ForwardDiff.seed!( - duals::GPUArraysCore.AbstractGPUArray{ForwardDiff.Dual{T, V, N}}, x, index, - seeds::NTuple{N, ForwardDiff.Partials{N, V}}, chunksize = N - ) where {T, V, N} - offset = index - 1 - dual_inds = (1 + offset):(offset + chunksize) - duals[dual_inds] .= ForwardDiff.Dual{T, V, N}.(view(x, dual_inds), seeds[1:chunksize]) - return duals -end - ### ------ Type Piracy Ends ------ ### # NOTE: LinearSolve.defaultalg for AbstractSciMLOperator + AbstractGPUArray was upstreamed diff --git a/test/Project.toml b/test/Project.toml index 7c51ea3..32a3d1b 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -8,7 +8,6 @@ DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7" FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" -JLArrays = "27aeb0d3-9eb9-45fb-866b-73c2ecf80fcb" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" Optimisers = "3bd65402-5787-11e9-1adc-39752487f4e2" @@ -36,7 +35,6 @@ DiffEqBase = "6, 7.5" ExplicitImports = "1" FiniteDifferences = "0.12" ForwardDiff = "0.10, 1" -JLArrays = "0.1, 0.2, 0.3" NonlinearSolve = "1, 2, 3, 4" Optimisers = "0.3, 0.4" Optimization = "3, 4, 5.6" diff --git a/test/core/cpu/forwarddiff_seed.jl b/test/core/cpu/forwarddiff_seed.jl deleted file mode 100644 index b722cb9..0000000 --- a/test/core/cpu/forwarddiff_seed.jl +++ /dev/null @@ -1,22 +0,0 @@ -using ComplementaritySolve, ForwardDiff, JLArrays, Test - -# ForwardDiff 1.x seeds dual arrays with scalar indexing, which errors on GPU -# arrays. ComplementaritySolve pirates broadcast-based `seed!` methods for -# `AbstractGPUArray`; JLArrays emulates GPU array semantics on the CPU so the -# overloads can be exercised without a GPU. -JLArrays.allowscalar(false) - -@testset "ForwardDiff seeding on GPU arrays" begin - f(x) = x .^ 2 .+ 2 .* x - - @testset "vector mode (length $(n))" for n in (4, 8) - x = collect(Float64, 1:n) - @test Array(ForwardDiff.jacobian(f, JLArray(x))) == ForwardDiff.jacobian(f, x) - end - - # lengths above the default chunk size exercise the chunked `seed!` methods - @testset "chunk mode (length $(n))" for n in (16, 20, 27) - x = collect(Float64, 1:n) - @test Array(ForwardDiff.jacobian(f, JLArray(x))) == ForwardDiff.jacobian(f, x) - end -end diff --git a/test/runtests.jl b/test/runtests.jl index 9b31163..383c36d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -25,9 +25,6 @@ end @safetestset "Mixed Complementarity Problems" begin include("core/cpu/mcp.jl") end - @safetestset "ForwardDiff Seeding on GPU Arrays" begin - include("core/cpu/forwarddiff_seed.jl") - end end @testif BACKEND_GROUP "CUDA" begin From b329b5f7f80cf4ddfa66215a7cabc6369d4aaeeb Mon Sep 17 00:00:00 2001 From: "Chris Rackauckas (Claude)" Date: Wed, 10 Jun 2026 05:12:24 -0400 Subject: [PATCH 6/7] Document ForwardDiff GPU seeding breakage in an issue; mark CUDA tests broken Per review feedback: move the ForwardDiff >= 1 GPU limitation documentation out of the README/source and into https://github.com/SciML/ComplementaritySolve.jl/issues/65, and mark the affected CUDA tests with @test_broken so the GPU CI reflects the known upstream breakage instead of erroring. Only InteriorPointMethod (which does not use ForwardDiff jacobians) keeps regular assertions on CUDA; the NonlinearReformulation/BokhovenIterativeAlgorithm solve and adjoint paths are @test_broken until ForwardDiff restores broadcast-compatible seeding. Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Fable 5 --- README.md | 6 +- src/ComplementaritySolve.jl | 8 --- test/core/cuda/lcp.jl | 133 ++++++++++++++++++++++-------------- test/core/cuda/mcp.jl | 20 ++++-- 4 files changed, 97 insertions(+), 70 deletions(-) diff --git a/README.md b/README.md index 094e5ab..794ea78 100644 --- a/README.md +++ b/README.md @@ -42,11 +42,7 @@ All `LCP`s, `MLCP`s, and `NCP`s can be converted to `MCP`s, and these solvers ca | MixedComplementarityAdjoint | MCP | | :heavy_check_mark: | :heavy_check_mark: | | [^1]: Solvers internally using `NonlinearSolve.jl` need to use a CUDA compatible solver -(like `SimpleNewtonRaphson()`). Additionally, with ForwardDiff 1.x these solvers -currently error on CUDA arrays: ForwardDiff >= 1 seeds dual arrays with scalar -`setindex!` loops (ForwardDiff 0.10 used broadcast), which is disallowed on GPU -arrays, so ForwardDiff-based jacobians (the default AD) fail until this is fixed -upstream in ForwardDiff. +(like `SimpleNewtonRaphson()`). ## Usage diff --git a/src/ComplementaritySolve.jl b/src/ComplementaritySolve.jl index e5b3905..3042a3c 100644 --- a/src/ComplementaritySolve.jl +++ b/src/ComplementaritySolve.jl @@ -44,14 +44,6 @@ const AV = AbstractVector const AM = AbstractMatrix const AA3 = AbstractArray{T, 3} where {T} -# KNOWN LIMITATION: with ForwardDiff >= 1, ForwardDiff-based jacobians (the -# default AD for SimpleNewtonRaphson and friends) seed dual arrays with scalar -# `setindex!` loops (via `structural_eachindex`), which errors on GPU arrays -# with scalar indexing disallowed. ForwardDiff 0.10 used broadcast-based -# seeding and worked on GPU arrays. Until this is fixed upstream in -# ForwardDiff, solvers that compute jacobians via ForwardDiff (e.g. -# NonlinearReformulation with the default nonlinear solver) error on CUDA -# arrays when run with ForwardDiff >= 1. const DEFAULT_NLSOLVER = SimpleNewtonRaphson() ### ----- Type Piracy Starts ----- ### diff --git a/test/core/cuda/lcp.jl b/test/core/cuda/lcp.jl index 896f88b..a95b9f9 100644 --- a/test/core/cuda/lcp.jl +++ b/test/core/cuda/lcp.jl @@ -20,12 +20,23 @@ include("../../test_utils.jl") InteriorPointMethod(), NonlinearReformulation(), BokhovenIterativeAlgorithm(), ] - sol = solve(prob, solver) - - u = Array(sol.u) - @test u ≈ [4.0 / 3, 7.0 / 3] rtol = 1.0e-3 - w = Array(A * sol.u .+ q) - @test w ≈ [0.0, 0.0] atol = 1.0e-3 + if solver isa InteriorPointMethod + sol = solve(prob, solver) + + u = Array(sol.u) + @test u ≈ [4.0 / 3, 7.0 / 3] rtol = 1.0e-3 + w = Array(A * sol.u .+ q) + @test w ≈ [0.0, 0.0] atol = 1.0e-3 + else + # ForwardDiff >= 1 scalar-indexes GPU dual arrays during jacobian + # seeding, breaking ForwardDiff-jacobian-based solvers on CUDA. See + # https://github.com/SciML/ComplementaritySolve.jl/issues/65 + @test_broken begin + sol = solve(prob, solver) + isapprox(Array(sol.u), [4.0 / 3, 7.0 / 3]; rtol = 1.0e-3) && + isapprox(Array(A * sol.u .+ q), [0.0, 0.0]; atol = 1.0e-3) + end + end end @testset "Batched Version" begin @@ -35,10 +46,18 @@ include("../../test_utils.jl") BokhovenIterativeAlgorithm(), NonlinearReformulation(), ] - sol = solve(prob, solver; maxiters = 10000) - - @test all(z -> ≈(Array(z), [4.0 / 3, 7.0 / 3]; rtol = 1.0e-3), eachcol(sol.u)) - @test all(z -> ≈(Array(A * z .+ q), [0.0, 0.0]; atol = 1.0e-3), eachcol(sol.u)) + # Broken by ForwardDiff >= 1 GPU dual seeding, see + # https://github.com/SciML/ComplementaritySolve.jl/issues/65 + @test_broken begin + sol = solve(prob, solver; maxiters = 10000) + all( + z -> ≈(Array(z), [4.0 / 3, 7.0 / 3]; rtol = 1.0e-3), + eachcol(sol.u) + ) && all( + z -> ≈(Array(A * z .+ q), [0.0, 0.0]; atol = 1.0e-3), + eachcol(sol.u) + ) + end end end @@ -48,28 +67,35 @@ include("../../test_utils.jl") solver = NonlinearReformulation() for loss_function in (sum, Base.Fix1(sum, abs2)) - ∂A, ∂q = Zygote.gradient(A, q) do A, q - prob = LinearComplementarityProblem{false}(A, q, u0) - sol = solve(prob, solver) - return loss_function(sol.u) + # Broken by ForwardDiff >= 1 GPU dual seeding, see + # https://github.com/SciML/ComplementaritySolve.jl/issues/65 + @test_broken begin + ∂A, ∂q = Zygote.gradient(A, q) do A, q + prob = LinearComplementarityProblem{false}(A, q, u0) + sol = solve(prob, solver) + return loss_function(sol.u) + end + + θ = ComponentArray((; A = Array(A), q = Array(q))) + ∂θ_fd = ForwardDiff.gradient(θ) do θ + prob = LinearComplementarityProblem{false}(θ.A, θ.q, Array(u0)) + sol = solve(prob, solver) + return loss_function(sol.u) + end + + (∂A !== nothing && !iszero(∂A) && size(∂A) == size(A)) && + (∂q !== nothing && !iszero(∂q) && size(∂q) == size(q)) && + isapprox(Array(∂A), ∂θ_fd.A; atol = 1.0e-3, rtol = 1.0e-3) && + isapprox(Array(∂q), ∂θ_fd.q; atol = 1.0e-3, rtol = 1.0e-3) end - θ = ComponentArray((; A = Array(A), q = Array(q))) - ∂θ_fd = ForwardDiff.gradient(θ) do θ - prob = LinearComplementarityProblem{false}(θ.A, θ.q, Array(u0)) - sol = solve(prob, solver) - return loss_function(sol.u) - end - - @test ∂A !== nothing && !iszero(∂A) && size(∂A) == size(A) - @test ∂q !== nothing && !iszero(∂q) && size(∂q) == size(q) - @test Array(∂A) ≈ ∂θ_fd.A atol = 1.0e-3 rtol = 1.0e-3 - @test Array(∂q) ≈ ∂θ_fd.q atol = 1.0e-3 rtol = 1.0e-3 - - @test_nowarn Zygote.gradient(A, q) do M, q - prob = LinearComplementarityProblem{true}(M, q, u0) - sol = solve(prob) - return loss_function(sol.u) + @test_broken begin + Zygote.gradient(A, q) do M, q + prob = LinearComplementarityProblem{true}(M, q, u0) + sol = solve(prob) + return loss_function(sol.u) + end + true end end end @@ -86,31 +112,38 @@ include("../../test_utils.jl") q_ = randn(rng, Float32, szq...) |> cu for loss_function in (sum, Base.Fix1(sum, abs2)) - ∂M, ∂q = Zygote.gradient(M_, q_) do M, q - prob = LinearComplementarityProblem{false}(M, q) - sol = solve(prob) - return loss_function(sol.u) + # Broken by ForwardDiff >= 1 GPU dual seeding, see + # https://github.com/SciML/ComplementaritySolve.jl/issues/65 + @test_broken begin + ∂M, ∂q = Zygote.gradient(M_, q_) do M, q + prob = LinearComplementarityProblem{false}(M, q) + sol = solve(prob) + return loss_function(sol.u) + end + + θ = ComponentArray((; M = Array(M_), q = Array(q_))) + ∂θ_fd = ForwardDiff.gradient(θ) do θ + prob = LinearComplementarityProblem{false}(θ.M, θ.q) + sol = solve(prob) + return loss_function(sol.u) + end + + (∂M !== nothing && size(∂M) == size(M_)) && + (∂q !== nothing && size(∂q) == size(q_)) && + isapprox(Array(∂M), ∂θ_fd.M; atol = 1.0e-2, rtol = 1.0e-2) && + isapprox(Array(∂q), ∂θ_fd.q; atol = 1.0e-2, rtol = 1.0e-2) end - θ = ComponentArray((; M = Array(M_), q = Array(q_))) - ∂θ_fd = ForwardDiff.gradient(θ) do θ - prob = LinearComplementarityProblem{false}(θ.M, θ.q) - sol = solve(prob) - return loss_function(sol.u) - end - - @test ∂M !== nothing && size(∂M) == size(M_) - @test ∂q !== nothing && size(∂q) == size(q_) - @test Array(∂M) ≈ ∂θ_fd.M atol = 1.0e-2 rtol = 1.0e-2 - @test Array(∂q) ≈ ∂θ_fd.q atol = 1.0e-2 rtol = 1.0e-2 - # We can't check for correctness with FwdDiff & FiniteDifferences # for inplace problems since the in-place batched solvers are not as # accurate as the non-inplace ones. - @test_nowarn Zygote.gradient(M_, q_) do M, q - prob = LinearComplementarityProblem{true}(M, q) - sol = solve(prob) - return loss_function(sol.u) + @test_broken begin + Zygote.gradient(M_, q_) do M, q + prob = LinearComplementarityProblem{true}(M, q) + sol = solve(prob) + return loss_function(sol.u) + end + true end end end diff --git a/test/core/cuda/mcp.jl b/test/core/cuda/mcp.jl index 56ee035..8163827 100644 --- a/test/core/cuda/mcp.jl +++ b/test/core/cuda/mcp.jl @@ -35,9 +35,12 @@ rng = StableRNG(0) NonlinearReformulation(:smooth), NonlinearReformulation(:minmax), ) - sol = solve(prob, solver; verbose = false) - - @test sol.u[1:2] ≈ θ atol = 1.0e-4 rtol = 1.0e-4 + # Broken by ForwardDiff >= 1 GPU dual seeding, see + # https://github.com/SciML/ComplementaritySolve.jl/issues/65 + @test_broken begin + sol = solve(prob, solver; verbose = false) + isapprox(sol.u[1:2], θ; atol = 1.0e-4, rtol = 1.0e-4) + end end end end @@ -59,10 +62,13 @@ rng = StableRNG(0) end θ = [2.0, -3.0] |> cu - ∂θ_zygote = only(Zygote.gradient(loss_function, θ)) - ∂θ_forwarddiff = ForwardDiff.gradient(loss_function_cpu, Array(θ)) - - @test Array(∂θ_zygote) ≈ ∂θ_forwarddiff atol = 1.0e-6 rtol = 1.0e-6 + # Broken by ForwardDiff >= 1 GPU dual seeding, see + # https://github.com/SciML/ComplementaritySolve.jl/issues/65 + @test_broken begin + ∂θ_zygote = only(Zygote.gradient(loss_function, θ)) + ∂θ_forwarddiff = ForwardDiff.gradient(loss_function_cpu, Array(θ)) + isapprox(Array(∂θ_zygote), ∂θ_forwarddiff; atol = 1.0e-6, rtol = 1.0e-6) + end end end end From ec0932d3bc979f37d054c443b8f5b2eccf747d60 Mon Sep 17 00:00:00 2001 From: "Chris Rackauckas (Claude)" Date: Wed, 10 Jun 2026 06:15:13 -0400 Subject: [PATCH 7/7] Keep regular assertions for the size(q) = (2, 1) batched CUDA adjoint CI reported Unexpected Pass for this combination: it takes the batched code path that does not go through ForwardDiff jacobian seeding, so it still works on CUDA with ForwardDiff >= 1 and must not be @test_broken. Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Fable 5 --- test/core/cuda/lcp.jl | 59 +++++++++++++++++++++++++++++++++---------- 1 file changed, 45 insertions(+), 14 deletions(-) diff --git a/test/core/cuda/lcp.jl b/test/core/cuda/lcp.jl index a95b9f9..51bf4e7 100644 --- a/test/core/cuda/lcp.jl +++ b/test/core/cuda/lcp.jl @@ -112,9 +112,9 @@ include("../../test_utils.jl") q_ = randn(rng, Float32, szq...) |> cu for loss_function in (sum, Base.Fix1(sum, abs2)) - # Broken by ForwardDiff >= 1 GPU dual seeding, see - # https://github.com/SciML/ComplementaritySolve.jl/issues/65 - @test_broken begin + if szq == (2, 1) + # this combination takes the batched code path that does not + # go through ForwardDiff jacobian seeding, so it works on CUDA ∂M, ∂q = Zygote.gradient(M_, q_) do M, q prob = LinearComplementarityProblem{false}(M, q) sol = solve(prob) @@ -128,22 +128,53 @@ include("../../test_utils.jl") return loss_function(sol.u) end - (∂M !== nothing && size(∂M) == size(M_)) && - (∂q !== nothing && size(∂q) == size(q_)) && - isapprox(Array(∂M), ∂θ_fd.M; atol = 1.0e-2, rtol = 1.0e-2) && - isapprox(Array(∂q), ∂θ_fd.q; atol = 1.0e-2, rtol = 1.0e-2) - end + @test ∂M !== nothing && size(∂M) == size(M_) + @test ∂q !== nothing && size(∂q) == size(q_) + @test Array(∂M) ≈ ∂θ_fd.M atol = 1.0e-2 rtol = 1.0e-2 + @test Array(∂q) ≈ ∂θ_fd.q atol = 1.0e-2 rtol = 1.0e-2 - # We can't check for correctness with FwdDiff & FiniteDifferences - # for inplace problems since the in-place batched solvers are not as - # accurate as the non-inplace ones. - @test_broken begin - Zygote.gradient(M_, q_) do M, q + # We can't check for correctness with FwdDiff & FiniteDifferences + # for inplace problems since the in-place batched solvers are not + # as accurate as the non-inplace ones. + @test_nowarn Zygote.gradient(M_, q_) do M, q prob = LinearComplementarityProblem{true}(M, q) sol = solve(prob) return loss_function(sol.u) end - true + else + # Broken by ForwardDiff >= 1 GPU dual seeding, see + # https://github.com/SciML/ComplementaritySolve.jl/issues/65 + @test_broken begin + ∂M, ∂q = Zygote.gradient(M_, q_) do M, q + prob = LinearComplementarityProblem{false}(M, q) + sol = solve(prob) + return loss_function(sol.u) + end + + θ = ComponentArray((; M = Array(M_), q = Array(q_))) + ∂θ_fd = ForwardDiff.gradient(θ) do θ + prob = LinearComplementarityProblem{false}(θ.M, θ.q) + sol = solve(prob) + return loss_function(sol.u) + end + + (∂M !== nothing && size(∂M) == size(M_)) && + (∂q !== nothing && size(∂q) == size(q_)) && + isapprox(Array(∂M), ∂θ_fd.M; atol = 1.0e-2, rtol = 1.0e-2) && + isapprox(Array(∂q), ∂θ_fd.q; atol = 1.0e-2, rtol = 1.0e-2) + end + + # We can't check for correctness with FwdDiff & FiniteDifferences + # for inplace problems since the in-place batched solvers are not + # as accurate as the non-inplace ones. + @test_broken begin + Zygote.gradient(M_, q_) do M, q + prob = LinearComplementarityProblem{true}(M, q) + sol = solve(prob) + return loss_function(sol.u) + end + true + end end end end