Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Make sparse operations less dependent on inference #19611

Closed
wants to merge 2 commits into from
Closed

Make sparse operations less dependent on inference #19611

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
e1ad722
Make sparse operations less dependent on inference
pabloferz Dec 15, 2016
2a6a83d
Tests for #19611
pabloferz Dec 15, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
95 changes: 53 additions & 42 deletions base/sparse/sparsematrix.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1401,52 +1401,65 @@ sparse(S::UniformScaling, m::Integer, n::Integer=m) = speye_scaled(S.λ, m, n)
## map/map! and broadcast/broadcast! over sparse matrices

# map/map! entry points
function map!{Tf,N}(f::Tf, C::SparseMatrixCSC, A::SparseMatrixCSC, Bs::Vararg{SparseMatrixCSC,N})
function map!{F,N}(f::F, C::SparseMatrixCSC, A::SparseMatrixCSC, Bs::Vararg{SparseMatrixCSC,N})
_checksameshape(C, A, Bs...)
return _map!(f, C, A, Bs...)
end
@inline function _map!(f, C, A, Bs::Vararg)
fofzeros = f(_zeros_eltypes(A, Bs...)...)
fpreszeros = fofzeros == zero(fofzeros)
return fpreszeros ? _map_zeropres!(f, C, A, Bs...) :
_map_notzeropres!(f, fofzeros, C, A, Bs...)
end
function map{Tf,N}(f::Tf, A::SparseMatrixCSC, Bs::Vararg{SparseMatrixCSC,N})
function map{F,N}(f::F, A::SparseMatrixCSC, Bs::Vararg{SparseMatrixCSC,N})
_checksameshape(A, Bs...)
fofzeros = f(_zeros_eltypes(A, Bs...)...)
fpreszeros = fofzeros == zero(fofzeros)
maxnnzC = fpreszeros ? min(length(A), _sumnnzs(A, Bs...)) : length(A)
return _map(f, A, Bs...)
end
@inline function _map(f, A, Bs::Vararg)
entrytypeC = _broadcast_type(f, A, Bs...)
indextypeC = _promote_indtype(A, Bs...)
Ccolptr = Vector{indextypeC}(A.n + 1)
Crowval = Vector{indextypeC}(maxnnzC)
Cnzval = Vector{entrytypeC}(maxnnzC)
C = SparseMatrixCSC(A.m, A.n, Ccolptr, Crowval, Cnzval)
return fpreszeros ? _map_zeropres!(f, C, A, Bs...) :
_map_notzeropres!(f, fofzeros, C, A, Bs...)
if isleaftype(entrytypeC)
indextypeC = _promote_indtype(A, Bs...)
fofzeros = f(_zeros_eltypes(A, Bs...)...)
fpreszeros = fofzeros == zero(fofzeros)
maxnnzC = fpreszeros ? min(length(A), _sumnnzs(A, Bs...)) : length(A)
Ccolptr = Vector{indextypeC}(A.n + 1)
Crowval = Vector{indextypeC}(maxnnzC)
Cnzval = Vector{entrytypeC}(maxnnzC)
C = SparseMatrixCSC(A.m, A.n, Ccolptr, Crowval, Cnzval)
return fpreszeros ? _map_zeropres!(f, C, A, Bs...) :
_map_notzeropres!(f, fofzeros, C, A, Bs...)
end
return sparse(collect(Base.Generator(f, A, Bs...)))
end
# broadcast/broadcast! entry points
broadcast{Tf}(f::Tf, A::SparseMatrixCSC) = map(f, A)
broadcast!{Tf}(f::Tf, C::SparseMatrixCSC, A::SparseMatrixCSC) = map!(f, C, A)
function broadcast!{Tf,N}(f::Tf, C::SparseMatrixCSC, A::SparseMatrixCSC, Bs::Vararg{SparseMatrixCSC,N})
_aresameshape(C, A, Bs...) && return map!(f, C, A, Bs...) # could avoid a second dims check in map
broadcast{F}(f::F, A::SparseMatrixCSC) = map(f, A)
broadcast!{F}(f::F, C::SparseMatrixCSC, A::SparseMatrixCSC) = map!(f, C, A)
function broadcast!{F,N}(f::F, C::SparseMatrixCSC, A::SparseMatrixCSC, Bs::Vararg{SparseMatrixCSC,N})
_aresameshape(C, A, Bs...) && return _map!(f, C, A, Bs...)
Base.Broadcast.check_broadcast_indices(indices(C), A, Bs...)
fofzeros = f(_zeros_eltypes(A, Bs...)...)
fpreszeros = fofzeros == zero(fofzeros)
return fpreszeros ? _broadcast_zeropres!(f, C, A, Bs...) :
_broadcast_notzeropres!(f, fofzeros, C, A, Bs...)
end
function broadcast{Tf,N}(f::Tf, A::SparseMatrixCSC, Bs::Vararg{SparseMatrixCSC,N})
_aresameshape(A, Bs...) && return map(f, A, Bs...) # could avoid a second dims check in map
fofzeros = f(_zeros_eltypes(A, Bs...)...)
fpreszeros = fofzeros == zero(fofzeros)
indextypeC = _promote_indtype(A, Bs...)
function broadcast{F,N}(f::F, A::SparseMatrixCSC, Bs::Vararg{SparseMatrixCSC,N})
_aresameshape(A, Bs...) && return _map(f, A, Bs...)
entrytypeC = _broadcast_type(f, A, Bs...)
Cm, Cn = Base.to_shape(Base.Broadcast.broadcast_indices(A, Bs...))
maxnnzC = fpreszeros ? _checked_maxnnzbcres(Cm, Cn, A, Bs...) : (Cm * Cn)
Ccolptr = Vector{indextypeC}(Cn + 1)
Crowval = Vector{indextypeC}(maxnnzC)
Cnzval = Vector{entrytypeC}(maxnnzC)
C = SparseMatrixCSC(Cm, Cn, Ccolptr, Crowval, Cnzval)
return fpreszeros ? _broadcast_zeropres!(f, C, A, Bs...) :
_broadcast_notzeropres!(f, fofzeros, C, A, Bs...)
shape = Base.Broadcast.broadcast_indices(A, Bs...)
if isleaftype(entrytypeC)
indextypeC = _promote_indtype(A, Bs...)
fofzeros = f(_zeros_eltypes(A, Bs...)...)
fpreszeros = fofzeros == zero(fofzeros)
Cm, Cn = Base.to_shape(shape)
maxnnzC = fpreszeros ? _checked_maxnnzbcres(Cm, Cn, A, Bs...) : (Cm * Cn)
Ccolptr = Vector{indextypeC}(Cn + 1)
Crowval = Vector{indextypeC}(maxnnzC)
Cnzval = Vector{entrytypeC}(maxnnzC)
C = SparseMatrixCSC(Cm, Cn, Ccolptr, Crowval, Cnzval)
return fpreszeros ? _broadcast_zeropres!(f, C, A, Bs...) :
_broadcast_notzeropres!(f, fofzeros, C, A, Bs...)
end
return sparse(Base.Broadcast.broadcast_t(f, Any, shape, CartesianRange(shape), A, Bs...))
end
# map/map! and broadcast/broadcast! entry point helper functions
@inline _sumnnzs(A) = nnz(A)
Expand All @@ -1468,7 +1481,7 @@ _broadcast_type(f, As...) = Base._promote_op(f, Base.Broadcast.typestuple(As...)

# _map_zeropres!/_map_notzeropres! specialized for a single sparse matrix
"Stores only the nonzero entries of `map(f, Matrix(A))` in `C`."
function _map_zeropres!{Tf}(f::Tf, C::SparseMatrixCSC, A::SparseMatrixCSC)
function _map_zeropres!(f, C::SparseMatrixCSC, A::SparseMatrixCSC)
spaceC = min(length(C.rowval), length(C.nzval))
Ck = 1
@inbounds for j in 1:C.n
Expand All @@ -1491,7 +1504,7 @@ end
Densifies `C`, storing `fillvalue` in place of each unstored entry in `A` and
`f(A[i,j])` in place of each stored entry `A[i,j]` in `A`.
"""
function _map_notzeropres!{Tf}(f::Tf, fillvalue, C::SparseMatrixCSC, A::SparseMatrixCSC)
function _map_notzeropres!(f, fillvalue, C::SparseMatrixCSC, A::SparseMatrixCSC)
# Build dense matrix structure in C, expanding storage if necessary
_densestructure!(C)
# Populate values
Expand Down Expand Up @@ -2278,17 +2291,15 @@ round{To}(::Type{To}, A::SparseMatrixCSC) = round.(To, A)

## Binary arithmetic and boolean operators

# TODO: These seven functions should probably be reimplemented in terms of sparse map
# when a better sparse map exists. (And vectorized min, max, &, |, and xor should be
# deprecated in favor of compact-broadcast syntax.)
_checksameshape(A, B) = size(A) == size(B) || throw(DimensionMismatch("size(A) must match size(B)"))
(+)(A::SparseMatrixCSC, B::SparseMatrixCSC) = (_checksameshape(A, B); broadcast(+, A, B))
(-)(A::SparseMatrixCSC, B::SparseMatrixCSC) = (_checksameshape(A, B); broadcast(-, A, B))
min(A::SparseMatrixCSC, B::SparseMatrixCSC) = (_checksameshape(A, B); broadcast(min, A, B))
max(A::SparseMatrixCSC, B::SparseMatrixCSC) = (_checksameshape(A, B); broadcast(max, A, B))
(&)(A::SparseMatrixCSC, B::SparseMatrixCSC) = (_checksameshape(A, B); broadcast(&, A, B))
(|)(A::SparseMatrixCSC, B::SparseMatrixCSC) = (_checksameshape(A, B); broadcast(|, A, B))
xor(A::SparseMatrixCSC, B::SparseMatrixCSC) = (_checksameshape(A, B); broadcast(xor, A, B))
# TODO: Vectorized min, max, &, |, and xor should be deprecated in favor of
# compact-broadcast syntax.
(+)(A::SparseMatrixCSC, B::SparseMatrixCSC) = map(+, A, B)
(-)(A::SparseMatrixCSC, B::SparseMatrixCSC) = map(-, A, B)
min(A::SparseMatrixCSC, B::SparseMatrixCSC) = map(min, A, B)
max(A::SparseMatrixCSC, B::SparseMatrixCSC) = map(max, A, B)
(&)(A::SparseMatrixCSC, B::SparseMatrixCSC) = map(&, A, B)
(|)(A::SparseMatrixCSC, B::SparseMatrixCSC) = map(|, A, B)
xor(A::SparseMatrixCSC, B::SparseMatrixCSC) = map(xor, A, B)

(.+)(A::SparseMatrixCSC, B::Number) = Array(A) .+ B
( +)(A::SparseMatrixCSC, B::Array ) = Array(A) + B
Expand Down
5 changes: 5 additions & 0 deletions test/sparse/sparse.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1828,3 +1828,8 @@ let
@test_throws DimensionMismatch broadcast(+, A, B, speye(N))
@test_throws DimensionMismatch broadcast!(+, X, A, B, speye(N))
end

let A = sparse(Real[1 1])
@test (A + A)::SparseMatrixCSC{Int,Int} == [2 2] #19595
Copy link
Contributor

Choose a reason for hiding this comment

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

the test will be good to port over

Copy link
Member

Choose a reason for hiding this comment

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

Would you say we should verify the return type to have eltype Int? Presently

julia> Real[1 1] + Real[1 1]
1×2 Array{Real,2}:
 2  2

so maybe Real would be an acceptable return type, too?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Personally, I'd prefer the dense case to return an Array{Int,2}, but we'd need a revision of the behavior of arithmetic operation on dense arrays. So, checking for eltype(A + A) === Int here seems better to me.

Copy link
Member

Choose a reason for hiding this comment

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

See #19669.

@test_throws DimensionMismatch A + A'
end