# What is the state-of-the-art in parallel sparse matrix and dense vector multiplication?

For sure, there has been many highly optimized library on this. But I am working on a matrix-free context since the problem size does not allow explicit storage of sparse matrix elements. I'd love have a general knowledge on how people handle this sort of problem, especially on how to alleviate the memory bottleneck.

=========================update========================

As suggested by other commentators, I shall further elaborate on my situation.

My current implementation is below, which is just the basic one showing up in various textbooks.

void SpMVMtply(vector<det>& basis, ColVec& x, ColVec& b)
{

/*
"basis" is the orthonormal basis that spans the physical space
"x" ,"b" are just those in Ax = b
*/
int len = x.size();

#pragma omp parallel
{
#pragma omp for nowait
for (int i=0; i<len; i++)
{
vector<det> pool_det;
vector<cdouble> pool_amp;

/*
this function generates the non-zero off diagonal elements
with respect to the given basis vector. pool_det is the
collection of generated basis vectors, and pool_map collection
of the amplitude.
*/

OffDiagGen(basis[i], pool_det, pool_amp);

// this function generates the diagonal element
cdouble sum = x(i)* DiagCal(basis[i]);

/*
Below, basis_pos is an unordered_map with basis vector as key
and that basis vector's position in "basis" array as value.

*/
for (int j=0; j< pool_det.size(); j++)
sum += x(basis_pos[pool_det[j]]) * conj(pool_amp[j]);

b(i) = sum ;
}
}
}


The size of vector could reach 1 billion, so the explicit storage of sparse matrix is impossible. A rough profiling shows that the non-sequential accessing of memory in x(basis_pos[pool_det[j]]) is the most time-consuming part. On a Xeon Gold 6242 with 12 cores in use, it takes 6.2 seconds to such part for a vector with 3 million elements, while it takes only 1.5 seconds when that part is replaced by a constant. That's why I call it a memory bottleneck. My main concern is that if such bottleneck is inevitable, since even when the sparse matrix is explicitly stored, such non-sequential memory access also shows up. Based on my knowledge I haven't figured out a way to avoid it.

• It is a bit unclear what aspects you are interested in. By memory bottleneck, do you imply memory access or general "not enough memory" aspect? Any details on your matrix-free context and why do you still have memory bottleneck in the matrix-free solver? What paradigm of parallel computing are you interested in (distributed, shared, hybrid, GPU, FPGA)? Have you looked into any libraries that implement this functionality? Currently the question is extremely broad, I want to help you narrow it down. – Anton Menshov Apr 1 '20 at 20:01
• You want to read the papers by Martin Kronbichler! – Wolfgang Bangerth Apr 2 '20 at 1:13
• So you have a sparse matrix but not enough storage to hold the sparse matrix- you can only compute elements of the sparse matrix as needed, right? Can you iterate over all the nonzero entries in parallel in a reasonable way? Do you have enough storage to hold the vector input and the vector output? – Brian Borchers Apr 2 '20 at 3:55
• I have updated my question in order to address some ambiguities. – Izzy Vang Apr 2 '20 at 22:02
• @BrianBorchers Yes, only input and output vectors are stored explicitly, and there is a way to iterate over all non-zero elements. – Izzy Vang Apr 2 '20 at 22:12