Fix bug in bareiss algorithm

We need to maintain the invariant that the rows of the sparse
matrix are sorted, otherwise the bareiss algorithm will be incorrect.
Without this, for appropriate choice of pivots, our bareiss
implementation turned

```
[1 1 0
 1 0 1
 0 1 1]
```

into

```
[1 1 0
 0 0 -1
 0 0 0]
```

rather than the correct

```
[1 1 0
 0 1 -1
 0 0 0]
```

Fix that and refactor a bit to make this testable.

Author: Keno

src/systems/alias_elimination.jl

@@ -2,15 +2,30 @@ using SymbolicUtils: Rewriters
 
 const KEEP = typemin(Int)
 
-function alias_elimination(sys)
+function alias_eliminate_graph(sys::AbstractSystem)
     sys = initialize_system_structure(sys; quick_cancel=true)
     s = structure(sys)
 
     mm = linear_subsys_adjmat(sys)
-    size(mm, 1) == 0 && return sys # No linear subsystems
+    size(mm, 1) == 0 && return sys, nothing, mm # No linear subsystems
 
     ag, mm = alias_eliminate_graph!(s.graph, complete(s.var_to_diff), mm)
+    return sys, ag, mm
+end
+
+# For debug purposes
+function aag_bareiss(sys::AbstractSystem)
+    sys = initialize_system_structure(sys; quick_cancel=true)
+    s = structure(sys)
+    mm = linear_subsys_adjmat(sys)
+    return aag_bareiss!(s.graph, complete(s.var_to_diff), mm)
+end
+
+function alias_elimination(sys)
+    sys, ag, mm = alias_eliminate_graph(sys)
+    ag === nothing && return sys
 
+    s = structure(sys)
     @unpack fullvars, graph = s
 
     subs = OrderedDict()
@@ -193,9 +208,7 @@ count_nonzeros(a::AbstractArray) = count(!iszero, a)
 # Here we have a guarantee that they won't, so we can make this identification
 count_nonzeros(a::SparseVector) = nnz(a)
 
-function alias_eliminate_graph!(graph, var_to_diff, mm_orig::SparseMatrixCLIL)
-    diff_to_var = invview(var_to_diff)
-
+function aag_bareiss!(graph, var_to_diff, mm_orig::SparseMatrixCLIL)
     mm = copy(mm_orig)
     is_linear_equations = falses(size(AsSubMatrix(mm_orig), 1))
     for e in mm_orig.nzrows
@@ -245,10 +258,10 @@ function alias_eliminate_graph!(graph, var_to_diff, mm_orig::SparseMatrixCLIL)
         (rank1, rank2, rank3, pivots)
     end
 
-    # mm2 = Array(copy(mm))
-    # @show do_bareiss!(mm2)
-    # display(mm2)
+    return mm, solvable_variables, do_bareiss!(mm, mm_orig)
+end
 
+function alias_eliminate_graph!(graph, var_to_diff, mm_orig::SparseMatrixCLIL)
     # Step 1: Perform bareiss factorization on the adjacency matrix of the linear
     #         subsystem of the system we're interested in.
     #
@@ -261,7 +274,8 @@ function alias_eliminate_graph!(graph, var_to_diff, mm_orig::SparseMatrixCLIL)
     # -------------------|------------------------
     # rank3 | [ 0    0   | M₃₃  M₃₄ ]   [v₃] = [0]
     #         [ 0    0   | 0    0   ]   [v₄] = [0]
-    (rank1, rank2, rank3, pivots) = do_bareiss!(mm, mm_orig)
+    mm, solvable_variables, (rank1, rank2, rank3, pivots) =
+        aag_bareiss!(graph, var_to_diff, mm_orig)
 
     # Step 2: Simplify the system using the Bareiss factorization
     ag = AliasGraph(size(mm, 2))
@@ -272,6 +286,7 @@ function alias_eliminate_graph!(graph, var_to_diff, mm_orig::SparseMatrixCLIL)
     # Kind of like the backward substitution, but we don't actually rely on it
     # being lower triangular. We eliminate a variable if there are at most 2
     # variables left after the substitution.
+    diff_to_var = invview(var_to_diff)
     function lss!(ei::Integer)
         vi = pivots[ei]
         # the lowest differentiated variable can be eliminated

src/systems/sparsematrixclil.jl

@@ -32,6 +32,13 @@ function swaprows!(S::SparseMatrixCLIL, i, j)
     swap!(S.row_vals, i, j)
 end
 
+function SparseMatrixCLIL(mm::AbstractMatrix)
+    nrows, ncols = size(mm)
+    row_cols = [findall(!iszero, row) for row in eachrow(mm)]
+    row_vals = [row[cols] for (row, cols) in zip(eachrow(mm), row_cols)]
+    SparseMatrixCLIL(nrows, ncols, Int[1:length(row_cols);], row_cols, row_vals)
+end
+
 struct CLILVector{T, Ti} <: AbstractSparseVector{T, Ti}
     vec::SparseVector{T, Ti}
 end
@@ -166,7 +173,9 @@ function bareiss_update_virtual_colswap_mtk!(zero!, M::SparseMatrixCLIL, k, swap
 
         tmp_incidence = similar(eadj[ei], 0)
         tmp_coeffs = similar(old_cadj[ei], 0)
-        vars = union(ivars, kvars)
+        # TODO: We know both ivars and kvars are sorted, we could just write
+        # a quick iterator here that does this without allocation/faster.
+        vars = sort(union(ivars, kvars))
 
         for v in vars
             v == vpivot && continue
@@ -197,7 +206,7 @@ struct AsSubMatrix{T, Ti<:Integer} <: AbstractSparseMatrix{T, Ti}
 end
 Base.size(S::AsSubMatrix) = (S.M.nparentrows, S.M.ncols)
 
-function Base.getindex(S::SparseMatrixCLIL{T}, i1, i2) where {T}
+function Base.getindex(S::SparseMatrixCLIL{T}, i1::Integer, i2::Integer) where {T}
     checkbounds(S, i1, i2)
 
     col = S.row_cols[i1]
@@ -207,7 +216,7 @@ function Base.getindex(S::SparseMatrixCLIL{T}, i1, i2) where {T}
     return S.row_vals[i1][nncol]
 end
 
-function Base.getindex(S::AsSubMatrix{T}, i1, i2) where {T}
+function Base.getindex(S::AsSubMatrix{T}, i1::Integer, i2::Integer) where {T}
     checkbounds(S, i1, i2)
     S = S.M
 

test/structural_transformation/tearing.jl

@@ -139,6 +139,14 @@ eqs = [
        0 ~ x + z,
       ]
 @named nlsys = NonlinearSystem(eqs, [x, y, z], [])
+let (mm, _, _) = ModelingToolkit.aag_bareiss(nlsys)
+    @test mm == [
+        -1  1  0;
+         0 -1 -1;
+         0  0  0
+    ]
+end
+
 newsys = tearing(nlsys)
 @test length(equations(newsys)) == 1