A fast way to perform sparse matrix-free product with a vector

Dear all,

I was wondering if it is possible to increase the performance of a matrix-free product of a sparse matrix defined by 3 vectors (rows, columns and values) with another vector {B}. That is:

A(l,m)=v and performing {C}=[A]*{B}

So far I have the following strategy:

function C = mfree_times(l,m,v,B)
     aux = B(m);
     prod = v.*aux;
     C = accumarray(l,prod);
end

My intentions are to increase this performance as much as possible. The bottleneck is majorly given by accumarray function. I am already using l, m and v as gpuArrays. Please share your thoughts. Thanks!

[Joss Knight]

Can you not use the 3-argument sparse constructor?

[Paulo Ribeiro]

I think that is a quick option. However this a very large sparse matrix and the sparse constructor demands a considerable processing time. The matrix must be reassembled constantly and I’m trying to avoid this step. Any thoughts? Thanks!

[James Tursa]

Is the vector B full or sparse? What is changing at each iteration? All of l, m, v, and B? Are the column indexes m really 1-based as you seem to be using them?

[Paulo Ribeiro]

B is set as full and can be used as a constant during the iterations (conjugate gradient method). This is also true for l,m and v. The first two are the nonzero positions (lines and columns) of values v. These are not unique positions and duplicates are expected. As proposed by Joss the sparse constructor may be used, but it is generally time consuming for these large matrices.

Once converged, these quantities must be recalculated. My intentions are to provide a matrix-free approach to the conventional pcg solver.

[James Tursa]

Also, are the l & m & v vectors really 3-tuples of row & column & values rather than CSR format vectors? Are they sorted in a particular order?

[Paulo Ribeiro]

Yes. These are column vectors and not sorted by any order. Like this:

l=[1 1 2 3 1 2]'
m=[1 2 2 3 1 2]'
v=[0.5 1.3 4.5 2.1 -2.3 1.0]'

defining an array of order 3.

[Bruno Luong]

The only reason accumarray is slow because it creates full matrix, might be quite large depending on max(l) even for small-length vector.

You might try using sparse storage (using SPARSE or ACCUMARRAY with SPARSE option)

[Paulo Ribeiro]

Thanks Bruno. I tried with:

   C = accumarray(l,prod,[],[],[],true);

and with the "false" flag also. Since {C} is a full vector the performance using "true" was worse. And it makes no difference with the "false" flag or no flag at all.

The accumarray performs quite fast on the GPU but over so many iterations it becomes a bottleneck.

[James Tursa]

@Paulo: So, based on your replies, your variables appear to have nothing to do with CSR format sparse matrices ... you simply have 3-tuples of row,column,value to use. My answer below is based on that assumption.

[Paulo Ribeiro]

James, Thanks. That’s right. The question is now updated.

Answer 1

[James Tursa]

Assuming the vectors are 3-tuples of row, column, and value, and that all of the variables involved are full double vectors, here is a naive mex routine to do the multiplication and produce a full vector result. If you are working with R2018a and compiling with the -R2018a option then replace the mxGetPr calls with mxGetDoubles calls.

/*   C = mfree_timesx(r,c,v,B) does the equivalent of:
 *
 *   C = sparse(r,c,v,numel(B),numel(B)) * B
 *
 * where
 *
 *   r,c,v are 3-tuples of values from an MxM matrix:
 *     r = row index (1-based)
 *     c = column index (1-based)
 *     v = values
 *   B = column vector of size Mx1
 *   C = column vector of size Mx1
 *
 * It is assumed that the indexing contained in the r and c vectors are
 * all within an MxM square matrix (i.e., are integers between 1 and M).
 * (The code does *not* check for this. A violation will bomb this routine)
 *
 * Duplicate indexing within the r and c vectors is allowed, and is handled
 * by simply adding the results together (same as what would happen if you
 * did C = sparse(r,c,v,numel(B),numel(B)) * B).
 *
 * All inputs must be full real double vectors.
 *
 * The output is a full real double column vector.
 *
 * Programmer: James Tursa
 *
 */

/* Includes ----------------------------------------------------------- */

#include "mex.h"

/* Gateway ------------------------------------------------------------ */

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    mwSize i, n, m;
    double *Ar, *Ac, *Av, *B, *C;

      if( nrhs != 4 ) {
          mexErrMsgTxt("Need exactly four double input vectors");
      }
      if( nlhs > 1 ) {
          mexErrMsgTxt("Too many outputs");
      }
      if( !(mxIsDouble(prhs[0]) && mxIsDouble(prhs[1]) && mxIsDouble(prhs[2]) && mxIsDouble(prhs[3])) ||
            mxIsComplex(prhs[0]) || mxIsComplex(prhs[1]) || mxIsComplex(prhs[2]) || mxIsComplex(prhs[3]) ||
            mxIsSparse(prhs[0]) || mxIsSparse(prhs[1]) || mxIsSparse(prhs[2]) || mxIsSparse(prhs[3]) ) {
          mexErrMsgTxt("All inputs must be full real double vectors");
      }
      n = mxGetNumberOfElements(prhs[0]);
      if( mxGetNumberOfElements(prhs[1]) != n || mxGetNumberOfElements(prhs[2]) != n ) {
          mexErrMsgTxt("The first three inputs must be the same size");
      }
      m = mxGetNumberOfElements(prhs[3]);
      if( mxGetM(prhs[3]) != m ) {
          mexErrMsgTxt("The last input must be a column vector");
      }
      Ar = mxGetPr(prhs[0]);
      Ac = mxGetPr(prhs[1]);
      Av = mxGetPr(prhs[2]);
      B = mxGetPr(prhs[3]);
      plhs[0] = mxCreateDoubleMatrix(m,1,mxREAL);
      C = mxGetPr(plhs[0]);
      for( i=0; i<n; i++ ) {
          C[(mwSize)*Ar++ - 1] += *Av++ * B[(mwSize)*Ac++ - 1];
      }
  }

