/* File:    omp_simpson.c
 * Purpose: Estimate area using Simpson's rule.
 *
 * Input:   a, b, n
 * Output:  estimate of area between x = a, x = b, x-axis and graph of f(x)
 *          using n subdivisions with Simpson's rule.
 *
 * Compile: gcc -g -Wall -fopenmp -o omp_simpson omp_simpson.c
 * Usage:   ./omp_simpson <number of threads>
 *
 * Notes:   
 *   1.  The function f(x) is hardwired.
 *   2.  In this version, each thread explicitly estimates the area
 *       over its assigned subinterval, a critical directive is used
 *       for the global sum.
 *   3.  This version assumes that n is evenly divisble by 
 *       2*thread_count
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <omp.h>
#include "timer.h"


void Usage(char* prog_name);
double f(double x);    /* Function we're integrating */
void Simpson(double a, double b, int n, double* area_p);

int main(int argc, char* argv[]) {
   double  area;       /* Store result in area          */
   double  a, b;       /* Left and right endpoints      */
   int     n;          /* Total number of trapezoids    */
   int     thread_count;
   double  start, finish;

   if (argc != 2) Usage(argv[0]);
   thread_count = strtol(argv[1], NULL, 10);

   printf("Enter a, b, and n\n");
   scanf("%lf", &a);
   scanf("%lf", &b);
   scanf("%d", &n);
   
   GET_TIME(start);
   area = 0.0;
#  pragma omp parallel num_threads(thread_count)
   Simpson(a, b, n, &area);
   GET_TIME(finish);

   printf("With n = %d our estimate\n", n);
   printf("of the area from %f to %f = %.15f\n",
      a, b, area);
   printf("Elapsed time = %e seconds\n", finish-start);

   return 0;
}  /* main */

/*--------------------------------------------------------------------
 * Function:    Usage
 * Purpose:     Print command line for function and terminate
 * In arg:      prog_name
 */
void Usage(char* prog_name) {

   fprintf(stderr, "usage: %s <number of threads>\n", prog_name);
   exit(0);
}  /* Usage */


/*------------------------------------------------------------------
 * Function:    f
 * Purpose:     Compute value of function to be integrated
 * Input arg:   x
 * Return val:  f(x)
 */
double f(double x) {
   double return_val;

   return_val = x*x;
// return_val = sin(x);
   return return_val;
}  /* f */

/*------------------------------------------------------------------
 * Function:    Simpson
 * Purpose:     Use Simpson's rule to estimate area
 * Input args:  
 *    a: left endpoint
 *    b: right endpoint
 *    n: number of trapezoids
 * 
 * Output arg:
 *    area_p:  estimate of area
 */
void Simpson(double a, double b, int n, double* area_p) {
   double  h, x, my_area, temp;
   double  local_a, local_b;
   int  i, local_n;
   int my_rank = omp_get_thread_num();
   int thread_count = omp_get_num_threads();

   h = (b-a)/n;
   local_n = n/thread_count;
   local_a = a + my_rank*local_n*h;
   local_b = local_a + local_n*h;

   my_area = f(local_a) + f(local_b);

   temp = 0.0;
   for (i = 1; i <= local_n-1; i += 2) {
      x = local_a + i*h;
      temp += f(x);
   }
   my_area += 4*temp;
   
   temp = 0.0;
   for (i = 2; i <= local_n-2; i += 2) {
      x = local_a + i*h;
      temp += f(x);
   }
   my_area += 2*temp;

   my_area *= h;
   my_area /= 3;
   
#  pragma omp critical
   *area_p += my_area;
}  /* Simpson */