title: Sparse matrix-matrix multiplication in CUDA using cuSPARSE
body: <p>I'm benchmarking the sparse matrix-matrix multiplication on Nvidia K40 using <code>cuSPARSE</code> library.
  I'm creating my own sparse matrix in <code>CSR</code> format and I'm using the <code>cusparseXcsrgemmNnz</code> routine of the <code>cuSPARSE</code> library.
  However, as I increase the data size, an error occurs when calling <code>cusparseXcsrgemmNnz</code>, i.e <code>CUSPARSE_STATUS_SUCCESS</code> is not returned.
  Also a <code>cudaMemcpy</code> fails and <code>CUDA_SUCCESS</code> is not returned.
  The code works fine for <code>8 x 8</code> and <code>16 x 16</code> matrices. However, from <code>32 x 32</code> on, this error is observed.</p>
  
  <p>Edit: I'm receiving <code>CUSPARSE_STATUS_EXECUTION_FAILED</code> from <code>cusparseXcsrgemmNnz</code> for the third matrix size. The execution happens correctly for the first two matrix sizes. </p>
  
  <pre><code>#include &lt;cusparse_v2.h&gt;
  #include &lt;stdio.h&gt;
  #include &lt;time.h&gt;
  #include &lt;sys/time.h&gt;
  
  // error check macros
  #define CUSPARSE_CHECK(x) {cusparseStatus_t _c=x; if (_c != CUSPARSE_STATUS_SUCCESS) {printf("cusparse fail: %d, line: %d\n", (int)_c, __LINE__); exit(-1);}}
  
  #define cudaCheckErrors(msg) \
  do { \
      cudaError_t __err = cudaGetLastError(); \
      if (__err != cudaSuccess) { \
          fprintf(stderr, "Fatal error: %s (%s at %s:%d)\n", \
              msg, cudaGetErrorString(__err), \
              __FILE__, __LINE__); \
          fprintf(stderr, "*** FAILED - ABORTING\n"); \
          exit(1); \
      } \
  } while (0)
  
  
  double timerval()
  {
      struct timeval st;
      gettimeofday(&amp;st, NULL);
      return (st.tv_sec+st.tv_usec*1e-6);
  }
  
  // perform sparse-matrix multiplication C=AxB
  int main(){ 
  double avg_time = 0, s_time, e_time;
  
  cusparseStatus_t stat;
  cusparseHandle_t hndl;
  cusparseMatDescr_t descrA, descrB, descrC;
  int *csrRowPtrA, *csrRowPtrB, *csrRowPtrC, *csrColIndA, *csrColIndB, *csrColIndC;
  int *h_csrRowPtrA, *h_csrRowPtrB, *h_csrRowPtrC, *h_csrColIndA, *h_csrColIndB, *h_csrColIndC,*pos;
  float *csrValA, *csrValB, *csrValC, *h_csrValA, *h_csrValB, *h_csrValC;
  int nnzA, nnzB, nnzC;
  int m=4,n,k,loop;
  int i,j;
  int iterations;
  for (iterations=0;iterations&lt;10;iterations++)
  {
      m *=2;
      n = m;
      k = m;
      //density of the sparse matrix to be created. Assume 5% density.
      double dense_const = 0.05;
      int temp5, temp6,temp3,temp4;
      int density=(m*n)*(dense_const);
      nnzA = density;
      nnzB = density;
      h_csrRowPtrA = (int *)malloc((m+1)*sizeof(int));
      h_csrRowPtrB = (int *)malloc((n+1)*sizeof(int));
      h_csrColIndA = (int *)malloc(density*sizeof(int));
      h_csrColIndB = (int *)malloc(density*sizeof(int));
      h_csrValA  = (float *)malloc(density*sizeof(float));
      h_csrValB  = (float *)malloc(density*sizeof(float));
      if ((h_csrRowPtrA == NULL) || (h_csrRowPtrB == NULL) || (h_csrColIndA == NULL) || (h_csrColIndB == NULL) || (h_csrValA == NULL) || (h_csrValB == NULL))
      {printf("malloc fail\n"); return -1;}
  
      //position array for random initialisation of positions in input matrix
      pos= (int *)calloc((m*n), sizeof(int));
      int temp,temp1;
  
      //  printf("the density is %d\n",density);
      //  printf("check 1:\n");
  
      //randomly initialise positions
      for(i=0;i&lt;density;i++)
      {
          temp1=rand()%(m*n);
          pos[i]=temp1;
  
      }
      //  printf("check 2:\n");
  
      //sort the 'pos' array
      for (i = 0 ; i &lt; density; i++) {
          int d = i;
          int t;
  
          while ( d &gt; 0 &amp;&amp; pos[d] &lt; pos[d-1]) {
              t          = pos[d];
              pos[d]   = pos[d-1];
              pos[d-1] = t;        
              d--;
          }
      }
      // initialise with non zero elements and extract column and row ptr vector
      j=1;
      //ja[0]=1;
  
      int p=0;
      int f=0;
  
      for(i = 0; i &lt; density; i++)
      {
          temp=pos[i];
           h_csrValA[f] = rand();
           h_csrValB[f] = rand();
           h_csrColIndA[f] = temp%m;
           h_csrColIndB[f] = temp%m;
          f++;
          p++;
          temp5= pos[i];
          temp6=pos[i+1];
          temp3=temp5-(temp5%m);
          temp4=temp6-(temp6%m);
  
          if(!(temp3== temp4))
          {   
              if((temp3+m==temp6))
              {}
              else    
              {   
                  h_csrRowPtrA[j]=p;
                  h_csrRowPtrB[j]=p;
                  j++;
              }
  
          }       
      }
  
      // transfer data to device
  
      cudaMalloc(&amp;csrRowPtrA, (m+1)*sizeof(int));
      cudaMalloc(&amp;csrRowPtrB, (n+1)*sizeof(int));
      cudaMalloc(&amp;csrColIndA, density*sizeof(int));
      cudaMalloc(&amp;csrColIndB, density*sizeof(int));
      cudaMalloc(&amp;csrValA, density*sizeof(float));
      cudaMalloc(&amp;csrValB, density*sizeof(float));
      cudaCheckErrors("cudaMalloc fail");
      cudaMemcpy(csrRowPtrA, h_csrRowPtrA, (m+1)*sizeof(int), cudaMemcpyHostToDevice);
      cudaMemcpy(csrRowPtrB, h_csrRowPtrB, (n+1)*sizeof(int), cudaMemcpyHostToDevice);
      cudaMemcpy(csrColIndA, h_csrColIndA, density*sizeof(int), cudaMemcpyHostToDevice);
      cudaMemcpy(csrColIndB, h_csrColIndB, density*sizeof(int), cudaMemcpyHostToDevice);
      cudaMemcpy(csrValA, h_csrValA, density*sizeof(float), cudaMemcpyHostToDevice);
      cudaMemcpy(csrValB, h_csrValB, density*sizeof(float), cudaMemcpyHostToDevice);
      cudaCheckErrors("cudaMemcpy fail");
  
      // set cusparse matrix types
      CUSPARSE_CHECK(cusparseCreate(&amp;hndl));
      stat = cusparseCreateMatDescr(&amp;descrA);
      CUSPARSE_CHECK(stat);
      stat = cusparseCreateMatDescr(&amp;descrB);
      CUSPARSE_CHECK(stat);
      stat = cusparseCreateMatDescr(&amp;descrC);
      CUSPARSE_CHECK(stat);
      stat = cusparseSetMatType(descrA, CUSPARSE_MATRIX_TYPE_GENERAL);
      CUSPARSE_CHECK(stat);
      stat = cusparseSetMatType(descrB, CUSPARSE_MATRIX_TYPE_GENERAL);
      CUSPARSE_CHECK(stat);
      stat = cusparseSetMatType(descrC, CUSPARSE_MATRIX_TYPE_GENERAL);
      CUSPARSE_CHECK(stat);
      stat = cusparseSetMatIndexBase(descrA, CUSPARSE_INDEX_BASE_ZERO);
      CUSPARSE_CHECK(stat);
      stat = cusparseSetMatIndexBase(descrB, CUSPARSE_INDEX_BASE_ZERO);
      CUSPARSE_CHECK(stat);
      stat = cusparseSetMatIndexBase(descrC, CUSPARSE_INDEX_BASE_ZERO);
      CUSPARSE_CHECK(stat);
      cusparseOperation_t transA = CUSPARSE_OPERATION_NON_TRANSPOSE;
      cusparseOperation_t transB = CUSPARSE_OPERATION_NON_TRANSPOSE;
  
      // figure out size of C
      int baseC;
      // nnzTotalDevHostPtr points to host memory
      int *nnzTotalDevHostPtr = &amp;nnzC;
      stat = cusparseSetPointerMode(hndl, CUSPARSE_POINTER_MODE_HOST);
      CUSPARSE_CHECK(stat);
      cudaMalloc((void**)&amp;csrRowPtrC, sizeof(int)*(m+1));
      cudaCheckErrors("cudaMalloc fail");
  
      s_time=timerval();
  
      stat = cusparseXcsrgemmNnz(hndl, transA, transB, m, n, k,
      descrA, nnzA, csrRowPtrA, csrColIndA,
      descrB, nnzB, csrRowPtrB, csrColIndB,
      descrC, csrRowPtrC, nnzTotalDevHostPtr );
      CUSPARSE_CHECK(stat);
      if (NULL != nnzTotalDevHostPtr){
      nnzC = *nnzTotalDevHostPtr;}
      else{
      cudaMemcpy(&amp;nnzC, csrRowPtrC+m, sizeof(int), cudaMemcpyDeviceToHost);
      cudaMemcpy(&amp;baseC, csrRowPtrC, sizeof(int), cudaMemcpyDeviceToHost);
      cudaCheckErrors("cudaMemcpy fail");
      nnzC -= baseC;}
      cudaMalloc((void**)&amp;csrColIndC, sizeof(int)*nnzC);
      cudaMalloc((void**)&amp;csrValC, sizeof(float)*nnzC);
      cudaCheckErrors("cudaMalloc fail");
      // perform multiplication C = A*B
  
      for(loop=0;loop&lt;1000;loop++)
      {
          stat = cusparseScsrgemm(hndl, transA, transB, m, n, k,
          descrA, nnzA,
          csrValA, csrRowPtrA, csrColIndA,
          descrB, nnzB,
          csrValB, csrRowPtrB, csrColIndB,
          descrC,
          csrValC, csrRowPtrC, csrColIndC);
          CUSPARSE_CHECK(stat);
      }
  
      e_time=timerval();
  
      avg_time=avg_time/1000;
      // copy result (C) back to host
      h_csrRowPtrC = (int *)malloc((m+1)*sizeof(int));
      h_csrColIndC = (int *)malloc(nnzC *sizeof(int));
      h_csrValC  = (float *)malloc(nnzC *sizeof(float));
      if ((h_csrRowPtrC == NULL) || (h_csrColIndC == NULL) || (h_csrValC == NULL))
      {printf("malloc fail\n"); return -1;}
      cudaMemcpy(h_csrRowPtrC, csrRowPtrC, (m+1)*sizeof(int), cudaMemcpyDeviceToHost);
      cudaMemcpy(h_csrColIndC, csrColIndC,  nnzC*sizeof(int), cudaMemcpyDeviceToHost);
      cudaMemcpy(h_csrValC, csrValC, nnzC*sizeof(float), cudaMemcpyDeviceToHost);
      cudaCheckErrors("cudaMemcpy fail");
  
      printf ("\n Input size: %d x %d ,Time: %lf and density is %d \n", m,n, avg_time, density); 
  
      cudaFree(csrRowPtrC);
      cudaFree(csrColIndC);
      cudaFree(csrValC);
  
      cudaFree(csrRowPtrA);
      cudaFree(csrColIndA);
      cudaFree(csrValA);
  
      cudaFree(csrRowPtrB);
      cudaFree(csrColIndB);
      cudaFree(csrValB);
  
      free(h_csrRowPtrC);
      free(h_csrColIndC);
      free(h_csrValC);
  
      free(h_csrRowPtrA);
      free(h_csrColIndA);
      free(h_csrValA);
  
      free(h_csrRowPtrB);
      free(h_csrColIndB);
      free(h_csrValB);
  }
  return 0;
  </code></pre>
  
  <p>}</p>