/* L-13 MCS 572 Wednesday 8 Feb 2006: the composite trapezoidal rule,
 * compile as "mpicc -o compute_pi compute_pi.c -lm",
 * and run as "mpirun -np 4 compute_pi" for 4 processors. */

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

#define v 0      /* 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,p;
   int n = 1000;
   double a,b,h,y,my_pi,pi,error,startwtime,endwtime;

   MPI_Init(&argc,&argv);
   MPI_Comm_size(MPI_COMM_WORLD,&p);
   MPI_Comm_rank(MPI_COMM_WORLD,&i);

   if(i==0) startwtime = MPI_Wtime();

   h = 1.0/p; a = i*h; b = (i+1)*h;

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

   if(v>0)
   {
      printf("Node %d computes %.15e as approximation.\n",i,y); fflush(stdout);
   }
   MPI_Reduce(&y,&my_pi,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD); 

   if(i==0)
   {
      endwtime = MPI_Wtime();
      if(v>0)
      {
         printf("Node 0 receives %.15e as sum.\n",my_pi); fflush(stdout);
      }
      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);
      printf("Total wall time : %lf seconds.\n",endwtime-startwtime);
   }

   MPI_Finalize();
   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;
}