Skip to content

Commit

Permalink
Merge pull request #17265 from Sacha0/unarybcastspmat
Browse files Browse the repository at this point in the history 
De-eval-ify vectorized unary functions over SparseMatrixCSCs, and then transition them to compact broadcast syntax
  • Loading branch information
@tkelman
tkelman authored Aug 27, 2016
2 parents c20a8bc + bd7da66 commit 379c248
Show file tree
Hide file tree
Showing 3 changed files with 155 additions and 145 deletions.
13 changes: 13 additions & 0 deletions base/deprecated.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -852,4 +852,17 @@ function convert(::Type{UpperTriangular}, A::Bidiagonal)
end
end

# Deprecate vectorized unary functions over sparse matrices in favor of compact broadcast syntax (#17265).
for f in (:sin, :sinh, :sind, :asin, :asinh, :asind,
:tan, :tanh, :tand, :atan, :atanh, :atand,
:sinpi, :cosc, :ceil, :floor, :trunc, :round, :real, :imag,
:log1p, :expm1, :abs, :abs2, :conj,
:log, :log2, :log10, :exp, :exp2, :exp10, :sinc, :cospi,
:cos, :cosh, :cosd, :acos, :acosd,
:cot, :coth, :cotd, :acot, :acotd,
:sec, :sech, :secd, :asech,
:csc, :csch, :cscd, :acsch)
@eval @deprecate $f(A::SparseMatrixCSC) $f.(A)
end

# End deprecations scheduled for 0.6
220 changes: 102 additions & 118 deletions base/sparse/sparsematrix.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1305,134 +1305,118 @@ end

## Unary arithmetic and boolean operators

macro _unary_op_nz2z_z2z(op,A,Tv,Ti)
esc(quote
nfilledA = nnz($A)
colptrB = Array{$Ti}($A.n+1)
rowvalB = Array{$Ti}(nfilledA)
nzvalB = Array{$Tv}(nfilledA)

nzvalA = $A.nzval
colptrA = $A.colptr
rowvalA = $A.rowval

k = 0 # number of additional zeros introduced by op(A)
@inbounds for i = 1 : $A.n
colptrB[i] = colptrA[i] - k
for j = colptrA[i] : colptrA[i+1]-1
opAj = $(op)(nzvalA[j])
if opAj == 0
k += 1
else
rowvalB[j - k] = rowvalA[j]
nzvalB[j - k] = opAj
end
end
end
colptrB[end] = $A.colptr[end] - k
deleteat!(rowvalB, colptrB[end]:nfilledA)
deleteat!(nzvalB, colptrB[end]:nfilledA)
return SparseMatrixCSC($A.m, $A.n, colptrB, rowvalB, nzvalB)
end) # quote
end

# Operations that may map nonzeros to zero, and zero to zero
# Result is sparse
for op in (:ceil, :floor, :trunc, :round,
:sin, :tan, :asin, :atan,
:sinh, :tanh, :asinh, :atanh,
:sinpi, :cosc,
:sind, :tand, :asind, :atand)
@eval begin
$(op){Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}) = @_unary_op_nz2z_z2z($op,A,Tv,Ti)
end # quote
end # macro

for op in (:real, :imag)
@eval begin
($op){Tv<:Complex,Ti}(A::SparseMatrixCSC{Tv,Ti}) = @_unary_op_nz2z_z2z($op,A,Tv.parameters[1],Ti)
end # quote
end # macro
real{Tv<:Number,Ti}(A::SparseMatrixCSC{Tv,Ti}) = copy(A)
imag{Tv<:Number,Ti}(A::SparseMatrixCSC{Tv,Ti}) = spzeros(Tv, Ti, A.m, A.n)

for op in (:ceil, :floor, :trunc, :round)
@eval begin
($op){T,Tv,Ti}(::Type{T},A::SparseMatrixCSC{Tv,Ti}) = @_unary_op_nz2z_z2z($op,A,T,Ti)
end # quote
end # macro


# Operations that map nonzeros to nonzeros, and zeros to zeros
# Result is sparse
for op in (:-, :log1p, :expm1)
@eval begin

function ($op)(A::SparseMatrixCSC)
B = similar(A)
nzvalB = B.nzval
nzvalA = A.nzval
@simd for i=1:length(nzvalB)
@inbounds nzvalB[i] = ($op)(nzvalA[i])
end
return B
end

end
end

