Newer
Older
* A first-touch/affinity frendly OpenMP implementation of the Matrix-Vector multiplication
*
* Author: Petros Anastasiadis(panastas@cslab.ece.ntua.gr)
*
* For more info about OpenMP thread affinity see http://pages.tacc.utexas.edu/~eijkhout/pcse/html/omp-affinity.html
*/
#include <errno.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <omp.h>
/* Need to include External_Functions for these */
#include "matrix_op.h"
#include "util.h"
#include "input.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");
double *x = (double *) malloc(m * sizeof(*x));
double *y = (double *) malloc(n * sizeof(*y));
double *M = (double *) malloc(n * m * sizeof(*M));
#pragma omp parallel for schedule(static) /* Initialize data for each thread in corresponding socket/cache with first-touch policy */
for( i=0 ; i<n ; ++i){
for ( j=0 ; j<m ; ++j) M[i*m+j]=0.0;
//printf( "Initialize data Thread=%d i=%d\n", omp_get_thread_num(), i);
}
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);
/* OpenMP Affinity Kernel */
printf("OpenMP_aff Version(N=%d, M=%d, Threads=%s): ", n, m, getenv("OMP_NUM_THREADS"));
timer = csecond();
for (i = 0; i < NR_ITER; ++i){
register double yi = 0;
#pragma omp parallel for private(j,yi) shared(n,m,M,y) schedule(static) /* Each thread computes n/thread_num contiguous elements of y */
for (k = 0; k < n; ++k) {
//printf( "Compute Thread=%d i=%d\n", omp_get_thread_num(), k);
yi = 0.0;
for (j = 0; j < m; ++j) yi += M[k*m+j]*x[j];
y[k] = yi;
}
}
timer = csecond() - timer;
#ifdef _DEBUG_
/* Output y vector to a file for debugging */
FILE * fp;
char * filename = "OpenMP_aff.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
report_results(timer);
return 0;
}