/*  --------------------------------------------------- !
! The current OpenMP code has three purposes:           !
! (1) how to parallalize nested loops.                  !
! (2) matrix-matrix product.                            !
! (3) which loop should or should not be parallelized.  !
!-------------------------------------------------------!
! 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,id,nothread,nosize,index,inda,indb;
   const int nsize = 6000;
   double temp,time_start,time_end;
   double a[nsize*nsize-1],b[nsize*nsize-1],c[nsize*nsize-1];


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


time_start = omp_get_wtime();

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


time_end = omp_get_wtime();

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


}
