Skip to content

Commit ae78443

Browse files
committed
dropzeros the mm
1 parent 214d022 commit ae78443

4 files changed

Lines changed: 39 additions & 7 deletions

File tree

lib/ModelingToolkitTearing/src/reassemble.jl

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -787,7 +787,8 @@ function __reduce_linear_system!(A::StateSelection.CLIL.SparseMatrixCLIL{Num, In
787787
cst += coeff * other_cst
788788
end
789789

790-
aliases[var] = SparseArrays.sparsevec(new_I, new_V, length(coeffs))
790+
# `sparsevec` sums duplicate indices but keeps explicit zeros; drop them.
791+
aliases[var] = SparseArrays.dropzeros!(SparseArrays.sparsevec(new_I, new_V, length(coeffs)))
791792
constants[var] = cst
792793
end
793794

src/math/sparsematrixclil.jl

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -91,6 +91,24 @@ zero!(a::SparseVector) = (empty!(a.nzind); empty!(a.nzval))
9191
zero!(a::CLILVector) = zero!(a.vec)
9292
SparseArrays.dropzeros!(a::CLILVector) = SparseArrays.dropzeros!(a.vec)
9393

94+
# Remove explicitly-stored zeros from each row, in place.
95+
function SparseArrays.dropzeros!(S::SparseMatrixCLIL)
96+
for r in eachindex(S.row_vals)
97+
cols = S.row_cols[r]
98+
vals = S.row_vals[r]
99+
j = 0
100+
for k in eachindex(vals)
101+
iszero(vals[k]) && continue
102+
j += 1
103+
cols[j] = cols[k]
104+
vals[j] = vals[k]
105+
end
106+
resize!(cols, j)
107+
resize!(vals, j)
108+
end
109+
return S
110+
end
111+
94112
struct NonZeros{T <: AbstractArray}
95113
v::T
96114
end

src/utils.jl

Lines changed: 7 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -341,10 +341,8 @@ function get_new_mm(
341341
push!(final_row_vals, new_row_val_i[indices[1]])
342342
for i in Iterators.drop(eachindex(indices), 1)
343343
col = new_row_col_i[indices[i]]
344-
# Compare against the last *retained* column rather than the previous
345-
# sorted entry: a prior cancellation may have `pop!`ed the matching
346-
# entry, in which case the remaining duplicates of `col` must start a
347-
# fresh accumulator instead of folding into an unrelated column.
344+
# Compare against the last retained column, not the previous sorted
345+
# entry: a prior cancellation may have `pop!`ed the matching entry.
348346
if !isempty(final_row_cols) && col == final_row_cols[end]
349347
final_row_vals[end] += new_row_val_i[indices[i]]
350348
if iszero(final_row_vals[end])
@@ -356,13 +354,16 @@ function get_new_mm(
356354
push!(final_row_vals, new_row_val_i[indices[i]])
357355
end
358356
end
359-
357+
360358
push!(new_row_cols, final_row_cols)
361359
push!(new_row_vals, final_row_vals)
362360
push!(new_nzrows, old_to_new_eq[eq])
363361
end
364362

365-
return typeof(mm)(new_nparentrows, count(!iszero, old_to_new_var), new_nzrows, new_row_cols, new_row_vals)
363+
# Drop explicit zeros (e.g. from an alias that cancelled to zero) so they
364+
# don't survive as phantom structural nonzeros.
365+
new_mm = typeof(mm)(new_nparentrows, count(!iszero, old_to_new_var), new_nzrows, new_row_cols, new_row_vals)
366+
return SparseArrays.dropzeros!(new_mm)
366367
end
367368

368369
struct BadMMAliasError <: Exception

test/runtests.jl

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -59,4 +59,16 @@ include("carpanzano_tearing.jl")
5959
@test mm2.row_cols == [[1, 2]]
6060
@test mm2.row_vals == [[7, 5]]
6161
@test mm2.ncols == 2
62+
63+
# An alias carrying an explicit stored zero (as `sparsevec` produces when
64+
# duplicate indices cancel) must not leak in as a phantom structural nonzero.
65+
# Row references var2 (removed, aliased to 0*var1) and var3 (retained); there
66+
# is no direct var1 term, so the only contribution to col 1 is the zero.
67+
mm = SSel.CLIL.SparseMatrixCLIL([0 1 1])
68+
old_to_new_eq = [1]
69+
old_to_new_var = [1, 0, 2]
70+
aliases = Dict(2 => sparsevec([1, 1], [2, -2], 3))
71+
mm2 = SSel.get_new_mm(aliases, old_to_new_eq, old_to_new_var, mm)
72+
@test mm2.row_cols == [[2]]
73+
@test mm2.row_vals == [[1]]
6274
end

0 commit comments

Comments
 (0)