Fix build_J_W cache init for structured jac_prototypes

[ChrisRackauckas-Claude]

Note

This PR should be ignored until reviewed by @ChrisRackauckas.

Problem

lib/DiffEqBase's [Downstream] test group fails on master in test/downstream/default_linsolve_structure.jl:

ArgumentError: array of type Diagonal{Float64, Vector{Float64}} and size (2, 2) can
not be filled with 1.0, since some of its entries are constrained.

thrown from build_J_W at lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl when solving an ODEProblem whose jac_prototype is a structured matrix (Diagonal, Tridiagonal, ...) with a default (factorization) linsolve, e.g.

fun = ODEFunction(f; jac = jac, jac_prototype = Diagonal(zeros(2)))
solve(ODEProblem(fun, ones(2), (1.0, 10.0)), Rosenbrock23())

Root cause

Bisected to d8cbb9e ("Initialize sparse Jacobian caches with stored values", #3833). That commit initializes the J/W caches allocated in build_J_W with ones so factorizations never see undef garbage:

if is_sparse(J)
    set_all_nzval!(J, one(eltype(J)))
    set_all_nzval!(W, one(eltype(W)))
else
    fill!(J, one(eltype(J)))
    fill!(W, one(eltype(W)))
end

But LinearAlgebra's fill! on banded structured matrices (Diagonal, Bidiagonal, Tridiagonal, SymTridiagonal) throws for any nonzero value because the off-band entries are constrained to zero. So every structured (non-sparse, non-dense) jac_prototype on the default-linsolve path now errors at cache construction. Verified by running the reproducer on d8cbb9e~1 with the same current dependency versions: it passes there and fails on master, so this is a code regression from #3833, not registry drift (the failure only surfaced later because the monorepo's own CI lanes for #3833 did not run the DiffEqBase Downstream group).

Fix

Add a fill_stored! helper that writes only the stored bands of the LinearAlgebra structured matrix types and falls back to fill! otherwise, and use it in build_J_W's non-sparse initialization. Semantics for dense and sparse prototypes are unchanged.

Tests

• New fill_stored! unit tests and a structured-jac_prototype init/solve regression test in lib/OrdinaryDiffEqDifferentiation/test/nzval_helpers_tests.jl (asserts the J/W caches are initialized with stored ones, mirroring the Initialize sparse Jacobian caches with stored values #3833 testset, and that the solve completes).
• The previously failing downstream file lib/DiffEqBase/test/downstream/default_linsolve_structure.jl passes locally with this change (output in the PR conversation).

This resolves the red lib/DiffEqBase [Downstream] lane seen on #3692 (which merely surfaced the pre-existing master failure).

🤖 Generated with Claude Code

[ChrisRackauckas-Claude]

Local verification (Julia 1.12, current registry, SciMLBase 3.36.0):

Bisect evidence — reproducer on d8cbb9e44~1 (b797bac) with identical dep versions:

PRE3833 RESULT: [10.0, 10.0]  nsteps=46

Same reproducer on unmodified master: ArgumentError: array of type Diagonal{Float64, Vector{Float64}} ... cannot be filled with 1.0 from build_J_W at derivative_utils.jl:1210.

With this PR — reproducer:

Diagonal => [10.0, 10.0] nsteps=46
Tridiagonal => [10.0, 10.0] nsteps=46

Previously failing downstream file lib/DiffEqBase/test/downstream/default_linsolve_structure.jl:

Test Summary:                   | Pass  Broken  Total   Time
Default linsolve with structure |    8       2     10  43.1s

ODEDIFFEQ_TEST_GROUP=Core julia --project -e 'using Pkg; Pkg.test()' (includes the new fill_stored! unit tests and the structured jac_prototype init/solve regression testset): passed, exit 0.

ODEDIFFEQ_TEST_GROUP=QA: Testing OrdinaryDiffEqDifferentiation tests passed, exit 0.

[ChrisRackauckas-Claude]

Status check from an independent triage run today (2026-08-08, master c11ff74bb): this is still live, and it is now breaking a downstream repo's CI.

Three open drafts fix the same fill! in build_J_W and are waiting on review: this one, #3913, and #4072. They need consolidating down to one.

Downstream impact. The Downstream test group of JuliaDiff/FiniteDiff.jl fails on FiniteDiff master because of this: test/downstream/ordinarydiffeq_tridiagonal_solve.jl errors at cache construction with OrdinaryDiffEqDifferentiation v3.7.0. Standalone reproducer (no FiniteDiff, default AutoForwardDiff/AutoSparse — the AD backend is irrelevant, this fires before any Jacobian is computed):

using OrdinaryDiffEq, OrdinaryDiffEqRosenbrock, LinearAlgebra

const n = 10
f(du, u, p, t) = (du .= 0;
    du[2:(end - 1)] .= p[1] .* (u[3:end] .- 2 .* u[2:(end - 1)] .+ u[1:(end - 2)]);
    nothing)

u0 = sin.(range(0, 1, length = n))
jp = Tridiagonal(similar(u0, n - 1), similar(u0), similar(u0, n - 1))
prob = ODEProblem(ODEFunction(f; jac_prototype = jp), u0, (0.0, 1.0), [0.42])
solve(prob, Rodas4P(), saveat = 0.1)

ERROR: ArgumentError: array of type Tridiagonal{Float64, Vector{Float64}} and size (10, 10) can
    not be filled with 1.0, since some of its entries are constrained.
  [1] fill!  @ .../LinearAlgebra/src/special.jl:358 [inlined]
  [2] build_J_W  @ .../OrdinaryDiffEqDifferentiation/G755Y/src/derivative_utils.jl:1242

Regression range confirms this PR's bisect: git log -S 'fill!(J, one(eltype(J)))' gives only d8cbb9e44 (#3833, 2026-07-08), first released in OrdinaryDiffEqDifferentiation v3.2.3 (#3838).

One measured difference between the three PRs, for whoever picks. BandedMatrix prototypes hit the same fill! and are also broken today:

ArrayInterface.has_sparsestruct(jp) = true
ArrayInterface.isstructured(jp) = true
fill! on BandedMatrix threw: BandError: attempt to access BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}} with bandwidths (1, 1) at band 7
solve threw: BandError: ... at band 7

This PR's fill_stored! keeps fill! as the fallback, so a BandedMatrix prototype still throws BandError; #3913 and #4072 dispatch through ArrayInterface.has_sparsestruct / isstructured, which are true for BandedMatrix, so they cover it. Caveat: with the fill! removed for banded, that path then fails on a separate pre-existing lu! ambiguity between ArrayLayouts and LinearAlgebra, so BandedMatrix is not a working configuration either way — it just argues for the ArrayInterface-based dispatch of #3913/#4072 over enumerating the LinearAlgebra band types.

Independent verification of an equivalent patch (band-writing fill_stored!, same as this PR) on today's master, Julia 1.12.6, lib/OrdinaryDiffEqDifferentiation:

• reproducer above: ArgumentError before, ReturnCode.Success after;
• FiniteDiff's ordinarydiffeq_tridiagonal_solve.jl: errors before, passes after;
• GROUP=Core: all 12 testsets green (Testing OrdinaryDiffEqDifferentiation tests passed), including a new solve-level Tridiagonal/Diagonal testset;
• GROUP=QA: JET 1/1, Aqua 18/18, passed.

No fourth PR opened. Suggest merging whichever of the three you prefer (my read: #4072 or #3913 for the ArrayInterface dispatch) and closing the other two, then cutting a patch release of OrdinaryDiffEqDifferentiation — the fix does not reach users, or FiniteDiff's CI, without one.