diff --git a/src/carpanzano_tearing.jl b/src/carpanzano_tearing.jl index 6b235e7..584eab1 100644 --- a/src/carpanzano_tearing.jl +++ b/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 """ diff --git a/src/interface.jl b/src/interface.jl index b26be7d..29cc201 100644 --- a/src/interface.jl +++ b/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) diff --git a/src/partial_state_selection.jl b/src/partial_state_selection.jl index 13c1054..f7ed304 100644 --- a/src/partial_state_selection.jl +++ b/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 diff --git a/src/singularity_removal.jl b/src/singularity_removal.jl index ce5924f..2ecd53e 100644 --- a/src/singularity_removal.jl +++ b/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]]