Make length(A.nzval)==nnz(A)

[andreasnoack]

We currently allow that length(A.nzval)>nnz(A) such that memory for future elements in the matrix can be preallocated. Instead nnz(A)=A.colptr[n+1]-1. However, as pointed out in JuliaLang/julia#30435 (comment), this could simply be handled by sizehint! and having two different mechanisms for ensuring extras allocations seems superfluous. Hence, it might be worth simplifying the sparse matrix code and enforce that nnz(A)==length(A.nzval).

[Sacha0]

Ref. JuliaLang/julia#20464 for previous discussion of this topic. Decoupling of buffer length and number of stored entries is a well established property of compressed sparse data structures, and is regularly exploited in sparse linear algebra libraries; ref. JuliaLang/julia#20464 (comment) particularly for expansion on use cases. Best!

[andreasnoack]

Decoupling of buffer length and number of stored entries is a well established property of compressed sparse data structures

The point is that this is already possible because the buffer size of a Vector is decoupled from the length of the Vector.

[StefanKarpinski]

The most that seems necessary is the ability to ask a Vector how much underlying storage it has, which we can surely do in a stdlib. We'll need an official API for that if we ever move SparseArrays out of stdlib, but until then, we can violate abstraction boundaries as necessary.

[Sacha0]

The point is that this is already possible because the buffer size of a Vector is decoupled from the length of the Vector.

Yes, had we control over the buffer size of a Vector, that would do the trick. But last I checked we do not have such control? Best!

[Sacha0]

To clarify, last I checked sizehint!'s contract was that sizehint! is just that --- a hint, not a guarantee. Is that no longer the case?

[abraunst]

To clarify, last I checked sizehint!'s contract was that sizehint! is just that --- a hint, not a guarantee. Is that no longer the case?

If I read the code correctly in jl_array_grow_at_end in src/array.c, you are right, it is not a guarantee: it seems to always allocate at least the requested amount, but sometimes more. But resize! calls the very same function so no guarantees either...

[Sacha0]

Perhaps we would benefit from an API that provides a stronger contract than that of sizehint!, implementation not being contract. (Tossing the some-time-ago-discussed buffersize[!] name into the ring.)

Another consideration mentioned in JuliaLang/julia#20464 (comment) (that I forgot when responding over the weekend) is symbolic stages of sparse matrix operations. There the traditional decoupling is used to save roughly a factor of two of storage by keeping nzval empty while working with a populated pattern (colptr and rowval ). Best!

[martinholters]

I just hit a bug where we blindly map an operation over all values of nzval, even though they might be invalid:

julia> using SparseArrays

julia> # create a sparse matrix with extra entries in nzval
       B = SparseMatrixCSC(Complex{BigInt}[1 2]')'[1:1, 2:2]
1×1 SparseMatrixCSC{Complex{BigInt},Int64} with 1 stored entry:
  [1, 1]  =  2+0im

julia> dump(B)
SparseMatrixCSC{Complex{BigInt},Int64}
  m: Int64 1
  n: Int64 1
  colptr: Array{Int64}((2,)) [1, 2]
  rowval: Array{Int64}((3,)) [1, 140548820313728, 9]
  nzval: Array{Complex{BigInt}}((3,)) Complex{BigInt}[2 + 0im, #undef, #undef]

julia> -B
ERROR: UndefRefError: access to undefined reference
Stacktrace:
 [1] getindex at ./array.jl:728 [inlined]
 [2] iterate at ./array.jl:704 [inlined]
 [3] iterate at ./generator.jl:44 [inlined]
 [4] collect_to!(::Array{Complex{BigInt},1}, ::Base.Generator{Array{Complex{BigInt},1},typeof(-)}, ::Int64, ::Int64) at ./array.jl:651
 [5] collect_to_with_first! at ./array.jl:630 [inlined]
 [6] _collect(::Array{Complex{BigInt},1}, ::Base.Generator{Array{Complex{BigInt},1},typeof(-)}, ::Base.EltypeUnknown, ::Base.HasShape{1}) at ./array.jl:624
 [7] map at ./array.jl:548 [inlined]
 [8] -(::SparseMatrixCSC{Complex{BigInt},Int64}) at ./stdlib/v1.2/SparseArrays/src/sparsematrix.jl:1507
 [9] top-level scope at REPL[4]:1

julia> conj(B)
ERROR: UndefRefError: access to undefined reference
Stacktrace:
 [1] getindex at ./array.jl:728 [inlined]
 [2] _broadcast_getindex at ./broadcast.jl:590 [inlined]
 [3] _getindex at ./broadcast.jl:621 [inlined]
 [4] _broadcast_getindex at ./broadcast.jl:596 [inlined]
 [5] getindex at ./broadcast.jl:557 [inlined]
 [6] macro expansion at ./broadcast.jl:887 [inlined]
 [7] macro expansion at ./simdloop.jl:73 [inlined]
 [8] copyto! at ./broadcast.jl:886 [inlined]
 [9] copyto! at ./broadcast.jl:841 [inlined]
 [10] copy at ./broadcast.jl:817 [inlined]
 [11] materialize at ./broadcast.jl:797 [inlined]
 [12] broadcast(::typeof(conj), ::Array{Complex{BigInt},1}) at ./broadcast.jl:751
 [13] conj at ./arraymath.jl:30 [inlined]
 [14] conj(::SparseMatrixCSC{Complex{BigInt},Int64}) at ./stdlib/v1.2/SparseArrays/src/sparsematrix.jl:1510
 [15] top-level scope at REPL[5]:1

Definitions at https://github.com/JuliaLang/julia/blob/de48421024439b873c7616c607e3310a8f24c2d3/stdlib/SparseArrays/src/sparsematrix.jl#L1506-L1511

The existence of these bugs makes me think that the general idea of not having extra entries in nzval is a good one. Short term, this particular bug is an easy fix, but I'm not 100% clear whether it is currently preferred to retain the extra entries in the result (just not applying any operation to them) or to truncate nzval to nnz. If someone clarifies, I'll prepare a PR.

[KristofferC]

Fixed by JuliaLang/julia#40523 I think