Add isequal for adjoint/transpose of sparse matrices

src/sparsematrix.jl

@@ -2378,41 +2378,52 @@ Base.isequal(A1::AbstractSparseMatrixCSC, A2::AbstractSparseMatrixCSC) = _iseq(i
 ## Explicit efficient comparisons with transposed arrays
 
 # Check whether all nonzero elements of A are equal to the respective elements in B
-function nzeq(A::AbstractSparseMatrixCSC, B::AbstractSparseMatrixCSCInclAdjointAndTranspose)
+# under the elementwise predicate `eq` (`==` or `isequal`)
+function nzeq(eq::F, A::AbstractSparseMatrixCSC, B::AbstractMatrix) where {F}
     @inbounds for j in axes(A,2)
         for k in nzrange(A, j)
             i = rowvals(A)[k]
             val = nonzeros(A)[k]
-            val ≠ B[i,j] && return false
+            eq(val, B[i,j]) || return false
         end
     end
     return true
 end
 # Peel off `Adjoint` and `Transpose` from first argument
-nzeq(A::Adjoint{<:Any,<:AbstractSparseMatrixCSCInclAdjointAndTranspose},
-     B::AbstractSparseMatrixCSCInclAdjointAndTranspose) =
-    nzeq(A', B')
-nzeq(A::Transpose{<:Any,<:AbstractSparseMatrixCSCInclAdjointAndTranspose},
-     B::AbstractSparseMatrixCSCInclAdjointAndTranspose) =
-    nzeq(transpose(A), transpose(B))
+# `B` may be a nested wrapper such as `Adjoint{<:Any,<:Transpose}` (from `A' == transpose(B)`),
+# hence the loose `AbstractMatrix` bound: `B` is only ever indexed
+nzeq(eq::F, A::Adjoint{<:Any,<:AbstractSparseMatrixCSCInclAdjointAndTranspose},
+     B::AbstractMatrix) where {F} =
+    nzeq(eq, A', B')
+nzeq(eq::F, A::Transpose{<:Any,<:AbstractSparseMatrixCSCInclAdjointAndTranspose},
+     B::AbstractMatrix) where {F} =
+    nzeq(eq, transpose(A), transpose(B))
 
 # Compare by walking both matrices
 # (We could further optimize the case `AbstractSparseMatrixCSC ==
 # Adjoint(Transpose(AbstractSparseMatrixCSC))` more efficiently, i.e.
 # the case where the RHS is both adjoint and transposed, i.e. where it
 # is in CSC format again.)
-function ==(A::AbstractSparseMatrixCSC,
-            B::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSCInclAdjointAndTranspose})
+function _iseq(eq::F, A::AbstractSparseMatrixCSC,
+               B::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSCInclAdjointAndTranspose}) where {F}
     # Different sizes are always different
     size(A) ≠ size(B) && return false
     # Compare nonzero elements
-    return nzeq(A, B) && nzeq(B, A)
+    return nzeq(eq, A, B) && nzeq(eq, B, A)
 end
+==(A::AbstractSparseMatrixCSC, B::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSCInclAdjointAndTranspose}) =
+    _iseq(==, A, B)
+Base.isequal(A::AbstractSparseMatrixCSC, B::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSCInclAdjointAndTranspose}) =
+    _iseq(isequal, A, B)
 # Peel off `Adjoint` and `Transpose` from first argument
 ==(A::Adjoint{<:Any,<:AbstractSparseMatrixCSCInclAdjointAndTranspose}, B::AbstractSparseMatrixCSCInclAdjointAndTranspose) =
     A' == B'
 ==(A::Transpose{<:Any,<:AbstractSparseMatrixCSCInclAdjointAndTranspose}, B::AbstractSparseMatrixCSCInclAdjointAndTranspose) =
     transpose(A) == transpose(B)
+Base.isequal(A::Adjoint{<:Any,<:AbstractSparseMatrixCSCInclAdjointAndTranspose}, B::AbstractSparseMatrixCSCInclAdjointAndTranspose) =
+    isequal(A', B')
+Base.isequal(A::Transpose{<:Any,<:AbstractSparseMatrixCSCInclAdjointAndTranspose}, B::AbstractSparseMatrixCSCInclAdjointAndTranspose) =
+    isequal(transpose(A), transpose(B))
 
 ## Reductions
 

test/sparsematrix_ops.jl

@@ -74,6 +74,41 @@ end
     end
 end
 
+@testset "isequal for adjoint/transpose of sparse matrices" begin
+    n = 10^5
+    A = spzeros(n, n); A[1, 1] = 1
+    B = copy(A)
+    for (L, R) in ((A', B'), (transpose(A), transpose(B)), (A, B'), (A', B),
+                   (A, transpose(B)), (transpose(A), B), (A', transpose(B)))
+        @test isequal(L, R)
+        @test @elapsed(isequal(L, R)) < 0.1
+    end
+    A[1, 2] = 1; B[2, 1] = 1
+    @test isequal(A, B') && isequal(A', B) && !isequal(A', B') && !isequal(A, B)
+    @test !isequal(spzeros(2, 3)', spzeros(2, 3))
+    # adjoint vs transpose of a complex matrix nests wrappers (`Adjoint{<:Any,<:Transpose}`)
+    C = sparse([1, 2], [2, 3], [1.0im, 2.0], 3, 3)
+    @test C' == transpose(conj(C)) && isequal(C', transpose(conj(C)))
+    @test C' != transpose(C) && !isequal(C', transpose(C))
+    # isequal semantics for NaN, signed zeros and conjugation must match dense arrays
+    X = sparse([1, 2, 3, 1], [1, 1, 2, 3], [NaN, -0.0, 2.0, 1.0im], 3, 3)
+    for Y in (sparse([1, 2, 3, 1], [1, 1, 2, 3], [NaN, -0.0, 2.0, 1.0im], 3, 3),
+              sparse([1, 2, 3, 1], [1, 1, 2, 3], [NaN, -0.0, 2.0, -1.0im], 3, 3),
+              sparse([1, 1, 2, 3], [1, 2, 3, 1], [NaN, -0.0, 2.0, 1.0im], 3, 3),
+              sparse([1, 1, 2, 3], [1, 2, 3, 1], [NaN, -0.0, 2.0, -1.0im], 3, 3),
+              sparse([1, 3, 1], [1, 2, 3], [NaN, 2.0, 1.0im], 3, 3),
+              sparse([1, 2, 3, 1], [1, 1, 2, 3], [NaN, 0.0, 2.0, 1.0im], 3, 3),
+              sparse([1, 2, 3, 1, 3], [1, 1, 2, 3, 3], [NaN, -0.0, 2.0, 1.0im, 0.0], 3, 3),
+              sparse([1, 2, 3, 1], [1, 1, 2, 3], [1.0, -0.0, 2.0, 1.0im], 3, 3))
+        for (L, R) in ((X', Y'), (transpose(X), transpose(Y)), (X, Y'), (X', Y),
+                       (X, transpose(Y)), (transpose(X), Y), (X', transpose(Y)))
+            @test isequal(L, R) == isequal(Matrix(L), Matrix(R))
+            @test isequal(R, L) == isequal(Matrix(R), Matrix(L))
+            @test (L == R) == (Matrix(L) == Matrix(R))
+        end
+    end
+end
+
 @testset "iszero specialization for SparseMatrixCSC" begin
     @test !iszero(sparse(I, 3, 3))                  # test failure
     @test iszero(spzeros(3, 3))                     # test success with no stored entries