/* ///////////////////////////////////////////////////////////////////// */
/*! 
  \file  
  \brief Compute \f$\pi\f$ in parallel.
 
  Compute \f$ \pi \f$ by solving the integral
  \f[
    \int_0^1 \frac{4}{1+x^2} \, dx = \pi
  \f]
  with a midpoint rule with N sub-interval.
  The number of processes gives the stride between one interval and the
  next. So, for example, with 4 processes we have

  proc #0 will sum intervals i = 1, 5,  9, ...
  proc #1 will sum intervals i = 2, 6, 10, ...
  proc #2 will sum intervals i = 3, 7, 11, ...
  proc #3 will sum intervals i = 4, 8, 12, ...
  
  \author A. Mignone (mignone@to.infn.it)
  \date   March 10, 2020
*/
/* ///////////////////////////////////////////////////////////////////// */
#include <mpi.h>
#include <stdio.h>
#include <math.h>

int main(int argc, char *argv[])
{
  int rank, size;
  long int i, n;
  double tbeg, tend;
  double mypi, pi, h, sum, x;

/* --------------------------------------------------------
   0. Initialize the MPI execution environment
   -------------------------------------------------------- */

  MPI_Init(&argc,&argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  MPI_Comm_size(MPI_COMM_WORLD, &size);

/* --------------------------------------------------------
   1. Start computation.
   -------------------------------------------------------- */

  while (1) {    /* Continue until n = 0 is input */

  /* -- 1a. Read n from input & Broadcast (n = 0 to quit) -- */

    if (rank == 0) {
      printf("> Number of intervals: (0 quits) ");
      scanf("%ld",&n);
    }
    MPI_Bcast(&n, 1, MPI_LONG, 0, MPI_COMM_WORLD);
    if (n == 0) break;

  /* -- 1b. Begin local sum -- */

    tbeg = MPI_Wtime();
    h   = 1.0 / (double) n;  // 1.0/n will do as well 
    sum = 0.0;
    for (i = rank + 1; i <= n; i += size) {
      x    = h * ((double)i - 0.5);
      sum += 4.0 / (1.0 + x*x);
    }
    mypi = h*sum;
  
  /* -- 1c. Global reduction and print -- */

    MPI_Reduce(&mypi, &pi, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
  
    tend = MPI_Wtime();
    if (rank == 0){
      printf("> pi \approx  %12.6e, Error is %12.6e\n", pi, fabs(pi/M_PI-1.0));
      printf ("> Elapsed time = %f (s)\n", tend-tbeg);
    }
  }

  MPI_Finalize();
  return 0;
}