Fix GPU compat of sparse dae solvers

[hexaeder]

This PR attempts to make (most) of the sparse DAE solvers GPU compatible. The main changes are:

• special handling of Diagonal mass matrices (improves performance by order of magnitude for certain dae init calls)
• moves the find_algebraic_vars_eqs function to OrdinaryDiffEqCore
• fixing of (mostly scalar indexing) problems in the init or step code of several methods

This PR needs JuliaGPU/CUDA.jl#3032 to slice into the sparse jacobians.

The tests in the simple_dae.jl script in the GPU tests. Currently, those are more like a debug script. I define a custom test-system and execute it with all kinds of algorithms. The tests compare:

• the CPU solution of the different solvers against the "known good" CPU solution of Rodas5P
• compare the GPU solution of the different solvers to their CPU counterpart
• do this all for Diagonal mass matrices (non-diagonal are much harder to support)
• the system is testes without sparcity pattern, with sparcity pattern as cuda CSR and with sparsity pattern as coda CSC

I would like input by the maintainers on how the tests should look like. Currently I see several problems:

• it takes a long time since (too) many algorithms are tested
• some algos seem to perform badly on my toymodel, resulting in large discrepancy between the CPU solutions of one alg vs another alg
• some algos perform weirdly with CSC vs CSR jacobian prototoypes. CSC is always worse. The CUDA default is CSC, however CUDSS only supports LU factorization on CSR. So the tests currently compare CSC to CPU with Krylov and CSR to CPU with default linsolve.

Algorithm overview

Fully working

ROS2, ROS3, ROS3PRL, ROS3PRL2, Rodas3, Scholz4_7, ROS34PW3, RosShamp4, Veldd4, Velds4, GRK4T, GRK4A,
Ros4LStab, Rodas4, Rodas42, Rodas4P, Rodas4P2, ROK4a, Rodas5, Rodas5P, Rodas5Pe, Rodas5Pr, Rodas6P,
ImplicitEuler, SDIRK2, ABDF2, QNDF1, QBDF1, Trapezoid

CSC path with elevated errors (Krylov vs direct LU accuracy)

These pass but need relaxed csc_tol. The dense and CSR jacobian paths are fine:
Rosenbrock23, ROS2PR, ROS2S, ROS34PW2, ROS34PRw, Cash4, Hairer4, Hairer42, QNDF2, QBDF2

Working, but solver is a poor fit for this DAE problem

GPU correctly reproduces the CPU result, but the solver itself diverges from the Rodas5P reference:
ROS3PR, ROS3P, ROS34PW1a, ROS34PW1b, QBDF1

Special case: Rosenbrock32/CSC: GPU and CPU Krylov agree, but the combination gives catastrophic
errors vs the reference. Effectively unusable with CSC.

Not yet working (require code changes)

The following algorithms fail and I consider them out of scope of this PR since they'd require substantial work

• Rodas23W, Rodas3P Scalar indexing in calculate_interpoldiff!
• QNDF, QBDF DeviceMemory error in LinAlg (some mul! call within the step)
• FBDF Scalar indexing in reinitFBDF!
• RadauIIA3/5/9, AdaptiveRadau lots of problems

Checklist

• Appropriate tests were added
• Any code changes were done in a way that does not break public API
• All documentation related to code changes were updated
• The new code follows the
  contributor guidelines, in particular the SciML Style Guide and
  COLPRAC.
• Any new documentation only uses public API

Additional context

Add any other context about the problem here.

[ChrisRackauckas]

This should get rebased for the new GPU testing infra

[hexaeder]

I rebased. Still waiting for the merge in CUDA though. Is there anything else to do for the "new GPU testing infra" or did you just mean the different ENV variable based testset selection?

[ChrisRackauckas]

Rebase, I had to spend a bit of time fixing the machines. It should all be good now though.

[hexaeder]

CUDA PR was merged. Now we still need:

• new CUDA tag
• merge and release of CUDSS with CUDA 6 compat (compat with CUDA@6.0.0 exanauts/CUDSS.jl#141)
• tag of ArrayInterface (CUDA@6 compat allready merged there)
• PR to fix LinearSolveCUDSSExt + tag once the CUDSS release is there.

Even with those fixes locally applied I see new problems with this PR (forward diff incompatible f? weird). Will investigate parallel to the listed todos...

[ChrisRackauckas]

ArrayInterface tagged.

[ChrisRackauckas-Claude]

Heads up: the GPU tests this PR adds/edits live in the root test/gpu/ directory, which no longer runs in CI. The root GPU lane was removed during the 2026-03-18 sublibrary split (f803bdf53 deleted the group: GPU block from the root CI.yml), and the centralized-matrix migration (#3730) reproduced that no-GPU matrix. So test/runtests.jl has no GPU group and nothing includes test/gpu/*.jl — as written, dae_tests.jl here will not execute.

I opened #3813 to restore a top-level [GPU] group (self-hosted gpu runner) with a gpu_group() that globs test/gpu/*.jl. Once both land, the DAE GPU tests here will run automatically. Two options for sequencing:

1. Merge Restore top-level GPU test lane; drop vestigial AlgConvergence_II group #3813 first, then rebase this PR — its dae_tests.jl is picked up by the glob with no further wiring.
2. Or fold the Restore top-level GPU test lane; drop vestigial AlgConvergence_II group #3813 wiring into this PR.

Either way, since #3813 turns the lane back on, please double-check the test/gpu/ files are current (several had gone stale/diverged from their sublibrary copies while the lane was off).

[ChrisRackauckas]

Rebase onto latest master for #3813

[hexaeder]

Full testset runs through now, other GPU tests are failing though. My debuggin tests took ~20 Minutes (mainly compiling of many different algs). Is that fine for CI or should i thin it down to 1-2 solvers per class?

[hexaeder]

A few things to note: this PR moves find_algebraic_vars_eqs into OrdinaryDiffEqCore. I'm not sure about the release strategy, but it means the subsolver packages that use it (OrdinaryDiffEqNonlinearSolve, OrdinaryDiffEqBDF, OrdinaryDiffEqRosenbrock) need to bump their compat to the new OrdinaryDiffEqCore.

That requirement is what led me to add [sources] blocks to the ODEInterfaceRegression test env. I'm not sure that's in the spirit of the tests though, maybe those should stay red to signal during the release process that the compat bound needs bumping?

There are also a few failing CIs, all of which look unrelated / pre-existing:

• test/gpu group: a pre-existing failure from an upstream GPUArrays bug (csr_type(::CuSparseMatrixCSC) vs. csr_type(::Type{<:CuSparseMatrixCSC})).
• sublibrary-ci / lib/OrdinaryDiffEqLinear [GPU]: the same GPUArrays bug.
• sublibrary-ci / lib/OrdinaryDiffEqNonlinearSolve: pre-existing -- the Dependent Mass Matrix Tests OOP case; fails identically on unrelated PRs
• StochasticDiffEq [QA]: pre-existing Aqua/JET failures; fail identically on unrelated PRs
• downgrade-sublibraries / *: red on master
• Sublibrary CI / sublibrary-ci / lib/OrdinaryDiffEqNonlinearSolve [ModelingToolkit] -- preexisting fail

[ChrisRackauckas]

This needs a rebase for the explicit imports QA changes. I can handle it on Friday if you don't get around to it.

[ChrisRackauckas]

The test seemed to work, it just needs the new docker so the machine won't segfault, but I'll merge given it was working.