/*
  A simplification of the Stream benchmark for OpenMP
  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 <petscsys.h>

//#define N 2*4*20000000
#define N 80000000
//#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[N + OFFSET], b[N + OFFSET], c[N + OFFSET];
static double mintime = FLT_MAX;
static double bytes   = 3 * sizeof(double) * N;

int main(int argc, char **argv)
{
  MPI_Init(&argc, &argv);
  const static double scalar = 3.0;
#pragma omp threadprivate(scalar)
  double times[NTIMES], rate;
  int    size;
  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 < N; j++) {
    a[j] = 1.0;
    b[j] = 2.0;
    c[j] = 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 (register int j = 0; j < N; j++) a[j] = b[j] + scalar * c[j];
    times[k] = MPI_Wtime() - times[k];
  }
  for (int k = 1; k < NTIMES; k++) { /* note -- skip first iteration */
    mintime = MIN(mintime, times[k]);
  }

  if (size == 65) printf("Never printed %g\n", a[11]);
  rate = 1.0E-06 * bytes / 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);
  }
  MPI_Finalize();
  return 0;
}