!******************************************************************************
!    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 2002 (c) Critical Software SA
!    . http://www.criticalsoftware.com
!    . http://www.criticalsoftware.com/hpc
!    . csWMPI@criticalsoftware.com
!******************************************************************************

      include "mpi_init.f"

!******************************************************************************
!    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
      program main
      implicit none
      include "mpif.h"

      integer nCommRank, nCommSize, ierr, nIntervals, i, master
      real dPi, dIntervalArea

      master = 0

!     Initialize csWMPI II:
      call MPI_INIT(ierr)

!     Determine what the world looks like and our own position in it:
      call MPI_COMM_RANK(MPI_COMM_WORLD, nCommRank, ierr)
      call MPI_COMM_SIZE(MPI_COMM_WORLD, nCommSize, ierr)

      print *, 'I am rank',nCommRank,' of ',nCommSize,
     +         ' in MPI_COMM_WORLD'

!     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 .eq. master ) then
            print *,'Calculation of PI by Numerical Integration'
            print *,'Input of intervals: '
            read  *, nIntervals
            print *,'Master: Sending # of intervals to MPI-Processes '
            do i=1,nCommSize-1 !Send for all the number of intervals 
                  call MPI_SEND(nintervals, 1, MPI_INTEGER, 
     +                  i,0,MPI_COMM_WORLD,ierr)
            enddo
      else
!     If I am a Slave just wait for the Master to tell me how many      
!     intervals we are going compute:
            call MPI_RECV(nintervals, 1, MPI_INTEGER, master, 0,
     +            MPI_COMM_WORLD, MPI_STATUS_IGNORE,ierr)
      endif
!     Compute my contribution:
      call ComputeInterval(nCommRank,nCommSize,nintervals,
     +                     dIntervalArea)

!     If I am the Master, I should collect results, compute the sum, and    
!     display the result:
      if ( nCommRank .eq. master ) then !recolher os resultados
            dPi = dIntervalArea
            do i=1,nCommSize-1
                  call MPI_RECV(dIntervalArea, 1, MPI_REAL,
     +                  MPI_ANY_SOURCE,0, MPI_COMM_WORLD, 
     +                  MPI_STATUS_IGNORE, ierr)
                  dPi = dPi + dIntervalArea
            enddo

            print *,'Master: Has collected results from MPI processes '
            print *,'PI approximation:',dPi
      else
!     If I am a Slave send my result to the Master                      
            print *,'Slave: Sending my contribution to Master'
            call MPI_SEND(dIntervalArea, 1, MPI_REAL, master, 0,
     +                  MPI_COMM_WORLD,ierr)
      end if

      call MPI_FINALIZE(ierr)

!     Pause rank 0 so that the output can be verified:
      if (nCommRank .eq. master) then
            pause 'Press ENTER to exit...'
      end if
      stop
      end