SparseMatrixCSC as a general sparse container: Setindex!(...) not working for SparseMatrixCSC with vector elements (v0.7)

HTSykora · September 4, 2018, 10:28am

In my application I have a sparse matrix, but every element of this matrix is a discretised function. If I try to use SparseMatrixCSC{Array{Float64,1},Int64} as the container, but setindex!(...) in 0.7 throws a stackoverflow:

julia> sa=sparse([1],[1],[Vector{Float64}(undef,3)],3,3)
3×3 SparseMatrixCSC{Array{Float64,1},Int64} with 1 stored entry:
  [1, 1]  =  [6.93279e-310, 6.93279e-310, 6.93279e-310]

julia> sa[1,1] = [1.0,2.0,3.0]
StackOverflowError:
Stacktrace:
 [1] setindex!(::SparseMatrixCSC{Array{Int64,1},Int64}, ::Array{Int64,1}, ::Int64, ::Int64) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v0.7/SparseArrays/src/sparsematrix.jl:2540 (repeats 80000 times)

This can be solved by using broadcasted setindex!(...), but it does not solve the problem adding new elements to the array:

julia> sa[1,1] .= [1.0,2.0,3.0]
3-element Array{Float64,1}:
 1.0
 2.0
 3.0

julia> sa[2,2] = [1.0,2.0,3.0]
StackOverflowError:
Stacktrace:
 [1] setindex!(::SparseMatrixCSC{Array{Float64,1},Int64}, ::Array{Float64,1}, ::Int64, ::Int64) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v0.7/SparseArrays/src/sparsematrix.jl:2540 (repeats 80000 times)

I can work around it using custom structures, but this has worked in 0.6.x. Is it intentional, or a bug?

pablosanjose · September 4, 2018, 10:44am

Wow, good catch. I’m pretty sure this is not intentional. I would wait a bit for some pro feedback, and if it’s not clarified I’d definitely file an issue

Sukera · September 4, 2018, 10:58am

I don’t think it’s intentional.

github.com

JuliaLang/julia/blob/c8450d862f0e2653011f68118daecfe12b398c90/stdlib/SparseArrays/src/sparsematrix.jl#L2538


      
                  deleteat!(rowvalA, rowidx:nnzA)
                  deleteat!(nzvalA, rowidx:nnzA)
              end
              return A
          end
          
          
# Nonscalar A[I,J] = B: Convert B to a SparseMatrixCSC of the appropriate shape first
          _to_same_csc(::SparseMatrixCSC{Tv, Ti}, V::AbstractMatrix, I...) where {Tv,Ti} = convert(SparseMatrixCSC{Tv,Ti}, V)
          _to_same_csc(::SparseMatrixCSC{Tv, Ti}, V::AbstractVector, I...) where {Tv,Ti} = convert(SparseMatrixCSC{Tv,Ti}, reshape(V, map(length, I)))
          
          
setindex!(A::SparseMatrixCSC{Tv}, B::AbstractVecOrMat, I::Integer, J::Integer) where {Tv} = setindex!(A, convert(Tv, B), I, J)
          function setindex!(A::SparseMatrixCSC{Tv,Ti}, V::AbstractVecOrMat, Ix::Union{Integer, AbstractVector{<:Integer}, Colon}, Jx::Union{Integer, AbstractVector{<:Integer}, Colon}) where {Tv,Ti<:Integer}
              @assert !has_offset_axes(A, V, Ix, Jx)
              (I, J) = Base.ensure_indexable(to_indices(A, (Ix, Jx)))
              checkbounds(A, I, J)
              B = _to_same_csc(A, V, I, J)
          
          
    issortedI = issorted(I)
              issortedJ = issorted(J)
          
          
    if !issortedI && !issortedJ

That line is probably calling itself again even after a conversion if B is already a Tv, since what you put in is in my opinion always going to be <:AbstractVecOrMat, which in itself is just const AbstractVecOrMat{T} = Union{AbstractVector{T}, AbstractMatrix{T}}. I think an issue is appropriate.

Sukera · September 4, 2018, 11:20am

@HTSykora please also reference this discourse thread in the issue - I don’t have access to my github account right now, so I can’t add that reference. In particular the snippet from the source I linked is probably going to be useful.

HTSykora · September 4, 2018, 2:39pm

Of course, I’m on it

danspielman · November 25, 2018, 10:01pm

Is there any word on when a fix for this will appear in Julia?
I am facing the same problem, but with a sparse array of StaticArrays, so the broadcast solution will not work for me.