"""
    RosenbrockMutableCache <: OrdinaryDiffEqMutableCache

Abstract supertype for the in-place caches of the Rosenbrock (and Rosenbrock-W)
methods. Concrete Rosenbrock caches subtype this; the shared integrator interface
dispatches on it to access the stage buffers, Jacobian/`W` matrices, and
differentiation configs common to the Rosenbrock family. Declared public so
cross-sublibrary references to the Rosenbrock cache hierarchy are recognized as a
supported extension point.
"""
abstract type RosenbrockMutableCache <: OrdinaryDiffEqMutableCache end
abstract type RosenbrockConstantCache <: OrdinaryDiffEqConstantCache end

"""
    JacReuseState{T, J, D, W}

Lightweight mutable state for tracking Jacobian reuse in Rosenbrock-W methods.
W-methods guarantee correctness with a stale Jacobian, so we can skip expensive
Jacobian recomputations when conditions allow it.

`T` is the element type of `dtgamma` at solve time — usually `Float64`, but
`ForwardDiff.Dual` (including nested Duals under `ForwardDiff.hessian`) when
the solve is being differentiated. Callers must construct this with a zero
value of the right type (typically `zero(dt * gamma)`) so the same cache can hold
Dual-valued dtgammas without any `ForwardDiff.value` unwrapping — unconditional
unwrapping is unsafe under nested AD because it collapses a still-active
inner derivative into a primal.

`J`, `D`, and `W` are the concrete types of the cached Jacobian, time
derivative, and factorized W matrix respectively. They are seeded with the
matching real values built in `alg_cache` (rather than `nothing`
placeholders) so the field types are concrete — required for AD frameworks
such as Mooncake whose `_copy_output` does an invariant `typeassert` against
the inferred field types and rejects `Any` slots. The "first use" guard
inside `_rosenbrock_jac_reuse_decision` (`iszero(last_dtgamma)` →
`(true, true)`) ensures these initial seed values are never *read* during
solve: step 1 always recomputes fresh and overwrites them.

Fields:
- `last_dtgamma`: The dtgamma value from the last Jacobian computation
- `pending_dtgamma`: The dtgamma from the current step (committed on accept)
- `last_naccept`: The naccept count at last Jacobian computation
- `max_jac_age`: Maximum number of accepted steps between Jacobian updates
  (populated from `alg.max_jac_age`, default 20)
- `cached_J::J`: Cached Jacobian for OOP reuse
- `cached_dT::D`: Cached time derivative for OOP reuse
- `cached_W::W`: Cached factorized W for OOP reuse
- `last_step_iter`: The integrator.iter at the last Rosenbrock step
- `last_u_length`: The length of u at the last Jacobian computation
"""
mutable struct JacReuseState{T, J, D, W}
    last_dtgamma::T
    pending_dtgamma::T
    last_naccept::Int
    max_jac_age::Int
    cached_J::J
    cached_dT::D
    cached_W::W
    last_step_iter::Int
    last_u_length::Int
end

"""
    _make_jac_reuse_state(dtgamma, max_jac_age, J, dT, W)

Allocate a `JacReuseState` seeded with the real `J`, `dT`, `W` values from
`alg_cache`. Always allocates so cache types stay concrete and `@inferred`
tests pass; when `max_jac_age ≤ 1` the age check in
`_rosenbrock_jac_reuse_decision` triggers every step anyway, so the overhead
is negligible.

`J`, `dT`, `W` are forwarded from the call site so each `cached_*` slot is
concretely typed (no `Union{Nothing, ...}` placeholders). AD frameworks such
as Mooncake require this — `Any` or wrapper-typed slots break `_copy_output`'s
invariant typeassert.
"""
@inline function _make_jac_reuse_state(
        dtgamma::T, max_jac_age::Int, J::JT, dT::DT, W::WT
    ) where {T, JT, DT, WT}
    return JacReuseState{T, JT, DT, WT}(
        dtgamma, dtgamma, 0, max_jac_age, J, dT, W, 0, 0
    )
end

# Fake values since non-FSAL
get_fsalfirstlast(cache::RosenbrockMutableCache, u) = (nothing, nothing)

tabtype(::HybridExplicitImplicitRK) = Tsit5DATableau

################################################################################

# Shampine's Low-order Rosenbrocks

