/* L-30 MCS 572 Friday 30 March 2006: the composite trapezoidal rule,
 * adapted from "compute_pi.c" in L-13, save this file as "comptrap1.c"
 * and compile as "cc_r -qsmp=omp -o comptrap1 comptrap1.c -lm". */

#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[8];
   double a,b,h,y,pi,error;

   omp_set_num_threads(p);

   h = 1.0/p;

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

      if(v>0)
      {
         printf("Thread %d integrates from %.2e to %.2e\n",i,a,b);
         fflush(stdout);
      }

      my_pi[i] = traprule(integrand,a,b,n);

      if(v>0)
      {
         printf("Thread %d computes %.15e as approximation.\n",i,my_pi[i]);
         fflush(stdout);
      }
   }
   for(i=1; i<p; i++) my_pi[0] += my_pi[i];

   my_pi[0] = 4.0*my_pi[0]; pi = 2.0*asin(1.0); error = my_pi[0]-pi;
   printf("Approximation for pi = %.15e with error = %.3e\n",my_pi[0],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;
}