Fix Barreiss algorithm (#2089)

This has a number of changes.
The first is that we were using ÷ instead of exact_div in bareiss_update_virtual_colswap_mtk
which led to undetected overflows. We also had a number of bugs that arise when using BigInt precsion in bareiss (e.g. checking `!==0` vs `!=0`.
This PR makes it so we first try to bareiss with the regular Int64 matrix, but if that gets an overflow, we fallback to
the BigInt version.

Author: Oscar Smith

src/systems/alias_elimination.jl

@@ -1,8 +1,6 @@
 using SymbolicUtils: Rewriters
 using Graphs.Experimental.Traversals
 
-const KEEP = typemin(Int)
-
 function alias_eliminate_graph!(state::TransformationState; kwargs...)
     mm = linear_subsys_adjmat!(state; kwargs...)
     if size(mm, 1) == 0
@@ -225,8 +223,9 @@ the `constraint`.
         vertices = eadj[i]
         if constraint(length(vertices))
             for (j, v) in enumerate(vertices)
-                (mask === nothing || mask[v]) &&
+                if (mask === nothing || mask[v])
                     return (CartesianIndex(i, v), M.row_vals[i][j])
+                end
             end
         end
     end
@@ -241,7 +240,6 @@ end
         row = @view M[i, :]
         if constraint(count(!iszero, row))
             for (v, val) in enumerate(row)
-                iszero(val) && continue
                 if mask === nothing || mask[v]
                     return CartesianIndex(i, v), val
                 end
@@ -325,7 +323,8 @@ function Base.setindex!(ag::AliasGraph, v::Integer, i::Integer)
     return 0 => 0
 end
 
-function Base.setindex!(ag::AliasGraph, p::Union{Pair{Int, Int}, Tuple{Int, Int}},
+function Base.setindex!(ag::AliasGraph,
+                        p::Union{Pair{<:Integer, Int}, Tuple{<:Integer, Int}},
                         i::Integer)
     (c, v) = p
     if c == 0 || v == 0
@@ -530,7 +529,7 @@ function find_linear_variables(graph, linear_equations, var_to_diff, irreducible
     return linear_variables
 end
 
-function aag_bareiss!(graph, var_to_diff, mm_orig::SparseMatrixCLIL)
+function aag_bareiss!(graph, var_to_diff, mm_orig::SparseMatrixCLIL{T, Ti}) where {T, Ti}
     mm = copy(mm_orig)
     linear_equations_set = BitSet(mm_orig.nzrows)
 
@@ -554,7 +553,16 @@ function aag_bareiss!(graph, var_to_diff, mm_orig::SparseMatrixCLIL)
     end
     solvable_variables = findall(is_linear_variables)
 
-    return mm, solvable_variables, do_bareiss!(mm, mm_orig, is_linear_variables)
+    local bar
+    try
+        bar = do_bareiss!(mm, mm_orig, is_linear_variables)
+    catch e
+        e isa OverflowError || rethrow(e)
+        mm = convert(SparseMatrixCLIL{BigInt, Ti}, mm_orig)
+        bar = do_bareiss!(mm, mm_orig, is_linear_variables)
+    end
+
+    return mm, solvable_variables, bar
 end
 
 function do_bareiss!(M, Mold, is_linear_variables)
@@ -589,6 +597,7 @@ function do_bareiss!(M, Mold, is_linear_variables)
     end
     bareiss_ops = ((M, i, j) -> nothing, myswaprows!,
                    bareiss_update_virtual_colswap_mtk!, bareiss_zero!)
+
     rank2, = bareiss!(M, bareiss_ops; find_pivot = find_and_record_pivot)
     rank1 = something(rank1r[], rank2)
     (rank1, rank2, pivots)
@@ -983,7 +992,7 @@ function locally_structure_simplify!(adj_row, pivot_var, ag)
     nirreducible = 0
     # When this row only as the pivot element, the pivot is zero by homogeneity
     # of the linear system.
-    alias_candidate::Union{Int, Pair{Int, Int}} = 0
+    alias_candidate::Union{Int, Pair{eltype(adj_row), Int}} = 0
 
     # N.B.: Assumes that the non-zeros iterator is robust to modification
     # of the underlying array datastructure.

src/systems/sparsematrixclil.jl

@@ -29,11 +29,20 @@ function Base.copy(S::SparseMatrixCLIL{T, Ti}) where {T, Ti}
                      map(copy, S.row_vals))
 end
 function swaprows!(S::SparseMatrixCLIL, i, j)
+    i == j && return
     swap!(S.nzrows, i, j)
     swap!(S.row_cols, i, j)
     swap!(S.row_vals, i, j)
 end
 
+function Base.convert(::Type{SparseMatrixCLIL{T, Ti}}, S::SparseMatrixCLIL) where {T, Ti}
+    return SparseMatrixCLIL(S.nparentrows,
+                            S.ncols,
+                            copy.(S.nzrows),
+                            copy.(S.row_cols),
+                            [T.(row) for row in S.row_vals])
+end
+
 function SparseMatrixCLIL(mm::AbstractMatrix)
     nrows, ncols = size(mm)
     row_cols = [findall(!iszero, row) for row in eachrow(mm)]
@@ -59,6 +68,7 @@ function Base.setindex!(S::SparseMatrixCLIL, v::CLILVector, i::Integer, c::Colon
     if v.vec.n != S.ncols
         throw(BoundsError(v, 1:(S.ncols)))
     end
+    any(iszero, v.vec.nzval) && error("setindex failed")
     S.row_cols[i] = copy(v.vec.nzind)
     S.row_vals[i] = copy(v.vec.nzval)
     return v
@@ -91,7 +101,7 @@ function Base.iterate(nzp::NonZerosPairs{<:CLILVector}, (idx, col))
         idx = length(col)
     end
     oldcol = nzind[idx]
-    if col !== oldcol
+    if col != oldcol
         # The vector was changed since the last iteration. Find our
         # place in the vector again.
         tail = col > oldcol ? (@view nzind[(idx + 1):end]) : (@view nzind[1:idx])
@@ -189,13 +199,14 @@ function bareiss_update_virtual_colswap_mtk!(zero!, M::SparseMatrixCLIL, k, swap
             v == vpivot && continue
             ck = getcoeff(kvars, kcoeffs, v)
             ci = getcoeff(ivars, icoeffs, v)
-            ci = (pivot * ci - coeff * ck) ÷ last_pivot
-            if ci !== 0
+            p1 = Base.Checked.checked_mul(pivot, ci)
+            p2 = Base.Checked.checked_mul(coeff, ck)
+            ci = exactdiv(Base.Checked.checked_sub(p1, p2), last_pivot)
+            if !iszero(ci)
                 push!(tmp_incidence, v)
                 push!(tmp_coeffs, ci)
             end
         end
-
         eadj[ei] = tmp_incidence
         old_cadj[ei] = tmp_coeffs
     end