Auto nonstructural zeros: decide on filled values, deactivate when nothing is droppable (#1110)

Two fixes to the `NonstructuralZeros.Auto` state machine, both aimed at the same
failure: the reduction latching into per-solve `dropzeros` that drops nothing.

1. Defer the activation decision when the starting matrix is entirely zero.
   A cache is commonly initialized from a *prototype* — the intended sparsity
   pattern with values not filled in yet — and an all-zero sample is
   indistinguishable from "100% nonstructural zeros". `Auto` activated on it with
   an all-`false` mask, and the first filled solve then flipped every entry, so
   `activated == nstart_zeros` unconditionally tripped the non-persistence guard
   and set `cache_union = false` permanently. Every later solve then allocated a
   fresh `dropzeros` copy while dropping nothing. Deferring to the first matrix
   that carries a nonzero decides on real values instead.

2. Deactivate once the union covers every stored entry. The union only grows, so
   at that point no future matrix can have anything dropped — the reduction is
   provably dead weight. Switching it off is also structure-preserving, since
   `reduced` equals the full pattern exactly when this triggers, so the operand
   handed to the factorization keeps the structure it just had.

Measured on a prototype-initialized cache (n=200, nnz=598, KLU): 17.7 kB per
solve with `nnz(reduced) == nnz(A)` before, 0 B after. The reduction itself is
unchanged where it pays: a matrix with genuine nonstructural zeros still reduces
(359/598 stored entries kept) at 0 B per solve, and prototype-initialized callers
now *gain* that reduction on the first filled solve instead of latching into the
degraded path.

This also removes a 30x sparse allocation gap in OrdinaryDiffEq's
`NonlinearSolveAlg`, whose inner Jacobian buffer is exactly such a prototype.



Claude-Session: https://claude.ai/code/session_01NB3oFTNzGW79UtDt8AyjsM

Co-authored-by: ChrisRackauckas-Claude <accounts@chrisrackauckas.com>
Co-authored-by: Claude Fable 5 <noreply@anthropic.com>

Author: 3 people

Project.toml

@@ -1,6 +1,6 @@
 name = "LinearSolve"
 uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
-version = "5.1.1"
+version = "5.1.2"
 authors = ["SciML"]
 
 [deps]

ext/LinearSolveSparseArraysExt.jl

@@ -1200,6 +1200,13 @@ end
 # so `init_sparse_reduction` returns a concrete type either way — type-stable.
 mutable struct SparseReduction{Tv, Ti}
     active::Bool
+    # `Auto` only: the activation decision is deferred until the first solve whose matrix
+    # carries a nonzero. A cache is commonly initialized from a *prototype* — the intended
+    # sparsity pattern with values not filled in yet (all zero) — which says nothing about
+    # nonstructural zeros. Taking it at face value activates the reduction on an all-zero
+    # mask, and the first real solve then bloats the union and degrades to per-solve
+    # dropzeros for the rest of the run. Deciding on the first filled matrix sees real values.
+    pending::Bool
     cache_union::Bool             # true: union caching; false: per-solve dropzeros
     auto::Bool                    # may switch union -> per-solve on bloat
     nstart_zeros::Int             # # stored entries zero on the starting matrix
@@ -1255,15 +1262,18 @@ function LinearSolve.init_sparse_reduction(
     nsz = LinearSolve.__nonstructural_zeros(assumptions)
     NZ = LinearSolve.NonstructuralZeros
     nz = nonzeros(A)
+    auto = nsz == NZ.Auto
+    # An all-zero starting matrix is a prototype whose values are not filled in yet: its zero
+    # fraction carries no information, so defer the `Auto` decision to the first filled solve.
+    pending = auto && !isempty(nz) && all(iszero, nz)
     active = if nsz == NZ.Persistent || nsz == NZ.Present
         true
-    elseif nsz == NZ.None
+    elseif nsz == NZ.None || pending
         false
-    else  # Auto
+    else  # Auto, on a matrix that carries values
         !isempty(nz) &&
             count(iszero, nz) / length(nz) >= LinearSolve.PERSISTENT_ZERO_FRACTION_THRESHOLD
     end
-    auto = nsz == NZ.Auto
     cache_union = nsz != NZ.Present     # Present => per-solve dropzeros from the start
     colptr = copy(getcolptr(A))
     rowval = copy(rowvals(A))
@@ -1272,19 +1282,43 @@ function LinearSolve.init_sparse_reduction(
         nstart_zeros = count(iszero, nz)
         keep, reduced = _persistent_reduced(A, colptr, rowval, mask, nz)
         return SparseReduction{Tv, Ti}(
-            true, cache_union, auto, nstart_zeros, colptr, rowval, mask, keep, reduced,
-            1, 0, nnz(reduced), false
+            true, false, cache_union, auto, nstart_zeros, colptr, rowval, mask, keep,
+            reduced, 1, 0, nnz(reduced), false
         )
     else
         reduced = LinearSolve.make_SparseMatrixCSC(A)
         return SparseReduction{Tv, Ti}(
-            false, cache_union, auto, 0, colptr, rowval, Bool[], Int[], reduced,
+            false, pending, cache_union, auto, 0, colptr, rowval, Bool[], Int[], reduced,
             0, 0, 0, false
         )
     end
 end
 
+# Resolve a deferred `Auto` activation (see `SparseReduction.pending`) against the first
+# matrix that carries values. A still-all-zero matrix remains uninformative, so keep waiting;
+# there is nothing worth reducing in it either way.
+function _resolve_pending!(red::SparseReduction, A)
+    nz = nonzeros(A)
+    (isempty(nz) || all(iszero, nz)) && return red
+    red.pending = false
+    count(iszero, nz) / length(nz) >=
+        LinearSolve.PERSISTENT_ZERO_FRACTION_THRESHOLD || return red
+    red.active = true
+    red.nstart_zeros = count(iszero, nz)
+    red.colptr = copy(getcolptr(A))
+    red.rowval = copy(rowvals(A))
+    red.mask = Bool[!iszero(v) for v in nz]
+    red.keep, red.reduced = _persistent_reduced(A, red.colptr, red.rowval, red.mask, nz)
+    red.nanalyze += 1
+    red.reduced_nnz = nnz(red.reduced)
+    # The operand handed to the factorization now has a different structure than the one the
+    # cache was initialized with, so cached symbolic factorizations must be redone.
+    red.structure_changed = true
+    return red
+end
+
 function LinearSolve.reduce_operand!(red::SparseReduction, A)
+    red.pending && _resolve_pending!(red, A)
     red.active || return A
     nz = nonzeros(A)
     # A change in the stored nnz breaks union caching: an explicit `Persistent`
@@ -1319,6 +1353,15 @@ function LinearSolve.reduce_operand!(red::SparseReduction, A)
             A, red.colptr, red.rowval, red.mask, nz
         )
         red.nanalyze += 1
+        # The union only ever grows, so once it covers every stored entry no future matrix
+        # can have anything dropped: the reduction is dead weight from here on. Deactivating
+        # is safe rather than merely cheap — `reduced` currently equals the full pattern, so
+        # the operand handed to the factorization keeps the same structure it just had.
+        if length(red.keep) == length(red.mask)
+            red.active = false
+            red.structure_changed = false
+            return A
+        end
         # auto mode: if more than NONPERSISTENT_ZERO_FRACTION of the starting zeros
         # have activated, the zeros are not persistent — stop caching the union and
         # drop per solve instead (this matrix's own zeros only). `activated` is

test/Core/nonstructural_zeros.jl

@@ -282,4 +282,64 @@ reduction(cache) = cache.cacheval.sparse_reduction
         @test sol2.u ≈ Matrix(A_new2) \ b rtol = 1.0e-8
         @test !reduction(cache).cache_union
     end
+
+    @testset "auto defers the decision on an all-zero starting matrix" begin
+        # A cache is often initialized from a *prototype*: the intended sparsity pattern
+        # with values not filled in yet. Its zero fraction says nothing about
+        # nonstructural zeros, so `auto` must not decide from it -- taking it at face
+        # value activates on an all-zero mask, and the first filled solve then bloats the
+        # union and latches per-solve dropzeros while dropping nothing.
+        proto = copy(mats[1])
+        nonzeros(proto) .= 0.0
+
+        cache = init(LinearProblem(copy(proto), copy(b)))   # auto
+        @test !reduction(cache).active
+        @test reduction(cache).pending
+
+        # filled with a matrix that genuinely carries nonstructural zeros: activate
+        filled = copy(mats[1])
+        cache.A = copy(filled)
+        sol = solve!(cache)
+        @test sol.retcode == ReturnCode.Success
+        @test sol.u ≈ Matrix(filled) \ b rtol = 1.0e-8
+        @test !reduction(cache).pending
+        @test reduction(cache).active && reduction(cache).cache_union
+        @test reduction(cache).reduced_nnz < nnz(filled)
+
+        # filled with a matrix that has no droppable zeros: stay inactive
+        dense_valued = copy(mats[1])
+        nonzeros(dense_valued) .= 1.5
+        c2 = init(LinearProblem(copy(proto), copy(b)))
+        @test c2.cacheval.sparse_reduction.pending
+        c2.A = copy(dense_valued)
+        sol2 = solve!(c2)
+        @test sol2.retcode == ReturnCode.Success
+        @test sol2.u ≈ Matrix(dense_valued) \ b rtol = 1.0e-8
+        @test !reduction(c2).pending
+        @test !reduction(c2).active
+    end
+
+    @testset "auto deactivates once the union covers every stored entry" begin
+        # The union only grows, so when it covers the whole stored pattern nothing can
+        # ever be dropped again: the reduction is dead weight and must switch off rather
+        # than fall back to allocating a per-solve dropzeros copy that drops nothing.
+        start = copy(mats[1])                      # has stored zeros -> auto activates
+        cache = init(LinearProblem(copy(start), copy(b)))
+        @test reduction(cache).active
+        full = copy(mats[1])
+        nonzeros(full) .= 2.0                      # every stored entry becomes nonzero
+        cache.A = copy(full)
+        sol = solve!(cache)
+        @test sol.retcode == ReturnCode.Success
+        @test sol.u ≈ Matrix(full) \ b rtol = 1.0e-8
+        @test !reduction(cache).active
+        # and it keeps solving correctly afterwards, with no reduction overhead
+        full2 = copy(full)
+        nonzeros(full2) .= 3.0
+        cache.A = copy(full2)
+        sol2 = solve!(cache)
+        @test sol2.retcode == ReturnCode.Success
+        @test sol2.u ≈ Matrix(full2) \ b rtol = 1.0e-8
+        @test !reduction(cache).active
+    end
 end