// omp_matmul.c  number of threads set using	export OMP_NUM_THREADS=4

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

#define DIM 2048  // why limit ? C handles 2048 4096

int main(int argc, char *argv[])
{
  int NCPU, tid, NPR, NTHR;
  int i, j, k;
  double a[DIM][DIM];
  double b[DIM][DIM];
  double c[DIM][DIM];
  double sum;
  double twall, tcpu;

  // initialize before going parallel
  printf("omp_matmu.c running, DIM=%d \n", DIM);
  for(i=0; i<DIM; i++)
  {
    for(j=0; j<DIM; j++)
    {
      a[i][j] = (double)(i+j);
      b[i][j] = (double)(i+j+1);
    } // end j
  } // end i
  printf("initialization finished \n");

  #pragma omp parallel
  {
    NCPU = omp_get_num_procs();
    tid = omp_get_thread_num();
    NPR = omp_get_num_threads();
    NTHR = omp_get_max_threads();
    if (tid == 0) 
    {
      printf("omp_matmu.c running \n");
      // tid is current thread running
      printf("tid=%i : available CPU     NCPU = %i\n",tid,NCPU);
      printf("tid=%i : available threads NTHR = %i\n",tid,NTHR);
      printf("tid=%i : threads requested NPR  = %i\n",tid,NPR);
    } // end if

    twall = (double)time(NULL);
    tcpu = (double)clock()/(double)CLOCKS_PER_SEC;

    //  firstprivate(DIM) did not work at end of next line
#pragma omp parallel for private(i,j,k,sum) shared(a,b,c)
    for(i=0; i<DIM; i++)
    {
      for(j=0; j<DIM; j++)
      {
        sum=0.0;
        for(k=0; k<DIM; k++)
        {
          sum += a[i][k]*b[k][j];
        } // end k
        c[i][j]=sum;
      } // end j
    } // end i
  } // end parallel

  printf("c[1][1]=%f, c[2][2]=%f \n", c[1][1], c[2][2]);  

  printf("took %f wall sec, %f cpu sec \n",
         (double)time(NULL)-twall,
         (double)clock()/(double)CLOCKS_PER_SEC-tcpu);
  printf("average number of parallel threads %f\n",
         ((double)clock()/(double)CLOCKS_PER_SEC-tcpu)/
         ((double)time(NULL)-twall));

  return(0);
} // end omp_matmul.c