Some sample test code:

N = 5000;
r = ceil(N*rand(N*N,1));
c = ceil(N*rand(N*N,1));
v = rand(N*N,1);
B = rand(N,1);

disp(' ')
disp('mex routine')
tic
x = mfree_timesx(r,c,v,B);
toc

disp('sparse')
tic
s = sparse(r,c,v,numel(B),numel(B))*B;
toc

disp('accumarray')
tic
aux = B(c);
prod = v.*aux;
C = accumarray(r,prod);
toc

disp(' ')

disp('mex vs sparse')
max(abs(x-s))

disp('mex vs accumarray')
max(abs(x-C))

disp('sparse vs accumarray')
max(abs(s-C))

disp(' ')

and results:

mex routine
Elapsed time is 0.150372 seconds.
sparse
Elapsed time is 5.806035 seconds.
accumarray
Elapsed time is 0.666030 seconds.

mex vs sparse
ans =
  1.2960e-011
mex vs accumarray
ans =
     0
sparse vs accumarray
ans =
  1.2960e-011

[Paulo Ribeiro]

Dear James,

Fantastic! Thanks for this very detailed answer. I have just compiled your code (in Linux, with gcc 4.8) and tried your example. It's running :)

> mex mfree_timesx.c

I think many others will probably follow this thread. The core of your code is given by the final for loop? Outstanding!

I will send some results.Thanks!!!

[James Tursa]

Yes, that small loop does all of the real work. Everything else is just overhead. If you want I can add in code to check the individual r and c indexes to make sure they are legit.

Also, if these variables are really mxGPUArray variables, then you will need to replace all of the mxETC( ) functions with their mxGPU( ) counterparts. E.g., use (double *)mxGPUGetData(etc) instead of mxGetPr(etc).

[Paulo Ribeiro]

Dear James,

A GPU version for the 3 tuples and B will be more than welcome. Do you think it is manageable? Many thanks for your attention and work on this thread. Glad to find an expert :)

Best regards,

Paulo

[Paulo Ribeiro]

Here are some values after @James contributions:

Problem: finite element, 2D elasticity, 116.000 nodes. Setup: i7 7700HQ, NVIDIA GTX 1050, Linux (gcc compiler); Solver: matrix-free conjugate gradient, as described on this thread

Solver time for CPU with MEX compilation (mfree_timesx) = 42s;

Solver time for CPU with mfree_times (no MEX) = 317s;

Solver time for GPU with mfree_times (no MEX) = 82s.

Really impressive results with MEX on CPU

Answer 2

[Bruno Luong]

Can you give the detail about how you iterate ?

There seems no room to improve the basic block mfree_times

If you output is sparse and you add them, and internally it requires new memory allocation and merge sorted array.

Better building all collection triplets then call ACCUMARRAY once.

[Paulo Ribeiro]

Bruno,

This basic block mfree_times is the only difference between my code and the conjugate gradient function:

 function [ x ] = conjgrad_mfree(l,m,v, b, x,maxit) 
    r = b - mfree_times(l,m,v,x);
    p = r;
    rsold = r' * r;

    for i = 1:maxit
        Ap = mfree_times(l,m,v,p);
        alpha = rsold / (p' * Ap);
        x = x + alpha * p;
        r = r - alpha * Ap;
        rsnew = r' * r;
        if sqrt(rsnew) < 1e-5
              break;
        end
        p = r + (rsnew / rsold) * p;
        rsold = rsnew;
    end 
 end

Please notice that mfree_times is being called twice.

Thanks!

PS: I will send some updates.

[Bruno Luong]

It seems index complement of "l" has no true impact on CG algorithm. If it is a small vector with large value, you might want to squeeze as consecutive vector (using UNIQUE command) and work on the subspace of the new value.

That prevent CG work unnecessary large vector/matrix.

Then after converging, you can recover original index of the solution if that is matter to use the result.

[Paulo Ribeiro]

Bruno,

After several attempts, I am really out of options. Vectors l,m and v are very large, up to 1e9. I imagine two possible scenarios:

1. Fast convergence and reduced number of iterations; 2. Translate this code to Fortran or C/C++.

Thanks!

[Bruno Luong]

It seems that your problem is memory and not inneficient calculation.

If the matrix is sparse (a lot of zero) you should consider my suggestion if squeezing subspace.

By do try to reduce memory by using storage the indexes with INT32, matrix values with SINGLE, or using TALL ARRAY for more efficient HD cache.