diff --git a/src/abstractsparse.jl b/src/abstractsparse.jl index d23eebca..962759a0 100644 --- a/src/abstractsparse.jl +++ b/src/abstractsparse.jl @@ -114,6 +114,8 @@ end Return a tuple `(I, J, V)` where `I` and `J` are the row and column indices of the stored ("structurally non-zero") values in sparse matrix `A`, and `V` is a vector of the values. +`A` may also be the adjoint or transpose of a sparse matrix or vector, in which case the +values in `V` are correspondingly adjointed or transposed. # Examples ```jldoctest diff --git a/src/sparsematrix.jl b/src/sparsematrix.jl index 56ee216a..3f927991 100644 --- a/src/sparsematrix.jl +++ b/src/sparsematrix.jl @@ -1958,6 +1958,11 @@ function findnz(S::AbstractSparseMatrixCSC{Tv,Ti}) where {Tv,Ti} return (I, J, V) end +# Materializing the (conjugate) transpose is linear in `nnz` and yields the +# indices in column-major order of the wrapped matrix, consistent with the +# `AbstractSparseMatrixCSC` method above. +findnz(S::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC}) = findnz(copy(S)) + function _sparse_findnextnz(m::AbstractSparseMatrixCSC, ij::CartesianIndex{2}) row, col = Tuple(ij) col > size(m, 2) && return nothing diff --git a/src/sparsevector.jl b/src/sparsevector.jl index d6f83fd6..7f522d7e 100644 --- a/src/sparsevector.jl +++ b/src/sparsevector.jl @@ -924,6 +924,25 @@ function findnz(x::SparseVectorUnion) return (I, V) end +function findnz(x::AdjOrTransSparseVectorUnion) + p = parent(x) + numnz = nnz(p) + I = ones(indtype(p), numnz) + J = Vector{indtype(p)}(undef, numnz) + V = Vector{eltype(x)}(undef, numnz) + + nzind = nonzeroinds(p) + nzval = nonzeros(p) + f = wrapperop(x) + + @inbounds for i = 1 : numnz + J[i] = nzind[i] + V[i] = f(nzval[i]) + end + + return (I, J, V) +end + function _sparse_findnextnz(v::AbstractCompressedVector, i::Integer) n = searchsortedfirst(nonzeroinds(v), i) if n > length(nonzeroinds(v)) diff --git a/test/sparsematrix_ops.jl b/test/sparsematrix_ops.jl index ab2065cb..3db2acb8 100644 --- a/test/sparsematrix_ops.jl +++ b/test/sparsematrix_ops.jl @@ -29,6 +29,27 @@ end @test nnz(zero(sparse(fill(1,5,5)))) == 0 end +@testset "findnz for adjoint/transpose (issue #632)" begin + A = sparse([1, 1, 2, 3], [1, 2, 3, 2], [1.0+2.0im, 3.0, 4.0-1.0im, 0.0], 3, 4) + for T in (Float64, ComplexF64), op in (adjoint, transpose) + B = op(T == Float64 ? real(A) : A) + I, J, V = findnz(B) + @test (I, J, V) == findnz(SparseMatrixCSC(B)) + @test issorted(collect(zip(J, I))) # column-major order of the wrapper + @test all(B[i, j] == v for (i, j, v) in zip(I, J, V)) + @test length(I) == nnz(B) + @test typeof(I) == typeof(J) == Vector{Int} && eltype(V) == T + @test all(isempty, findnz(op(spzeros(T, 2, 3)))) + end + x = sparsevec([2, 4], [1.0+im, 0.0], 5) + for op in (adjoint, transpose) + I, J, V = findnz(op(x)) + @test I == [1, 1] && J == [2, 4] && V == op.([1.0+im, 0.0]) + @test (I, J, V) == findnz(sparse(op(x))) + @test all(isempty, findnz(op(spzeros(3)))) + 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