mutable struct RosenbrockCache{
        uType, rateType, tabType, uNoUnitsType, JType, WType, TabType,
        TFType, UFType, F, JCType, GCType, RTolType, A, StepLimiter, StageLimiter, JRType,
    } <:
    RosenbrockMutableCache
    u::uType
    uprev::uType
    dense::Vector{rateType}
    du::rateType
    du1::rateType
    du2::rateType
    dtC::Matrix{tabType}
    dtd::Vector{tabType}
    ks::Vector{rateType}
    fsalfirst::rateType
    fsallast::rateType
    dT::rateType
    J::JType
    W::WType
    tmp::rateType
    atmp::uNoUnitsType
    weight::uNoUnitsType
    tab::TabType
    tf::TFType
    uf::UFType
    linsolve_tmp::rateType
    linsolve::F
    jac_config::JCType
    grad_config::GCType
    reltol::RTolType
    alg::A
    step_limiter!::StepLimiter
    stage_limiter!::StageLimiter
    interp_order::Int
    jac_reuse::JRType
end
function full_cache(c::RosenbrockCache)
    return [
        c.u, c.uprev, c.dense..., c.du, c.du1, c.du2,
        c.ks..., c.fsalfirst, c.fsallast, c.dT, c.tmp, c.atmp, c.weight, c.linsolve_tmp,
    ]
end

struct RosenbrockCombinedConstantCache{TF, UF, Tab, JType, WType, F, AD, JRType} <:
    RosenbrockConstantCache
    tf::TF
    uf::UF
    tab::Tab
    J::JType
    W::WType
    linsolve::F
    autodiff::AD
    interp_order::Int
    jac_reuse::JRType
end

@cache mutable struct Rosenbrock23Cache{
        uType, rateType, uNoUnitsType, JType, WType,
        TabType, TFType, UFType, F, JCType, GCType,
        RTolType, A, AV, StepLimiter, StageLimiter, JRType,
    } <: RosenbrockMutableCache
    u::uType
    uprev::uType
    k₁::rateType
    k₂::rateType
    k₃::rateType
    du1::rateType
    du2::rateType
    f₁::rateType
    fsalfirst::rateType
    fsallast::rateType
    dT::rateType
    J::JType
    W::WType
    tmp::rateType
    atmp::uNoUnitsType
    weight::uNoUnitsType
    tab::TabType
    tf::TFType
    uf::UFType
    linsolve_tmp::rateType
    linsolve::F
    jac_config::JCType
    grad_config::GCType
    reltol::RTolType
    alg::A
    algebraic_vars::AV
    step_limiter!::StepLimiter
    stage_limiter!::StageLimiter
    jac_reuse::JRType
end

@cache mutable struct Rosenbrock32Cache{
        uType, rateType, uNoUnitsType, JType, WType,
        TabType, TFType, UFType, F, JCType, GCType,
        RTolType, A, AV, StepLimiter, StageLimiter, JRType,
    } <: RosenbrockMutableCache
    u::uType
    uprev::uType
    k₁::rateType
    k₂::rateType
    k₃::rateType
    du1::rateType
    du2::rateType
    f₁::rateType
    fsalfirst::rateType
    fsallast::rateType
    dT::rateType
    J::JType
    W::WType
    tmp::rateType
    atmp::uNoUnitsType
    weight::uNoUnitsType
    tab::TabType
    tf::TFType
    uf::UFType
    linsolve_tmp::rateType
    linsolve::F
    jac_config::JCType
    grad_config::GCType
    reltol::RTolType
    alg::A
    algebraic_vars::AV
    step_limiter!::StepLimiter
    stage_limiter!::StageLimiter
    jac_reuse::JRType
end

