Dividing sparse matrix by vector produces fill in

[gerw]

Code says more than thousand words 😄

using SparseArrays

function test(m = 3, n = 2)
	p = .2
	A = sprandn(m, n, p)
	display("text/plain", A)

	x = randn(m)
	B = A .* x
	display("text/plain", B)

	C = B ./ x
	display("text/plain", C)

	nothing
end

This produces

3×2 SparseMatrixCSC{Float64, Int64} with 1 stored entry:
  ⋅         ⋅ 
 0.186969   ⋅ 
  ⋅         ⋅ 
3×2 SparseMatrixCSC{Float64, Int64} with 1 stored entry:
  ⋅         ⋅ 
 0.409903   ⋅ 
  ⋅         ⋅ 
3×2 SparseMatrixCSC{Float64, Int64} with 6 stored entries:
  0.0        0.0
  0.186969   0.0
 -0.0       -0.0

Note that the multiplication with x is done correctly, but the division adds zero entries everywhere. I would expect that C has the same sparsity structure as A.

[ViralBShah]

Yes, we probably need to squeeze out the zeros after ./.

[gerw]

I think the problem lies in https://github.com/JuliaSparse/SparseArrays.jl/blob/main/src/higherorderfns.jl#L214. When computing A ./ x, this line evaluates 0.0 / 0.0 (which results in NaN). Since NaN != 0.0, the routine _broadcast_notzeropres! is called on line 221.

For A .* x, we get 0.0 * 0.0 == 0.0 and _broadcast_zeropres! is called.