[WIP] generic sparse indices

[matbesancon]

This tackles https://github.com/JuliaLang/julia/issues/31330
Linked issue: https://github.com/JuliaLang/julia/issues/25118

Will solve in a cleaner way: https://github.com/JuliaLang/julia/issues/39887

Design goal:
a generic index iterator for structural non-zeros of abstract arrays. This will allow efficient generic algorithms exploiting sparsity.

The function is defined in Base, so it can be implemented by the types in both LinearAlgebra and SparseArrays, I found it awkward to have it defined in LinearAlgebra itself

ping @dkarrasch

[dkarrasch]

I'm curious how you're going to implement the iterator for *Triangulars and *Diagonals. I was playing around with a little predicate-based draft. It turned out to improve on the current state, but wasn't nearly as fast as your latest Diagonal contribution.

[matbesancon]

I'm curious how you're going to implement the iterator for *Triangulars and *Diagonals.

I think some indices might require allocating? If you have some funky array where there non-zero pattern is well-defined but not trivially implementable with ranges or CartesianIndices.

The compromise would be:
--> non-allocating indices struct: more zeros
--> possibly arbitrary indices container: true structural non-zeros

It turned out to improve on the current state, but wasn't nearly as fast as your latest Diagonal contribution

I don't think we will do as fast as specialized methods for end-functions with specific types (like dot on Diagonal). However, this can still be a better base than having to rely on eachindex

[matbesancon]

@dkarrasch looking at SparseArrays.findnz, I think it could be ok for eachstoredindex to allocate a vector of non-zeros in the case where a smart type is not possible IMO

[matbesancon]

Not sure why SparseArrays does not depend on Printf, it seems like it should?

[dkarrasch]

@dkarrasch looking at SparseArrays.findnz, I think it could be ok for eachstoredindex to allocate a vector of non-zeros in the case where a smart type is not possible IMO

I guess it depends on where you want to use it. I think we rarely, if ever, call findnz internally, but rather work with the fields directly and thereby avoid allocations. In any case, for sparse matrices, I think you want to consider introducing a

getindex(A::AbstractSparseMatrixCSC, I::CartesianIndex{2}) = getindex(A, I[1], I[2])

method.

[matbesancon]

I would be ok with adding the method, but it already has the corresponding behaviour implemented with the current fallback:

julia> using SparseArrays

julia> m = spzeros(200,100);

julia> m[CartesianIndex(1,1)]
0.0

why do you have it in mind?

[dkarrasch]

why do you have it in mind?

Because I saw that we have such a method specifically for tuples of integers, so it seemed we don't rely on fallbacks there either. But otherwise no other reason.

[mbauman]

Do we want to define this as being tied to structural nonzeros? Or simply stored values? We could simultaneously introduce a unstoredvalue to give room for other sorts of structures.

In any case, I like this and think it's very needed.

[matbesancon]

Good question. What would be expected for a FillArrays.Ones? In the case of stored values, it would return an empty iterator? I am not sure we can make this useful in general, except for constructors?

[DilumAluthge]

@matbesancon It looks like this PR touches seven files. Two of those files (base/abstractarray.jl and NEWS.md) are in Base, one of those files (stdlib/LinearAlgebra/src/diagonal.jl) is in the LinearAlgebra stdlib, and the remaining four files are in the SparseArrays stdlib.

We have moved the SparseArrays stdlib to an external repository. Therefore, for the portion of this PR that modifies the SparseArrays stdlib, can you please:

1. Remove the SparseArrays changes from this PR.
2. Open a new PR against https://github.com/JuliaLang/SparseArrays.jl with the SparseArrays changes.

Thank you!