function alg_cache(
        alg::Rosenbrock23, u, rate_prototype, ::Type{uEltypeNoUnits},
        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
        dt, reltol, p, calck,
        ::Val{true}, verbose
    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
    k₁ = zero(rate_prototype)
    k₂ = zero(rate_prototype)
    k₃ = zero(rate_prototype)
    du1 = zero(rate_prototype)
    du2 = zero(rate_prototype)
    # f₀ = zero(u) fsalfirst
    f₁ = zero(rate_prototype)
    fsalfirst = zero(rate_prototype)
    fsallast = zero(rate_prototype)
    dT = zero(rate_prototype)
    tmp = zero(rate_prototype)
    atmp = similar(u, uEltypeNoUnits)
    recursivefill!(atmp, false)
    weight = similar(u, uEltypeNoUnits)
    recursivefill!(weight, false)
    tab = Rosenbrock23Tableau(constvalue(uBottomEltypeNoUnits))
    tf = TimeGradientWrapper(f, uprev, p)
    uf = UJacobianWrapper(f, t, p)
    linsolve_tmp = zero(rate_prototype)

    grad_config = build_grad_config(alg, f, tf, du1, t)
    jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, du2)

    J, W = build_J_W(alg, u, uprev, p, t, dt, f, jac_config, uEltypeNoUnits, Val(true))

    linprob = LinearProblem(W, _vec(linsolve_tmp), (nothing, u, p, t); u0 = _vec(tmp))
    linsolve = init(
        linprob, wrapprecs(alg.linsolve, W, weight),
        alias = LinearAliasSpecifier(alias_A = true, alias_b = true),
        abstol = reltol, reltol = reltol,
        assumptions = LinearSolve.OperatorAssumptions(true),
        verbose = verbose.linear_verbosity
    )

    algebraic_vars = f.mass_matrix === I ? nothing :
        find_algebraic_vars_eqs(f.mass_matrix)[1]

    return Rosenbrock23Cache(
        u, uprev, k₁, k₂, k₃, du1, du2, f₁,
        fsalfirst, fsallast, dT, J, W, tmp, atmp, weight, tab, tf, uf,
        linsolve_tmp,
        linsolve, jac_config, grad_config, reltol, alg, algebraic_vars, alg.step_limiter!,
        alg.stage_limiter!, _make_jac_reuse_state(
            zero(dt), alg.max_jac_age, J, zero(rate_prototype), W
        )
    )
end

function alg_cache(
        alg::Rosenbrock32, u, rate_prototype, ::Type{uEltypeNoUnits},
        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
        dt, reltol, p, calck,
        ::Val{true}, verbose
    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
    k₁ = zero(rate_prototype)
    k₂ = zero(rate_prototype)
    k₃ = zero(rate_prototype)
    du1 = zero(rate_prototype)
    du2 = zero(rate_prototype)
    # f₀ = zero(u) fsalfirst
    f₁ = zero(rate_prototype)
    fsalfirst = zero(rate_prototype)
    fsallast = zero(rate_prototype)
    dT = zero(rate_prototype)
    tmp = zero(rate_prototype)
    atmp = similar(u, uEltypeNoUnits)
    recursivefill!(atmp, false)
    weight = similar(u, uEltypeNoUnits)
    recursivefill!(weight, false)
    tab = Rosenbrock32Tableau(constvalue(uBottomEltypeNoUnits))

    tf = TimeGradientWrapper(f, uprev, p)
    uf = UJacobianWrapper(f, t, p)
    linsolve_tmp = zero(rate_prototype)

    grad_config = build_grad_config(alg, f, tf, du1, t)
    jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, du2)

    J, W = build_J_W(alg, u, uprev, p, t, dt, f, jac_config, uEltypeNoUnits, Val(true))

    linprob = LinearProblem(W, _vec(linsolve_tmp), (nothing, u, p, t); u0 = _vec(tmp))
    linsolve = init(
        linprob, wrapprecs(alg.linsolve, W, weight),
        alias = LinearAliasSpecifier(alias_A = true, alias_b = true),
        abstol = reltol, reltol = reltol,
        assumptions = LinearSolve.OperatorAssumptions(true),
        verbose = verbose.linear_verbosity
    )

    algebraic_vars = f.mass_matrix === I ? nothing :
        find_algebraic_vars_eqs(f.mass_matrix)[1]

    return Rosenbrock32Cache(
        u, uprev, k₁, k₂, k₃, du1, du2, f₁, fsalfirst, fsallast, dT, J, W,
        tmp, atmp, weight, tab, tf, uf, linsolve_tmp, linsolve, jac_config,
        grad_config, reltol, alg, algebraic_vars, alg.step_limiter!, alg.stage_limiter!,
        _make_jac_reuse_state(
            zero(dt), alg.max_jac_age, J, zero(rate_prototype), W
        )
    )
end

struct Rosenbrock23ConstantCache{T, TF, UF, JType, WType, F, AD, JRType} <:
    RosenbrockConstantCache
    c₃₂::T
    d::T
    tf::TF
    uf::UF
    J::JType
    W::WType
    linsolve::F
    autodiff::AD
    jac_reuse::JRType
