/*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);
}//}}}