calc_W! never marks an operator Jacobian as updated, so a caching linear solver cannot see J move

[ChrisRackauckas-Claude]

In calc_W! (lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl), the WOperator branch guards its bookkeeping behind

if W.J !== nothing && !(W.J isa AbstractSciMLOperator)
    ...
    new_jac && mark_jacobian_updated!(W)
    ...
end

So when W.J is an operator — the islin path, and any jac_prototype isa AbstractSciMLOperator — mark_jacobian_updated! is never called, while the update_coefficients!(W; gamma = dtgamma) above it moves the operators numbers in place. The object identity is unchanged and jac_stale stays false, so a linear solver caching a factorization of J across steps has no way to learn that J moved.

Why it matters

Any solver that caches something derived from J and re-uses it across gamma changes will silently answer with the previous Jacobian. Demonstrated at the LinearSolve level with LHLFactorization, which caches a Hessenberg reduction of J:

step 1 relres = 2.46e-16
step 2 relres = 1.34       <- J moved in place; wrong answer, nothing raised
   (vs OLD J) = 1.07

No exception, just a wrong number.

Current mitigation

SciML/LinearSolve.jl#1215 (merged) makes defaultalg decline the split form whenever W.J is an operator, so LHLFactorization is never selected on this path and the wrong answer above cannot occur today. That is a workaround on the consumer side, not a fix: it also means the islin path — a linear ODE with a dense operator, exactly the shape where reducing J once and re-shifting per step pays best — cannot use the O(n²) re-shift at all.

Suggested fix

Call mark_jacobian_updated!(W) on the operator path too, whenever update_coefficients! may have moved J. Once an operator J reliably reports staleness, _lhl_defaultable in LinearSolve can be widened back to accept it, and the islin path becomes the best case for the reduction rather than an excluded one.

Related: #4266, SciML/LinearSolve.jl#1215