end

function alg_cache(
        alg::Rosenbrock23, u, rate_prototype, ::Type{uEltypeNoUnits},
        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
        dt, reltol, p, calck,
        ::Val{false}, verbose
    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
    tf = TimeDerivativeWrapper(f, u, p)
    uf = UDerivativeWrapper(f, t, p)
    J, W = build_J_W(alg, u, uprev, p, t, dt, f, nothing, uEltypeNoUnits, Val(false))
    linprob = nothing #LinearProblem(W,copy(u); u0=copy(u))
    linsolve = nothing #init(linprob,alg.linsolve,alias_A=true,alias_b=true)
    tab = Rosenbrock23Tableau(constvalue(uBottomEltypeNoUnits))
    # Seed JacReuseState with `zero(dt * tab.d)` rather than a
    # `constvalue`-stripped type: this carries the exact type dtgamma has at
    # solve time, including Dual-of-Dual under nested ForwardDiff
    # (e.g. `hessian`), where a Float64 field would reject the assign.
    dtgamma_seed = zero(dt * tab.d)
    return Rosenbrock23ConstantCache(
        tab.c₃₂, tab.d, tf, uf, J, W, linsolve, alg_autodiff(alg),
        _make_jac_reuse_state(
            dtgamma_seed, alg.max_jac_age, J, zero(rate_prototype), W
        )
    )
end

struct Rosenbrock32ConstantCache{T, TF, UF, JType, WType, F, AD, JRType} <:
    RosenbrockConstantCache
    c₃₂::T
    d::T
    tf::TF
    uf::UF
    J::JType
    W::WType
    linsolve::F
    autodiff::AD
    jac_reuse::JRType
end

function alg_cache(
        alg::Rosenbrock32, u, rate_prototype, ::Type{uEltypeNoUnits},
        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
        dt, reltol, p, calck,
        ::Val{false}, verbose
    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
    tf = TimeDerivativeWrapper(f, u, p)
    uf = UDerivativeWrapper(f, t, p)
    J, W = build_J_W(alg, u, uprev, p, t, dt, f, nothing, uEltypeNoUnits, Val(false))
    linprob = nothing #LinearProblem(W,copy(u); u0=copy(u))
    linsolve = nothing #init(linprob,alg.linsolve,alias_A=true,alias_b=true)
    tab = Rosenbrock32Tableau(constvalue(uBottomEltypeNoUnits))
    # See the Rosenbrock23 OOP alg_cache above for why we pass `zero(dt * tab.d)`
    # here rather than a `constvalue`-stripped type.
    dtgamma_seed = zero(dt * tab.d)
    return Rosenbrock32ConstantCache(
        tab.c₃₂, tab.d, tf, uf, J, W, linsolve, alg_autodiff(alg),
        _make_jac_reuse_state(
            dtgamma_seed, alg.max_jac_age, J, zero(rate_prototype), W
        )
    )
end

### Rodas4+ methods and consolidated Rosenbrock methods (using RodasTableau)

# Helper accessors for step_limiter!/stage_limiter! — algorithms that have these fields
# return them directly; algorithms without return trivial_limiter!.
_get_step_limiter(alg) = trivial_limiter!
_get_stage_limiter(alg) = trivial_limiter!
for Alg in (
        :Rosenbrock23, :Rosenbrock32, :ROS3P, :Rodas3, :Rodas3d, :Rodas23W, :Rodas3P,
        :Rodas4, :Rodas42, :Rodas4P, :Rodas4P2, :Rodas4PW, :Rodas5,
        :Rodas5P, :Rodas5Pe, :Rodas5Pr, :Rodas6P,
    )
    @eval _get_step_limiter(alg::$Alg) = alg.step_limiter!
    @eval _get_stage_limiter(alg::$Alg) = alg.stage_limiter!
end

# Tableau type dispatch
tabtype(::Rodas4) = Rodas4Tableau
tabtype(::Rodas42) = Rodas42Tableau
tabtype(::Rodas4P) = Rodas4PTableau
tabtype(::Rodas4P2) = Rodas4P2Tableau
tabtype(::Rodas4PW) = Rodas4PWTableau
tabtype(::Rodas5) = Rodas5Tableau
tabtype(::Rodas5P) = Rodas5PTableau
tabtype(::Rodas5Pr) = Rodas5PTableau
tabtype(::Rodas5Pe) = Rodas5PeTableau
tabtype(::Rodas6P) = Rodas6PTableau

# Consolidated methods: tableau type dispatch
tabtype(::ROS3P) = ROS3PRodasTableau
tabtype(::Rodas3) = Rodas3RodasTableau
tabtype(::Rodas3d) = Rodas3dRodasTableau
tabtype(::Rodas3P) = Rodas3PRodasTableau
tabtype(::Rodas23W) = Rodas23WRodasTableau
tabtype(::ROS2) = ROS2RodasTableau
tabtype(::ROS2PR) = ROS2PRRodasTableau
tabtype(::ROS2S) = ROS2SRodasTableau
tabtype(::ROS3) = ROS3RodasTableau
tabtype(::ROS3PR) = ROS3PRRodasTableau
tabtype(::Scholz4_7) = Scholz4_7RodasTableau
tabtype(::ROS34PW1a) = ROS34PW1aRodasTableau
tabtype(::ROS34PW1b) = ROS34PW1bRodasTableau
tabtype(::ROS34PW2) = ROS34PW2RodasTableau
tabtype(::ROS34PW3) = ROS34PW3RodasTableau
tabtype(::ROS34PRw) = ROS34PRwRodasTableau
tabtype(::ROS3PRL) = ROS3PRLRodasTableau
tabtype(::ROS3PRL2) = ROS3PRL2RodasTableau
tabtype(::ROK4a) = ROK4aRodasTableau
tabtype(::RosShamp4) = RosShamp4RodasTableau
tabtype(::Veldd4) = Veldd4RodasTableau
tabtype(::Velds4) = Velds4RodasTableau
tabtype(::GRK4T) = GRK4TRodasTableau
tabtype(::GRK4A) = GRK4ARodasTableau
tabtype(::Ros4LStab) = Ros4LStabRodasTableau
tabtype(::RosenbrockW6S4OS) = RosenbrockW6S4OSRodasTableau

# Union of all algorithms using RodasTableau-based RosenbrockCache
const RodasTableauAlgorithms = Union{
    Rodas4, Rodas42, Rodas4P, Rodas4P2, Rodas4PW,
    Rodas5, Rodas5P, Rodas5Pe, Rodas5Pr, Rodas6P,
    ROS3P, Rodas3, Rodas3d, Rodas3P, Rodas23W,
    ROS2, ROS2PR, ROS2S, ROS3, ROS3PR, Scholz4_7,
    ROS34PW1a, ROS34PW1b, ROS34PW2, ROS34PW3,
    ROS34PRw, ROS3PRL, ROS3PRL2, ROK4a,
    RosShamp4, Veldd4, Velds4, GRK4T, GRK4A, Ros4LStab,
    RosenbrockW6S4OS,
}

function alg_cache(
        alg::RodasTableauAlgorithms,
        u, rate_prototype, ::Type{uEltypeNoUnits},
        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
        dt, reltol, p, calck,
        ::Val{false}, verbose
    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
    tf = TimeDerivativeWrapper(f, u, p)
    uf = UDerivativeWrapper(f, t, p)
    J, W = build_J_W(alg, u, uprev, p, t, dt, f, nothing, uEltypeNoUnits, Val(false))
    linprob = nothing #LinearProblem(W,copy(u); u0=copy(u))
    linsolve = nothing #init(linprob,alg.linsolve,alias_A=true,alias_b=true)
    tab = tabtype(alg)(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
    H_rows = size(tab.H, 1)
    if H_rows == 0
        interp_order = -1
    elseif alg isa Union{Rodas3P, Rodas23W}
        interp_order = 2
    else
        interp_order = H_rows
    end
    # Seed JacReuseState with `zero(dt * tab.gamma)` so its dtgamma fields carry
    # the full (possibly ForwardDiff.Dual) type dtgamma will have at solve time.
    dtgamma_seed = zero(dt * tab.gamma)
    return RosenbrockCombinedConstantCache(
        tf, uf,
        tab, J, W, linsolve,
        alg_autodiff(alg), interp_order,
        _make_jac_reuse_state(
            dtgamma_seed, alg.max_jac_age, J, zero(rate_prototype), W
        )
    )
end

function alg_cache(
        alg::RodasTableauAlgorithms,
        u, rate_prototype, ::Type{uEltypeNoUnits},
        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
        dt, reltol, p, calck,
        ::Val{true}, verbose
    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
    tab = tabtype(alg)(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
    # Initialize vectors: kshortsize depends on whether H has rows
    H_rows = size(tab.H, 1)
    kshortsize = H_rows > 0 ? H_rows : 2
    # Rodas3P/Rodas23W: H has 3 rows but only 2 are for interpolation;
    # the 3rd row is for interpoldiff error estimation.
    # Methods with empty H use generic Hermite (f₀, f₁), not Rosenbrock
    # dense output coefficients — signal this with interp_order = -1.
    if H_rows == 0
        interp_order = -1
    elseif alg isa Union{Rodas3P, Rodas23W}
        interp_order = 2
    else
        interp_order = kshortsize
    end
    dense = [zero(rate_prototype) for _ in 1:kshortsize]
    ks = [zero(rate_prototype) for _ in 1:size(tab.A, 1)]
    du = zero(rate_prototype)
    du1 = zero(rate_prototype)
    du2 = zero(rate_prototype)

    # Promote t-type for AD
    dtC = zero(tab.C) .* dt
    dtd = zero(tab.d) .* dt

    # Initialize other variables
    fsalfirst = zero(rate_prototype)
    fsallast = zero(rate_prototype)
    dT = zero(rate_prototype)

    # Temporary and helper variables
    tmp = zero(rate_prototype)
    atmp = similar(u, uEltypeNoUnits)
    recursivefill!(atmp, false)
    weight = similar(u, uEltypeNoUnits)
    recursivefill!(weight, false)

    tf = TimeGradientWrapper(f, uprev, p)
    uf = UJacobianWrapper(f, t, p)

    grad_config = build_grad_config(alg, f, tf, du1, t)
    jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, du2)

    J, W = build_J_W(alg, u, uprev, p, t, dt, f, jac_config, uEltypeNoUnits, Val(true))

    linsolve_tmp = zero(rate_prototype)
    linprob = LinearProblem(W, _vec(linsolve_tmp), (nothing, u, p, t); u0 = _vec(tmp))

    linsolve = init(
        linprob, wrapprecs(alg.linsolve, W, weight),
        alias = LinearAliasSpecifier(alias_A = true, alias_b = true),
        abstol = reltol, reltol = reltol,
        assumptions = LinearSolve.OperatorAssumptions(true),
        verbose = verbose.linear_verbosity
    )

    # Return the cache struct with vectors
    return RosenbrockCache(
        u, uprev, dense, du, du1, du2, dtC, dtd, ks, fsalfirst, fsallast,
        dT, J, W, tmp, atmp, weight, tab, tf, uf, linsolve_tmp,
        linsolve, jac_config, grad_config, reltol, alg,
        _get_step_limiter(alg), _get_stage_limiter(alg), interp_order,
        _make_jac_reuse_state(
            zero(dt), alg.max_jac_age, J, zero(rate_prototype), W
        )
    )
end

function get_fsalfirstlast(
        cache::Union{
            Rosenbrock23Cache, Rosenbrock32Cache,
            RosenbrockCache,
        },
        u
    )
    return (cache.fsalfirst, cache.fsallast)
end

################################################################################

### Tsit5DA - hybrid explicit/linear-implicit method for DAEs

struct HybridExplicitImplicitConstantCache{TF, UF, Tab, JType, WType, F, AD} <: RosenbrockConstantCache
    tf::TF
    uf::UF
    tab::Tab
    J::JType
    W::WType
    linsolve::F
    autodiff::AD
    interp_order::Int
end

mutable struct HybridExplicitImplicitCache{
        uType, rateType, uNoUnitsType, JType, WType, TabType,
        TFType, UFType, F, JCType, GCType, RTolType, A,
        StepLimiter, StageLimiter, DVType, AVType,
        GZType, GYType, WZType, FZ,
    } <: RosenbrockMutableCache
    u::uType
    uprev::uType
    dense::Vector{rateType}
    du::rateType
    du1::rateType
    du2::rateType
    ks::Vector{rateType}
    fsalfirst::rateType
    fsallast::rateType
    dT::rateType
    J::JType
    W::WType
    tmp::rateType
    atmp::uNoUnitsType
    weight::uNoUnitsType
    tab::TabType
    tf::TFType
    uf::UFType
    linsolve_tmp::rateType
    linsolve::F
    jac_config::JCType
    grad_config::GCType
    reltol::RTolType
    alg::A
    step_limiter!::StepLimiter
    stage_limiter!::StageLimiter
    interp_order::Int
    # DAE-specific fields
    diff_vars::DVType
    alg_vars::AVType
    g_z::GZType         # n_g x n_g algebraic Jacobian block
    g_y::GYType         # n_g x n_f coupling block
    W_z::WZType         # -gamma * g_z (used for linear solve)
    linsolve_tmp_z::FZ  # n_g-sized RHS for algebraic solve
end
function full_cache(c::HybridExplicitImplicitCache)
    return [
        c.u, c.uprev, c.dense..., c.du, c.du1, c.du2,
        c.ks..., c.fsalfirst, c.fsallast, c.dT, c.tmp, c.atmp, c.weight, c.linsolve_tmp,
    ]
end

function alg_cache(
        alg::HybridExplicitImplicitRK, u, rate_prototype, ::Type{uEltypeNoUnits},
        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
        dt, reltol, p, calck,
        ::Val{false}, verbose
    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
    tf = TimeDerivativeWrapper(f, u, p)
    uf = UDerivativeWrapper(f, t, p)
    J, W = build_J_W(alg, u, uprev, p, t, dt, f, nothing, uEltypeNoUnits, Val(false))
    tab = tabtype(alg)(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
    return HybridExplicitImplicitConstantCache(
        tf, uf, tab, J, W, nothing, alg_autodiff(alg), size(tab.H, 1)
    )
end

function alg_cache(
        alg::HybridExplicitImplicitRK, u, rate_prototype, ::Type{uEltypeNoUnits},
        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
        dt, reltol, p, calck,
        ::Val{true}, verbose
    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
    tab = tabtype(alg)(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
    num_stages = size(tab.A, 1)
    interp_order = size(tab.H, 1)

    # Initialize vectors
    dense = [zero(rate_prototype) for _ in 1:interp_order]
    ks = [zero(rate_prototype) for _ in 1:num_stages]
    du = zero(rate_prototype)
    du1 = zero(rate_prototype)
    du2 = zero(rate_prototype)
    fsalfirst = zero(rate_prototype)
    fsallast = zero(rate_prototype)
    dT = zero(rate_prototype)
    tmp = zero(rate_prototype)
    atmp = similar(u, uEltypeNoUnits)
    recursivefill!(atmp, false)
    weight = similar(u, uEltypeNoUnits)
    recursivefill!(weight, false)
    linsolve_tmp = zero(rate_prototype)

    tf = TimeGradientWrapper(f, uprev, p)
    uf = UJacobianWrapper(f, t, p)

    grad_config = build_grad_config(alg, f, tf, du1, t)
    jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, du2)
    J, W = build_J_W(alg, u, uprev, p, t, dt, f, jac_config, uEltypeNoUnits, Val(true))

    linprob = LinearProblem(W, _vec(linsolve_tmp), (nothing, u, p, t); u0 = _vec(tmp))
    linsolve = init(
        linprob, wrapprecs(alg.linsolve, W, weight),
        alias = LinearAliasSpecifier(alias_A = true, alias_b = true),
        abstol = reltol, reltol = reltol,
        assumptions = LinearSolve.OperatorAssumptions(true),
        verbose = verbose.linear_verbosity
    )

    # Detect algebraic variables from mass matrix
    mass_matrix = f.mass_matrix
    n = length(u)
    if mass_matrix === I
        diff_vars = collect(1:n)
        alg_vars = Int[]
    else
        diff_vars = findall(i -> mass_matrix[i, i] != 0, 1:n)
        alg_vars = findall(i -> mass_matrix[i, i] == 0, 1:n)
    end
    n_g = length(alg_vars)
    n_f = length(diff_vars)
    g_z = zeros(eltype(u), n_g, n_g)
    g_y = zeros(eltype(u), n_g, n_f)
    W_z = zeros(eltype(u), n_g, n_g)
    linsolve_tmp_z = zeros(eltype(u), n_g)

    return HybridExplicitImplicitCache(
        u, uprev, dense, du, du1, du2, ks,
        fsalfirst, fsallast, dT, J, W, tmp, atmp, weight, tab, tf, uf,
        linsolve_tmp, linsolve, jac_config, grad_config, reltol, alg,
        alg.step_limiter!, alg.stage_limiter!, interp_order,
        diff_vars, alg_vars, g_z, g_y, W_z, linsolve_tmp_z
    )
end
