Commit c4c72685 authored by tempy@maggie's avatar tempy@maggie
Browse files

SpMV using MKL for shared memory has been added

parent 30b1861a
# ==================================================================================================
# 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)
# ==================================================================================================
/*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);
}//}}}
/*
* 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))
{