/*

  C = A*B / row-column 
  
*/
#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define N 1024
#define N_THREADS 4

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

	/* set initial value */
  	#pragma omp parallel shared(A,B,C) private(j,k) num_threads(N_THREADS) 
  	{
	    #pragma omp for
            for(i=0; i<N; i++)
                for(j=0; j<N; j++)
                    A[i][j] = 1.0/(i+j+1);  	
	    #pragma omp for
            for(i=0; i<N; i++)
                for(j=0; j<N; j++)
                    B[i][j] = i+j+1;
    }
    
    // serial
  	t = (double)clock()/CLOCKS_PER_SEC;
  	for(i=0; i<N; i++)
        for(j=0; j<N; j++)
        {
            C[i][j]=0;
            for(k=0; k<N; 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 : row-column 
  	t = (double)clock()/CLOCKS_PER_SEC;
  	#pragma omp parallel shared(A,B,C) private(i,j,k) num_threads(N_THREADS) 
    {
        tid = omp_get_thread_num(); 
        for(i=tid; i<N; i+=N_THREADS)
            for(j=0; j<N; j++)
            {
                C[i][j]=0;
                for(k=0; k<N; 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;
}