!******************************************************************************
!    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 'mpif.f90'

!******************************************************************************
!    This function computes an interval for a single process.
!******************************************************************************
      subroutine ComputeInterval(n_my_rank,n_comm_size,n_intervals,&
                                   r_dPi)
      implicit none
      integer n_my_rank, n_comm_size, n_intervals
      real r_dPi,dWidth, dX, dLocalSum
      integer  nCounter

      dWidth = 1.0 / n_intervals
      dLocalSum = 0.0

      do nCounter=1,n_intervals,n_comm_size
            dX = (nCounter + 0.5) * dWidth
          dLocalSum = dLocalSum + (4 / (1 + (dX * dX)))
      enddo

      r_dPi = (dLocalSum * dWidth)

      end !subroutine ComputeInterval


!******************************************************************************
!    main:
!******************************************************************************
      program main
      use MPI
      implicit none
      integer nCommRank, nCommSize, nIntervals, ierr
      integer nMaxProcs, nScanfResult, mcInterComm
      integer, allocatable ::pnArrayOfErrorCodes(:)
      real dPI, dIntervalArea
      character*(*) COMMAND
      parameter (COMMAND = 'Dynamic_PI_Calc')

!     Initialize csWMPI II II:
      call MPI_INIT(ierr)

!     Determine what the world looks like and our own position in it:
      call MPI_COMM_SIZE(MPI_COMM_WORLD, nCommSize, ierr)
      call MPI_COMM_RANK(MPI_COMM_WORLD, nCommRank, ierr)
      print *,'I am rank',nCommRank,' of ',nCommSize,&
               ' in MPI_COMM_WORLD'

!     Get the parent of this process (if any):
      call MPI_COMM_GET_PARENT(mcInterComm, ierr)

!     If I do not have a parent that means that I am in the master group:
      if (mcInterComm .eq. MPI_COMM_NULL) then
            if (nCommRank .eq. 0) then

                  print *,'Calculation of PI by Numerical Integration'
                  nMaxProcs = 1
                  !cycle
10                      print *,'Input number of processes ', &
                                 'to spawn (0 = quit):'
                        read *, nMaxProcs

                        if (nMaxProcs .gt. 0) then
                            allocate(pnArrayOfErrorCodes(nMaxProcs))
!                             Do the spawn:
                              call MPI_COMM_SPAWN(COMMAND,MPI_ARGV_NULL,&
                                nMaxProcs,MPI_INFO_NULL,0,MPI_COMM_SELF,&
                                mcInterComm,pnArrayOfErrorCodes,ierr)
                              deallocate(pnArrayOfErrorCodes)

                              print *,'Input of intervals: '
                              read *, nIntervals

!                    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.
                              call MPI_Bcast(nIntervals, 1, MPI_INTEGER,&
                                   MPI_ROOT, mcInterComm, ierr)

                              call MPI_REDUCE(0,dPI,1,MPI_REAL,&
                                   MPI_SUM,MPI_ROOT,mcInterComm, ierr)

                              print *,'Master: Collected results from',&
                                    'MPI processes'
                              print *,'PI approximation: ',dPI

!                    Disconnect the inter-communicator since we are done
!                    using it:
                              call MPI_COMM_DISCONNECT(mcInterComm,ierr)
                              goto 10
                      else
                              goto 20
                        end if
!                 end cycle
            end if
      else
!           If I am a Slave just wait for the Master to tell me how many
!           intervals we are going compute:
            call MPI_BCAST(nIntervals,1,MPI_INTEGER,0, mcInterComm,&
                           ierr)

            call ComputeInterval(nCommRank,nCommSize,nIntervals,&
                                 dIntervalArea)

            call MPI_REDUCE(dIntervalArea, 0, 1, MPI_REAL,&
                            MPI_SUM, 0, mcInterComm, ierr)

!           Disconnect from the Master so we can call MPI_Finalize and exit
!           without waiting for the Master to call MPI_Finalize:
            call MPI_COMM_DISCONNECT(mcInterComm, ierr)
      end if

20    call MPI_FINALIZE(ierr)

      stop
      end