Simple example for sparse matrix-vector multiplication?

[albertomercurio]

Hello,

is it possible to have a very minimal working example of a matrix-vector multiplication? I know there is the example on the documentation here, but I was wondering if that example could be simplified.

For example, the initialization of the vector

$$ \begin{pmatrix} 1 \\ 0 \\ 0 \\ 0 \\ -1 \end{pmatrix} $$

is performed with

IV = map(row_partition) do row_indices
    I,V = Int[], Float64[]
    for global_row in local_to_global(row_indices)
        if global_row == 1
            v = 1.0
        elseif global_row == n
            v = -1.0
        else
            continue
        end
        push!(I,global_row)
        push!(V,v)
    end
    I,V
end
II,VV = tuple_of_arrays(IV)
b = pvector(II,VV,row_partition) |> fetch

But can this be simplified with the following code?

IV = ([1, 5], [1, -1])
II,VV = tuple_of_arrays(IV)
b = pvector(II,VV,row_partition) |> fetch

Same for the matrix initialization.

Moreover, does the application of the mul! operation allocate memory? In principle, if the vector is dense, there should not be any allocations, but I don't know in this case. Can the vector be considered dense if I initialize also the zero entries?

IV = ([1, 2, 3, 4, 5], [1, 0, 0, 0, -1])
II,VV = tuple_of_arrays(IV)
b = pvector(II,VV,row_partition) |> fetch

Finally, it would be great to have a tutorial on how to promote a SparseMatrixCSC into a PSparseMatrix

[fverdugo]

Hi @albertomercurio

IV = ([1, 5], [1, -1])
II,VV = tuple_of_arrays(IV)
b = pvector(II,VV,row_partition) |> fetch

This works only if you are partitioning b into 2 parts.

This is how you would implement the SpMV:

muladd!(b,A,x) = mul!(b,A,x,1,1)
function myspmv!(b,A,x)
    t = PA.consistent!(x)
    map(mul!,PA.own_values(b),PA.own_own_values(A),PA.own_values(x))
    wait(t)
    map(muladd!,PA.own_values(b),PA.own_ghost_values(A),PA.ghost_values(x))
    b
end

In recent versions, this is allocation free.

To convert a sequential sparse matrix into a PSparseMatrix you would need a graph partitioner (e.g., metis). I agree that it would be a good idea to have an example on how to do this. Feel free to propose one!

If you want you learn more about the library, I am giving this seminar soon:

https://www.linkedin.com/posts/bernatfontgarcia_hpc-mpi-julialang-activity-7308470412051406850-RPW7?utm_source=share&utm_medium=member_desktop&rcm=ACoAABi6k4YByX9-JF8E5UevQLVRaNzTpnwwhlc

[albertomercurio]

Thank you for your reply.

I still don't understand how it can be non-allocating, as the sparsity of the b vector is not known a priori. Moreover, the documentation uses the IterativeSolvers.cg! function, which calls mul! internally. Does mul! exactly implement the myspmv! you wrote here?

Moreover, I see that the documentation uses the MPI backend, but is there a way to make it work also for Distributed.jl? This would be useful in cases where the users want to use pure Julia codes or helpers like ClusterManagers.jl or SlurmClusterManager.jl.

Sorry for the many questions. I'm very curious to use this package, in order to apply it also to simulations of large quantum systems using QuantumToolbox.jl.

[fverdugo]

b is dense. We work with dense vectors and sparse matrices. It is possible to work with sparse vectors, but this is not implemented.

Indeed, mul! is implemented like myspmv!

For the moment only MPI for parallel computing. MPI is widely supported in supercomputers.

[albertomercurio]

The IterativeSolvers code works, but if I try doing simply

x = similar(b,axes(A,2))
mul!(x, A, b)

I get the error

ERROR: AssertionError: matching_ghost_indices(axes(a, 2), axes(b, 1))
Stacktrace:
 [1] mul!(c::PVector{…}, a::PSparseMatrix{…}, b::PVector{…})
   @ PartitionedArrays ~/.julia/packages/PartitionedArrays/Tq4ld/src/p_sparse_matrix.jl:2093
 [2] top-level scope
   @ ~/GitHub/Research/Undef/PartitionedArrays Benchmark/partitionedarrays.jl:71
Some type information was truncated. Use `show(err)` to see complete types.

I'm wondering how the iterative solver can work without a working implementation of mul!. What am I missing?