csrmultcsr confusing calling API, not matching documentation

[philvdm]

I am struggling to use CSRMULTCSR ( sparse matrix-matrix multiplication)

I have my own reference version of Sparse matrix-matrix multiplication for CSR running on CPU. I would like to test it against MKL and I have hard time to match MKL when m,n,k are not equal.

The documentation and the MKL header file do not even have the same variable convention and their m,n,k mapping do not follow the convention.

MKL Documentation:
mkl_dcsrmultcsr(&trans, &request, &sort, &m, &n, &k, a, ja, ia, b, jb, ib, c, jc, ic, &nzmax, &info);
m INTEGER. Number of rows of the matrix A.
n INTEGER. Number of columns of the matrix A.
k INTEGER. Number of columns of the matrix B.

Header file :
void mkl_dcsrmultcsr(char *transa, MKL_INT *job, MKL_INT *sort, MKL_INT *n, MKL_INT *k, MKL_INT *m, double *a, MKL_INT *ja, MKL_INT *ia, double *b, MKL_INT *jb, MKL_INT *ib, double *c, MKL_INT *jc, MKL_INT *ic, MKL_INT *nnzmax, MKL_INT *ierr);

But for example CSRMM follows the normal convention m,n,k in that order where m= rows of A, n = cols of B, k = cols of A,rows of B :

void mkl_dcsrmm(char *transa, MKL_INT *m, MKL_INT *n, MKL_INT *k, double *alpha, char *matdescra, double *val, MKL_INT *indx, MKL_INT *pntrb, MKL_INT *pntre, double *b, MKL_INT *ldb, double *beta, double *c, MKL_INT *ldc);

How should I call mkl_dcsrmultcsr ?
What do m,n,k represent?

[Ying_H_Intel]

Hello,

thanks for the asking. Yes, it seemsthat the documentaion and declaration in header file should be corrected. we will fix it in the future release.

If use it, you may read it as Amxn* Bnxk.call it based on it's discription

mkl_dcsrmultcsr(&trans, &request, &sort, &m, &n, &k, a, ja, ia, b, jb, ib, c, jc, ic, &nzmax, &info);
m INTEGER. Number of rows of the matrix A.
n INTEGER. Number of columns of the matrix A.
k INTEGER. Number of columns of the matrix B

And thereisa sample given in the directory /mkl/examples/spblas/source/dcsr_multiplication.f

Regards,
Ying H.

[philvdm]

When m,n,k are NOT equal, it seems that there is some problem.
The example given in /mkl/examples/spblas/source is not helping because n=m=k

With m=3, n=4,k=5, I tried all the 6 combinations. 4 seg faults, 2 do not match my reference code with 3 errors.

I have even a smaller example using 2 small dense matrix in csr format

matrix A
RowPtr :1 3 5
ColId :1 2 1 2
Dense representation :
1.00, 2.00,
3.00, 4.00,

matrix B
RowPtr :1 4 7
ColId :1 2 3 1 2 3
Dense repesentation
1.00, 2.00, 3.00,
4.00, 5.00, 6.00,

My reference code gives:
matrix C
RowPtr :1 4 7
ColId :1 2 3 1 2 3
9.00, 12.00, 15.00,
19.00, 26.00, 33.00,

MKL 12
matric C
RowPtr :1 3 5
ColId :1 2 1 2
57.00, 12.00, 0.00,
19.00, 26.00, 0.00,

[Ying_H_Intel]

Hi Philvdm,

The discriptionof m, n, kis not aligned with other routines (we will fix this),but it is correctifreadit asCmxk=Amxn* Bnxk.
mkl_dcsrmultcsr(&trans, &request, &sort, &m, &n, &k, a, ja, ia, b, jb, ib, c, jc, ic, &nzmax, &info);
m INTEGER. Number of rows of the matrix A.
n INTEGER. Number of columns of the matrix A.
k INTEGER. Number of columns of the matrix B

So for your case, you may try m=2, n=2, k=3.
for example,

