/*
* mpipi2.c
* Neil Gershenfeld  7/1/19
* use MPI to evaluate pi by summation, including overhead
*/

#include <stdio.h>
#include <mpi.h>
#include <sys/time.h>
#include <stdint.h>

#define NPTS 10000000000

void main(int argc, char** argv) {
   uint32_t rank,nproc;
   uint64_t i,istart,iend,init_time,start_time,end_time;
   struct timeval init,start,end;
   double sum,pi,mflops;

   gettimeofday(&init,NULL);
   MPI_Init(&argc,&argv);
   MPI_Comm_rank(MPI_COMM_WORLD,&rank);
   MPI_Comm_size(MPI_COMM_WORLD,&nproc);
   if (rank == 0) {
      MPI_Barrier(MPI_COMM_WORLD);
      gettimeofday(&start,NULL);
      istart = 1+rank*NPTS;
      iend = (rank+1)*NPTS;
      sum = 0.0;
      for (i = istart; i <= iend; ++i)
         sum += 0.5/((i-0.75)*(i-0.25));
      MPI_Reduce(&sum,&pi,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
      gettimeofday(&end,NULL);
      init_time = init.tv_sec*1e6+init.tv_usec;
      start_time = start.tv_sec*1e6+start.tv_usec;
      end_time = end.tv_sec*1e6+end.tv_usec;
      mflops = nproc*NPTS*(5.0/(end_time-start_time));
      printf("processes = %d, NPTS = %ld, pi = %f\n",nproc,NPTS,pi);
      printf("init = %f, time = %f, estimated MFlops = %f\n",(start_time-init_time)/1.0e6,(end_time-start_time)/1.0e6,mflops);
      }
   else {
      MPI_Barrier(MPI_COMM_WORLD);
      istart = 1+rank*NPTS;
      iend = (rank+1)*NPTS;
      sum = 0.0;
      for (i = istart; i <= iend; ++i)
         sum += 0.5/((i-0.75)*(i-0.25));
      MPI_Reduce(&sum,&pi,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
      }
   MPI_Finalize();
   }