feat: better leverage mm in tearing

src/carpanzano_tearing.jl

@@ -31,6 +31,16 @@ $TYPEDFIELDS
     filter do not participate in the maximal matching and subsequent SCC decomposition.
     """
     eqfilter::F3 = (_ -> true)
+    """
+    The integer-linear subsystem matrix (see `get_mm`), or `nothing`. When provided,
+    SCCs consisting entirely of integer-linear equations are matched exactly: a
+    fraction-free Bareiss factorization of the SCC's rows over its own variables
+    replaces the equations with their triangular reduced forms and matches each
+    equation to its pivot — a numerically proven assignment. A rank-deficient
+    factorization means the SCC is genuinely singular as scheduled and falls back
+    to the structural heuristics.
+    """
+    mm::Union{Nothing, SparseMatrixCLIL{Int, Int}} = nothing
 end
 
 """
@@ -75,8 +85,99 @@ function update_full_var_eq_matching!(
     end
 end
 
+"""
+    $TYPEDSIGNATURES
+
+Try to match the SCC given by `active_vars`/`active_eqs` exactly. Applicable when the
+SCC is square and every equation is a row of the integer-linear subsystem `mm` (with
+`mm_row_of` mapping equation index to `mm` row index, and each row in sync with the
+structural graph). Runs a fraction-free Bareiss factorization of those rows with
+pivoting restricted to `active_vars` (preferring variables satisfying `isder`):
+
+- Full rank: the SCC is proven nonsingular. The equations are replaced — in `mm`, the
+  structural graph and the solvable graph — by their reduced, triangular forms, and
+  each is matched to its pivot variable. The triangular structure makes the
+  solved-equation dependency graph acyclic (matching against the *original* row
+  contents could produce cycles, which `reassemble` cannot schedule), and every pivot
+  coefficient is proven nonzero, so downstream elimination can never manufacture an
+  exactly singular system from these rows. The `(equation, columns, coefficients)`
+  triples are appended to `rewrites` so the caller can update the symbolic equations
+  to match.
+- Rank deficient: the SCC's coefficient matrix over its own variables is exactly
+  singular — the block cannot determine its variables no matter the matching. Warn
+  and return `false` so the caller falls back to the structural heuristics.
+
+Returns `true` iff the SCC was matched exactly.
+"""
+function exact_scc_matching!(
+        structure::SystemStructure, mm::SparseMatrixCLIL{T, Ti}, mm_row_of::Dict{Int, Int},
+        var_eq_matching::MatchingT, active_vars::AbstractSet{Int},
+        active_eqs::AbstractSet{Int}, isder::F,
+        rewrites::Vector{Tuple{Int, Vector{Int}, Vector{Int}}}
+    ) where {T, Ti, F}
+    n = length(active_eqs)
+    (n >= 2 && n == length(active_vars)) || return false
+    (; graph, solvable_graph) = structure
+
+    # Every equation must be an integer-linear row whose cached coefficients are
+    # in sync with the structural graph.
+    rowids = Vector{Int}(undef, n)
+    for (j, e) in enumerate(active_eqs)
+        i = get(mm_row_of, e, nothing)
+        i === nothing && return false
+        cols = mm.row_cols[i]
+        nbrs = 𝑠neighbors(graph, e)
+        length(nbrs) == length(cols) || return false
+        sort(nbrs) == cols || return false
+        rowids[j] = i
+    end
+
+    nvars = ndsts(graph)
+    active_mask = falses(nvars)
+    tier1 = falses(nvars)
+    for v in active_vars
+        active_mask[v] = true
+        isder !== nothing && isder(v) && (tier1[v] = true)
+    end
+
+    sub = SparseMatrixCLIL{T, Ti}(
+        mm.nparentrows, mm.ncols,
+        Ti[mm.nzrows[i] for i in rowids],
+        Vector{Ti}[copy(mm.row_cols[i]) for i in rowids],
+        Vector{T}[copy(mm.row_vals[i]) for i in rowids])
+    ctx = RestrictedBareissContext(tier1, active_mask, nothing)
+    update! = RestrictedContextUpdate(ctx, bareiss_update_virtual_colswap_mtk!)
+    bareiss_ops = (noop_colswap, SyncedSwapRows(nothing), update!, bareiss_zero!)
+    bareiss!(sub, bareiss_ops; find_pivot = ctx)
+    rank = length(ctx.pivots)
+
+    if rank < n
+        @warn lazy"An SCC of $n integer-linear equations is exactly singular over its own \
+        variables (rank $rank): the block cannot determine the variables it is scheduled \
+        to solve. Falling back to structural tearing; expect a singular linear system \
+        downstream. Equations: $(sort!(collect(active_eqs)))." maxlog = 10
+        return false
+    end
+
+    haskey(ENV, "EXACT_SCC_DEBUG") &&
+        println("exact SCC match: n=$n eqs=$(sort!(collect(active_eqs)))")
+    for k in 1:n
+        e = sub.nzrows[k]
+        i = mm_row_of[e]
+        mm.row_cols[i] = sub.row_cols[k]
+        mm.row_vals[i] = sub.row_vals[k]
+        set_neighbors!(graph, e, sub.row_cols[k])
+        if solvable_graph isa BipartiteGraph{Int, Nothing}
+            set_neighbors!(solvable_graph, e, sub.row_cols[k])
+        end
+        var_eq_matching[ctx.pivots[k]] = e
+        push!(rewrites, (e, sub.row_cols[k], sub.row_vals[k]))
+    end
+    return true
+end
+
 function (alg::CarpanzanoTearing)(structure::SystemStructure)
-    (; isder, varfilter, eqfilter) = alg
+    (; isder, varfilter, eqfilter, mm) = alg
     (; graph, solvable_graph) = structure
 
     var_eq_matching = maximal_matching(
@@ -95,6 +196,15 @@ function (alg::CarpanzanoTearing)(structure::SystemStructure)
     full_var_eq_matching = copy(var_eq_matching)
     var_sccs = find_var_sccs(graph, var_eq_matching)
 
+    # Exact matching of pure integer-linear SCCs is only implemented for
+    # continuous systems.
+    mm_row_of = if mm !== nothing && !is_only_discrete(structure)
+        Dict{Int, Int}(e => i for (i, e) in enumerate(mm.nzrows))
+    else
+        nothing
+    end
+    rewrites = Tuple{Int, Vector{Int}, Vector{Int}}[]
+
     active_vars = OrderedSet{Int}()
     active_eqs = OrderedSet{Int}()
     remaining_eqs = OrderedSet{Int}()
@@ -116,11 +226,17 @@ function (alg::CarpanzanoTearing)(structure::SystemStructure)
             # Tearing will now determine the matching
             var_eq_matching[var] = unassigned
         end
-        carpanzano_tear_scc!(alg, structure, var_eq_matching, active_vars, active_eqs)
+        exact = mm_row_of !== nothing && exact_scc_matching!(
+            structure, mm, mm_row_of, var_eq_matching, active_vars, active_eqs,
+            isder, rewrites)
+        if !exact
+            carpanzano_tear_scc!(alg, structure, var_eq_matching, active_vars, active_eqs)
+        end
         update_full_var_eq_matching!(graph, full_var_eq_matching, var_eq_matching, vars, remaining_eqs; varfilter)
     end
 
-    return TearingResult(var_eq_matching, full_var_eq_matching, var_sccs), (;)
+    extra = (; linear_rewrite = rewrites)
+    return TearingResult(var_eq_matching, full_var_eq_matching, var_sccs), extra
 end
 
 """

src/interface.jl

@@ -106,6 +106,17 @@ function linear_subsys_adjmat! end
 function eq_derivative! end
 function var_derivative! end
 
+"""
+    get_mm(state::TransformationState)
+
+Return the cached integer-linear subsystem matrix (a `SparseMatrixCLIL` whose
+rows are kept in sync with the structural graph, including through
+`eq_derivative!`), or `nothing` when the state does not maintain one. Used by
+tearing algorithms to make exact (rank-aware) decisions about integer-linear
+equations.
+"""
+get_mm(::TransformationState) = nothing
+
 function eq_derivative_graph!(s::SystemStructure, eq::Int)
     add_vertex!(s.graph, SRC)
     s.solvable_graph === nothing || add_vertex!(s.solvable_graph, SRC)

src/partial_state_selection.jl

@@ -179,7 +179,10 @@ function dummy_derivative_graph!(state::TransformationState, jac = nothing;
     state.structure.solvable_graph === nothing && find_solvables!(state; kwargs...)
     complete!(state.structure)
     var_eq_matching = complete(pantelides!(state; kwargs...))
-    dummy_derivative_graph!(state.structure, var_eq_matching, jac, state_priority, log)
+    # NOTE: `get_mm` must be queried after `pantelides!`, which extends the
+    # linear subsystem matrix with differentiated rows (`eq_derivative!`).
+    dummy_derivative_graph!(
+        state.structure, var_eq_matching, jac, state_priority, log; mm = get_mm(state))
 end
 
 struct DummyDerivativeSummary
@@ -203,7 +206,8 @@ Perform the dummy derivatives algorithm.
 function dummy_derivative_graph!(
         structure::SystemStructure, var_eq_matching, jac = nothing,
         state_priority = nothing, ::Val{log} = Val(false);
-        tearing_alg::TearingAlgorithm = DummyDerivativeTearing(), kwargs...) where {log}
+        tearing_alg::TearingAlgorithm = DummyDerivativeTearing(),
+        mm = nothing, kwargs...) where {log}
     (; eq_to_diff, var_to_diff, graph) = structure
     diff_to_eq = invview(eq_to_diff)
     diff_to_var = invview(var_to_diff)
@@ -370,7 +374,11 @@ function dummy_derivative_graph!(
         @warn "The number of dummy derivatives ($n_dummys) does not match the number of differentiated equations ($n_diff_eqs)."
     end
 
-    tearing_result, extra = tearing_alg(structure, BitSet(dummy_derivatives))
+    tearing_result, extra = if tearing_alg isa DummyDerivativeTearing
+        tearing_alg(structure, BitSet(dummy_derivatives); mm)
+    else
+        tearing_alg(structure, BitSet(dummy_derivatives))
+    end
     extra = (; extra..., ddsummary = DummyDerivativeSummary(var_dummy_scc, var_state_priority))
     return tearing_result, extra
 end
@@ -428,7 +436,8 @@ struct DummyDerivativeTearing{T <: TearingAlgorithm} <: TearingAlgorithm end
 DummyDerivativeTearing() = DummyDerivativeTearing{CarpanzanoTearing}()
 
 function (::DummyDerivativeTearing{T})(
-        structure::SystemStructure, dummy_derivatives::Union{BitSet, Tuple{}} = ()
+        structure::SystemStructure, dummy_derivatives::Union{BitSet, Tuple{}} = ();
+        mm = nothing
     ) where {T}
     (; var_to_diff) = structure
     # We can eliminate variables that are not selected (differential
@@ -441,11 +450,19 @@ function (::DummyDerivativeTearing{T})(
             can_eliminate[v] = true
         end
     end
-    inner_tearing_alg = T(;
-        isder = Base.Fix1(isdiffed, (structure, dummy_derivatives)),
-        varfilter = Base.Fix1(getindex, can_eliminate)
-    )
-    tearing_result, _ = inner_tearing_alg(structure)
+    inner_tearing_alg = if :mm in fieldnames(T)
+        T(;
+            isder = Base.Fix1(isdiffed, (structure, dummy_derivatives)),
+            varfilter = Base.Fix1(getindex, can_eliminate),
+            mm
+        )
+    else
+        T(;
+            isder = Base.Fix1(isdiffed, (structure, dummy_derivatives)),
+            varfilter = Base.Fix1(getindex, can_eliminate)
+        )
+    end
+    tearing_result, inner_extra = inner_tearing_alg(structure)
 
     for v in 𝑑vertices(structure.graph)
         is_present(structure, v) || continue
@@ -454,5 +471,5 @@ function (::DummyDerivativeTearing{T})(
         tearing_result.var_eq_matching[v] = SelectedState()
     end
 
-    return tearing_result, (; can_eliminate)
+    return tearing_result, (; can_eliminate, inner_extra...)
 end

src/singularity_removal.jl

@@ -540,6 +540,53 @@ struct PivotInfo
     pivots::Vector{Int}
 end
 
+"""
+    $(TYPEDEF)
+
+Tiered pivot search used for exact Bareiss factorization of pure-integer SCC
+subsets during tearing. Unlike [`BareissContext`](@ref), there is no
+unrestricted final tier: once `tier2` is exhausted the factorization stops.
+This guarantees every chosen pivot column is an eligible (solvable-for) variable.
+"""
+mutable struct RestrictedBareissContext{V1 <: AbstractVector{Bool}, V2 <: AbstractVector{Bool}, P <: Union{Nothing, AbstractVector{Int}}}
+    pivots::Vector{Int}
+    tier1::V1
+    tier2::V2
+    valid_pivot_mask::BitVector
+    var_priorities::P
+    tier1_done::Bool
+end
+
+function RestrictedBareissContext(tier1, tier2, var_priorities = nothing)
+    return RestrictedBareissContext(
+        Int[], tier1, tier2, trues(length(tier1)), var_priorities, false)
+end
+
+function (ctx::RestrictedBareissContext)(M, k::Int)
+    if !ctx.tier1_done
+        r = find_masked_pivot(LazyMaskAnd(ctx.tier1, ctx.valid_pivot_mask), M, k, ctx.var_priorities)
+        if r !== nothing
+            push!(ctx.pivots, r[1][2])
+            return r
+        end
+        ctx.tier1_done = true
+    end
+    r = find_masked_pivot(LazyMaskAnd(ctx.tier2, ctx.valid_pivot_mask), M, k, ctx.var_priorities)
+    r === nothing && return nothing
+    push!(ctx.pivots, r[1][2])
+    return r
+end
+
+struct RestrictedContextUpdate{C <: RestrictedBareissContext, F}
+    context::C
+    inner_update::F
+end
+
+function (bcu::RestrictedContextUpdate)(zero!, M, k, swapto, pivot, last_pivot; kw...)
+    bcu.context.valid_pivot_mask[swapto[2]] = false
+    return bcu.inner_update(zero!, M, k, swapto, pivot, last_pivot; kw...)
+end
+
 function _uf_find!(parent::Vector{Int}, x::Int)
     while parent[x] != x
         parent[x] = parent[parent[x]]