/* L-30 MCS 572 Friday 30 March 2006: the composite trapezoidal rule,
 * adapted from comptrap2.c using a reduction, save as comptrap3.c
 * and compile as "cc_r -qsmp=omp -o comptrap3 comptrap3.c -lm".
 * Observe how close the code resembles a sequential program! */

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

#define v 1      /* verbose flag */

double traprule ( double (*f) ( double x ), double a, double b, int n );
/* applies the composite Trapezoidal rule to approximate the integral
 * of f over [a,b] using n+1 function evaluations, use only for n > 0 */

double integrand ( double x ); /* the function we integrate */

int main ( int argc, char *argv[] )
{
   int i;
   int p = 8;
   int n = 1000;
   double my_pi = 0.0;
   double a,b,h,y,pi,error;

   omp_set_num_threads(p);

   h = 1.0/p;

   #pragma omp parallel for private(a,b) /* each thread has its own a,b */
      for(i=0; i<p; i++)
      {
         a = i*h;
         b = (i+1)*h;

         #pragma omp critical        /* critical section to protect my_pi */
            my_pi += traprule(integrand,a,b,n);
      }
   my_pi = 4.0*my_pi; pi = 2.0*asin(1.0); error = my_pi-pi;

   printf("Approximation for pi = %.15e with error = %.3e\n",my_pi,error);

   return 0;
}

double integrand ( double x )
{
   return sqrt(1.0 - x*x);
}

double traprule ( double (*f) ( double x ), double a, double b, int n )
{
   int i;
   double h = (b-a)/n; 
   double y = (f(a) + f(b))/2.0;
   double x;

   for(i=1,x=a+h; i < n; i++,x+=h) y += f(x);

   return h*y;
}