3
$\begingroup$

I've written a C++ function that multiplies a sparse matrix (stored in CSR format) by a dense vector. Here's the code:

VectorXd csrMult(VectorXd x, vector<double> Adata, vector<int> Aindices, vector<int> Aindptr, int numRowsA) 
{
    VectorXd Ax = VectorXd::Zero(numRowsA); 
    for (int i = 0; i < numRowsA; i++)
    {
        for (int dataIdx = Aindptr[i]; dataIdx < Aindptr[i + 1]; dataIdx++)
        {
            Ax[i] += Adata[dataIdx] * x[Aindices[dataIdx]];
        }
    }

    return Ax;
}

Here VectorXd is a data type provided by the Eigen3 linear algebra library. The inputs Adata, Aindices, and Aindptr describe a matrix A in CSR format.

To be precise: Adata is a list of nonzero entries of A (stored in row major order), Aindices[i] tells us which column of A the nonzero entry Adata[i] belongs to, and Aindptr[i+1] - Aindptr[i] is the number of nonzero entries in the ith row of A.

I'm observing that Eigen3's sparse matrix-vector multiplication operation is about 5 times faster than my csrMult function, even when openMP is disabled. However, when I look at the source code SparseDenseProduct.h that I believe Eigen3 is using to compute this matrix-vector product, it's not clear to me how Eigen3 is faster.

Question: Do you have any suggestions to improve the speed of my code? Can you explain why Eigen3 is faster?


Edit 3: I noticed an important clue about what's going on. When I make the problem smaller by setting all but the top N entries of A equal to 0 (so that A becomes more sparse), my implementation compares better with Eigen3. In fact, if I set all but the top 10,000 nonzero entries of A equal to 0, my implementation beats Eigen3 by a factor of about 3. However, for the very large problem size I am interested in, Eigen3 beats my implementation by a factor of 2. I would still really like to tie with Eigen3 for large problem sizes.


Edit 2: I changed the inputs to pass by reference, as @TylerOlsen suggested, and now Eigen3 is only about twice as fast as my code (with openMP disabled). So that was a significant improvement. Here's the latest version of my code. I still need to figure out how to make my code twice as fast in order to tie with Eigen3.

void csrMult_v3(VectorXd& Ax, VectorXd& x, vector<double>& Adata, vector<int>& Aindices, vector<int>& Aindptr)
{
    // This code assumes that the size of Ax is numRowsA.
    for (int i = 0; i < Ax.size(); i++)
    {
        double Ax_i = 0.0;
        for (int dataIdx = Aindptr[i]; dataIdx < Aindptr[i + 1]; dataIdx++)
        {
            Ax_i += Adata[dataIdx] * x[Aindices[dataIdx]];
        }

        Ax[i] = Ax_i;
    }       
}

Edit 1: I also tried the following code, where the value Ax[i] is accumulated in a temporary variable, but the effect on runtime was negligible. I'm still observing Eigen3 is about 5 times faster (without openMP enabled).

