diff --git a/lib/OrdinaryDiffEqCore/src/misc_utils.jl b/lib/OrdinaryDiffEqCore/src/misc_utils.jl
index ed84e22c06..64e18db9f7 100644
--- a/lib/OrdinaryDiffEqCore/src/misc_utils.jl
+++ b/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)
diff --git a/lib/OrdinaryDiffEqCore/test/algebraic_vars_detection_tests.jl b/lib/OrdinaryDiffEqCore/test/algebraic_vars_detection_tests.jl
index 6685c7c067..671539d5be 100644
--- a/lib/OrdinaryDiffEqCore/test/algebraic_vars_detection_tests.jl
+++ b/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
diff --git a/lib/OrdinaryDiffEqNonlinearSolve/test/mass_matrix_tests.jl b/lib/OrdinaryDiffEqNonlinearSolve/test/mass_matrix_tests.jl
index fbd06cd3fb..71c48ef406 100644
--- a/lib/OrdinaryDiffEqNonlinearSolve/test/mass_matrix_tests.jl
+++ b/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
diff --git a/test/gpu/dae_tests.jl b/test/gpu/dae_tests.jl
index 21fba7e85a..ad4cb57e22 100644
--- a/test/gpu/dae_tests.jl
+++ b/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.