function abs{Tv<:Complex,Ti}(A::SparseMatrixCSC{Tv,Ti})
T = Tv.parameters[1]
(T <: Integer) && (T = (T <: BigInt) ? BigFloat : Float64)
@_unary_op_nz2z_z2z(abs,A,T,Ti)
end
abs2{Tv<:Complex,Ti}(A::SparseMatrixCSC{Tv,Ti}) = @_unary_op_nz2z_z2z(abs2,A,Tv.parameters[1],Ti)
for op in (:abs, :abs2)
@eval begin
function ($op){Tv<:Number,Ti}(A::SparseMatrixCSC{Tv,Ti})
B = similar(A)
nzvalB = B.nzval
nzvalA = A.nzval
@simd for i=1:length(nzvalB)
@inbounds nzvalB[i] = ($op)(nzvalA[i])
"""
Helper macro for the unary broadcast definitions below. Takes parent method `fp` and a set
of desired child methods `fcs`, and builds an expression defining each of the child methods
such that `broadcast(::typeof(fc), A::SparseMatrixCSC) = fp(fc, A)`.
"""
macro _enumerate_childmethods(fp, fcs...)
fcexps = Expr(:block)
for fc in fcs
push!(fcexps.args, :( broadcast(::typeof($(esc(fc))), A::SparseMatrixCSC) = $(esc(fp))($(esc(fc)), A) ) )
end
return fcexps
end

# Operations that map zeros to zeros and may map nonzeros to zeros, yielding a sparse matrix
"""
Takes unary function `f` that maps zeros to zeros and may map nonzeros to zeros, and returns
a new `SparseMatrixCSC{TiA,TvB}` `B` generated by applying `f` to each nonzero entry in
`A` and retaining only the resulting nonzeros.
"""
function _broadcast_unary_nz2z_z2z_T{TvA,TiA,TvB}(f::Function, A::SparseMatrixCSC{TvA,TiA}, ::Type{TvB})
Bcolptr = Array{TiA}(A.n + 1)
Browval = Array{TiA}(nnz(A))
Bnzval = Array{TvB}(nnz(A))
Bk = 1
@inbounds for j in 1:A.n
Bcolptr[j] = Bk
for Ak in nzrange(A, j)
x = f(A.nzval[Ak])
if x != 0
Browval[Bk] = A.rowval[Ak]
Bnzval[Bk] = x
Bk += 1
end
return B
end
end
end

Bcolptr[A.n + 1] = Bk
resize!(Browval, Bk - 1)
resize!(Bnzval, Bk - 1)
return SparseMatrixCSC(A.m, A.n, Bcolptr, Browval, Bnzval)
end
function _broadcast_unary_nz2z_z2z{Tv}(f::Function, A::SparseMatrixCSC{Tv})
_broadcast_unary_nz2z_z2z_T(f, A, Tv)
end
@_enumerate_childmethods(_broadcast_unary_nz2z_z2z,
sin, sinh, sind, asin, asinh, asind,
tan, tanh, tand, atan, atanh, atand,
sinpi, cosc, ceil, floor, trunc, round)
broadcast(::typeof(real), A::SparseMatrixCSC) = copy(A)
broadcast{Tv,Ti}(::typeof(imag), A::SparseMatrixCSC{Tv,Ti}) = spzeros(Tv, Ti, A.m, A.n)
broadcast{TTv}(::typeof(real), A::SparseMatrixCSC{Complex{TTv}}) = _broadcast_unary_nz2z_z2z_T(real, A, TTv)
broadcast{TTv}(::typeof(imag), A::SparseMatrixCSC{Complex{TTv}}) = _broadcast_unary_nz2z_z2z_T(imag, A, TTv)
ceil{To}(::Type{To}, A::SparseMatrixCSC) = _broadcast_unary_nz2z_z2z_T(ceil, A, To)
floor{To}(::Type{To}, A::SparseMatrixCSC) = _broadcast_unary_nz2z_z2z_T(floor, A, To)
trunc{To}(::Type{To}, A::SparseMatrixCSC) = _broadcast_unary_nz2z_z2z_T(trunc, A, To)
round{To}(::Type{To}, A::SparseMatrixCSC) = _broadcast_unary_nz2z_z2z_T(round, A, To)

