title: Writing a function to use Intel MKL Library for Sparse Matrix Multiplication
body: <p>I am trying to develop a function that multiplies matrix A and B which are in general format but sparse in nature. These matrices contain complex numbers. My issue is, when I do not use the function and I write everything in the main(), the multiplication works perfectly for any size of array. But when I am using my own function, result is corrupt and most of the time I get random error messages.</p>
  
  <p>Here is what this function does:</p>
  
  <ol>
  <li>Convert A and B to CSR format (mkl_zdnscsr).</li>
  <li>Create CSR handles for A and B using data in step 1 (mkl_sparse_z_create_csr).</li>
  <li>Multiply A and B using handles in step 2 (mkl_sparse_spmm) and store the output in "result".</li>
  </ol>
  
  <p>My guess is the way I return 'result'  from the function is somehow incorrect because I have checked the output of steps 1 &amp; 2 and they produce correct output.</p>
  
  <p>Any idea what the problem is? I am gonna include the summarized version of my code below for your reference.</p>
  
  <p>Thank you very much in advance.</p>
  
  <p>Afshin</p>
  
  <pre><code>/* ***************** Macro ********************* */
  #define ALIGN 128
  
  /* To avoid constantly repeating the part of code that checks different functions' status, using the below macros */
  #define CHECK_SPARSE(function) do { \
  if(function != SPARSE_STATUS_SUCCESS)             \
   {                                                 \
   status = 2;                                       \
   goto memory_free;                                 \
   }                                                 \
  } while(0)
  
  
  /* ****************** Main ******************** */
  int main()
  {
  
   &lt;&lt; matrices A and B are generated using some data &gt;&gt;
  
      MKL_INT stat = 0;
  
      // This part calls the function to multiply matrices as I discussed.
  
      // A x B --&gt; csrC
      sparse_matrix_t  csrC = NULL;
      stat = dnmm_sp_CSR_handle(CfPrime, Num_of_Buses, Num_of_Branches, CfPrime_nonzero, Yf, Num_of_Branches, Num_of_Buses, Yf_nonzero, &amp;csrC);
      printf("\nstat = %i", stat);
  
      // Now I convert the csrC to 4-array version of CSR.
      MKL_INT rows, cols;
      sparse_index_base_t indexing = 0;
      MKL_INT *columns_C = NULL, *pointerB_C = NULL, *pointerE_C = NULL;
      MKL_Complex16  *values_C = NULL;
      mkl_sparse_z_export_csr(csrC, &amp;indexing, &amp;rows, &amp;cols, &amp;pointerB_C, &amp;pointerE_C, &amp;columns_C, &amp;values_C);
  
      // Print the number of rows and columns of converted matrix (which are incorrect sizes)
      printf("\nrows = %i , cols = %i", rows, cols);
  }
  
  /* ****************** Function ******************** */
  // This function receives two dense matrice, convert them to sparse CSR format, multiply them, and returns the result in CSR handle
  int dnmm_sp_CSR_handle(MKL_Complex16 *A, MKL_INT A_rownum, MKL_INT A_colnum, MKL_INT A_nnz, MKL_Complex16 *B, MKL_INT B_rownum, MKL_INT B_colnum, MKL_INT B_nnz, sparse_matrix_t *result) {
  
      // A : Matrix A
      // A_rownum : Number of rows in matrix A
      // A_colnum : Number of columns in matrix A
      // A_nnz : Number of nonzero elements in matrix A
      // B : Matrix B
      // B_rownum : Number of rows in matrix B
      // B_colnum : Number of columns in matrix B
      // B_nnz : Number of nonzero elements in matrix B
      // result : return CSR handle for A x B
  
      MKL_INT job[8];
      job[0] = 0; // the rectangular matrix A is converted to the CSR format;
      job[1] = 0; // zero-based indexing for the rectangular matrix A is used;
      job[2] = 0; // zero-based indexing for the matrix in CSR format is used;
      job[3] = 2; // whole matrix
      //job[4] // maximum number of the non-zero elements allowed if job[0] = 0
      job[5] = 5; // If job[5]&gt;0, arrays acsr, ia, ja are generated for the output storage. If job[5]=0, only array ia is generated for the output storage.
      MKL_INT info = 0; // If info = 0, execution of mkl_zdnscsr was successful.
  
  
      MKL_INT status = 1; // return this value to check the execution status  
      //(1 : successfull, 2: error in sparse functions, 3: error in deallocating memory)
  
      MKL_Complex16 *A_val = (MKL_Complex16 *)mkl_malloc(A_nnz * sizeof(MKL_Complex16), ALIGN);
      MKL_INT *A_col = (MKL_INT *)mkl_malloc(A_nnz * sizeof(MKL_INT), ALIGN);
      MKL_INT *A_row = (MKL_INT *)mkl_malloc( (A_rownum + 1) * sizeof(MKL_INT), ALIGN); // +1 is because we are using 3-aaray variation
      job[4] = A_nnz;
      mkl_zdnscsr(job, &amp;A_rownum, &amp;A_colnum, A, &amp;A_colnum, A_val, A_col, A_row, &amp;info);
  
      MKL_Complex16 *B_val = (MKL_Complex16 *)mkl_malloc(B_nnz * sizeof(MKL_Complex16), ALIGN);
      MKL_INT *B_col = (MKL_INT *)mkl_malloc(B_nnz * sizeof(MKL_INT), ALIGN);
      MKL_INT *B_row = (MKL_INT *)mkl_malloc((B_rownum + 1) * sizeof(MKL_INT), ALIGN); // +1 is because we are using 3-aaray variation
      job[4] = B_nnz;
      mkl_zdnscsr(job, &amp;B_rownum, &amp;B_colnum, B, &amp;B_colnum, B_val, B_col, B_row, &amp;info);
  
      sparse_matrix_t   csrA = NULL, csrB = NULL;
      CHECK_SPARSE( mkl_sparse_z_create_csr(&amp;csrA, SPARSE_INDEX_BASE_ZERO, A_rownum, A_colnum, A_row, A_row + 1, A_col, A_val) );
      CHECK_SPARSE( mkl_sparse_z_create_csr(&amp;csrB, SPARSE_INDEX_BASE_ZERO, B_rownum, B_colnum, B_row, B_row + 1, B_col, B_val) );
      CHECK_SPARSE( mkl_sparse_spmm(SPARSE_OPERATION_NON_TRANSPOSE, csrA, csrB, &amp;result) );
  
  memory_free:
  
      //Release matrix handle and deallocate arrays for which we allocate memory ourselves.
      if (mkl_sparse_destroy(csrA) != SPARSE_STATUS_SUCCESS) status = 3;
      if (mkl_sparse_destroy(csrB) != SPARSE_STATUS_SUCCESS) status = 3;
  
      //Deallocate arrays for which we allocate memory ourselves.
      mkl_free(A_val); mkl_free(A_col); mkl_free(A_row);
      mkl_free(B_val); mkl_free(B_col); mkl_free(B_row);
  
      return status;
  }
  </code></pre>