/*
  A simplification of the Stream benchmark for OpenMP
  The array for each thread is a large distance from the array for the other threads
  Original code developed by John D. McCalpin
*/
#include <stdio.h>
#include <math.h>
#include <limits.h>
#include <float.h>
#include <sys/time.h>
#include <stdlib.h>
#include <omp.h>
#include <petscsys.h>

#define NTIMESINNER 1
#define N 2*4*20000000
//#define N 1200000
//#define N 120000
#define NTIMES       50
#define OFFSET       0

# if !defined(MIN)
# define MIN(x,y) ((x)<(y) ? (x) : (y))
# endif
# if !defined(MAX)
# define MAX(x,y) ((x)>(y) ? (x) : (y))
# endif

static double a[64][10000000],b[64][10000000],c[64][10000000];
static double  mintime = FLT_MAX;
static double bytes =  3 * sizeof(double) * N;

int main()
{
  const static double scalar = 3.0;
#pragma omp threadprivate(scalar)
  double       times[NTIMES],rate;
  int          size;
  static int   n;
#pragma omp threadprivate(n)
  char         *env;
  FILE         *fd;

  env = getenv("OMP_NUM_THREADS");
  if (!env) env = (char *) "1";
  sscanf(env,"%d",&size);

#pragma omp parallel for schedule(static)
  for (int j=0; j<size; j++) {
    n = (N / size + ((N % size) > omp_get_thread_num()));
    for (int i=0; i<n; i++){
      a[j][i] = 1.0;
      b[j][i] = 2.0;
      c[j][i] = 3.0;
    }
  }

  /*  --- MAIN LOOP --- repeat test cases NTIMES times --- */
  for (int k=0; k<NTIMES; k++)
  {
    times[k] = MPI_Wtime();
    // https://www.openmp.org/wp-content/uploads/OpenMP-API-Specification-5-2.pdf
// #pragma omp parallel for  (same performance as below)
// #pragma omp parallel for simd schedule(static)  (same performance as below)
#pragma omp parallel for schedule(static)
    for (int j=0; j<size; j++) {
      n = (N / size + ((N % size) > omp_get_thread_num()));
      double *aa = a[j];  // these don't change the timings
      const double *bb = b[j];
      const double *cc = c[j];
      for (int l=0; l<NTIMESINNER; l++) {
        for (register int i=0; i<n; i++) aa[i] = bb[i]+scalar*cc[i];
          if (size == 65) printf("never printed %g\n",a[0][11]);
      }
    }
    times[k] = MPI_Wtime() - times[k];
  }
  for (int k=1; k<NTIMES; k++) {  /* note -- skip first iteration */
      mintime = MIN(mintime, times[k]);
  }
  rate = 1.0E-06 * bytes*NTIMESINNER/mintime;

  if (size == 1) {
    printf("%d %11.4f   Rate (MB/s) 1\n",size, rate);
    fd = fopen("flops","w");
    fprintf(fd,"%g\n",rate);
    fclose(fd);
  } else {
    double prate;
    fd = fopen("flops","r");
    fscanf(fd,"%lg",&prate);
    fclose(fd);
    printf("%d %11.4f   Rate (MB/s) %g \n", size, rate,rate/prate);
  }
  return 0;
}