# Operations that map zeros to zeros and map nonzeros to nonzeros, yielding a sparse matrix
"""
Takes unary function `f` that maps zeros to zeros and nonzeros to nonzeros, and returns a
new `SparseMatrixCSC{TiA,TvB}` `B` generated by applying `f` to each nonzero entry in `A`.
"""
function _broadcast_unary_nz2nz_z2z_T{TvA,TiA,TvB}(f::Function, A::SparseMatrixCSC{TvA,TiA}, ::Type{TvB})
Bcolptr = Vector{TiA}(A.n + 1)
Browval = Vector{TiA}(nnz(A))
Bnzval = Vector{TvB}(nnz(A))
copy!(Bcolptr, 1, A.colptr, 1, A.n + 1)
copy!(Browval, 1, A.rowval, 1, nnz(A))
@inbounds @simd for k in 1:nnz(A)
Bnzval[k] = f(A.nzval[k])
end
return SparseMatrixCSC(A.m, A.n, Bcolptr, Browval, Bnzval)
end
function _broadcast_unary_nz2nz_z2z{Tv}(f::Function, A::SparseMatrixCSC{Tv})
_broadcast_unary_nz2nz_z2z_T(f, A, Tv)
end
@_enumerate_childmethods(_broadcast_unary_nz2nz_z2z,
log1p, expm1, abs, abs2, conj)
broadcast{TTv}(::typeof(abs2), A::SparseMatrixCSC{Complex{TTv}}) = _broadcast_unary_nz2nz_z2z_T(abs2, A, TTv)
broadcast{TTv}(::typeof(abs), A::SparseMatrixCSC{Complex{TTv}}) = _broadcast_unary_nz2nz_z2z_T(abs, A, TTv)
broadcast{TTv<:Integer}(::typeof(abs), A::SparseMatrixCSC{Complex{TTv}}) = _broadcast_unary_nz2nz_z2z_T(abs, A, Float64)
broadcast{TTv<:BigInt}(::typeof(abs), A::SparseMatrixCSC{Complex{TTv}}) = _broadcast_unary_nz2nz_z2z_T(abs, A, BigFloat)
function conj!(A::SparseMatrixCSC)
nzvalA = A.nzval
@simd for i=1:length(nzvalA)
@inbounds nzvalA[i] = conj(nzvalA[i])
@inbounds @simd for k in 1:nnz(A)
A.nzval[k] = conj(A.nzval[k])
end
return A
end

conj(A::SparseMatrixCSC) = conj!(copy(A))

# Operations that map nonzeros to nonzeros, and zeros to nonzeros
# Result is dense
for op in (:cos, :cosh, :acos, :sec, :csc, :cot, :acot, :sech,
:csch, :coth, :asech, :acsch, :cospi, :sinc, :cosd,
:cotd, :cscd, :secd, :acosd, :acotd, :log, :log2, :log10,
:exp, :exp2, :exp10)
@eval begin

function ($op){Tv}(A::SparseMatrixCSC{Tv})
B = fill($(op)(zero(Tv)), size(A))
@inbounds for col = 1 : A.n
for j = A.colptr[col] : A.colptr[col+1]-1
row = A.rowval[j]
nz = A.nzval[j]
B[row,col] = $(op)(nz)
end
end
return B
# Operations that map both zeros and nonzeros to zeros, yielding a dense matrix
"""
Takes unary function `f` that maps both zeros and nonzeros to nonzeros, constructs a new
`Matrix{TvB}` `B`, populates all entries of `B` with the result of a single `f(one(zero(Tv)))`
call, and then, for each stored entry in `A`, calls `f` on the entry's value and stores the
result in the corresponding location in `B`.
"""
function _broadcast_unary_nz2nz_z2nz{Tv}(f::Function, A::SparseMatrixCSC{Tv})
B = fill(f(zero(Tv)), size(A))
@inbounds for j in 1:A.n
for k in nzrange(A, j)
i = A.rowval[k]
x = A.nzval[k]
B[i,j] = f(x)
end

end
return B
end
@_enumerate_childmethods(_broadcast_unary_nz2nz_z2nz,
log, log2, log10, exp, exp2, exp10, sinc, cospi,
cos, cosh, cosd, acos, acosd,
cot, coth, cotd, acot, acotd,
sec, sech, secd, asech,
csc, csch, cscd, acsch)


## Broadcasting kernels specialized for returning a SparseMatrixCSC
Expand Down Expand Up @@ -1734,7 +1718,7 @@ end # macro
(.^)(A::SparseMatrixCSC, B::Number) =
B==0 ? sparse(ones(typeof(one(eltype(A)).^B), A.m, A.n)) :
SparseMatrixCSC(A.m, A.n, copy(A.colptr), copy(A.rowval), A.nzval .^ B)
(.^)(::Irrational{:e}, B::SparseMatrixCSC) = exp(B)
(.^)(::Irrational{:e}, B::SparseMatrixCSC) = exp.(B)
(.^)(A::Number, B::SparseMatrixCSC) = (.^)(A, full(B))
(.^)(A::SparseMatrixCSC, B::Array) = (.^)(full(A), B)
(.^)(A::Array, B::SparseMatrixCSC) = (.^)(A, full(B))
Expand Down
67 changes: 40 additions & 27 deletions test/sparsedir/sparse.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -286,7 +286,7 @@ end

