#include <omp.h>
#include <stdlib.h>
#include <stdio.h>
#include <time.h>

#define N 2048

void initMatrix(float* A) {
	
#pragma omp parallel for
	for (int ii=0; ii<N; ii++)
		for (int jj=0; jj<N; jj++) 
			A[ii*N+jj] = (float)rand()/(float)RAND_MAX;
	
}

void clearMatrix(float* A) {
	
#pragma omp parallel for
	for (int ii=0; ii<N; ii++)
		for (int jj=0; jj<N; jj++) 
			A[ii*N+jj] = 0.0;
	
}

void printMatrix(float* A) {
	for (int ii=0; ii<N; ii++)
		for (int jj=0; jj<N; jj++) 
			printf("%f\t",A[ii*N+jj]);
	
}

void thrMM(float* A, float* B, float* C) {
	// Note: OpenMP allows
	//    more robust workload distribution among threads
	//    assign threads to specific cores (affinity)
	
	#pragma omp parallel for
	for (int ii=0; ii<N; ii++) {
		for (int kk=0; kk<N; kk++) {
			float r =  A[ii*N+kk];
			for (int jj=0; jj<N; jj++) {
				C[ii*N+jj] += r * B[kk*N+jj]; 
			}
		}
	}
}


int main (int argc, const char * argv[]) {

	float *A = (float*)malloc(sizeof(float)*N*N);
	float *B = (float*)malloc(sizeof(float)*N*N);
	float *C = (float*)malloc(sizeof(float)*N*N);
	
	// The #threads should be dynamically set and commented
	omp_set_num_threads(2);
	
	initMatrix(A);
	initMatrix(B);
        clearMatrix(C);         
	
	//Be carefull with time measuring (avg?????)
	double start = omp_get_wtime();
	for (int i=0; i<5; i++) {
		thrMM(A, B, C);
		clearMatrix(C); 
	}
	double end = omp_get_wtime();
	
	//Use a more adequate metric to measure performance 
	printf("Avg. time: %f\n",(end-start)/5);
	
}