Advantages/disadvantages with PSparseMatrix <: AbstractMatrix

[fredrikekre]

Currently PSparseMatrix <: AbstractMatrix and PVector <: AbstractVector (

PartitionedArrays.jl/src/p_sparse_matrix.jl

Line 971 in b211949

struct PSparseMatrix{V,B,C,D,T} <: AbstractMatrix{T}

,

PartitionedArrays.jl/src/p_vector.jl

Line 324 in b211949

struct PVector{V,A,B,C,T} <: AbstractVector{T}

).

I was wondering if this was always the case and if there was any discussion around it, concretely:

Why not subtype AbstractSparseMatrix (https://github.com/JuliaSparse/SparseArrays.jl/blob/313a04f4a78bbc534f89b6b4d9c598453e2af17c/src/abstractsparse.jl#L33)?

What are the advantages? You get some things like eltype for free and such but anything "major"?

What are the disadvantages? Does the subtyping cause any issues with hitting generic methods (for example)? I guess since indexing is not supported generic code paths will error eventually anyway?

The reason I am asking is because currently in HYPRE I don't subtype, see e.g. https://github.com/fredrikekre/HYPRE.jl/blob/10bdaec7230a90666c0672f74b2587cca90e8797/src/HYPRE.jl#L52 but I am considering if it is actually useful to do so.

[fverdugo]

Pros:

• One can use PSparseMatrix and PVector in any existing functions for ::AbstractArray as long as these functions do not index the matrices and vectors and only consume them via linear algebra operations (for which PSparseMatrix and PVector define methods). This is the case for many solvers, e.g. solvers from IterativeSolvers.jl. We are able to use the conjugate gradient solver of this package, even though it is not aware of PartitionedArrays.

cons:

• PSparseMatrix and PVector do not fully implement the AbstractArray interface for performance reasons. However they could provide an inefficient implementation of getindex/setindex!. Similar to CuArrays.

[fredrikekre]

Thanks, I don't think the pros apply to the HYPREMatrix because there is no operations (other than solving) exposed from the library. The main advantage I see though is that it would allow HYPRE matrices to pass through methods with ::AbstractMatrix to e.g. be stored in a struct or something.

I think a slow getindex/setindex! could be implemented for HYPRE too (in fact, I do it in the tests) so perhaps that could be exposed in some debug mode.

[fverdugo]

Makes sense. In PetscCall.jl I give up exposing the petsc matrix and vector to the user. My solution was to only expose solver functions that take Julia standard arrays or PSparseMatrix and PVector objects as input.

[fredrikekre]

Okay, yea HYPRE supports PSparseMatrix, SparseMatrixCSC and SparseMatrixCSR as input, but also allows constructing the HYPREMatrix directly which can save a bit of memory (and some time).

[fverdugo]

BTW, when using HYPRE.jl with PAs do you also support DebugArray, or only MPIArray? I made the PetscCall.jl wrappers also to work with DebugArray (by collecting the pieces, building a sequential matrix and passing this one to petsc).

[fredrikekre]

I only use the generic map(...) routines so I believe it should work with both backends. The tests are not super extensive but at least it tests it https://github.com/fredrikekre/HYPRE.jl/blob/10bdaec7230a90666c0672f74b2587cca90e8797/test/runtests.jl#L96.