# conj
cA = sprandn(5,5,0.2) + im*sprandn(5,5,0.2)
@test full(conj(cA)) == conj(full(cA))
@test full(conj.(cA)) == conj(full(cA))

# Test SparseMatrixCSC [c]transpose[!] and permute[!] methods
let smalldim = 5, largedim = 10, nzprob = 0.4
Expand Down Expand Up @@ -461,22 +461,47 @@ end
@test maximum(sparse(-ones(3,3))) == -1
@test minimum(sparse(ones(3,3))) == 1

# Unary functions
a = sprand(5,15, 0.5)
afull = full(a)
for op in (:sin, :cos, :tan, :ceil, :floor, :abs, :abs2)
@eval begin
@test ($op)(afull) == full($(op)(a))
end
end

for op in (:ceil, :floor)
@eval begin
@test ($op)(Int,afull) == full($(op)(Int,a))
# Test unary functions with specialized broadcast over SparseMatrixCSCs
let
A = sprand(5, 15, 0.5)
C = A + im*A
Afull = full(A)
Cfull = full(C)
# Test representatives of [unary functions that map zeros to zeros and may map nonzeros to zeros]
@test sin.(Afull) == full(sin.(A))
@test tan.(Afull) == full(tan.(A)) # should be redundant with sin test
@test ceil.(Afull) == full(ceil.(A))
@test floor.(Afull) == full(floor.(A)) # should be redundant with ceil test
@test real.(Afull) == full(real.(A))
@test imag.(Afull) == full(imag.(A))
@test real.(Cfull) == full(real.(C))
@test imag.(Cfull) == full(imag.(C))
# Test representatives of [unary functions that map zeros to zeros and nonzeros to nonzeros]
@test expm1.(Afull) == full(expm1.(A))
@test abs.(Afull) == full(abs.(A))
@test abs2.(Afull) == full(abs2.(A))
@test abs.(Cfull) == full(abs.(C))
@test abs2.(Cfull) == full(abs2.(C))
# Test representatives of [unary functions that map both zeros and nonzeros to nonzeros]
@test cos.(Afull) == full(cos.(A))
# Test representatives of remaining vectorized-nonbroadcast unary functions
@test ceil(Int, Afull) == full(ceil(Int, A))
@test floor(Int, Afull) == full(floor(Int, A))
# Tests of real, imag, abs, and abs2 for SparseMatrixCSC{Int,X}s previously elsewhere
for T in (Int, Float16, Float32, Float64, BigInt, BigFloat)
R = rand(T[1:100;], 2, 2)
I = rand(T[1:100;], 2, 2)
D = R + I*im
S = sparse(D)
@test R == real.(S)
@test I == imag.(S)
@test real.(sparse(R)) == R
@test nnz(imag.(sparse(R))) == 0
@test abs.(S) == abs(D)
@test abs2.(S) == abs2(D)
end
end


# getindex tests
ni = 23
nj = 32
Expand Down Expand Up @@ -872,7 +897,7 @@ end
@test_throws ArgumentError sparsevec(Dict(-1=>1,1=>2))

# issue #8976
@test conj(sparse([1im])) == sparse(conj([1im]))
@test conj.(sparse([1im])) == sparse(conj([1im]))
@test conj!(sparse([1im])) == sparse(conj!([1im]))

# issue #9525
Expand Down Expand Up @@ -1038,18 +1063,6 @@ end
x = speye(100)
@test_throws BoundsError x[-10:10]

for T in (Int, Float16, Float32, Float64, BigInt, BigFloat)
let R=rand(T[1:100;],2,2), I=rand(T[1:100;],2,2)
D = R + I*im
S = sparse(D)
@test R == real(S)
@test I == imag(S)
@test real(sparse(R)) == R
@test nnz(imag(sparse(R))) == 0
@test abs(S) == abs(D)
@test abs2(S) == abs2(D)
end
end

# issue #10407
@test maximum(spzeros(5, 5)) == 0.0
Expand Down

2 comments on commit 379c248

@nanosoldier
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Executing the daily benchmark build, I will reply here when finished:

@nanosoldier runbenchmarks(ALL, isdaily = true)

@nanosoldier
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Your benchmark job has completed - possible performance regressions were detected. A full report can be found here. cc @jrevels

Please sign in to comment.