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

    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 <stdio.h>
#include <malloc.h>
#include <mpi.h>
#include <windows.h>


#define COMMAND         "Dynamic_PI_Calc.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);
}



/******************************************************************************
    main:
 *****************************************************************************/
int main(int argc, char ** argv)
{
    int         nCommRank;
    int         nCommSize;
    int         nIntervals;
    double      dPI;
    double      dIntervalArea;
    int         nMaxProcs;
    int*        pnArrayOfErrorCodes;
    int         nScanfResult;
    MPI_Comm    mcInterComm;


    /* Initialize csWMPI II 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);

    /* Get the parent of this process (if any):                              */
    MPI_Comm_get_parent(&mcInterComm);

    /* If I do not have a parent that means that I am in the master group:   */
    if (mcInterComm == MPI_COMM_NULL) {
        if (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:                                         */
                    MPI_Comm_spawn(COMMAND,
                                   MPI_ARGV_NULL,
                                   nMaxProcs,
                                   MPI_INFO_NULL,
                                   0,
                                   MPI_COMM_SELF,
                                   &mcInterComm,
                                   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.  */
                    MPI_Bcast(&nIntervals, 1, MPI_INT, MPI_ROOT, mcInterComm);

                    MPI_Reduce(NULL,
                               &dPI,
                               1,
                               MPI_DOUBLE,
                               MPI_SUM,
                               MPI_ROOT,
                               mcInterComm);

                    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:                                             */
                    MPI_Comm_disconnect(&mcInterComm);
                }
            } 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:                                   */
        MPI_Bcast(&nIntervals,
                  1,
                  MPI_INT,
                  0,
                  mcInterComm);

        dIntervalArea = ComputeInterval (nCommRank, nCommSize, nIntervals);
        MPI_Reduce(&dIntervalArea,
                   NULL,
                   1,
                   MPI_DOUBLE,
                   MPI_SUM,
                   0,
                   mcInterComm);

        /* Disconnect from the Master so we can call MPI_Finalize and exit   */
        /* without waiting for the Master to call MPI_Finalize:              */
        MPI_Comm_disconnect(&mcInterComm);
    }

    MPI_Finalize ();

    return 0;
}