/*  --------------------------------------------------- !
! The current OpenMP code has three purposes:           !
! (1) how to parallalize nested loops.                  !
! (2) matrix-matrix product.                            !
! (3) algorithm's importance.                           !
!-------------------------------------------------------!
! Xiaoxu Guan at LSU HPC.                               !
! March 07, 2016                                        !
!----------------------------------------------------- */

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


int main () {

   int i,j,k,nothread,nosize,index,inda,indb,jn,kn;
   const int nsize = 6000;
   double temp,time_start,time_end;

   double **a,**b,**c;
  

   nosize = nsize;
   
   a = malloc( nsize * sizeof(double *) );
   b = malloc( nsize * sizeof(double *) );
   c = malloc( nsize * sizeof(double *) );
   if(a == NULL) printf("Failure to allocate array a!\n");
   if(b == NULL) printf("Failure to allocate array b!\n");
   if(c == NULL) printf("Failure to allocate array c!\n");


   for(i = 0; i < nosize; i++) 
   {
      a[i] = malloc( nosize * sizeof(double) ); 
      if(a[i] == NULL) { printf("Failure to allocate array element a[i] %d\n",i); exit(0); }

      b[i] = malloc( nosize * sizeof(double) ); 
      if(b[i] == NULL) { printf("Failure to allocate array element b[i] %d\n",i); exit(0); }

      c[i] = malloc( nosize * sizeof(double) ); 
      if(c[i] == NULL) { printf("Failure to allocate array element c[i] %d\n",i); exit(0); }
   }


/* initialize arrays: a, b, and c:   */
   for(i = 0; i < nosize; i++)  {
   for(j = 0; j < nosize; j++)  {
       a[i][j] = sin( (double) (i-j) );
       b[i][j] = cos( (double) (i+j) );
       c[i][j] = 0.0;
     } }


time_start = omp_get_wtime();

/* only parallelize the outer loop with kji: */
#pragma omp parallel for default(none) private(i,j,k) shared(nosize,a,b,c)
   for(k = 0; k < nosize; k++)  { 
   for(j = 0; j < nosize; j++)  { 
       for(i = 0; i < nosize; i++)  c[k][i] +=  a[k][j] * b[j][k]; 
   } } 



time_end = omp_get_wtime();

printf("elapsed time (second) = %.5f\n", time_end - time_start);


}
