Store a zero in an sparse matrix

Hello,

do you know how can I store a zero in an sparse matrix? it seems that it automatically neglects the zero allocation: A[1,1]=0.0 is neglected. I know is far from efficient, but the need of it is to debug and replicate an old FORTAN code I am working with.

Thank you.

The point of sparse matrices is not to store zeros, so it’s surprising to me that you’d get a difference based on whether fortran sets an explicit zero - maybe fortran is storing explicitly set zeros? Can you maybe show some more info about what you’re trying to replicate?

1 Like

It turns out that getindex(::SparseMatrixCSC, v, i,j) first checks whether A[i,j] is structurally nonzero, and if it isn’t then it turns it into a structural nonzero only if v is nonzero, see https://github.com/JuliaLang/julia/blob/4f3a89e3abecb96bbe221d027278c55dccf23da7/stdlib/SparseArrays/src/sparsematrix.jl#L2634-L2665.

For debugging purposes, you can work around this by first making A[i,j] nonzero and then zeroing it out again:

julia> A = spzeros(3,3);

julia> A[2,2] = 0.0;  # This doesn't work

julia> A
3Γ—3 SparseMatrixCSC{Float64, Int64} with 0 stored entries:
  β‹…    β‹…    β‹… 
  β‹…    β‹…    β‹… 
  β‹…    β‹…    β‹… 

julia> A[2,2] = 1.0;  # This will make A[2,2] structurally nonzero

julia> A
3Γ—3 SparseMatrixCSC{Float64, Int64} with 1 stored entry:
  β‹…    β‹…    β‹… 
  β‹…   1.0   β‹… 
  β‹…    β‹…    β‹… 

julia> A[2,2] = 0.0;  # Now we get where we wanted

julia> A
3Γ—3 SparseMatrixCSC{Float64, Int64} with 1 stored entry:
  β‹…    β‹…    β‹… 
  β‹…   0.0   β‹… 
  β‹…    β‹…    β‹… 

To add my unsolicited five cents, I’m not sure whether special-casing v == 0 in setindex() is good design. But it is what it is and is highly unlikely to change anytime soon.

3 Likes

I am aware that is the whole point of a sparse matrix, but in this case after some coordinate transformations I have as a side product some zeros. My intention is to replicate and compare some sparsity maps as I have in the FORTRAN code.

Thanks for the idea!

You might use eps() as your structural nonzero.