Implementing Sparse-Matrix Dense-Vector Multiplication (SpMV) using a data-centric graph abstraction.

[neoblizz]

One cool, weird thing you can do with essentials/gunrock is implementing sparse-linear algebra using graphs. This is the opposite of what GraphBLAS offers. The core idea is that the advance operator, is essentially visiting all of your nonzeros of a sparse matrix. To formally define the sparse-matrix:

A sparse matrix or sparse array is a matrix in which most of the elements are zero.

In computation (and in fact in gunrock as well), graphs are actually internally represented using the same sparse formats as what you'd use to represent sparse matrices, in a sense, they are the sparse matrices. A lot of graph datasets also are available as sparse matrices in matrix-market format (or coordinate format). All of which is supported within essentials.

Therefore, to implement SpMV you simply have to advance from all nodes of the graph (all row indices), which will give you the following source (row index), destination (column index), weight (nonzero). The lambda-op will look something like this:

auto spmv = [=] __host__ __device__(
                             vertex_t const& source,     // ... source (row index)
                             vertex_t const& neighbor,   // neighbor (column index)
                             edge_t const& edge,         // edge (row ↦ column)
                             weight_t const& weight      // weight (nonzero).
                             ) -> bool {
      // y[row index] += weight * x[column index]
      math::atomic::add(&(y[source]), weight * x[neighbor]); 
      return false; // ignored.
    };

// Perform advance on the above lambda-op
operators::advance::execute<operators::load_balance_t::block_mapped,    // flavor of load-balancing schedule
                            operators::advance_direction_t::forward,    // direction (backward for transpose)
                            operators::advance_io_type_t::graph,        // entire graph as input
                            operators::advance_io_type_t::none>(        // no output frontier needed
      G, E, spmv, context);

Fin.

[neoblizz]

Avoiding the atomic

1. Use a load-balancing schedule that owns the entire row of the matrix (such as thread_mapped).
2. 2-step approach where we accumulate to temporary memory and then perform a key-value-reduction.

[neoblizz]

Complete code is available here, by @DanLoran and myself (@neoblizz).