Help with Sparse matrix vector multiplication using MKL DIA routine

[Hazem_A_]

I am using the MKL library to perform the sparse matrix vector multiplication using diagonal format, When I use the MKL mkl_sdiagemv function I get a "MKL ERROR: Parameter 4 was incorrect on entry to MKL_SDIAGEMV. " error.

void mm_read(char* filename, int *m, int *n, int *nnz, int **rowptrs, int **colinds, float **vals, float **adia, int **distance, int idiag, int ndiag) {
// open file 
FILE* mmfile = fopen(filename, "r");
assert(mmfile != NULL && "Read matrix file.");

// read MatrixMarket header
int status;
MM_typecode matcode;
status = mm_read_banner(mmfile, &matcode);
assert(status == 0 && "Parsed banner.");



status = mm_read_mtx_crd_size(mmfile, m, n, nnz);
assert(status == 0 && "Parsed matrix m, n, and nnz.");
printf("- matrix is %d-by-%d with %d nnz.\n", *m, *n, *nnz);


int *coo_rows = (int*) malloc(*nnz * sizeof(int));
int *coo_cols = (int*) malloc(*nnz * sizeof(int));
float *coo_vals = (float*) malloc(*nnz * sizeof(float));

// read COO values
int i = 0;
for ( i = 0; i < *nnz; i++)
    status = fscanf(mmfile, "%d %d %g\n",
        &coo_rows, &coo_cols, &coo_vals);


*rowptrs = (int*) malloc((*m+1)*sizeof(int));
*colinds = (int*) malloc(*nnz*sizeof(int));
*vals = (float*) malloc(*nnz*sizeof(int));

// convert to CSR matrix
int info;
int job[] = {
    2, // job(1)=2 (coo->csr with sorting)
    1, // job(2)=1 (one-based indexing for csr matrix)
    1, // job(3)=1 (one-based indexing for coo matrix)
    0, // empty
    *nnz, // job(5)=nnz (sets nnz for csr matrix)
    0  // job(6)=0 (all output arrays filled)
};

       mkl_scsrcoo(  job,   m, *vals,  *colinds, *rowptrs,  nnz,   coo_vals,   coo_rows,    coo_cols,  &info);
      assert(info == 0 && "Converted COO->CSR");



 // DIA matrix dimensions and values    

ndiag = 4;
idiag = 3;

*adia = (float*) malloc(*nnz * idiag * sizeof(int));
*distance = (int* ) malloc(idiag * sizeof(int));


       int job1[] = {
    0, // job(0)=2 
    1, // job(1)=1 
    1, // job(2)=1 
    2, // empty3
    *nnz, // empty4
    10, // job(5)=nnz 
};

mkl_scsrdia (job1, m, *vals, *colinds, *rowptrs,  *adia,  &ndiag, *distance, &idiag, *vals, *colinds, *rowptrs, &info);
  assert(info == 0 && "Converted CSR->DIA");


// free COO matrix
free(coo_rows);
free(coo_cols);
free(coo_vals);   }



 float * randvec(int n) {
float *v = (float*) malloc(n * sizeof(float));
int i = 0;
for (i = 0; i < n; i++)
    v = rand() / (float) RAND_MAX;
return v;
}


   int main(int argc, char* argv[]) {

// require filename for matrix to test
if (argc != 2) {
    fprintf(stderr, "Usage: %s MATRIX.mm\n", argv[0]);
    exit(1);
}

int m, n, nnz;
int *rowptrs, *colinds;
float *vals;

float *adia;
int  *distance;
int idiag,  ndiag;
// read matrix from file
mm_read(argv[1], &m, &n, &nnz, &rowptrs, &colinds, &vals, &adia, &distance, &idiag, &ndiag);

// allocate vectors for computation
float *v = randvec(n);
float *cpu_answer = (float*) malloc(m*sizeof(float));


struct timeval start, end;
printf (" Running Intel(R) MKL from 1 to %i threads \n\n", mkl_get_max_threads());


mkl_sdiagemv ((char*)"N", &m, adia, &idiag, distance, &ndiag, v, cpu_answer);


// release memory
free(rowptrs);
free(colinds);
free(vals);
free(cpu_answer);
free(v);
free(adia);
free(distance);

return 0; }

[mecej4]

I am not familiar with the Matrix Market conventions, and you have not shown the #include statements in your program, so it is not possible for me to verify if the arguments being passed to mkl_sdiagemv() are correct.

[Hazem_A_]

Sorry about that here are the included header files:

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <float.h>
#include <sys/time.h>
#include "mmio.h"
#include "mmio.c"
#include <mkl_spblas.h>

and you can find the "mmio.c" and "mmio.h" files from this site:

http://math.nist.gov/MatrixMarket/mmio-c.html

[mecej4]

Did you have a specific test matrix to input to the program -- one that gave you the MKL runtime error?

[Hazem_A_]

I used the following matrix https://docs.google.com/file/d/0B828bj6uvoMVcEdjaTFfNUdyYmM/edit

which has three diagonals and consists of 2,999,997 non zero elements in a 1,000,000*1,000,000 sparse matrix

[Sergey_P_Intel2]

Hi,

I can't build your reproducer:

test.cpp(102): error: argument of type "int *" is incompatible with parameter of type "int"
        mm_read(argv[1], &m, &n, &nnz, &rowptrs, &colinds, &vals, &adia, &distance, &idiag, &ndiag);
                                                                                                                           ^

test.cpp(102): error: argument of type "int *" is incompatible with parameter of type "int"
        mm_read(argv[1], &m, &n, &nnz, &rowptrs, &colinds, &vals, &adia, &distance, &idiag, &ndiag);
                                                                                                                                      ^
Parameters are declared in the function as

int idiag, int ndiag)

But in the call you try to use pointers:

 &idiag, &ndiag); 

Regards,

Sergey

[mecej4 · Accepted solution]

Rather than contend with a 10,000 X 10,000 matrix as the first test case, perhaps we should start with a simple, short example. Since I could not find an example in the files distributed in examples_core.zip in the MKL/examples directory, I made one up.

The diagonal-storage matrix from the MKL Ref. Man. (https://software.intel.com/en-us/node/522243#MKL_APPA_SMSF_6 ) is multiplied with a vector with the following code.

// MECEJ4 10/14/2014

#include <stdio.h>
#define N 5
#define NDIAG 5

int main(){

   int dist = { -3,-1, 0, 1, 2 };
   float u = 1e37;
   float val[N*NDIAG] =    { u, u, u,-4, 8,
                             u,-2, 0, 2, 0,
                             1, 5, 4, 7,-5,
                            -1, 0, 6, 0, u,
                            -3, 0, 4, u, u };
   float x = { 1, 2, 3, 4, 5 }, y;
   char transa = 'N';
   int n=N, ndiag=NDIAG;

   mkl_sdiagemv(&transa, &n, val, &n, dist, &ndiag, x, y);

   printf(" %12.4e %12.4e %12.4e %12.4e %12.4e\n",y[0],y[1],y[2],y[3],y[4]);
}

[Hazem_A_]

Thank you very much mecej4 as my error says I had a wrong value passed to the fourth parameter in the mkl_sdiagemv. Your example hinted me to what should I pass instead.