get_differential_vars: check Diagonal.diag, not the full matrix

all(!iszero, mm) visits every (i, j) pair of the mass matrix via
scalar getindex. For a GPU-backed diagonal mass matrix
(Diagonal{Float64, <:CuVector}) that indexing is disallowed by
default, so the check throws whenever it runs on a device array,
and it is O(n^2) for a structurally O(n) question. Diagonal
guarantees off-diagonal zeros, so check mm.diag directly; the
broadcast over the diagonal vector dispatches to a device kernel,
matching what find_algebraic_vars_eqs(::Diagonal) already does.

Adds a guard test: a Diagonal wrapping a vector that errors on
scalar getindex but supports whole-array broadcast, i.e. device
array semantics on the CPU, so a regression back to the
whole-matrix walk is caught without GPU hardware in CI.

The statically resolvable Diagonal branch also fixes inference of
solve for the two Rodas4 + ShampineCollocationInit cases in the
NonlinearSolve mass matrix tests. Their form
@test_broken sol = @inferred solve(...) errors with "expression
evaluated to non-Boolean" once @inferred succeeds, so those two
lines become @test (@inferred ...; true), still erroring if
inference regresses.

Also wires up show_errors(results) in the GPU DAE test entry
point, which was defined but never called, so a real failure
prints its actual exception instead of just a status glyph.

Author: singhharsh1708

lib/OrdinaryDiffEqCore/src/misc_utils.jl

@@ -213,6 +213,12 @@ function get_differential_vars(f, u)
         elseif mm isa SciMLOperators.ScalarOperator
             # λ·I: every variable is algebraic when λ is zero.
             return iszero(mm.val) ? falses(size(u)) : nothing
+        elseif mm isa Diagonal
+            # Check only `mm.diag`: the generic paths below visit every (i, j)
+            # pair via scalar `getindex`, which is disallowed when `mm.diag` is
+            # GPU-backed (e.g. `Diagonal{Float64, <:CuVector}`), while a
+            # broadcast over `mm.diag` dispatches to a proper GPU kernel.
+            return reshape(mm.diag .!= 0, size(u))
         elseif all(!iszero, mm)
             return trues(size(mm, 1))
         elseif !(mm isa SciMLOperators.AbstractSciMLOperator) && _isdiag(mm)

lib/OrdinaryDiffEqCore/test/algebraic_vars_detection_tests.jl

@@ -113,3 +113,37 @@ end
     ff_singular = ODEFunction(f!, mass_matrix = ScalarOperator(0.0))
     @test get_differential_vars(ff_singular, u) == falses(3)
 end
+
+@testset "get_differential_vars: Diagonal mass matrix" begin
+    f!(du, u, p, t) = (du .= -u; nothing)
+    u = zeros(4)
+    ff_mixed = ODEFunction(f!, mass_matrix = Diagonal([1.0, 1.0, 0.0, 0.0]))
+    @test get_differential_vars(ff_mixed, u) == Bool[1, 1, 0, 0]
+    ff_alldiff = ODEFunction(f!, mass_matrix = Diagonal([1.0, 2.0, 3.0, 4.0]))
+    @test get_differential_vars(ff_alldiff, u) == Bool[1, 1, 1, 1]
+    ff_allalg = ODEFunction(f!, mass_matrix = Diagonal(zeros(4)))
+    @test get_differential_vars(ff_allalg, u) == Bool[0, 0, 0, 0]
+end
+
+# Mimics GPU device-array semantics on the CPU: elementwise scalar `getindex`
+# is forbidden (like CUDA's default `allowscalar(false)`), while whole-array
+# broadcast still works. Any code path that walks a `Diagonal` of this via
+# `getindex` -- `all(!iszero, mm)` or a `diag(mm)` copy -- errors, so the test
+# below catches a regression back to the whole-matrix check without needing a
+# GPU. Must be top-level: structs cannot be defined inside a `@testset`.
+struct NoScalarIndexVector{T} <: AbstractVector{T}
+    data::Vector{T}
+end
+Base.size(v::NoScalarIndexVector) = size(v.data)
+function Base.getindex(::NoScalarIndexVector, ::Int)
+    error("scalar getindex disallowed on NoScalarIndexVector")
+end
+Base.Broadcast.broadcastable(v::NoScalarIndexVector) = v.data
+
+@testset "get_differential_vars: no scalar indexing of Diagonal mass matrix" begin
+    f!(du, u, p, t) = (du .= -u; nothing)
+    u = zeros(4)
+    mm = Diagonal(NoScalarIndexVector([1.0, 1.0, 0.0, 0.0]))
+    ff = ODEFunction(f!, mass_matrix = mm)
+    @test get_differential_vars(ff, u) == Bool[1, 1, 0, 0]
+end

lib/OrdinaryDiffEqNonlinearSolve/test/mass_matrix_tests.jl

@@ -382,9 +382,12 @@ prob = ODEProblem(f, x0, tspan)
 foop = ODEFunction{false, SciMLBase.AutoSpecialize}(dynamics, mass_matrix = M)
 proboop = ODEProblem(f, x0, tspan)
 @test_broken sol = @inferred solve(prob, Rosenbrock23(autodiff = adalg))
-@test_broken sol = @inferred solve(prob, Rodas4(autodiff = adalg), initializealg = ShampineCollocationInit())
+# The Rodas4 + ShampineCollocationInit cases infer since get_differential_vars
+# got a statically resolvable Diagonal branch; keep @inferred so an inference
+# regression errors, but @test_broken would error ("non-Boolean") on success.
+@test (@inferred solve(prob, Rodas4(autodiff = adalg), initializealg = ShampineCollocationInit()); true)
 @test_broken sol = @inferred solve(proboop, Rodas5())
-@test_broken sol = @inferred solve(proboop, Rodas4(), initializealg = ShampineCollocationInit())
+@test (@inferred solve(proboop, Rodas4(), initializealg = ShampineCollocationInit()); true)
 
 # InitialFailure with chunksize matching ODE size (> alg vars), #3157
 @testset "Mass Matrix: less alg vars than ODE AD chunksize" begin

test/gpu/dae_tests.jl

@@ -458,6 +458,11 @@ end
 @testset "GPU DAE solver compatibility" begin
     global results = run_all()
     show_results(results)
+    # Print the actual caught exceptions on failure -- `show_results` only prints
+    # a per-case status glyph, which is not enough on its own to diagnose what
+    # broke from a CI log.
+    any(r -> r.status ∈ (:gpu_error, :gpu_mismatch, :gpu_failed), results) &&
+        show_errors(results)
     @testset "$(r.case)" for r in results
         # :cpu_skip = the CPU couldn't solve this case, so there is no reference
         # to compare against — record as skipped, not pass/fail.