/*
 * COMP 633:  C/OpenMP example program
 *
 */
#include <stdio.h>
#include <math.h>
#include <sys/time.h>
#include <omp.h>
#include <unistd.h>
#define N_MAX 50000000 /* 50,000,000 */

double a[N_MAX],b[N_MAX];


/*
 * wall-clock timer in seconds for POSIX-compliant clocks
 */
double wctime()
{
  struct timeval tv;
  gettimeofday(&tv, NULL);
  return (tv.tv_sec + 1E-6 * tv.tv_usec);
}


main (int argc, char *argv[])
{
  double t1, t2, td, target_work;
  int i, t, num_threads, niter, n;

  /*
   * set iterations to keep parallel running time approximately constant
   * independent of problem size and number of threads
   */
  n = N_MAX;
  if (argc > 1) {
    n = atoi(argv[1]);
  }
  n = ( n > 0 && n <= N_MAX) ? n : N_MAX;
  num_threads = omp_get_max_threads();  
  target_work = 2.0 * N_MAX * num_threads;
  niter = (int) ceil( target_work / n );
  niter = 2 * ((niter + 1) / 2);

  /*
   * initialize arrays
   */
  /*  printf("System page size = %d bytes\n", getpagesize()); */

  printf("Init:");
  fflush(stdout);

#pragma omp parallel for shared(a,b) private(i) schedule(runtime)
  for (i = 0; i < n; i++){
    a[i] = 0.0;
    b[i] = 0.0;
  }

  a[0] = 1.0;
  b[0] = 1.0;

  /* 
   * time iterations
   */
  printf("\b\b\b\b\bTest:  n = %8d, p = %2d, niter = %8d, ...", 
	 n,num_threads,niter);
  fflush(stdout);
  t1 = wctime();

#pragma omp parallel shared(a,b) private(i,t) 
  for (t = 0; t < niter; t = t + 2){

#pragma omp for schedule(runtime)
      for (i = 1; i < n-1; i++)
	b[i] = (a[i-1] + a[i] + a[i+1]) / 3.0;

#pragma omp for schedule(runtime)
      for (i = 1; i < n-1; i++)
	a[i] = (b[i-1] + b[i] + b[i+1]) / 3.0;
  }

  t2 = wctime();
  td = t2 - t1;
  printf("\b\b\btime per elt = %6.1f ns\n", (td * 1E9) / (niter*n));
}