Merge pull request #97 from JuliaComputing/os/fix-inline-linear-solve

fix edge case in inline linear solve

Author: AayushSabharwal

lib/ModelingToolkitTearing/src/reassemble.jl

@@ -787,7 +787,8 @@ function __reduce_linear_system!(A::StateSelection.CLIL.SparseMatrixCLIL{Num, In
             cst += coeff * other_cst
         end
 
-        aliases[var] = SparseArrays.sparsevec(new_I, new_V, length(coeffs))
+        # `sparsevec` sums duplicate indices but keeps explicit zeros; drop them.
+        aliases[var] = SparseArrays.dropzeros!(SparseArrays.sparsevec(new_I, new_V, length(coeffs)))
         constants[var] = cst
     end
 

src/math/sparsematrixclil.jl

@@ -91,6 +91,24 @@ zero!(a::SparseVector) = (empty!(a.nzind); empty!(a.nzval))
 zero!(a::CLILVector) = zero!(a.vec)
 SparseArrays.dropzeros!(a::CLILVector) = SparseArrays.dropzeros!(a.vec)
 
+# Remove explicitly-stored zeros from each row, in place.
+function SparseArrays.dropzeros!(S::SparseMatrixCLIL)
+    for r in eachindex(S.row_vals)
+        cols = S.row_cols[r]
+        vals = S.row_vals[r]
+        j = 0
+        for k in eachindex(vals)
+            iszero(vals[k]) && continue
+            j += 1
+            cols[j] = cols[k]
+            vals[j] = vals[k]
+        end
+        resize!(cols, j)
+        resize!(vals, j)
+    end
+    return S
+end
+
 struct NonZeros{T <: AbstractArray}
     v::T
 end

src/utils.jl

@@ -340,24 +340,30 @@ function get_new_mm(
         push!(final_row_cols, new_row_col_i[indices[1]])
         push!(final_row_vals, new_row_val_i[indices[1]])
         for i in Iterators.drop(eachindex(indices), 1)
-            if new_row_col_i[indices[i]] == new_row_col_i[indices[i - 1]]
+            col = new_row_col_i[indices[i]]
+            # Compare against the last retained column, not the previous sorted
+            # entry: a prior cancellation may have `pop!`ed the matching entry.
+            if !isempty(final_row_cols) && col == final_row_cols[end]
                 final_row_vals[end] += new_row_val_i[indices[i]]
                 if iszero(final_row_vals[end])
                     pop!(final_row_cols)
                     pop!(final_row_vals)
                 end
             else
-                push!(final_row_cols, new_row_col_i[indices[i]])
+                push!(final_row_cols, col)
                 push!(final_row_vals, new_row_val_i[indices[i]])
             end
         end
-        
+
         push!(new_row_cols, final_row_cols)
         push!(new_row_vals, final_row_vals)
         push!(new_nzrows, old_to_new_eq[eq])
     end
 
-    return typeof(mm)(new_nparentrows, count(!iszero, old_to_new_var), new_nzrows, new_row_cols, new_row_vals)
+    # Drop explicit zeros (e.g. from an alias that cancelled to zero) so they
+    # don't survive as phantom structural nonzeros.
+    new_mm = typeof(mm)(new_nparentrows, count(!iszero, old_to_new_var), new_nzrows, new_row_cols, new_row_vals)
+    return SparseArrays.dropzeros!(new_mm)
 end
 
 struct BadMMAliasError <: Exception

test/runtests.jl

@@ -29,4 +29,46 @@ include("carpanzano_tearing.jl")
     @test mm2.ncols == 2
     aliases = Dict(1 => sparse([0, 0, 2, 1]))
     @test_throws SSel.BadMMAliasError SSel.get_new_mm(aliases, old_to_new_eq, old_to_new_var, mm)
+
+    # Three entries collapse onto the same column, and the first two cancel to
+    # zero. The remaining duplicate must survive as a fresh entry rather than
+    # being folded into the (now popped) accumulator.
+    # Row `[2 1 1]`: var1 retained, var2 -> -2*var1, var3 -> 3*var1.
+    # Contributions to new col 1: 2 (direct) + 1*(-2) + 1*3 = 3.
+    mm = SSel.CLIL.SparseMatrixCLIL([2 1 1])
+    old_to_new_eq = [1]
+    old_to_new_var = [1, 0, 0]
+    aliases = Dict(2 => sparse([-2, 0, 0]), 3 => sparse([3, 0, 0]))
+    mm2 = SSel.get_new_mm(aliases, old_to_new_eq, old_to_new_var, mm)
+    @test mm2.nzrows == [1]
+    @test mm2.row_cols == [[1]]
+    @test mm2.row_vals == [[3]]
+    @test mm2.ncols == 1
+
+    # Same cancellation, but with a distinct retained column already present
+    # before the cancelling group. The surviving duplicate must not leak into
+    # the earlier column.
+    # Row `[7 2 1 1]`: var1 -> col1, var2 -> col2, var3 -> -2*var2, var4 -> 5*var2.
+    # col1 = 7, col2 = 2 + 1*(-2) + 1*5 = 5.
+    mm = SSel.CLIL.SparseMatrixCLIL([7 2 1 1])
+    old_to_new_eq = [1]
+    old_to_new_var = [1, 2, 0, 0]
+    aliases = Dict(3 => sparse([0, -2, 0, 0]), 4 => sparse([0, 5, 0, 0]))
+    mm2 = SSel.get_new_mm(aliases, old_to_new_eq, old_to_new_var, mm)
+    @test mm2.nzrows == [1]
+    @test mm2.row_cols == [[1, 2]]
+    @test mm2.row_vals == [[7, 5]]
+    @test mm2.ncols == 2
+
+    # An alias carrying an explicit stored zero (as `sparsevec` produces when
+    # duplicate indices cancel) must not leak in as a phantom structural nonzero.
+    # Row references var2 (removed, aliased to 0*var1) and var3 (retained); there
+    # is no direct var1 term, so the only contribution to col 1 is the zero.
+    mm = SSel.CLIL.SparseMatrixCLIL([0 1 1])
+    old_to_new_eq = [1]
+    old_to_new_var = [1, 0, 2]
+    aliases = Dict(2 => sparsevec([1, 1], [2, -2], 3))
+    mm2 = SSel.get_new_mm(aliases, old_to_new_eq, old_to_new_var, mm)
+    @test mm2.row_cols == [[2]]
+    @test mm2.row_vals == [[1]]
 end