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

    PI Calc example:
    
    In this example PI approximated using numerical integration. By dividing 
    the interval [0..1] on the x-axis into n slices and calculating the area of
    the unit circle (radius = 1) in the first quadrant, the approximation is 
    computed. The computation is distributed by letting each of the m MPI 
    processes compute n/m slices.
    
    More specifically:
   --------------------------------
    - User is prompted for n (nIntervals) by the Master (rank 0)
    - The Master sends n to the Slaves.
    - Each process computes its slice and sends the result to the Master
    - The Master collects the result and write the PI approximation to stdout.
   --------------------------------

    Copyright 2003 (c) Critical Software SA
    . http://www.criticalsoftware.com
    . http://www.criticalsoftware.com/hpc
    . csWMPI@criticalsoftware.com

 *****************************************************************************/
#include <stdio.h>
#include <mpi.h>


/******************************************************************************
    This function computes an interval for a single process.
 *****************************************************************************/
double ComputeInterval (int n_my_rank, int n_comm_size, int n_intervals)
{
    double dWidth;
    double dX;
    double dLocalSum;
    int nCounter;

    dWidth = 1.0 / n_intervals;
    dLocalSum = 0.0;
    for (nCounter = n_my_rank; nCounter < n_intervals; nCounter += n_comm_size)
    {
      dX = (nCounter + 0.5) * dWidth;
      dLocalSum += 4 / (1 + dX * dX);
    }
    return (dLocalSum * dWidth);
}

/******************************************************************************
    main:
 *****************************************************************************/
int main(int argc, char ** argv)
{
    int         nCommRank;
    int         nCommSize;
    int         nIntervals;
    double      dPI;
    double      dIntervalArea;
    int         nCounter;

    /* Initialize csWMPI II:                                                   */
    MPI_Init(&argc, &argv);

    /* Determine what the world looks like and our own position in it:       */
    MPI_Comm_size(MPI_COMM_WORLD, &nCommSize);
    MPI_Comm_rank(MPI_COMM_WORLD, &nCommRank);
    printf ("I am rank %d of %d in MPI_COMM_WORLD\n", nCommRank, nCommSize);
    fflush(stdout);


    /* If I am the Master (rank 0) then prompt the user for the number of    */
    /* intervals to compute and send it to the Slaves.                       */
    if (nCommRank == 0) {
        printf ("Calculation of PI by Numerical Integration\n");
        printf ("Input of intervals:\n");
        fflush(stdout);
        scanf ("%ld", &nIntervals);
        getchar();

        printf ("Master: Sending # of intervals to MPI-Processes \n");
        for (nCounter = 1; nCounter < nCommSize; nCounter++)
            MPI_Send (&nIntervals, 1, MPI_LONG, nCounter, 98, MPI_COMM_WORLD);
    } else {
        /* If I am a Slave just wait for the Master to tell me how many      */
        /* intervals we are going compute:                                   */
        MPI_Recv (&nIntervals,
                  1,
                  MPI_LONG,
                  0,
                  98,
                  MPI_COMM_WORLD,
                  MPI_STATUS_IGNORE);
    }

    /* Compute my contribution:                                              */
    dIntervalArea = ComputeInterval (nCommRank, nCommSize, nIntervals);

    /* If I am the Master, I should collect results, compute the sum, and    */
    /* display the result:                                                   */
    if (nCommRank == 0) {
        dPI = dIntervalArea;
        for (nCounter = 1; nCounter < nCommSize; nCounter++) {
            MPI_Recv (&dIntervalArea,
                      1,
                      MPI_DOUBLE,
                      nCounter,
                      99,
                      MPI_COMM_WORLD,
                      MPI_STATUS_IGNORE);

            dPI += dIntervalArea;
        }

        printf ("Master: Has collected results from MPI processes \n");
        printf ("\nPI approximation: %.20lf\n", dPI);
        fflush(stdout);

    } else {
        /* If I am a Slave send my result to the Master                      */
        printf ("Slave: Sending my contribution to Master\n");
        fflush(stdout);
        MPI_Send (&dIntervalArea, 1, MPI_DOUBLE, 0, 99, MPI_COMM_WORLD);
    }

    MPI_Finalize ();

    /* Pause rank 0 so that the output can be verified:                      */
    if (nCommRank == 0) {
        printf("\nPress ENTER to exit...\n");
      fflush(stdout);
    }
    return 0;
}