/*
 * A Serial implementation of the Matrix-Vector multiplication
 * 
 * Author: Petros Anastasiadis(panastas@cslab.ece.ntua.gr) 
 */

#include <errno.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
/* Need to include External_Functions for these */
#include "matrix_op.h"
#include "util.h"
#include "input.h"
#include "mkl.h"
#include "mkl_blas.h"

int main(int argc, char **argv)
{
	/* Initializations */
	int i, j, k, n, m;
	double timer;

	if (argc < 3) error("Usage: ./Program N M");
	else if ( argc == 3) { /*./Program N M */
		n = atoi(argv[1]);
		m = atoi(argv[2]);		
	}
	else error("Too many Arguments");
	
	/* Allocate space */
	double *x 	= (double *) malloc(m * sizeof(*x));
	double *y	= (double *) malloc(n * sizeof(*y));
	double *M 			= (double *) malloc(n * m * sizeof(*M));
	if( !y || !x || !M ) error("memory allocation failed");

	/* Initialize matrices */
	ser_matrix_init_rand(M,n,m,1.0); /* Normal matrices generated randomly */

	/* Initialize vectors */
	vec_init_rand(x, m, 1.0);
	vec_init(y, n, 0.0);


	/* Serial Kernel */
	printf("Serial Version(N=%d, M=%d): ", n, m);
	timer = csecond();
	for (i = 0; i < NR_ITER; ++i){
		register double yi;
		for (k = 0; k < n; ++k) {
        	yi = 0.0 ;
        	for (j = 0; j < m; ++j) yi += M[n*k+j]*x[j];
        	y[k] = yi;
    	}
	}
	timer = csecond() - timer ;
	report_results(timer);

	/* BLAS Kernel */
	printf("BLAS dgemv Version(N=%d, M=%d): ", n, m);
	const double  a=1.0,b=0.0;
	const char trans='N';
	const int inc=1;
	timer = csecond();
	for (i = 0; i < NR_ITER; ++i){
		dgemv_(&trans, &n, &m, &a, M, &n, x, &inc, &b, y, &inc);
	}
	timer = csecond() - timer ;
	report_results(timer);

#ifdef _DEBUG_
	/* Output y vector to a file for debugging */
	FILE * fp;
	char filename[] = "Serial.debug" ;
	if(( fp = fopen( filename, "w")) == NULL)  error("Output file creation failed\n");
	for (k = 0; k < n; ++k) fprintf(fp, "%lf ", y[k]) ;
	fclose(fp) ;
#endif

	return 0;
}