Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.

mkl_sparse_d_mm problem

di_luca__daniele1
1,960 Views

Hello everybody,

I ran into problems trying to use both mkl_sparse_d_mm() and mkl_dcsrmm() 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 3 A's CSR index arrays are all 0-based.

Both results I obtained by the 2 routines are different from what I calculated by myself (or by using other software), and they're both different from each other.

In particular, the parameters I gave to the mkl_sparse_d_mm() routine, which is not deprecated as mkl_dcsrmm() is, were (you can find all arrays I used for calculations in "test.7z"):

A matrix:

pointerE= pointer to integers array [3, 8, 12, 14, 17, 20, 22, 25, 27, 28, 29, 30];

pointerB= pointer to integers array [0, 3, 8, 12, 14, 17, 20, 22, 25, 27, 28, 29];

columns= pointer to integers array [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];

values= pointer to doubles array [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];

B matrix:

B= pointer to doubles array, column wise, but in matrix form it looks like:

[0.447212, 3.77898e-21, 0;
0, 0, 0.447214;
-0.00115627, -0.447214 0;
0, 0, -0.00653261;
0.00466197, 5.05844e-22, 0;
0, 0, -3.29415e-20;
0.447212, 6.67055e-22, 0;
0, 0, 0.447214;
0.00115627, -0.447214 0;
0, 0, -0.00653261;
0.00466197, -4.18026e-22, 0;
0, 0, -3.44084e-21];

I got the A matrix handle from:

mkl_sparse_d_create_csr(&A,SPARSE_INDEX_BASE_ZERO, 12, 12, pointerB, pointerE, columns, values);

then:

struct matrix_descr Atype;

Atype.type=SPARSE_MATRIX_TYPE_TRIANGULAR;

Atype.mode=SPARSE_FILL_MODE_LOWER;

Atype.diag=SPARSE_DIAG_NON_UNIT;

mkl_sparse_d_mm(SPARSE_OPERATION_NON_TRASPOSE,1.0,A,Atype,SPARSE_LAYOUT_COLUMN_MAJOR, B, 3, 12, 0.0, C, 12);

results in matrix form looks like:

C=

[1.61172e+06,1.36192e-14, 0;
0, 0, 206947;
-14573.5, -5.63662e+06, 0;
0, 0, -1.06898e+07;
9.20715e+06, 0, 0;
0, 0, 0;
1.61172e+06, 2.40403-15, 0;
0, 0, 206947;
14573.5, -5.63662e+06, 0;
0, 0, -1.06898e+07;
9.20715e+06, -8.25579-13, 0;
0, 0, -2.52744e-12];

while by using other solvers I obtained something like (in matrix form):

C=

[97406.0182688134, 3.35091630430399e-14, 0;
0, 0, 54856.5214615356;
2.18278728425503e-11, -5633120.33746129, 0;
0, 0, 1.80443748831749e-09;
1.04773789644241e-09, -4.06713794809850e-12, 0;
0, 0, 9.00096050061251e-11;
97406.0182688133, 1.13296949003555e-12, 0;
0, 0, 54856.5214615356;
1.63709046319127e-11, -5633120.33746129, 0;
0, 0, 1.86264514923096e-09;
1.86264514923096e-09, -3.86485664837091e-11, 0;
0, 0, 2.30894555998986e-11].

as you can see, there are many differences, as an example, the first element:

mkl_sparse_d_mm: C(1,1)=1.61172e+06;

other routine: C(1,1)= 97406.0182688134

 

...what did I do wrong?

Thanks in advance,

Daniele

 

NOTE:

In the attached file "test.7z" I inserted the 4 CRS requested arrays for the A sparse matrix (low triangular, column wise):

A_ptrB.bin : 12 integers array (pointerB);

A_ptrE.bin: 12 integers array (pointerE);

A_idx.bin: 30 integers array (columns);

A_val.bin: 30 doubles array (values);

as well as the full dense matrix B  [12x3]:

B_dense_matrix.bin: 36 doubles array, ordered by column wise layout

0 Kudos
7 Replies
di_luca__daniele1
1,925 Views

here some code to prove what I wrote above.

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; 
}

 

0 Kudos
di_luca__daniele1
1,909 Views

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).

0 Kudos
mecej4
Honored Contributor III
1,887 Views

Daniela,

The manner in which you have presented your questions makes it difficult to answer without spending a lot of time. Here are some hints for how to pose the questions such that more forum members can respond, perhaps some to one question only, others to the second, etc.:

  • If the question is about how to use some MKL routines, how to organize the matrices and call the routines with the correct arguments, it is best do so verbally. Next best is with an example containing as few numbers as possible without making it trivial. For example, take matrices that are 2 X 2 or 3 X 3, get the test problem to work, then modify the code to 12 X 12 or whatever.
  • If the question is related to code that you have already written and tested, but is not working yet, make the complete code available. Fragments of code are next to useless in this regard. Even if all the fragments that you show are completely correct, the fragments that you did not show may contain errors. Secondly, we cannot run the fragments through the compiler and linker.
  • You presented some matrix data using *.bin files, but you did not tell us what the files contain and in what order and form. Your code fragments do not include code to read or write these files. There may be dozens of ways of doing this, and we are unlikely to try all or even one of them.

Among the example files provided with MKL is file sparse_csr.c  in the examples/spblasc/source directory, in which are included examples of multiplying a sparse matrix and a full matrix. Perhaps you will find it easier to modify that file to run your problem.

 

di_luca__daniele1
1,880 Views

Dear @mecej4, thanks a lot  for your kind reply and the helpful tips,

luckly I finally found my mistake (and I'm pretty embaressed about it): since I defined my symmetric sparse_matrix CSR row-indeces by column major layout for the lower triangle only, I had to to set the sparse_matrix description type on SPARSE_FILL_MODE_UPPER instead of LOWER to get the right CSR row-indeces. 

It's a little confusing, and probably I have to change other things too, but row-indexing was surely wrong.
Thank you a lot anyway,
Daniele

0 Kudos
di_luca__daniele1
1,862 Views
The main problem wasn't about indexing, as I thought, but about inserting the triangular part of the sparse matrix only(just like in Paradiso). May you confirm the routine needs the full Matrix in CSR form and not the upper/lower triangle only, even though the matrix is symmetric, or am I wrong again?
0 Kudos
mecej4
Honored Contributor III
1,847 Views

A lower triangular matrix is not a symmetric matrix unless it is also a diagonal matrix.

It is convenient to be able to enter only the elements of the lower triangle of a symmetric matrix and have the MKL routine fill in the strictly upper triangular part. On the other hand, you should not specify that the matrix itself is triangular (because that would require the strictly upper triangle to contain zeroes, instead of being a mirror image of the strictly lower triangle).

Is this the error in what you have done? Check by zeroing out the upper triangle and repeating your non-MKL calculation. Do the results match now?

di_luca__daniele1
1,788 Views

Thanks @mecej4 , that was the problem, I though I had to declare my symmetric matrix as "triangular," since I stored one triangle only, but I was wrong, I just needed to declare it as "symmetric". 

0 Kudos
Reply