VectorXd csrMult_v2(VectorXd x, vector<double> Adata, vector<int> Aindices, vector<int> Aindptr, int numRowsA)
{
    VectorXd Ax = VectorXd::Zero(numRowsA);
    for (int i = 0; i < numRowsA; i++)
    {
        double Ax_i = 0.0;
        for (int dataIdx = Aindptr[i]; dataIdx < Aindptr[i + 1]; dataIdx++)
        {
            Ax_i += Adata[dataIdx] * x[Aindices[dataIdx]];
        }

        Ax[i] = Ax_i;
    }

    return Ax;
}
$\endgroup$
15
  • $\begingroup$ Offhand, the immediate difference is that they accumulate their result for Ax[i] into a temporary variable (line 66 in their code). This lets you move the write to memory from the inner loop to a single write afterward and may let the compiler do more aggressive optimizations with loop unrolling, vectorization, etc. Since floating point arithmetic is not associative, the compiler is not allowed to do this optimization for you, in general. Try that and let us know how it goes. $\endgroup$ Oct 4, 2017 at 2:50
  • $\begingroup$ @TylerOlsen Thanks for the suggestion. I just tried that, but the effect on runtime was negligible. I edited the question to show the exact code that I tried. $\endgroup$
    – littleO
    Oct 4, 2017 at 3:37
  • 1
    $\begingroup$ @TylerOlsen I'm pretty sure the compiler is allowed that optimization, because floating point arithmetic doesn't need to be associative for it to be valid: it just needs to know that Ax[i] doesn't alias any other location that might be read from. The actual sequence of fp operations ends up being just the same. $\endgroup$
    – Kirill
    Oct 4, 2017 at 3:39
  • 1
    $\begingroup$ Is it possible the compiler manages to vectorize Eigen's code, but not yours? $\endgroup$
    – Kirill
    Oct 4, 2017 at 3:50
  • 6
    $\begingroup$ Can we first change all of those pass-by-value vector arguments to pass-by-reference first? Eating the copy overhead on the function call is potentially significant. No need to use an output parameter for Ax since we can rely on RVO to take care of that copy. In addition, make sure that you're passing "-march=native" as a command line argument. This will let the compiler do as much vectorizing as possible on your machine. $\endgroup$ Oct 4, 2017 at 3:53

1 Answer 1

1
$\begingroup$

After your second edit, you are already here:

for (int i = 0; i < Ax.size(); i++)
{
  double Ax_i = 0.0;
  for (int dataIdx = Aindptr[i]; dataIdx < Aindptr[i + 1]; dataIdx++)
  {
     Ax_i += Adata[dataIdx] * x[Aindices[dataIdx]];
  }

  Ax[i] = Ax_i;
}

But notice how you in row i+1 you are initializing dataIdx to the value Aindptr[i+1] it already has because you stopped the loop in row i when it had that value. So you can transform this to the following:

int dataIdx = Aindptr[0];
for (int i = 0; i < Ax.size(); i++)
{
  double Ax_i = 0.0;
  for (; dataIdx < Aindptr[i + 1]; dataIdx++)
  {
     Ax_i += Adata[dataIdx] * x[Aindices[dataIdx]];
  }

  Ax[i] = Ax_i;
}

The next piece of information you need to know is that you are walking the Aindices array linearly, from left to right because your indices dataIdx are incremented one by one. So you can save yourself the subscripting and do this:

int dataIdx = Aindptr[0];
int *Aindex = &Aindices[dataIdx];
for (int i = 0; i < Ax.size(); i++)
{
  double Ax_i = 0.0;
  for (; dataIdx < Aindptr[i + 1]; ++dataIdx,++Aindex)
  {
     Ax_i += Adata[dataIdx] * x[*Aindex];
  }

  Ax[i] = Ax_i;
}

Give this a try.

$\endgroup$
8
  • $\begingroup$ I've been searching for the "default" implementation of SparseMatrix< number >::vmult_add in the deal.II source, but I can't seem to find it. it should contain a loop similar to what you suggest above, right? could you point me into the right direction? $\endgroup$
    – GoHokies
    Oct 4, 2017 at 18:00
  • $\begingroup$ Thanks. I tried these optimizations but I only got a negligible improvement in runtime on smaller problems, and I didn't get an improvement on the larger problems that are most important in my application. $\endgroup$
    – littleO
    Oct 5, 2017 at 8:51
  • $\begingroup$ Actually, sometimes I observe a slight improvement with these optimizations for the large matrix size I've been testing with, but still Eigen3 is twice as fast as my code (when using the large sparse matrix). $\endgroup$
    – littleO
    Oct 5, 2017 at 11:07
  • $\begingroup$ The implementation in deal.II is here: github.com/dealii/dealii/blob/master/include/deal.II/lac/… $\endgroup$ Oct 5, 2017 at 13:39
  • $\begingroup$ I don't know what else to suggest. If the Eigen code is so much faster, why don't you use their implementation? $\endgroup$ Oct 5, 2017 at 13:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.