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

    Dynamic PI Calc example:
    
    This PI calculation example is similar to the simple PI Calc example, 
    where PI is approximated by numerical integration. However, in this example
    one of the new features of MPI-2 is used, namely Process Creation and 
    Management. In this example the user is prompted for the number of 
    processes to spawn and then for the number of intervals to use in the 
    computation. 

    When the new processes are spawned they receive the number of intervals 
    from the Master (rank 0 of the first group of processes). They then 
    compute their contribution and using an MPI::Reduce the result is computed 
    and send to the Master machine. The spawned processes exit, and the 
    procedure start all over.


    More specifically:
   --------------------------------
    - User is prompted for the number of processes to spawn
    - The Slaves are spawned
    - User is prompted for n (nIntervals) by the Master (rank 0)
    - The Master sends n to the Slaves.
    - Each Slave computes its slice and does a MPI::Reduce to compute and send
      the result to the Master
    - The Master collects the result and write the PI approximation to stdout.
    - The Slaves exits and the procedure restarts.
   --------------------------------

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

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

#include "Dynamic_PI_Calc.h"
#include <windows.h>
#include <stdio.h>

#include <mpi.h>

#define COMMAND     "Dynamic_PI_Calc_cpp.exe"
#define PATH_SIZE   512

/******************************************************************************
    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
 *****************************************************************************/
CDynamic_PI_Calc::CDynamic_PI_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
 *****************************************************************************/
CDynamic_PI_Calc::~CDynamic_PI_Calc()
{
}

/******************************************************************************
    main
 *****************************************************************************/
int main(int argc, char *argv[])
{
    int                  nIntervals;
    int                  nMaxProcs;
    int                  nScanfResult;
    int*                 pnArrayOfErrorCodes;
    double               dPI;
    double               dIntervalArea;
    MPI::Intercomm       mcInterComm;
    CDynamic_PI_Calc*    pcDynamicPICalc = new CDynamic_PI_Calc(argc, argv);

    /* Get the parent of this process (if any):                              */
    mcInterComm = MPI::Comm::Get_parent();

    /* If I do not have a parent that means that I am in the master group:   */
    if (mcInterComm == MPI::COMM_NULL) {
        if (pcDynamicPICalc->m_nCommRank == 0) {

            printf ("Calculation of PI by Numerical Integration\n");
            do {
                printf ("Input number of processes to spawn (0 = quit):\n");
                fflush(stdout);
                nScanfResult = scanf ("%ld", &nMaxProcs);
                getchar();

                if (nMaxProcs > 0 && nScanfResult > 0) {
                    pnArrayOfErrorCodes = (int*) malloc(nMaxProcs*sizeof(int));
                    /* Do the spawn:                                         */
                    mcInterComm = MPI::COMM_SELF.Spawn(COMMAND,
                                                       MPI::ARGV_NULL,
                                                       nMaxProcs,
                                                       MPI::INFO_NULL,
                                                       0,
                                                       pnArrayOfErrorCodes);

                    free(pnArrayOfErrorCodes);

                    printf ("Input of intervals:\n");
                    fflush(stdout);
                    scanf ("%ld", &nIntervals);
                    getchar();

                    /* Broadcast the number of intervals to the Slaves.      */
                    /* Notice that since we are using an inter-communicator  */
                    /* and WE are the root process (we have the data), the   */
                    /* constant MPI::ROOT is passed as opposed to a rank in a*/
                    /* intra-communicator broadcast. The same is true for    */
                    /* the MPI::Reduce, since we use the inter-communicator. */
                    mcInterComm.Bcast(&nIntervals, 1, MPI::INT, MPI::ROOT);

                    mcInterComm.Reduce(NULL,
                                       &dPI,
                                       1,
                                       MPI::DOUBLE,
                                       MPI::SUM,
                                       MPI::ROOT);

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

                    /* Disconnect the inter-communicator since we are done   */
                    /* using it:                                             */
                    mcInterComm.Disconnect();
                }
            } while (nMaxProcs > 0 && nScanfResult > 0);

        }
    } else {
        /* If I am a Slave just wait for the Master to tell me how many      */
        /* intervals we are going compute:                                   */
        mcInterComm.Bcast(&nIntervals,
                          1,
                          MPI::INT,
                          0);

        dIntervalArea = ComputeInterval (pcDynamicPICalc->m_nCommRank,
                                         pcDynamicPICalc->m_nCommSize,
                                         nIntervals);
        mcInterComm.Reduce(&dIntervalArea,
                           NULL,
                           1,
                           MPI::DOUBLE,
                           MPI::SUM,
                           0);

        /* Disconnect from the Master so we can call MPI::Finalize and exit  */
        /* without waiting for the Master to call MPI::Finalize:             */
        mcInterComm.Disconnect();
      }

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

    delete pcDynamicPICalc;
    return 0;
}