v3.85 default sparse nonstructural-zeros reduction regresses BoundaryValueDiffEq: precompile ArgumentError (3.85.0) and silent Stalled BVP solves (3.85.1)

[ChrisRackauckas-Claude]

Summary

LinearSolve v3.85.0/v3.85.1 (the persistent_nonstructural_zeros operator assumption + default sparse persistent-zero reduction, #1017, plus the Auto fall-back #1020) regresses downstream BoundaryValueDiffEq.jl:

• v3.85.0 (persistent_nonstructural_zeros operator assumption + default sparse persistent-zero reduction #1017, cdb6f77): hard ArgumentError during precompilation of BoundaryValueDiffEqShooting / BoundaryValueDiffEq (Float64 sparse Jacobians in the precompile workload):

  ArgumentError: nonstructural_zeros reduction requires a constant stored sparsity pattern across solves (stored nnz changed)
    [1] reduce_operand!(red::LinearSolveSparseArraysExt.SparseReduction{Float64, Int64}, A::SparseArrays.SparseMatrixCSC{Float64, Int64})
      @ LinearSolveSparseArraysExt ext/LinearSolveSparseArraysExt.jl:1091
    [2] macro expansion @ src/default.jl:971
    [3] solve!(::LinearSolve.LinearCache{SparseMatrixCSC{Float64, Int64}, ...}, ::LinearSolve.DefaultLinearSolver; ...)
  

  This is what turned BoundaryValueDiffEq master CI red on 2026-06-08 (run on commit 074eeb8e).

• v3.85.1 (Auto-fall-back nonstructural_zeros reduction on stored-pattern changes #1020): the hard error is gone (Auto falls back to per-solve dropzeros), but BVP solves that converged with v3.84.0 now silently return ReturnCode.Stalled:

  • RadauIIa5() on a BigFloat two-point BVP: Success on v3.84.0 → Stalled on v3.85.1.
  • MIRKN4() on a BigFloat second-order BVP: threw on v3.84.0 (known-broken test) → now returns a Stalled solution.

  This is the currently-failing Misc test group on BoundaryValueDiffEq master CI (all of lts/1/pre), test/misc/bigfloat_test.jl.

Reproduction (Julia 1.12.6)

using Pkg; Pkg.activate(temp = true)
Pkg.add(["BoundaryValueDiffEq", "Sparspak", "SciMLBase"])  # resolves LinearSolve v3.85.1
using BoundaryValueDiffEq, Sparspak, SciMLBase

f!(du, u, p, t) = (du[1] = u[2]; du[2] = u[1])
bca!(ra, ua, p) = (ra[1] = ua[1] - 1)
bcb!(rb, ub, p) = (rb[1] = ub[1])
fn = BVPFunction(f!, (bca!, bcb!), bcresid_prototype = (zeros(1), zeros(1)), twopoint = Val(true))
prob = BVProblem(fn, BigFloat.([1.0, 0.0]), (0.0, 1.0))
sol = solve(prob, RadauIIa5(), dt = 0.05)
@show sol.retcode   # Stalled on LinearSolve v3.85.1, Success on v3.84.0

Version bisect (BoundaryValueDiffEq master e625130a, Julia 1.12.6):

LinearSolve	result
v3.84.0 (pinned)	Success for MIRK4/RadauIIa5/LobattoIIIa4 BigFloat two-point
cdb6f77 (#1017, = v3.85.0)	ArgumentError above, already during BoundaryValueDiffEqShooting precompile
v3.85.1 (registry)	no error, but RadauIIa5 → Stalled

Only #1017 and #1020 are between v3.84.0 and v3.85.1.

Notes from attempted minimization

Direct LinearCache usage with nothing (default) solver did not reproduce wrong solutions on v3.85.1 for either Float64 or BigFloat in these scenarios:

• reassigning cache.A to sparse matrices with different stored nnz (per-solve fallback worked, results correct);
• constant stored pattern with stored zeros that later activate (results correct).

So the corruption/stall appears to need the actual NonlinearSolve/BVP iteration pattern (the failing case goes through NonlinearSolveBase.LinearSolveParameters with the default solver on sparse BigFloat (Sparspak) and Float64 banded/sparse BVP Jacobians). Happy to provide more instrumentation if useful.

Failing CI for reference: https://github.com/SciML/BoundaryValueDiffEq.jl/actions/runs/27135552430 (Misc group, lts/1/pre).

Related downstream issue: SciML/BoundaryValueDiffEq.jl (Misc group red on master).