get_differential_vars: check Diagonal.diag, not the full matrix (#3922)

all(!iszero, mm) visits every (i, j) pair, including scalar getindex
on the underlying diagonal vector. For a GPU-backed diagonal mass
matrix (Diagonal{Float64, <:CuVector}), that's disallowed and throws
on every DAE init, which is why every case in test/gpu/dae_tests.jl
fails regardless of solver or jac_prototype. Diagonal structurally
guarantees off-diagonal zeros, so check mm.diag directly.

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,14 @@ 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
+            # `Diagonal` structurally guarantees off-diagonal zeros, so check
+            # only `mm.diag` rather than the generic `all(!iszero, ::AbstractMatrix)`
+            # / `_isdiag` paths below, both of which visit every (i, j) pair via
+            # scalar `getindex` -- disallowed when `mm.diag` is GPU-backed
+            # (e.g. `Diagonal{Float64, <:CuVector}`), unlike a reduction over
+            # `mm.diag` itself, which 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,21 @@ 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
+    # `get_differential_vars` must only ever check `mm.diag`, not the full matrix
+    # via `all(!iszero, mm)` -- that visits every (i, j) pair, which is wrong (both
+    # slow: O(n^2) instead of O(n)) and, for a GPU-backed diagonal
+    # (`Diagonal{Float64, <:CuVector}`), unsafe: it calls `getindex` on the
+    # underlying vector one element at a time, which CUDA disallows by default.
+    # Can't construct a CuVector here (no GPU in CI); this pins the CPU behavior
+    # so a regression back to the whole-matrix check is caught.
+    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

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.