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

    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 II@criticalsoftware.com

 *****************************************************************************/
#include "PI_Calc.h"
#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);
}

/******************************************************************************
        Constructor
 *****************************************************************************/
CPI_Calc::CPI_Calc(int argc, char * argv[])
{
    /* Initialize csWMPI II II:                                                   */
    MPI::Init (argc, argv);

    /* Determine what the world looks like and our own position in it:       */
    m_nCommSize = MPI::COMM_WORLD.Get_size ();
    m_nCommRank = MPI::COMM_WORLD.Get_rank ();
    printf("I am rank %d of %d in MPI_COMM_WORLD\n", m_nCommRank, m_nCommSize);
    fflush(stdout);
}

/******************************************************************************
        Destructor
 *****************************************************************************/
CPI_Calc::~CPI_Calc()
{
}

/******************************************************************************
    main
 *****************************************************************************/
int main(int argc, char *argv[])
{
    int       nIntervals;
    double    dPI;
    double    dIntervalArea;
    int       nCounter;
    CPI_Calc* pcCPICalc = new CPI_Calc (argc, argv);

    /* 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 (pcCPICalc->m_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");
        fflush(stdout);
        for (nCounter = 1; nCounter < pcCPICalc->m_nCommSize; nCounter++)
            MPI::COMM_WORLD.Send (&nIntervals, 1, MPI::LONG, nCounter, 98);
    } else {
        /* If I am a Slave just wait for the Master to tell me how many      */
        /* intervals we are going compute:                                   */
        MPI::COMM_WORLD.Recv (&nIntervals,
                              1,
                              MPI::LONG,
                              0,
                              98);
    }

    /* Compute my contribution:                                              */
    dIntervalArea = ComputeInterval (pcCPICalc->m_nCommRank,
                                     pcCPICalc->m_nCommSize,
                                     nIntervals);

    /* If I am the Master, I should collect results, compute the sum, and    */
    /* display the result:                                                   */
    if (pcCPICalc->m_nCommRank == 0) {
        dPI = dIntervalArea;
        for (nCounter = 1; nCounter < pcCPICalc->m_nCommSize; nCounter++) {
            MPI::COMM_WORLD.Recv (&dIntervalArea,
                                  1,
                                  MPI::DOUBLE,
                                  nCounter,
                                  99);
            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::COMM_WORLD.Send (&dIntervalArea, 1, MPI::DOUBLE, 0, 99);
    }

    /* Finalize csWMPI II II:                                                     */
    MPI::Finalize ();

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

    delete pcCPICalc;
    return 0;
}