Bug when comparing and affecting SparseMatrixCSC

[Lecrapouille]

First bug

in sparsematrix.jl concerning

function ==(A1::AbstractSparseMatrixCSC, A2::AbstractSparseMatrixCSC)

please replace all valsX[jY]!=0 by !iszero Indeed when redefining zero() the wrong zero is not correctly take into account.

Here the fix:

function ==(A1::AbstractSparseMatrixCSC, A2::AbstractSparseMatrixCSC)
    size(A1) != size(A2) && return false
    vals1, vals2 = nonzeros(A1), nonzeros(A2)
    rows1, rows2 = rowvals(A1), rowvals(A2)
    m, n = size(A1)
    @inbounds for i = 1:n
        nz1,nz2 = nzrange(A1,i), nzrange(A2,i)
        j1,j2 = first(nz1), first(nz2)
        # step through the rows of both matrices at once:
        while j1 <= last(nz1) && j2 <= last(nz2)
            r1,r2 = rows1[j1], rows2[j2]
            if r1==r2
                vals1[j1]!=vals2[j2] && return false
                j1+=1
                j2+=1
            else
                if r1<r2
                    !iszero(vals1[j1]) && return false
                    j1+=1
                else
                    !iszero(vals2[j2]) && return false
                    j2+=1
                end
            end
        end
        # finish off any left-overs:
        for j = j1:last(nz1)
            !iszero(vals1[j]) && return false
        end
        for j = j2:last(nz2)
            !iszero(vals2[j]) && return false
        end
    end
    return true
end

Here some checks (only made on Julia 1.0.3):

using LinearAlgebra, SparseArrays, Printf

# function ==(A1::AbstractSparseMatrixCSC, A2::AbstractSparseMatrixCSC)
function FIX(A1::SparseMatrixCSC, A2::SparseMatrixCSC)
           size(A1) != size(A2) && return false
           vals1, vals2 = nonzeros(A1), nonzeros(A2)
           rows1, rows2 = rowvals(A1), rowvals(A2)
           m, n = size(A1)

           @inbounds for i = 1:n
               nz1,nz2 = nzrange(A1,i), nzrange(A2,i)
               j1,j2 = first(nz1), first(nz2)

               # step through the rows of both matrices at once:
               while j1 <= last(nz1) && j2 <= last(nz2)
                   r1,r2 = rows1[j1], rows2[j2]
                   if r1==r2
                       vals1[j1]!=vals2[j2] && return false
                       j1+=1
                       j2+=1
                   else
                       if r1<r2
                           !iszero(vals1[j1]) && return false
                           j1+=1
                       else
                           !iszero(vals2[j2]) && return false
                           j2+=1
                       end
                   end
               end
               # finish off any left-overs:
               for j = j1:last(nz1)
                   !iszero(vals1[j]) && return false
               end
               for j = j2:last(nz2)
                   !iszero(vals2[j]) && return false
               end
           end
           return true
       end


# New algebra redefining zero() !!
struct MP{T} <: Number λ::T end
Base.zero(::Type{MP{T}}) where T = MP(typemin(T))
Base.zero(x::MP{T})      where T = zero(typeof(x))
Base.one(::Type{MP{T}})  where T = MP(zero(T))
Base.one(x::MP{T})       where T = one(typeof(x))
Base.promote_rule(::Type{MP{T}}, ::Type{U}) where {T, U} = MP{T}
Base.convert(::MP{T}, x) where T = MP(T(x))
MP(x::MP) = MP(x.λ)
MP(S::SparseMatrixCSC{T,U}; preserve=true) where {T, U} = convert(SparseMatrixCSC{MP{T},U}, S)
mpzeros(::Type{T}, m::Int64, n::Int64) where T = spzeros(MP{T}, m, n)


# New algebra with false zeros
A = MP(sparse([1, 2], [1, 2], [0.0, 0.0]))
B = mpzeros(Float64, 2,2)
# A contains two MP(0.0) and B is empty
A == B # returns true which is incorrect
FIX(A, B) # returns false which is expected

# New algebra with real zeros
mp0 = zero(MP{Float64})
AA = sparse([1, 2], [1, 2], [mp0, mp0])
BB = mpzeros(Float64, 2,2)
# A contains two MP(-Inf) and B is empty
AA == BB # returns false which is incorrect
FIX(AA, BB) # returns true which is expected

# Normal algebra
AAA = sparse([1, 2], [1, 2], [0.0, 0.0])
BBB = spzeros(Float64, 2,2)
AAA == BBB # returns true which is correct
FIX(AAA, BBB) # returns true which is expected

Second bug

The = of Sparse is also impacted (but I could not locate the function).

# Inserting fake zero
B[1,1] = MP(0.0)
# B is still empty. This is an error because shall have MP(0.0) which is not equivalent to zero()

 # Inserting real zero
B[1,1] = mp0
# B is not empty. This is an error because mp0 is the zero()

[dkarrasch]

Would you mind preparing a PR? The function for the second issue is setindex!. Interestingly, it does call !iszero on the value to be set, to decide whether it should do anything or not:

https://github.com/JuliaLang/julia/blob/c5cc3f2811804b9d5ef8cd0b584ba2e7ad8f67c7/stdlib/SparseArrays/src/sparsematrix.jl#L2443-L2476

[abraunst]

I think this is (at least partially) fixed in master, see JuliaLang/julia#30580

[dkarrasch]

Nice, and I see it's actually included in the v1.3 release candidate.

[Lecrapouille]

Yes I can create a PR and I'll check with version 1.3
I also relate to this ticket older JuliaLang/julia#24790

[Lecrapouille]

@abraunst @dkarrasch when I opened this ticket I only tested on Julia 1.0 and 1.2 and the setindex! was effectively fixed on the master. I've created a PR for the operator==

[Lecrapouille]

I up this ticket since this bug is old and still present in Julia 1.6 and could be quickly fixed. There is confusion in sparse matrix between stored literal 0 and stored zero() so vals1[j1]!=0 shall be replaced by !iszero(vals1[j1]) idem for vals2.

At this time I made a PR https://github.com/Lecrapouille/julia/commit/af7430fd45393cc3051b5dc53ca95dde635b7201 but now the function is in function ==(A1::AbstractSparseMatrixCSC, A2::AbstractSparseMatrixCSC) in file julia/stdlib/SparseArrays/src/sparsematrix.jl

[Lecrapouille]

I did not program with Julia for a long time. I did not check back, but it seems to me the current the code is now correct by this commit #227 So I close this ticket.