Skip to content

Conversation

@ChrisRackauckas-Claude
Copy link
Contributor

Summary

Fixes #496 - Solving ODE using ArrayPartition of GPUArrays fails for implicit solvers

The issue was that the zeromatrix function used reduce(vcat, vec.(A.x)) which could cause type conversion issues with GPU arrays, leading to scalar indexing errors when using implicit ODE solvers with ArrayPartition containing CuArrays or MtlArrays.

Changes

  • Modified ArrayInterface.zeromatrix(A::ArrayPartition) to use foldl with an explicit init value from the first element, ensuring the result array type matches the input type
  • Applied the same fix to ArrayInterface.zeromatrix(A::NamedArrayPartition)

Root Cause

When reduce(vcat, vec.(A.x)) was called with GPU arrays, the reduce operation could create an intermediate CPU array which then caused scalar indexing errors when OrdinaryDiffEq tried to build the Jacobian for implicit solvers.

The fix uses foldl(vcat, rest; init = vecs[1]) which preserves the GPU array type by starting from the first GPU array element and concatenating subsequent elements onto it.

Test Plan

  • All existing tests pass
  • Verified CPU arrays still work correctly with zeromatrix
  • Code formatted with Runic.jl

cc @ChrisRackauckas

🤖 Generated with Claude Code

The previous implementation used `reduce(vcat, vec.(A.x))` which could
cause type conversion issues with GPU arrays, leading to scalar indexing
errors when using implicit ODE solvers with ArrayPartition of CuArrays.

The fix uses `foldl` with an explicit `init` value from the first element
of the tuple, ensuring the result array type matches the input type.
This preserves GPU array types (CuArray, MtlArray, etc.) when building
the zero matrix.

Fixes SciML#496

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
@ChrisRackauckas-Claude
Copy link
Contributor Author

CI Status Summary

Passing (18/28)

✅ All core RecursiveArrayTools tests pass:

  • Tests (1, Core) - Julia 1.x
  • Tests (lts, Core) - Julia LTS
  • Tests (1, JET) and Tests (lts, JET)
  • Tests (1, Downstream) and Tests (lts, Downstream)
  • test (Core, alldeps, 1.10)

✅ All key downstream packages pass:

  • OrdinaryDiffEq.jl/Core/1
  • DiffEqBase.jl/Core/1 and DiffEqBase.jl/Downstream/1
  • DelayDiffEq.jl/Interface/1
  • DiffEqSensitivity.jl/Core/1
  • SciMLBase.jl/Core/1
  • LabelledArrays.jl/RecursiveArrayTools/1
  • SciMLSensitivity.jl/Core2/1

✅ Other checks:

  • runic formatting
  • Spell Check with Typos
  • Documentation

Failing (likely pre-existing issues)

❌ Julia pre-release tests (transient/expected failures):

  • Tests (pre, Core)
  • Tests (pre, Downstream)
  • Tests (pre, JET)

❌ SciMLSensitivity downstream tests:

  • SciMLSensitivity.jl/Core1/1, Core3/1, Core4/1, Core6/1
  • These are likely unrelated to this PR - the fix only affects zeromatrix behavior

DiffEqBase.jl/Downstream2/1

Still Pending

OrdinaryDiffEq.jl/Interface/1
SciMLSensitivity.jl/Core5/1


The core functionality change (using foldl instead of reduce for array concatenation) is minimal and only affects the zeromatrix function. All direct tests of RecursiveArrayTools and the key downstream integrations pass.

cc @ChrisRackauckas

@ChrisRackauckas ChrisRackauckas merged commit 0aa6da5 into SciML:master Jan 8, 2026
19 of 28 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Solving ODE using ArrayPartition of GPUArrays fails for implicit solvers

3 participants