From c4c72685a3d05c6dfa225e57e2b7ded4b5b77d0d Mon Sep 17 00:00:00 2001
From: "tempy@maggie" <kadir.cs@gmail.com>
Date: Wed, 20 Jan 2016 11:18:45 +0200
Subject: [PATCH] SpMV using MKL for shared memory has been added

---
 spmv/mkl_shmem/CMakeLists.txt |  60 ++++
 spmv/mkl_shmem/mklspmv.c      | 229 +++++++++++++++
 spmv/mkl_shmem/mmio.c         | 511 ++++++++++++++++++++++++++++++++++
 spmv/mkl_shmem/mmio.h         | 133 +++++++++
 spmv/mkl_shmem/test.mtx       |  13 +
 5 files changed, 946 insertions(+)
 create mode 100644 spmv/mkl_shmem/CMakeLists.txt
 create mode 100644 spmv/mkl_shmem/mklspmv.c
 create mode 100644 spmv/mkl_shmem/mmio.c
 create mode 100644 spmv/mkl_shmem/mmio.h
 create mode 100644 spmv/mkl_shmem/test.mtx

diff --git a/spmv/mkl_shmem/CMakeLists.txt b/spmv/mkl_shmem/CMakeLists.txt
new file mode 100644
index 0000000..eb1efbd
--- /dev/null
+++ b/spmv/mkl_shmem/CMakeLists.txt
@@ -0,0 +1,60 @@
+
+# ==================================================================================================
+# This file is part of the CodeVault project. The project is licensed under Apache Version 2.0.
+# CodeVault is part of the EU-project PRACE-4IP (WP7.3.C).
+#
+# Author(s):
+#   Valeriu Codreanu <valeriu.codreanu@surfsara.nl>
+#
+# ==================================================================================================
+cmake_minimum_required(VERSION 2.8.10 FATAL_ERROR)
+list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/../../../../cmake/Modules")
+set(CMAKE_VERBOSE_MAKEFILE ON)
+
+find_package(MKL)
+
+
+include(${CMAKE_CURRENT_SOURCE_DIR}/../../../../cmake/common.cmake)
+
+# ==================================================================================================
+
+if ("${DWARF_PREFIX}" STREQUAL "")
+  set(DWARF_PREFIX 2_sparse)
+endif()
+set(NAME ${DWARF_PREFIX}_spmv_mkl_shmem)
+
+# ==================================================================================================
+# C++ compiler settings
+
+find_package(Common)
+
+select_compiler_flags(cxx_flags
+  GNU "-march=native"   # I suggest remove "-O3" as this is controlled by the CMAKE_BUILD_TYPE
+  CLANG "-march=native" # same here
+  Intel "-axavx2,avx")
+set(CXX_FLAGS ${cxx_flags})
+if("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU")
+  set(CXX_FLAGS "${CXX_FLAGS} -Wall -Wno-comment")
+  if(APPLE)
+    set(CXX_FLAGS "${CXX_FLAGS} -Wa,-q")
+  endif()
+endif()
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${CXX_FLAGS}")
+
+# ==================================================================================================
+
+# LUD with the MKL library
+if (MKL_FOUND)
+	include_directories(${MKL_INCLUDE_DIR})
+	link_directories(${MKL_LIBRARY_DIR})
+	add_executable(${NAME} mklspmv.c mmio.c)
+	target_link_libraries(${NAME} mkl_intel_lp64 mkl_sequential mkl_core pthread m)
+	install(TARGETS ${NAME} DESTINATION bin)
+else ()
+	message("## Skipping '${NAME}': no MKL support found")
+  install(CODE "MESSAGE(\"${NAME} can only be built with MKL.\")")
+endif()
+
+unset(NAME)
+
+# ==================================================================================================
diff --git a/spmv/mkl_shmem/mklspmv.c b/spmv/mkl_shmem/mklspmv.c
new file mode 100644
index 0000000..0c26b86
--- /dev/null
+++ b/spmv/mkl_shmem/mklspmv.c
@@ -0,0 +1,229 @@
+/*includes {{{*/
+#include <stdio.h>
+#include <assert.h>
+#include <math.h>
+#include <mkl.h>
+#include "mmio.h"
+/*}}}*/
+
+/* global varibles and defines {{{*/
+int nrepeat = 100;
+int option_print_matrices = 0;
+#define OPTION_NOPRINT_MATRICES 0
+#define OPTION_PRINT_MATRICES 1
+char transa = 'n';
+/*}}}*/
+
+/* method declarations {{{*/
+int  main(int argc, char* argv[]);
+void printmm_one(int m, double* Aval, int* AJ, int* AI);
+void printfilemm_one(char* file, int m, int n, double* Aval, int* AJ, int* AI);
+void printmm_zero(int m, double* Aval, int* AJ, int* AI);
+/*}}}*/
+
+/** Multiply two sparse matrices which are stored in CSR format. MKL is used */
+void mkl_cpu_spmv(const MKL_INT Am, const MKL_INT An,  double* Aval, MKL_INT* AJ, MKL_INT* AI, double* xvec, double* yvec, double* time) {  /*{{{*/
+
+MKL_INT ierr;
+MKL_INT request = 1;    // symbolic multiplication
+double alpha = 1.0;
+double beta = 0.0;
+char* matdescra = "G NF";
+
+double time_st = dsecnd();
+int i;
+for(i = 0; i < nrepeat; i++) {
+    /** y := alpha*A*x + beta*y */
+    mkl_dcsrmv(&transa,  &Am, &An, &alpha, matdescra, Aval, AJ, AI, (AI+1), xvec, &beta, yvec);
+}
+double time_end = dsecnd();
+*time = (time_end - time_st)/nrepeat;
+} /*}}}*/ 
+
+/** Read Matrix Market file into COO matrix */
+void read_mm(char* strpath, int* pM, int* pN, int* prealnnz, int** pI, int** pJ, double** pval){ /*{{{*/
+/* 
+*   taken from Matrix Market I/O library for ANSI C
+*
+*   See http://math.nist.gov/MatrixMarket for details.
+*
+*
+*/
+
+int i, M, N, nz, *I, *J;
+double* val;
+int ret_code;
+MM_typecode matcode;
+FILE* f;
+if ((f = fopen(strpath, "r")) == NULL) {fprintf(stderr, "Input matrix file %s cannot be opened to read.", strpath);exit(1);}
+/* READ MATRIX */
+if (mm_read_banner(f, &matcode) != 0) {
+    printf("Could not process Matrix Market banner.\n");    exit(1);
+}
+/*  This is how one can screen matrix types if their application */
+/*  only supports a subset of the Matrix Market data types.      */
+if (mm_is_complex(matcode) && mm_is_matrix(matcode) && mm_is_sparse(matcode) ) {
+    printf("Sorry, this application does not support ");
+    printf("Market Market type: [%s]\n", mm_typecode_to_str(matcode));exit(1);
+}
+/* find out size of sparse matrix .... */
+if ((ret_code = mm_read_mtx_crd_size(f, &M, &N, &nz)) !=0) exit(1);
+/* reseve memory for matrices */
+I = (int *) malloc((mm_is_symmetric(matcode) ? 2*nz : nz) * sizeof(int));
+J = (int *) malloc((mm_is_symmetric(matcode) ? 2*nz : nz) * sizeof(int));
+val = (double *) malloc((mm_is_symmetric(matcode) ? 2*nz : nz) * sizeof(double));
+*pI = I;
+*pJ = J;
+*pval = val;
+/* NOTE: when reading in doubles, ANSI C requires the use of the "l"  */
+/*   specifier as in "%lg", "%lf", "%le", otherwise errors will occur */
+/*  (ANSI C X3.159-1989, Sec. 4.9.6.2, p. 136 lines 13-15)            */
+int realnnz = 0;
+for (i=0; i<nz; i++) {
+    if(mm_is_pattern(matcode)) {
+        fscanf(f, "%d %d\n", &I[realnnz], &J[realnnz]);
+        val[realnnz] = 1.0;
+    }
+    else
+        fscanf(f, "%d %d %lg\n", &I[realnnz], &J[realnnz], &val[realnnz]);
+    I[realnnz]--;  /* adjust from 1-based to 0-based */
+        J[realnnz]--;
+    if(mm_is_symmetric(matcode) && I[realnnz] != J[realnnz]) {
+        I[realnnz+1] = J[realnnz];
+        J[realnnz+1] = I[realnnz];
+        val[realnnz+1] = val[realnnz];
+        realnnz++;
+    }
+    realnnz++;
+}
+if (f !=stdin) fclose(f);
+*pM = M;
+*pN = N;
+*prealnnz = realnnz;
+} /* ENDOF read_mm }}}*/
+
+/** Converts COO matrix to CSR matrix */
+void coo_to_csr(int m, int nnz, int* I, int* J, double* val, MKL_INT* AI, MKL_INT* AJ, double* Aval) { /*{{{*/
+
+MKL_INT info = 0;
+MKL_INT job[8];   
+//job[1]=0; // zero based indexing in csr
+job[1]=1; // one based indexing in csr
+
+job[2]=0; // zero based indexing in coo
+job[3]=2; // I don't know
+job[4]=nnz; // nnz
+
+job[0]=1;  // coo to csr
+job[5]=0;  // Acsr and AJR allocated by user
+//void mkl_dcsrcoo (MKL_INT * job, MKL_INT * n, double *Acsr, MKL_INT * AJR, MKL_INT *AIR, MKL_INT * nnz, double *Acoo, MKL_INT * ir, MKL_INT * jc, MKL_INT * info);
+mkl_dcsrcoo (job,&m, Aval, AJ, AI, &nnz, val, I, J, &info);
+} /* ENDOF coo_to_csr }}}*/
+
+int main(int argc, char* argv[]) { /*{{{*/
+/** usage */
+int nrequired_args = 4;
+if (argc != nrequired_args){
+    fprintf(stderr, "NAME:\n\tmkl_spmv - multiply a sparse matrix with a dense vector\n");
+    fprintf(stderr, "\nSYNOPSIS:\n");
+    fprintf(stderr, "\tmkl_spmv MATRIX_A NUMBER_OF_THREADS PRINT_MATRIX\n");
+    fprintf(stderr, "\nDESCRIPTION:\n");
+    fprintf(stderr, "\tNUMBER_OF_THREADS: {0,1,2,...}\n");
+    fprintf(stderr, "\t\t0: Use number of threads determined by MKL\n");
+    fprintf(stderr, "\tPRINT_MATRIX: PRINT_YES, PRINT_NO\n");
+    fprintf(stderr, "\nSAMPLE EXECUTION:\n");
+    fprintf(stderr, "\t%s test.mtx 2 PRINT_YES\n", argv[0]);
+    exit(1);
+}
+/** parse arguments */
+int iarg = 1;
+char* strpathA = argv[iarg];    iarg++;
+int nthreads = atoi(argv[iarg]);    iarg++;
+if (nthreads > 0) {
+    mkl_set_num_threads(nthreads);
+} else {
+    nthreads = mkl_get_max_threads();
+}
+option_print_matrices = strcmp(argv[iarg], "PRINT_YES")==0?OPTION_PRINT_MATRICES:OPTION_NOPRINT_MATRICES;    iarg++;
+assert(nrequired_args == iarg);
+
+/** read matrix market file for A matrix */
+MKL_INT Am, An, Annz;
+MKL_INT *Ax, *Ay;  
+double *Anz;
+read_mm(strpathA, &Am, &An, &Annz, &Ax, &Ay, &Anz);
+
+/** construct csr storage for A matrix */
+MKL_INT* AJ = (MKL_INT*) mkl_malloc( Annz * sizeof( MKL_INT ), 64 );
+MKL_INT* AI = (MKL_INT*) mkl_malloc( (Am+1) * sizeof( MKL_INT ), 64 );
+double* Aval = (double*) mkl_malloc( Annz * sizeof( double ),  64 );
+coo_to_csr(Am, Annz, Ax, Ay, Anz, AI, AJ, Aval);    
+
+double* xvec = (double*) mkl_malloc( An * sizeof( double ), 64 ); 
+double* yvec = (double*) mkl_malloc( Am * sizeof( double ), 64 );
+int i;
+for(i=0;i<An;i++)xvec[i]=1.0;
+for(i=0;i<Am;i++)yvec[i]=0.0;
+
+double time; 
+mkl_cpu_spmv(Am, An, Aval, AJ, AI, xvec, yvec, &time);
+
+printmm_one(Am, Aval, AJ, AI);
+
+/** number of multiply-and-add operations in terms of giga flops*/
+double gflop = 2 * (double) Annz / 1e9;
+
+/** print gflop per second and time */
+printf("%d\t", nthreads);
+printf("%g\t", (gflop/time));
+printf("%g\t", time);
+printf("%s\t", strpathA);
+printf("\n");
+
+printf("Output vector:\n");
+for(i=0;i<Am;i++)printf("%g ", yvec[i]);
+printf("\n");
+/** free allocated space */
+mkl_free(AI);
+mkl_free(AJ);
+mkl_free(Aval);
+return 0;
+} /* ENDOF main }}}*/
+
+/** Prints matrix in CSR format */
+void printmm_one(int m, double* Aval, int* AJ, int* AI){ //{{{
+
+    if (option_print_matrices == OPTION_NOPRINT_MATRICES)
+        return;
+    int i;
+    for(i = 0; i < m; i++) {
+        printf("%d: ", i+1);
+        int j;
+        for(j = AI[i]-1; j < AI[i+1]-1; j++) {
+            printf("%d:%g  ", AJ[j], Aval[j]);
+        }
+        printf("\n");
+    }
+    printf("\n");
+}//}}}
+
+/** Writes matrix in CSR format in to a file using Matrix Market format */
+void printfilemm_one(char* file, int m, int n, double* Aval, int* AJ, int* AI){//{{{
+
+    FILE* f = fopen(file, "w");
+    if(f == NULL){
+        printf("%s %s %d: %s cannot be opened to write matrix\n", __FILE__, __PRETTY_FUNCTION__, __LINE__, file);
+        exit(1); 
+    }
+    int i;
+    fprintf(f, "%%%%MatrixMarket matrix coordinate real general\n");
+    fprintf(f, "%d %d %d\n", m, n, AI[m]-1);
+    for(i = 0; i < m; i++) {
+        int j;
+        for(j = AI[i]-1; j < AI[i+1]-1; j++) {
+            fprintf(f, "%d %d %g\n", i+1, AJ[j], Aval[j]);
+        }
+    }
+    fclose(f);
+}//}}}
+
diff --git a/spmv/mkl_shmem/mmio.c b/spmv/mkl_shmem/mmio.c
new file mode 100644
index 0000000..c250ff2
--- /dev/null
+++ b/spmv/mkl_shmem/mmio.c
@@ -0,0 +1,511 @@
+/* 
+*   Matrix Market I/O library for ANSI C
+*
+*   See http://math.nist.gov/MatrixMarket for details.
+*
+*
+*/
+
+
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+#include <ctype.h>
+
+#include "mmio.h"
+
+int mm_read_unsymmetric_sparse(const char *fname, int *M_, int *N_, int *nz_,
+                double **val_, int **I_, int **J_)
+{
+    FILE *f;
+    MM_typecode matcode;
+    int M, N, nz;
+    int i;
+    double *val;
+    int *I, *J;
+ 
+    if ((f = fopen(fname, "r")) == NULL)
+            return -1;
+ 
+ 
+    if (mm_read_banner(f, &matcode) != 0)
+    {
+        printf("mm_read_unsymetric: Could not process Matrix Market banner ");
+        printf(" in file [%s]\n", fname);
+        return -1;
+    }
+ 
+ 
+ 
+    if ( !(mm_is_real(matcode) && mm_is_matrix(matcode) &&
+            mm_is_sparse(matcode)))
+    {
+        fprintf(stderr, "Sorry, this application does not support ");
+        fprintf(stderr, "Market Market type: [%s]\n",
+                mm_typecode_to_str(matcode));
+        return -1;
+    }
+ 
+    /* find out size of sparse matrix: M, N, nz .... */
+ 
+    if (mm_read_mtx_crd_size(f, &M, &N, &nz) !=0)
+    {
+        fprintf(stderr, "read_unsymmetric_sparse(): could not parse matrix size.\n");
+        return -1;
+    }
+ 
+    *M_ = M;
+    *N_ = N;
+    *nz_ = nz;
+ 
+    /* reseve memory for matrices */
+ 
+    I = (int *) malloc(nz * sizeof(int));
+    J = (int *) malloc(nz * sizeof(int));
+    val = (double *) malloc(nz * sizeof(double));
+ 
+    *val_ = val;
+    *I_ = I;
+    *J_ = J;
+ 
+    /* NOTE: when reading in doubles, ANSI C requires the use of the "l"  */
+    /*   specifier as in "%lg", "%lf", "%le", otherwise errors will occur */
+    /*  (ANSI C X3.159-1989, Sec. 4.9.6.2, p. 136 lines 13-15)            */
+ 
+    for (i=0; i<nz; i++)
+    {
+        fscanf(f, "%d %d %lg\n", &I[i], &J[i], &val[i]);
+        I[i]--;  /* adjust from 1-based to 0-based */
+        J[i]--;
+    }
+    fclose(f);
+ 
+    return 0;
+}
+
+int mm_is_valid(MM_typecode matcode)
+{
+    if (!mm_is_matrix(matcode)) return 0;
+    if (mm_is_dense(matcode) && mm_is_pattern(matcode)) return 0;
+    if (mm_is_real(matcode) && mm_is_hermitian(matcode)) return 0;
+    if (mm_is_pattern(matcode) && (mm_is_hermitian(matcode) || 
+                mm_is_skew(matcode))) return 0;
+    return 1;
+}
+
+int mm_read_banner(FILE *f, MM_typecode *matcode)
+{
+    char line[MM_MAX_LINE_LENGTH];
+    char banner[MM_MAX_TOKEN_LENGTH];
+    char mtx[MM_MAX_TOKEN_LENGTH]; 
+    char crd[MM_MAX_TOKEN_LENGTH];
+    char data_type[MM_MAX_TOKEN_LENGTH];
+    char storage_scheme[MM_MAX_TOKEN_LENGTH];
+    char *p;
+
+
+    mm_clear_typecode(matcode);  
+
+    if (fgets(line, MM_MAX_LINE_LENGTH, f) == NULL) 
+        return MM_PREMATURE_EOF;
+
+    if (sscanf(line, "%s %s %s %s %s", banner, mtx, crd, data_type, 
+        storage_scheme) != 5)
+        return MM_PREMATURE_EOF;
+
+    for (p=mtx; *p!='\0'; *p=tolower(*p),p++);  /* convert to lower case */
+    for (p=crd; *p!='\0'; *p=tolower(*p),p++);  
+    for (p=data_type; *p!='\0'; *p=tolower(*p),p++);
+    for (p=storage_scheme; *p!='\0'; *p=tolower(*p),p++);
+
+    /* check for banner */
+    if (strncmp(banner, MatrixMarketBanner, strlen(MatrixMarketBanner)) != 0)
+        return MM_NO_HEADER;
+
+    /* first field should be "mtx" */
+    if (strcmp(mtx, MM_MTX_STR) != 0)
+        return  MM_UNSUPPORTED_TYPE;
+    mm_set_matrix(matcode);
+
+
+    /* second field describes whether this is a sparse matrix (in coordinate
+            storgae) or a dense array */
+
+
+    if (strcmp(crd, MM_SPARSE_STR) == 0)
+        mm_set_sparse(matcode);
+    else
+    if (strcmp(crd, MM_DENSE_STR) == 0)
+            mm_set_dense(matcode);
+    else
+        return MM_UNSUPPORTED_TYPE;
+    
+
+    /* third field */
+
+    if (strcmp(data_type, MM_REAL_STR) == 0)
+        mm_set_real(matcode);
+    else
+    if (strcmp(data_type, MM_COMPLEX_STR) == 0)
+        mm_set_complex(matcode);
+    else
+    if (strcmp(data_type, MM_PATTERN_STR) == 0)
+        mm_set_pattern(matcode);
+    else
+    if (strcmp(data_type, MM_INT_STR) == 0)
+        mm_set_integer(matcode);
+    else
+        return MM_UNSUPPORTED_TYPE;
+    
+
+    /* fourth field */
+
+    if (strcmp(storage_scheme, MM_GENERAL_STR) == 0)
+        mm_set_general(matcode);
+    else
+    if (strcmp(storage_scheme, MM_SYMM_STR) == 0)
+        mm_set_symmetric(matcode);
+    else
+    if (strcmp(storage_scheme, MM_HERM_STR) == 0)
+        mm_set_hermitian(matcode);
+    else
+    if (strcmp(storage_scheme, MM_SKEW_STR) == 0)
+        mm_set_skew(matcode);
+    else
+        return MM_UNSUPPORTED_TYPE;
+        
+
+    return 0;
+}
+
+int mm_write_mtx_crd_size(FILE *f, int M, int N, int nz)
+{
+    if (fprintf(f, "%d %d %d\n", M, N, nz) != 3)
+        return MM_COULD_NOT_WRITE_FILE;
+    else 
+        return 0;
+}
+
+int mm_read_mtx_crd_size(FILE *f, int *M, int *N, int *nz )
+{
+    char line[MM_MAX_LINE_LENGTH];
+    int num_items_read;
+
+    /* set return null parameter values, in case we exit with errors */
+    *M = *N = *nz = 0;
+
+    /* now continue scanning until you reach the end-of-comments */
+    do 
+    {
+        if (fgets(line,MM_MAX_LINE_LENGTH,f) == NULL) 
+            return MM_PREMATURE_EOF;
+    }while (line[0] == '%');
+
+    /* line[] is either blank or has M,N, nz */
+    if (sscanf(line, "%d %d %d", M, N, nz) == 3)
+        return 0;
+        
+    else
+    do
+    { 
+        num_items_read = fscanf(f, "%d %d %d", M, N, nz); 
+        if (num_items_read == EOF) return MM_PREMATURE_EOF;
+    }
+    while (num_items_read != 3);
+
+    return 0;
+}
+
+
+int mm_read_mtx_array_size(FILE *f, int *M, int *N)
+{
+    char line[MM_MAX_LINE_LENGTH];
+    int num_items_read;
+    /* set return null parameter values, in case we exit with errors */
+    *M = *N = 0;
+	
+    /* now continue scanning until you reach the end-of-comments */
+    do 
+    {
+        if (fgets(line,MM_MAX_LINE_LENGTH,f) == NULL) 
+            return MM_PREMATURE_EOF;
+    }while (line[0] == '%');
+
+    /* line[] is either blank or has M,N, nz */
+    if (sscanf(line, "%d %d", M, N) == 2)
+        return 0;
+        
+    else /* we have a blank line */
+    do
+    { 
+        num_items_read = fscanf(f, "%d %d", M, N); 
+        if (num_items_read == EOF) return MM_PREMATURE_EOF;
+    }
+    while (num_items_read != 2);
+
+    return 0;
+}
+
+int mm_write_mtx_array_size(FILE *f, int M, int N)
+{
+    if (fprintf(f, "%d %d\n", M, N) != 2)
+        return MM_COULD_NOT_WRITE_FILE;
+    else 
+        return 0;
+}
+
+
+
+/*-------------------------------------------------------------------------*/
+
+/******************************************************************/
+/* use when I[], J[], and val[]J, and val[] are already allocated */
+/******************************************************************/
+
+int mm_read_mtx_crd_data(FILE *f, int M, int N, int nz, int I[], int J[],
+        double val[], MM_typecode matcode)
+{
+    int i;
+    if (mm_is_complex(matcode))
+    {
+        for (i=0; i<nz; i++)
+            if (fscanf(f, "%d %d %lg %lg", &I[i], &J[i], &val[2*i], &val[2*i+1])
+                != 4) return MM_PREMATURE_EOF;
+    }
+    else if (mm_is_real(matcode))
+    {
+        for (i=0; i<nz; i++)
+        {
+            if (fscanf(f, "%d %d %lg\n", &I[i], &J[i], &val[i])
+                != 3) return MM_PREMATURE_EOF;
+
+        }
+    }
+
+    else if (mm_is_pattern(matcode))
+    {
+        for (i=0; i<nz; i++)
+            if (fscanf(f, "%d %d", &I[i], &J[i])
+                != 2) return MM_PREMATURE_EOF;
+    }
+    else
+        return MM_UNSUPPORTED_TYPE;
+
+    return 0;
+        
+}
+
+int mm_read_mtx_crd_entry(FILE *f, int *I, int *J,
+        double *real, double *imag, MM_typecode matcode)
+{
+    if (mm_is_complex(matcode))
+    {
+            if (fscanf(f, "%d %d %lg %lg", I, J, real, imag)
+                != 4) return MM_PREMATURE_EOF;
+    }
+    else if (mm_is_real(matcode))
+    {
+            if (fscanf(f, "%d %d %lg\n", I, J, real)
+                != 3) return MM_PREMATURE_EOF;
+
+    }
+
+    else if (mm_is_pattern(matcode))
+    {
+            if (fscanf(f, "%d %d", I, J) != 2) return MM_PREMATURE_EOF;
+    }
+    else
+        return MM_UNSUPPORTED_TYPE;
+
+    return 0;
+        
+}
+
+
+/************************************************************************
+    mm_read_mtx_crd()  fills M, N, nz, array of values, and return
+                        type code, e.g. 'MCRS'
+
+                        if matrix is complex, values[] is of size 2*nz,
+                            (nz pairs of real/imaginary values)
+************************************************************************/
+
+int mm_read_mtx_crd(char *fname, int *M, int *N, int *nz, int **I, int **J, 
+        double **val, MM_typecode *matcode)
+{
+    int ret_code;
+    FILE *f;
+
+    if (strcmp(fname, "stdin") == 0) f=stdin;
+    else
+    if ((f = fopen(fname, "r")) == NULL)
+        return MM_COULD_NOT_READ_FILE;
+
+
+    if ((ret_code = mm_read_banner(f, matcode)) != 0)
+        return ret_code;
+
+    if (!(mm_is_valid(*matcode) && mm_is_sparse(*matcode) && 
+            mm_is_matrix(*matcode)))
+        return MM_UNSUPPORTED_TYPE;
+
+    if ((ret_code = mm_read_mtx_crd_size(f, M, N, nz)) != 0)
+        return ret_code;
+
+
+    *I = (int *)  malloc(*nz * sizeof(int));
+    *J = (int *)  malloc(*nz * sizeof(int));
+    *val = NULL;
+
+    if (mm_is_complex(*matcode))
+    {
+        *val = (double *) malloc(*nz * 2 * sizeof(double));
+        ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val, 
+                *matcode);
+        if (ret_code != 0) return ret_code;
+    }
+    else if (mm_is_real(*matcode))
+    {
+        *val = (double *) malloc(*nz * sizeof(double));
+        ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val, 
+                *matcode);
+        if (ret_code != 0) return ret_code;
+    }
+
+    else if (mm_is_pattern(*matcode))
+    {
+        ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val, 
+                *matcode);
+        if (ret_code != 0) return ret_code;
+    }
+
+    if (f != stdin) fclose(f);
+    return 0;
+}
+
+int mm_write_banner(FILE *f, MM_typecode matcode)
+{
+    char *str = mm_typecode_to_str(matcode);
+    int ret_code;
+
+    ret_code = fprintf(f, "%s %s\n", MatrixMarketBanner, str);
+    free(str);
+    if (ret_code !=2 )
+        return MM_COULD_NOT_WRITE_FILE;
+    else
+        return 0;
+}
+
+int mm_write_mtx_crd(char fname[], int M, int N, int nz, int I[], int J[],
+        double val[], MM_typecode matcode)
+{
+    FILE *f;
+    int i;
+
+    if (strcmp(fname, "stdout") == 0) 
+        f = stdout;
+    else
+    if ((f = fopen(fname, "w")) == NULL)
+        return MM_COULD_NOT_WRITE_FILE;
+    
+    /* print banner followed by typecode */
+    fprintf(f, "%s ", MatrixMarketBanner);
+    fprintf(f, "%s\n", mm_typecode_to_str(matcode));
+
+    /* print matrix sizes and nonzeros */
+    fprintf(f, "%d %d %d\n", M, N, nz);
+
+    /* print values */
+    if (mm_is_pattern(matcode))
+        for (i=0; i<nz; i++)
+            fprintf(f, "%d %d\n", I[i], J[i]);
+    else
+    if (mm_is_real(matcode))
+        for (i=0; i<nz; i++)
+            fprintf(f, "%d %d %20.16g\n", I[i], J[i], val[i]);
+    else
+    if (mm_is_complex(matcode))
+        for (i=0; i<nz; i++)
+            fprintf(f, "%d %d %20.16g %20.16g\n", I[i], J[i], val[2*i], 
+                        val[2*i+1]);
+    else
+    {
+        if (f != stdout) fclose(f);
+        return MM_UNSUPPORTED_TYPE;
+    }
+
+    if (f !=stdout) fclose(f);
+
+    return 0;
+}
+  
+
+/**
+*  Create a new copy of a string s.  mm_strdup() is a common routine, but
+*  not part of ANSI C, so it is included here.  Used by mm_typecode_to_str().
+*
+*/
+char *mm_strdup(const char *s)
+{
+	int len = strlen(s);
+	char *s2 = (char *) malloc((len+1)*sizeof(char));
+	return strcpy(s2, s);
+}
+
+char  *mm_typecode_to_str(MM_typecode matcode)
+{
+    char buffer[MM_MAX_LINE_LENGTH];
+    char *types[4];
+	char *mm_strdup(const char *);
+    int error =0;
+
+    /* check for MTX type */
+    if (mm_is_matrix(matcode)) 
+        types[0] = MM_MTX_STR;
+    else
+        error=1;
+
+    /* check for CRD or ARR matrix */
+    if (mm_is_sparse(matcode))
+        types[1] = MM_SPARSE_STR;
+    else
+    if (mm_is_dense(matcode))
+        types[1] = MM_DENSE_STR;
+    else
+        return NULL;
+
+    /* check for element data type */
+    if (mm_is_real(matcode))
+        types[2] = MM_REAL_STR;
+    else
+    if (mm_is_complex(matcode))
+        types[2] = MM_COMPLEX_STR;
+    else
+    if (mm_is_pattern(matcode))
+        types[2] = MM_PATTERN_STR;
+    else
+    if (mm_is_integer(matcode))
+        types[2] = MM_INT_STR;
+    else
+        return NULL;
+
+
+    /* check for symmetry type */
+    if (mm_is_general(matcode))
+        types[3] = MM_GENERAL_STR;
+    else
+    if (mm_is_symmetric(matcode))
+        types[3] = MM_SYMM_STR;
+    else 
+    if (mm_is_hermitian(matcode))
+        types[3] = MM_HERM_STR;
+    else 
+    if (mm_is_skew(matcode))
+        types[3] = MM_SKEW_STR;
+    else
+        return NULL;
+
+    sprintf(buffer,"%s %s %s %s", types[0], types[1], types[2], types[3]);
+    return mm_strdup(buffer);
+
+}
diff --git a/spmv/mkl_shmem/mmio.h b/spmv/mkl_shmem/mmio.h
new file mode 100644
index 0000000..7cfd0a1
--- /dev/null
+++ b/spmv/mkl_shmem/mmio.h
@@ -0,0 +1,133 @@
+/* 
+*   Matrix Market I/O library for ANSI C
+*
+*   See http://math.nist.gov/MatrixMarket for details.
+*
+*
+*/
+
+#ifndef MM_IO_H
+#define MM_IO_H
+
+#define MM_MAX_LINE_LENGTH 1025
+#define MatrixMarketBanner "%%MatrixMarket"
+#define MM_MAX_TOKEN_LENGTH 64
+
+typedef char MM_typecode[4];
+
+char *mm_typecode_to_str(MM_typecode matcode);
+
+int mm_read_banner(FILE *f, MM_typecode *matcode);
+int mm_read_mtx_crd_size(FILE *f, int *M, int *N, int *nz);
+int mm_read_mtx_array_size(FILE *f, int *M, int *N);
+
+int mm_write_banner(FILE *f, MM_typecode matcode);
+int mm_write_mtx_crd_size(FILE *f, int M, int N, int nz);
+int mm_write_mtx_array_size(FILE *f, int M, int N);
+
+
+/********************* MM_typecode query fucntions ***************************/
+
+#define mm_is_matrix(typecode)	((typecode)[0]=='M')
+
+#define mm_is_sparse(typecode)	((typecode)[1]=='C')
+#define mm_is_coordinate(typecode)((typecode)[1]=='C')
+#define mm_is_dense(typecode)	((typecode)[1]=='A')
+#define mm_is_array(typecode)	((typecode)[1]=='A')
+
+#define mm_is_complex(typecode)	((typecode)[2]=='C')
+#define mm_is_real(typecode)		((typecode)[2]=='R')
+#define mm_is_pattern(typecode)	((typecode)[2]=='P')
+#define mm_is_integer(typecode) ((typecode)[2]=='I')
+
+#define mm_is_symmetric(typecode)((typecode)[3]=='S')
+#define mm_is_general(typecode)	((typecode)[3]=='G')
+#define mm_is_skew(typecode)	((typecode)[3]=='K')
+#define mm_is_hermitian(typecode)((typecode)[3]=='H')
+
+int mm_is_valid(MM_typecode matcode);		/* too complex for a macro */
+
+
+/********************* MM_typecode modify fucntions ***************************/
+
+#define mm_set_matrix(typecode)	((*typecode)[0]='M')
+#define mm_set_coordinate(typecode)	((*typecode)[1]='C')
+#define mm_set_array(typecode)	((*typecode)[1]='A')
+#define mm_set_dense(typecode)	mm_set_array(typecode)
+#define mm_set_sparse(typecode)	mm_set_coordinate(typecode)
+
+#define mm_set_complex(typecode)((*typecode)[2]='C')
+#define mm_set_real(typecode)	((*typecode)[2]='R')
+#define mm_set_pattern(typecode)((*typecode)[2]='P')
+#define mm_set_integer(typecode)((*typecode)[2]='I')
+
+
+#define mm_set_symmetric(typecode)((*typecode)[3]='S')
+#define mm_set_general(typecode)((*typecode)[3]='G')
+#define mm_set_skew(typecode)	((*typecode)[3]='K')
+#define mm_set_hermitian(typecode)((*typecode)[3]='H')
+
+#define mm_clear_typecode(typecode) ((*typecode)[0]=(*typecode)[1]= \
+									(*typecode)[2]=' ',(*typecode)[3]='G')
+
+#define mm_initialize_typecode(typecode) mm_clear_typecode(typecode)
+
+
+/********************* Matrix Market error codes ***************************/
+
+
+#define MM_COULD_NOT_READ_FILE	11
+#define MM_PREMATURE_EOF		12
+#define MM_NOT_MTX				13
+#define MM_NO_HEADER			14
+#define MM_UNSUPPORTED_TYPE		15
+#define MM_LINE_TOO_LONG		16
+#define MM_COULD_NOT_WRITE_FILE	17
+
+
+/******************** Matrix Market internal definitions ********************
+
+   MM_matrix_typecode: 4-character sequence
+
+				    ojbect 		sparse/   	data        storage 
+						  		dense     	type        scheme
+
+   string position:	 [0]        [1]			[2]         [3]
+
+   Matrix typecode:  M(atrix)  C(oord)		R(eal)   	G(eneral)
+						        A(array)	C(omplex)   H(ermitian)
+											P(attern)   S(ymmetric)
+								    		I(nteger)	K(kew)
+
+ ***********************************************************************/
+
+#define MM_MTX_STR		"matrix"
+#define MM_ARRAY_STR	"array"
+#define MM_DENSE_STR	"array"
+#define MM_COORDINATE_STR "coordinate" 
+#define MM_SPARSE_STR	"coordinate"
+#define MM_COMPLEX_STR	"complex"
+#define MM_REAL_STR		"real"
+#define MM_INT_STR		"integer"
+#define MM_GENERAL_STR  "general"
+#define MM_SYMM_STR		"symmetric"
+#define MM_HERM_STR		"hermitian"
+#define MM_SKEW_STR		"skew-symmetric"
+#define MM_PATTERN_STR  "pattern"
+
+
+/*  high level routines */
+
+int mm_write_mtx_crd(char fname[], int M, int N, int nz, int I[], int J[],
+		 double val[], MM_typecode matcode);
+int mm_read_mtx_crd_data(FILE *f, int M, int N, int nz, int I[], int J[],
+		double val[], MM_typecode matcode);
+int mm_read_mtx_crd_entry(FILE *f, int *I, int *J, double *real, double *img,
+			MM_typecode matcode);
+
+int mm_read_unsymmetric_sparse(const char *fname, int *M_, int *N_, int *nz_,
+                double **val_, int **I_, int **J_);
+
+
+
+#endif
diff --git a/spmv/mkl_shmem/test.mtx b/spmv/mkl_shmem/test.mtx
new file mode 100644
index 0000000..638b954
--- /dev/null
+++ b/spmv/mkl_shmem/test.mtx
@@ -0,0 +1,13 @@
+%%MatrixMarket matrix coordinate real general
+5 5 11
+1 2 1
+1 3 2
+1 5 3
+2 3 4
+2 4 5
+2 5 6
+3 1 7
+3 2 8
+3 4 9
+3 5 10
+4 4 11
-- 
GitLab