mkl_sparse_d_mm() wrong results

[di_luca__daniele1]

Hello everybody,

I ran into wrong results when I was trying to use mkl_sparse_d_mm() to multiply a symmetric sparse matrix A [12 x 12]  (CSR packed, lower triangular) by a dense rectangular matrix B [12 x 3] (full), that is:

C = A * B

Both matrices A and B have column wise layout, while the A matrix's CSR indeces arrays are all 0-based.

Here the example code:

int main() 
{ 
   // sparse matrix A[12x12] (lower triangular, column wise layout, 0-based indexing) 
   struct matrix_descr A_descr; 
   A_descr.type = SPARSE_MATRIX_TYPE_TRIANGULAR; 
   A_descr.mode = SPARSE_FILL_MODE_LOWER; 
   A_descr.diag = SPARSE_DIAG_NON_UNIT; 

   // A rows_start 
   MKL_INT A_ptrb[12] = { 0, 3, 8, 12, 14, 17, 20, 22, 25, 27, 28, 29 }; 
   // A rows_end 
   MKL_INT A_ptre[12] = { 3, 8, 12, 14, 17, 20, 22, 25, 27, 28, 29, 30 }; 

   // A column_indeces 
   MKL_INT A_col_index[30] = { 0, 4, 6, 1, 3, 5, 7, 11, 2, 4, 8, 10, 3, 9, 4, 8, 10, 5, 7, 11, 6, 10, 7, 9, 11, 8, 10, 9, 10, 11 }; 

   // A values
   double A_values[30] = { 3.60394e+06, -2.27463e+07, -3.14901e+06, 462747, 2.27463e+07, -1.56399e+06, -7819.96, -1.56399e+06, 1.26039e+07, 1.56399e+06, -7819.96, 1.56399e+06, 1.63638e+09, -7.91937e+07, 1.97495e+09, -1.56399e+06, 2.07831e+08, 7.3454e+08, 1.56399e+06, 2.07831e+08, 3.60394e+06, -2.27463e+07, 462747, 2.27463e+07, 1.56399e+06, 1.26039e+07, -1.56399e+06, 1.63638e+09, 1.97495e+09, 7.3454e+08 };

   // handle to A matrix
   sparse_matrix_t A_handle;
   sparse_status_t outcome = mkl_sparse_d_create_csr(&A_handle, SPARSE_INDEX_BASE_ZERO, 12, 12, A_ptrb, A_ptre, A_col_index, A_values); 

   // dense matrix B [12x3] values (full, column wise layout) 
   double B[36] = { 0.447212, 0, -0.00115627, 0, 0.00466197, 0, 0.447212, 0, 0.00115627, 0, 0.00466197, 0, 3.77898e-21, 0, -0.447214, 0, 5.05844e-22, 0, 6.67055e-22, 0, -0.447214, 0, -4.18026e-22, 0, 0, 0.447214, 0, -0.00653261, 0, -3.29415e-20, 0, 0.447214, 0, -0.00653261, 0, -3.44084e-21 }; 

   // C = A * B double C[36]; 
   mkl_sparse_d_mm(SPARSE_OPERATION_NON_TRANSPOSE, 1.0, A_handle, A_descr, SPARSE_LAYOUT_COLUMN_MAJOR, B, 3, 12, 0.0, C, 12); 

   // print C 
   for (size_t r = 0; r != 12; r++) 
   { 
      for (size_t c = 0; c != 3; c++) 
         std::cout << C[r + 12*c] << " "; 
      std::cout << std::endl; 
   } 
   return 0; 
}

Could it be a problem with the triangular form assigned to the CSR matrix?  Does the routine need the full CSR matrix maybe? The results seem to be correct if the sparse matrix has no off-diagonal terms (i.e. it it's diagonal).