// calculate integral using client-server model
// compile:  mpicc integral2_mpi.c -lm
// usage: mpirun -n process# ./a.out trial#

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include <mpi.h>


#define CHUNK_SIZE 1000

// w is weight function
double w(double x, double y)
{
  return exp( - ( x*x + y*y ) );
}

//f is the function to be integrated
double f(double x, double y)
{

   return M_PI *  ( x*x + y*y ) ;

}


void generator(double * x, double * y)
{
   double xt, yt, ratio, tmp ;
   double delta = 0.5 ;


   tmp = (double)rand() / (double)RAND_MAX  ;
   xt = (*x) + delta * (2 * tmp - 1) ; 

   tmp = (double)rand() / (double)RAND_MAX  ;
   yt = (*y) + delta * (2 * tmp - 1) ; 

   ratio = w(xt, yt ) / w(*x, *y) ;
   tmp = (double)rand() / (double)RAND_MAX  ;

   if(ratio > tmp )
   {

      *x = xt ;
      *y = yt ;
   }


}


int main(int argc, char* argv[])
{
   double  x, y ;

   double sum ;
   double total_sum ;
   double integral ;
   int trials, trials_per_p;
   int i, j ;
   int pnumber, myrank;
   int x_array_tag = 1111 ;
   int y_array_tag = 2222 ;
   int demand_tag = 3333 ;
   int all_done_tag = 4444 ;
   int sum_tag = 5555 ;


   MPI_Init(&argc,&argv);
   MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
   MPI_Comm_size(MPI_COMM_WORLD, &pnumber);




   x = 0.0 ;
   y = 0.0 ;


   trials = atoi(argv[1]);


   if(myrank == (pnumber - 1 ))
   {

      int sent = 0 ;
      int demander ;
      int destination = 0;
      double * x_array ;
      double * y_array ;
      int flag=0;
      MPI_Status status;
      MPI_Request request;
      x_array = malloc( trials * sizeof(double)) ;
      y_array = malloc( trials * sizeof(double)) ;

      int generated = 0 ;

      for(i=0 ;  i < trials ; i++)
      {
         generator(&x, &y);
         x_array[generated] = x ;
         y_array[generated] = y ;
         generated++ ;
      }


      while(sent < trials)
      {
     
         MPI_Recv(&demander, 1, MPI_INT, MPI_ANY_SOURCE,  demand_tag, MPI_COMM_WORLD, &status) ;

            MPI_Send(&x_array[sent  ], CHUNK_SIZE , MPI_DOUBLE, demander, x_array_tag, MPI_COMM_WORLD);
            MPI_Send(&y_array[sent  ], CHUNK_SIZE , MPI_DOUBLE, demander, y_array_tag, MPI_COMM_WORLD);
            sent += CHUNK_SIZE ;

      }




      for (i=0 ; i < (pnumber-2) ; i++)
      {
         MPI_Recv(&demander, 1, MPI_INT, MPI_ANY_SOURCE,  demand_tag, MPI_COMM_WORLD, &status) ;

      }

      double all_done = 1.0 ;
     


      MPI_Send(&all_done, 1 , MPI_DOUBLE, 0, all_done_tag, MPI_COMM_WORLD);


   }


   if(myrank == 0)
   {
      int all_done = 0;
      double sum = 0.0 ;
      double temp;
      MPI_Status status;
      MPI_Request request;

      while(1)
      {
         MPI_Recv(&temp, 1, MPI_DOUBLE, MPI_ANY_SOURCE,  MPI_ANY_TAG , MPI_COMM_WORLD, &status) ;
         if(status.MPI_SOURCE == (pnumber - 1))
         {
            printf("all done received\n");
            break ;
         }
         sum += temp ;
      }

      printf("%f\n", sum / trials);
      MPI_Finalize();
   }



   
   if((myrank != 0) && (myrank != pnumber-1) )
   {

      MPI_Status status;
      MPI_Request request;

      int flag;
      double * local_x_array ;
      double * local_y_array ;

      local_x_array = malloc( CHUNK_SIZE * sizeof(double)) ;
      local_y_array = malloc( CHUNK_SIZE * sizeof(double)) ;

      while(1)
      {
         MPI_Send(&myrank, 1 , MPI_INT, (pnumber-1), demand_tag, MPI_COMM_WORLD);
         MPI_Recv(local_x_array, CHUNK_SIZE , MPI_DOUBLE, (pnumber - 1 ),  x_array_tag, MPI_COMM_WORLD, &status) ;
         MPI_Recv(local_y_array, CHUNK_SIZE , MPI_DOUBLE, (pnumber - 1 ),  y_array_tag, MPI_COMM_WORLD, &status) ;

     
            for(i=0 ; i < CHUNK_SIZE ; i++)
            {
               sum += f(local_x_array[i], local_y_array[i]) ;

            }
         MPI_Send(&sum, 1 , MPI_DOUBLE, 0, sum_tag, MPI_COMM_WORLD);
         sum = 0.0 ;

       }

         

   }


   return 0;

}