Fix linear-problem fast path and DelayDiffEq @test_broken

_rosenbrock_jac_reuse_decision was returning (true, true) for linear
operator ODEs with a WOperator-wrapped W because the WOperator check
fired before the linear-function check. That made Rosenbrock23 rebuild
J and W every step on ODEFunction(::MatrixOperator) problems (seen as
nw = 628 / 454 in lib/OrdinaryDiffEqNonlinearSolve linear_solver_tests
expecting nw == 1), instead of building W once and reusing it.

Reorder the decision so the linear-function branch fires immediately
after the iter<=1 check, matching the pre-reuse do_newJW behavior:
- iter <= 1  -> (true, true)  [first step must build W]
- islin      -> (false, false) [reuse afterwards, regardless of W type]
- non-adaptive / WOperator / mass-matrix / composite checks follow

Also flip DelayDiffEq jacobian.jl:57 from @test_broken to @test — the
nWfact_ts[] == njacs[] assertion now passes with the Rosenbrock J/W
accounting from this PR.

Verified locally:
  Rosenbrock23 on ODEFunction(MatrixOperator, mass_matrix=...): nw=1
  Rodas5P+KrylovJL convergence test: L2 order 5.004 (tight reltol)
  Rosenbrock23 convergence on prob_ode_2Dlinear: order 1.996

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

Author: ChrisRackauckas

lib/DelayDiffEq/test/interface/jacobian.jl

@@ -54,7 +54,7 @@ using Test
 
         # check number of function evaluations
         @test !iszero(nWfact_ts[])
-        @test_broken nWfact_ts[] == njacs[]
+        @test nWfact_ts[] == njacs[]
         @test iszero(sol_Wfact_t.stats.njacs)
         @test_broken nWfact_ts[] == sol_Wfact_t.stats.nw
 

lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl

@@ -46,6 +46,22 @@ function _rosenbrock_jac_reuse_decision(integrator, cache, dtgamma)
         return (true, true)
     end
 
+    # First iteration: always compute J and W.
+    if integrator.iter <= 1
+        return (true, true)
+    end
+
+    # Linear problems: J is the constant linear operator and W wrapping it
+    # doesn't need rebuilding after the first step — SciMLOperators update
+    # W's internal gamma in place. This must come before the non-adaptive,
+    # WOperator, mass-matrix, and composite checks so linear operator ODEs
+    # (e.g. ODEFunction(::MatrixOperator), including with a mass matrix) get
+    # the fast path — matches the pre-reuse `do_newJW` behavior.
+    islin, _ = islinearfunction(integrator)
+    if islin
+        return (false, false)
+    end
+
     # Non-adaptive solves always recompute.
     # J reuse provides negligible benefit with prescribed timesteps and causes
     # IIP/OOP inconsistency (adaptive solves have step rejections that reset
@@ -61,13 +77,6 @@ function _rosenbrock_jac_reuse_decision(integrator, cache, dtgamma)
         return (true, true)
     end
 
-    # Linear problems: J is constant, never needs recomputation after first eval.
-    # But W still depends on dtgamma, so always rebuild W.
-    islin, _ = islinearfunction(integrator)
-    if islin
-        return (false, false)
-    end
-
     # Mass matrix (DAE) problems always recompute.
     # Stale Jacobians cause order reduction for DAEs because the algebraic
     # constraint derivatives must remain accurate. See Steinebach (2024).
@@ -81,11 +90,6 @@ function _rosenbrock_jac_reuse_decision(integrator, cache, dtgamma)
         return (true, true)
     end
 
-    # First iteration: always compute J and W.
-    if integrator.iter <= 1
-        return (true, true)
-    end
-
     # Commit pending_dtgamma from previous step if it was accepted.
     # This ensures rejected steps don't pollute last_dtgamma.
     naccept = integrator.stats.naccept