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

#define NRA 500              /* number of rows in matrix A */
#define NCA 500               /* number of columns in matrix A */
#define NCB 500                /* number of columns in matrix B */

void seq_matrix(double *a, double *b, double *d)
{
  int i, j, k;

  for (i=0; i<NRA; i++)    
    {
    
    for(j=0; j<NCB; j++)       
      for (k=0; k<NCA; k++)
        //d[i][j] += a[i][k] * b[k][j];
        *((d + i * NCB) + j) = *((a + i * NCA) + k) * *((b + k * NCB) + j);

    }
}



int main (int argc, char *argv[]) 
{
int	tid, nthreads, i, j, k, chunk;
double	a[NRA][NCA],           /* matrix A to be multiplied */
	b[NCA][NCB],           /* matrix B to be multiplied */
	c[NRA][NCB] , wtime, d[NRA][NCB];           /* result matrix C */

chunk = 10;                    /* set loop iteration chunk size */

/*** Spawn a parallel region explicitly scoping all variables ***/
#pragma omp parallel shared(a,b,c,nthreads,chunk) private(tid,i,j,k)
  {
  tid = omp_get_thread_num();
  if (tid == 0)
    {
    nthreads = omp_get_num_threads();
    printf("Starting matrix multiple example with %d threads\n",nthreads);
    printf("Initializing matrices...\n");
    }
  

  /*** Initialize matrices ***/
  #pragma omp for schedule (static, chunk) 
  for (i=0; i<NRA; i++)
    for (j=0; j<NCA; j++)
      a[i][j]= 1;
  #pragma omp for schedule (static, chunk)
  for (i=0; i<NCA; i++)
    for (j=0; j<NCB; j++)
      b[i][j]= 1;
  #pragma omp for schedule (static, chunk)
  for (i=0; i<NRA; i++)
    for (j=0; j<NCB; j++)
      c[i][j]= 0;
  #pragma omp for schedule (static, chunk)
  for (i=0; i<NRA; i++)
    for (j=0; j<NCB; j++)
      d[i][j]= 0;

  

  /*** Do matrix multiply sharing iterations on outer loop ***/
  /*** Display who does which iterations for demonstration purposes ***/
  printf("Thread %d starting matrix multiply...\n",tid);
  


  wtime = omp_get_wtime ( );
  #pragma omp for schedule (static, chunk)
  for (i=0; i<NRA; i++)    
    {
   //printf("Thread=%d did row=%d\n",tid,i);
    for(j=0; j<NCB; j++)       
      for (k=0; k<NCA; k++)
        c[i][j] += a[i][k] * b[k][j];
    }
  }   /*** End of parallel region ***/

  wtime = omp_get_wtime ( ) - wtime;



/*** Print results ***/
if(NCB < 16 && NCA < 16 && NCB < 16) {
printf("******************************************************\n");
printf("Result Matrix:\n");
for (i=0; i<NRA; i++)
  {
  for (j=0; j<NCB; j++) 
    printf("%6.2f   ", c[i][j]);
  printf("\n"); 
  }
printf("******************************************************\n");
printf ("Done.\n");
}

printf(" \ntime taken for parallel running %lf\n", wtime);

//sequential part
wtime = omp_get_wtime();
seq_matrix((double *)a,(double *)b,(double *)d);
wtime = omp_get_wtime ( ) - wtime;
if(NCB < 16 && NCA < 16 && NCB < 16) {
for (i=0; i<NRA; i++)
  {
  for (j=0; j<NCB; j++) 
    printf("%6.2f   ", c[i][j]);
  printf("\n"); 
  }
}
  printf("\ntime taken for sequential running %6.8lf\n", wtime);

}
