Merge pull request #1380 from Keno/kf/bareissbug

Fix bug in bareiss algorithm

Author: YingboMa

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