/* =====================================================
   MPI example: find min in array
   ===================================================== */


#include <iostream.h>
#include <stdlib.h>
#include <sys/time.h>
#include <math.h>
#include <mpi.h>

struct timeval start_time, stop_time;
int elapsed;
int num_procs;


void myAdd( void *a, void *b, int *len, MPI_Datatype *datatype)
{
   int i;

   if ( *datatype == MPI_INT )
   {
      int *x = (int *)a;    // Turn the (void *) into an (int *) 
      int *y = (int *)b;    // Turn the (void *) into an (int *)

      for (i = 0; i < *len; i++)
      {
         *y = *x + *y;
	 x++;
	 y++;
      }
   }
   else if ( *datatype == MPI_DOUBLE )
   {
      double *x = (double *)a;    // Turn the (void *) into an (double *)
      double *y = (double *)b;    // Turn the (void *) into an (double *)

      for (i = 0; i < *len; i++)
      {
         *y = *x + *y;
         x++;
         y++;
      }
   }
}


double f(double a)
{
   return( 2.0 / sqrt(1 - a*a) );
}

/* =======================
   MAIN
   ======================= */

int main(int argc, char *argv[])
{
   int    N;
   double w, x;
   int    i, myid;
   double mypi, final_pi;
   MPI_Op myOp;


   MPI_Init(&argc,&argv);                      // Initialize

   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);  // Get # processors
   MPI_Comm_rank(MPI_COMM_WORLD, &myid);       // Get my rank (id)


   /* -----
      Check command line
      ----- */
   if ( argc != 2 )
   {
      if ( myid == 0 )
         cout << "Usage: " << argv[0] << " Num_Intervals" << endl;

      MPI_Finalize();
      exit(1);
   }

   if ( myid == 0 )
      gettimeofday(&start_time, NULL);

   if ( myid == 0 )
      N = atoi(argv[1]);

   MPI_Bcast (&N, 1, MPI_INT, 0, MPI_COMM_WORLD);

   w = 1.0/(double) N;

   /* ******************************************************************* */

   mypi = 0.0;

   for (i = myid; i < N; i = i + num_procs)
   {
      x = w*(i + 0.5);         
      mypi = mypi +  w*f(x);  
   }

   /* ******************************************************************* */

   MPI_Op_create( myAdd, 1, &myOp);

   MPI_Reduce ( &mypi, &final_pi, 1, MPI_DOUBLE, myOp, 0, MPI_COMM_WORLD);

   if ( myid == 0 )
   {
      gettimeofday(&stop_time, NULL);

      elapsed = (stop_time.tv_sec*1000000 + stop_time.tv_usec) -
                (start_time.tv_sec*1000000 + start_time.tv_usec);

      cout << "Pi = " << final_pi << endl << endl;
      cout << "Elapsed time = " << elapsed << " microseconds" << endl;
   }

   MPI_Finalize();
}