/*

  C = A*B / parallel for 
  
*/
#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define M 4*1024/2
#define N 3*1024/2 // 1024*1024*1024/8
#define L 2*1024/2
#define N_THREADS 4

int main()
{
  	int i, j, k;  	
  	static double A[M][L], B[L][N], C[M][N]; 
  	double t;
  	

	/* set initial value */
	for(i=0; i<M; i++)
        for(j=0; j<L; j++)
            A[i][j] = 1.0/(i+j+1);  	
	for(i=0; i<L; i++)
        for(j=0; j<N; j++)
            B[i][j] = i+j+1;
            
    // serial
  	t=(double)clock()/CLOCKS_PER_SEC;
  	for(i=0; i<M; i++)
        for(j=0; j<N; j++)
        {
            C[i][j]=0;
            for(k=0; k<L; k++)
                C[i][j] += A[i][k]*B[k][j];
        }
  	t=(double)clock()/CLOCKS_PER_SEC - t;
  	printf("Serial: C[0][0]=%.2f, elapsed time %.2e\n", C[0][0], t);
  	
  	// parallel
  	t=(double)clock()/CLOCKS_PER_SEC;
  	#pragma omp parallel for shared(A,B,C) private(j,k) num_threads(N_THREADS) 
        for(i=0; i<M; i++)
            for(j=0; j<N; j++)
            {
                C[i][j]=0;
                for(k=0; k<L; k++)
                    C[i][j] += A[i][k]*B[k][j];
            }
  	t=(double)clock()/CLOCKS_PER_SEC - t;
  	printf("Parallel(%d): C[0][0]=%.2f, elapsed time %.2e\n", N_THREADS, C[0][0], t);
  	
  	return 0;
}