Sparse matrix interfaces

[j-fu]

Discuss (missing) API exports / public labels and possible further functionality.

Currently, we have an implementation of SparseMatrixCSC and a de facto public API for AbstractSparseMatrixCSC which is not documented as such but used e.g. by LinearSolve.jl. #395 proposes to fix this situation, see also #52, #265.

Moreover, there are various discussions on SparseMatrixCSR implementations and related packages. To prevent rank growth it may make sense to have at least a baseline for this in SparseArrays, see #546, #538. May be this could go into some other package.

There seems to be a wish for more general interfaces, #167 goes into this direction.

Interfaces

AbstractSparseMatrixCSC <: AbstractSparseMatrix

• Exported now: nnz, rowvals, nonzeros, nzrange, dropzeros
• To be exported: AbstractSparseMatrixCSC, getcolptr
• Also to be mandatory: SparseMatrixCSC(::AbstractSparseMatrixCSC)

AbstractSparseMatrixCSR <: AbstractSparseMatrix (See #538, #546)

• To be exported: SparseMatrixCSR, AbstractSparseMatrixCSR, getrowptr, colvals, (+ nnz, nonzeros, nzrange, dropzeros)
• Also to be mandatory: SparseMatrixCSR(::AbstractSparseMatrixCSR)

AbstractSparseMatrix <: AbstractMatrix

• Exported now: AbstractSparseMatrix
• iteratenz (?, See add the iteratenz API #167)

AbstractMatrix:

• size, getindex

[j-fu]

• Rename fields of SparseMatrixCSC #52 (comment) argues that
  SparseMatrixCSC content should be addressed by accessor functions. So getcolptr/getrowptr would have
  to be exported. That discussion got stuck on the names of the accessor functions.

• IMHO the current naming of fields and accessors is not very consistent, but I also think it is
  too late for bikeshedding here.

• sparse: should return "something sparse" is argued in
  AbstractSparseMatrixCSC not exported #265 (comment) so we might keep the semantics that it has now - returning a SparseMatrixCSC. I think it is not needed as part of an AbstractSparseMatrixCSC/CSR interface.

[ChrisRackauckas]

For an AbstractSparseMatrix, what are the overloads that are required to be defined, and what are the rules on those functions?

[ViralBShah]

Are we discussing AbstractSparseMatrix or AbstractSparseMatrixCSC?

[j-fu]

In that sense, both - just to not to miss the larger picture.

For a general, implementation independent AbstractSparseMatrix interface my feeling is that more thinking is needed which might be worth to pursue, but needs time.

OTHOH a de-facto AbstractSparseMatrixCSC interface exists, and is used by other packages including LinearSolve.jl without being "official API" due to AbstractSparseMatrixCSC and getcolptr not exported. Other packages (e.g. for preconditioners) may also take advantage of this, but I would feel inconvenient to advertise this.

Disclaimer: I maintain "ExtendableSparse.jl" (which i am currently overhauling) which allows for easy sparse matrix build just using setindex! via an internal linked list structure. It exports an ExtendableSparseMatrixCSC <:AbstractSparseMatrixCSC, and in the moment of accessing it via the AbstractSparseMatrixCSC interface, the internal linked list structure is transformed to a SparseMatrixCSC structure.

[j-fu]

• Base.size, Base.getindex (from AbstractMatrix)
• nnz, rowvals, nonzeros, nzrange
• AbstractSparseMatrixCSC, getcolptr
• SparseMatrixCSC(A::AbstractSparseMatrixCSC) = SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), nonzeros(A))

Providing these, #395 is essentially about, in the moment, the export of
AbstractSparseMatrixCSC, getcolptr is missing, and SparseMatrixCSC(A::AbstractSparseMatrixCSC)=copy(A)

As far as I can tell, implementations of linear algebra in the package is already realized on the
level of the AbstractSparseMatrixCSC interface, so we would have default overloads already.
This appears to be also true for dropzeros! and findnz.
I would assume that this is comprehensive, but may be this should be checked.

I think some additional unit test with a mock MyMatrix<: AbstractMatrixCSC could make sense.

[ChrisRackauckas]

How well does this work for CUDA sparse matrices?

[SobhanMP]

cusparse supports lots of other format but for CuSparseArrayCSC, CUDA.jl already implement these

SparseArrays.nnz(g::AbstractCuSparseArray) = g.nnz
SparseArrays.nonzeros(g::AbstractCuSparseArray) = g.nzVal

SparseArrays.nonzeroinds(g::AbstractCuSparseVector) = g.iPtr

SparseArrays.rowvals(g::CuSparseMatrixCSC) = g.rowVal
SparseArrays.getcolptr(S::CuSparseMatrixCSC) = S.colPtr
function SparseArrays.findnz(S::MT) where {MT <: AbstractCuSparseMatrix} ...

dropzeros! does not seem to be supported by the julia package nor cuSparse.

[j-fu]

Yes, checked this, and SparseMatrixCSC(::CuSparseMatrixCSC) works and returns a SparseMatrixCSC.

[ChrisRackauckas]

We need to make it a bit more definite as an interface, like:

SparseArrays.nnz(::AbstractSparseMatrixCSC)::Int
SparseArrays.nonzeros(::AbstractSparseMatrixCSC{T})::AbstractVector{T}

Etc., we should be clear on some of the types we're requiring too on the inputs and outputs.

[ViralBShah]

@SobhanMP has some nice work on sparse iterators in https://github.com/SobhanMP/SparseExtra.jl/blob/main/src/iternz.jl