[bash]!*******************************************************************************
*   Copyright(C) 2004-2011 Intel Corporation. All Rights Reserved.
*   
*   The source code, information  and  material ("Material") contained herein is
*   owned  by Intel Corporation or its suppliers or licensors, and title to such
*   Material remains  with Intel Corporation  or its suppliers or licensors. The
*   Material  contains proprietary information  of  Intel or  its  suppliers and
*   licensors. The  Material is protected by worldwide copyright laws and treaty
*   provisions. No  part  of  the  Material  may  be  used,  copied, reproduced,
*   modified, published, uploaded, posted, transmitted, distributed or disclosed
*   in any way  without Intel's  prior  express written  permission. No  license
*   under  any patent, copyright  or  other intellectual property rights  in the
*   Material  is  granted  to  or  conferred  upon  you,  either  expressly,  by
*   implication, inducement,  estoppel or  otherwise.  Any  license  under  such
*   intellectual  property  rights must  be express  and  approved  by  Intel in
*   writing.
*   
*   *Third Party trademarks are the property of their respective owners.
*   
*   Unless otherwise  agreed  by Intel  in writing, you may not remove  or alter
*   this  notice or  any other notice embedded  in Materials by Intel or Intel's
*   suppliers or licensors in any way.
*
********************************************************************************
*   Content : MKL Sparse BLAS Fortran example
*
********************************************************************************
C----------------------------------------------------------------------
C Example program to show the use of the MKL Sparse BLAS routines
C for multiplication of two compressed sparse row format matrices
C---------------------------------------------------------------------
      PROGRAM multiplication_test
      IMPLICIT NONE
C..  Description of all variables
      CHARACTER trans
      INTEGER m,n,k, sort, job, ierr, nzmax
      INTEGER ia(9), ib(9)
      INTEGER ja(18), jb(18)
      REAL*8 a(18),b(18), cden(8, 8), answer(8, 8) 
      REAL*8 zero
      REAL*8 ddum
      REAL*8 tolerance, normA, normres
      INTEGER ic, jc
      REAL*8 c
      POINTER (IC_PTR, IC(1)), (JC_PTR, JC(1)), (C_PTR, C(1))
      INTEGER i, j ,js, idum, nnz, sizeint 
      REAL*8  DNRM2
c#ifdef _IA32
      INTEGER*4   MKL_MALLOC
      INTEGER*4 alloc_size
c#else
C      INTEGER*8 MKL_MALLOC
C      INTEGER*8 alloc_size
c#endif  
c ..  
c .. Fill all arrays containing matrix data.
c ..
      DATA m/2/, n /2/, k/3/,sort/0/, sizeint/8/ 
      DATA zero/0.D0/
      DATA ia /1,3,5/
      DATA ja
     1 /1, 2, 
     2  1, 2/

    
      DATA a
     1 /1.d0,  2.d0,  
     2   3.d0, 4.d0/
    
      
      DATA ib /1,4,7/
      DATA jb
     1 /1,  2,  3,
     2  1,  2,  3/
    
     
      DATA b
     1 /1.d0,      2.d0,     3.d0,
     2    4.d0 ,   5.d0,     6.d0/
    
     


C.. This test computes A*B 
C..
      trans='N'
           
c ..
      alloc_size = sizeint*(n+1)
      ic_ptr=mkl_malloc(alloc_size, 128)  
      if(ic_ptr.eq.0) then
         write(*,*) 'Cannot allocate pointer array of the length ',n+1
         stop 1
      endif
      nzmax=n*n
      alloc_size =sizeint*nzmax
      jc_ptr=mkl_malloc(alloc_size, 128)  
      if(jc_ptr.eq.0) then
         write(*,*) 'Cannot allocate column array of the length ',nnz
         stop 1
      endif
      alloc_size =8*nzmax
      c_ptr=mkl_malloc(alloc_size, 128)  
      if(c_ptr.eq.0) then
         write(*,*) 'Cannot allocate value array of the length ', nnz
         stop 1
      endif
      job=0
      call mkl_dcsrmultcsr(trans, job, sort, m, n, k, a,ja,ia,
     *        b, jb, ib, c, jc, ic, nzmax, ierr)
c  
      if(ierr.ne.0)then 
         print  *,' FIRST TEST FAILED  '
         stop 1
      endif
      

      write(*,*) 'ic RowPtr'
      write(*,*) (ic(i),i=1,m+1)
      write(*,*) 'jc colld'
      write(*,*) (jc(i),i=1,ic(m+1)-1)
      write(*,*) 'c array=A*B'
      write(*,*) (c(i),i=1,ic(m+1)-1)
      
      
      end[/bash]

it will give you expected result.
ic RowPtr
1 4 7
jc colld
1 2 3 1 2 3
c array=A*B
9.00000000000000 12.0000000000000 15.0000000000000
19.0000000000000 26.0000000000000 33.0000000000000
Press any key to continue . . .
Regards,
Ying