//
// row-block 
//
#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#define m 8192
#define n 8192

const int num_test=10; // times for testing

int main()
{
    int i, j, k;
	clock_t t0, t1, t2;
    static double A[m][n], x[n], y1[m], y2[m];

    srand(time(0)); 
    for(i=0; i<m; i++)
        for(j=0; j<n; j++)
            A[i][j] = (double)rand()/RAND_MAX;

	for(j=0; j<n; j++)
	   x[j] = (double)rand()/RAND_MAX;
	               
	// serial
	t0 = clock();
	for(k=0; k<num_test; k++)  
    for(i=0; i<m; i++)
    {
        y1[i] = 0.0;
        for(j=0; j<n; j++)
            y1[i] += A[i][j] * x[j]; 
	}
    t1 = clock() - t0;
    printf("Serial: %.2f\n", (double)t1/CLOCKS_PER_SEC);	    

	// parallel 
	int tid, nthreads;
	t0 = clock();
	for(k=0; k<num_test; k++) 
	{		
	#pragma omp parallel shared(A,x,y2) private(i,j,tid) num_threads(4) 
	{
		tid = omp_get_thread_num();
		nthreads = omp_get_num_threads();
    	for(i=tid; i<m; i+=nthreads)
    	{
        	y2[i] = 0.0;
        	for(j=0; j<n; j++)
            	y2[i] += A[i][j] * x[j];
    	}
	}
	}
	t2 = clock() - t0;
    printf("Parallel: %.2f\n", (double)t2/CLOCKS_PER_SEC);
    printf("Speedup: %.2f\n", (double)t1/t2);

	// check
	double err=0.0;
	for(i=0; i<m; i++)
		err += fabs(y1[i]-y2[i]);
	printf("Error: %.2e\n", err);

    return 0;
}