include "mpiSim.f90"
program main
include "useSim.f90"

!**********************************************************************
!  http://www.csit.fsu.edu/~burkardt/f_src/mpi/mpi.html
!
!  MONTE_CARLO computes PI using Monte Carlo techniques.
!   The MONTE_CARLO program computes PI by the Monte Carlo method,
!  testing whether random points in the unit square are in the unit
!  circle.
!  Discussion:
!
!    Generate N random points in the unit square.  Count M, the number of
!    points that are in the quarter circle.  Then PI is equal to the
!    ratio 4 * M / N.
!
!    The work will be divided as follows:
!
!      Process 0 is in charge.
!
!      Processes 0 through NUM_PROCS - 2 take sets of random points and count
!      the number of points in the quarter circle.
!
!      Process NUM_PROCS - 1 computes the sets of random points.
!
!    The communicator MPI_WORLD_COMM comprises all processes.
!
!    The communicator WORKER_COMM is created to comprise the processes
!    whose world ID's are 1 through NUM_PROCS-2, that is, all but the
!    random number generating process.
!
!    Message tag 0 is a signal from the master process to exit the loop.
!    Message tag 1 is a request for a set of random numbers.
!    Message tag 2 is a reply from the server with a set of random numbers.
!
!  Modified:
!
!    15 September 2002
!
!  Reference:
!
!    Gropp, Lusk, Skjellum,
!    Using MPI,
!    Portable Parallel Programming with the Message-Passing Interface,
!    MIT Press, 1997, pages 47-53.
!

!
!  Fortran77 include file:
!
! include 'mpif.h'
!
!  Fortran90 module:
!
!  use mpi
!
! implicit none
!
  integer, parameter :: chunk_size = 1000

  integer done
  integer :: dummy = 0
  double precision error
  integer i
  integer ierr
  integer in_local
  integer in_total
  integer, parameter :: master = 0
  integer num_procs
  integer out_local
  integer out_total
  double precision pi_est
  double precision, parameter :: pi_true = 3.141592653589793238462643D+00
  integer point_local

  integer :: point_max = 1000000
  integer point_total
  double precision rands(chunk_size)
  integer ranks(1)
  integer requestor
  integer server
  integer stat(MPI_STATUS_SIZE)
  integer tag
  integer, parameter :: tag_exit = 0
  integer, parameter :: tag_random_send = 1
  integer, parameter :: tag_random_sent = 2
  double precision tolerance
  integer worker_comm
  integer worker_group
  integer worker_id
  integer world_group
  integer world_id
  double precision x
  double precision y
!
!  Initialize MPI.
!
  call MPI_Init ( ierr )

  if ( ierr /= MPI_SUCCESS ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'MONTE_CARLO: Warning!'
    write ( *, '(a,i6)' ) '  MPI_INIT returns IERR = ', ierr
    call MPI_Finalize ( ierr )
    stop
  end if
!
!  Get the number of processors.
!
  call MPI_Comm_size ( MPI_COMM_WORLD, num_procs, ierr )
!
!  Get the rank of this processor.
!
  call MPI_Comm_rank ( MPI_COMM_WORLD, world_id, ierr )
!
!  The master process gets the value of the tolerance...
!
  if ( world_id == master ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'MONTE_CARLO - Master process:'
    write ( *, '(a)' ) '  A FORTRAN90 program, using MPI to'
    write ( *, '(a)' ) '  estimate PI by the Monte Carlo method.'
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) '  The number of processes is ', num_procs
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  Points in the unit square will be tested'
    write ( *, '(a)' ) '  to see if they lie in the unit quarter circle.'
  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a,i3,a)' ) '  Process ', world_id, ' is active.'

  if ( world_id == master ) then
    tolerance = 0.001D+00
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  The method continues computing until:'
    write ( *, '(a,g14.6)' ) '  PI is computed to a tolerance of ', tolerance
    write ( *, '(a,i9)' ) '  or the number of points used reaches ', point_max
  end if
!
!  ...and broadcasts it to all other processes.
!
  call MPI_Bcast ( tolerance, 1, MPI_DOUBLE_PRECISION, master, MPI_COMM_WORLD, &
    ierr )
!
!  Now create the WORKER communication group, which excludes process the
!  SERVER process.
!
!
!  Start by getting the group corresponding to the world communicator.
!
  call MPI_Comm_group ( MPI_COMM_WORLD, world_group, ierr )
!
!  Put SERVER on the list of processes to exclude, and create the new
!  worker group.
!
  server = num_procs - 1

  if ( world_id == master ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'MONTE_CARLO - Master process:'
    write ( *, '(a,i6)' ) '  The random number server process ID is ', server
  end if

  ranks(1) = server

  call MPI_Group_excl ( world_group, 1, ranks, worker_group, ierr )
!
!  Use the worker group to create the new worker communicator.
!
  call MPI_Comm_create ( MPI_COMM_WORLD, worker_group, worker_comm, ierr )
!
!  Since we only needed the worker group to create the worker communicator,
!  we can free the worker group now.
!
  call MPI_Group_free ( worker_group, ierr )
!
!  Here is where the computation is carried out.
!
!  The SERVER process waits to receive a request from any other process.
!
  if ( world_id == server ) then

    do

      tag = mpi_any_tag

      call MPI_Recv ( dummy, 1, MPI_INTEGER, MPI_ANY_SOURCE, tag, &
        MPI_COMM_WORLD, stat, ierr )

      tag = stat(mpi_tag)

      if ( tag == tag_exit ) then
        exit
      end if

      requestor = stat(mpi_source)

      call random_number ( harvest = rands(1:chunk_size) )

      tag = tag_random_sent

      call MPI_Send ( rands, chunk_size, MPI_DOUBLE_PRECISION, requestor, &
        tag, MPI_COMM_WORLD, ierr )

    end do
!
!  Each worker process sends requests for numbers to the random number server.
!
  else

    in_local = 0
    out_local = 0
    point_local = 0

    do

      tag = tag_random_send

      call MPI_Send ( dummy, 1, MPI_INTEGER, server, tag, MPI_COMM_WORLD, ierr )

      tag = mpi_any_tag

      call MPI_Recv ( rands, chunk_size, MPI_DOUBLE_PRECISION, server, &
        tag, MPI_COMM_WORLD, stat, ierr )

      do i = 1, chunk_size, 2

        x = rands(i)
        y = rands(i+1)

        point_local = point_local + 1

        if ( x**2 + y**2 <= 1.0D+00 ) then
          in_local = in_local + 1
        else
          out_local = out_local + 1
        end if

      end do

      call MPI_Reduce ( in_local, in_total, 1, MPI_INTEGER, mpi_sum, &
        master, worker_comm, ierr )

      call MPI_Reduce ( out_local, out_total, 1, MPI_INTEGER, mpi_sum, &
        master, worker_comm, ierr )

      call MPI_Reduce ( point_local, point_total, 1, MPI_INTEGER, mpi_sum, &
        master, worker_comm, ierr )
!
!  The Master process now checks the value of PI, and the size of POINT_TOTAL.
!
      if ( world_id == master ) then

        pi_est = 4.0D+00 * in_total / point_total

        error = abs ( pi_est - pi_true )

        done = 1
!
!  If it's time to stop, the Master process informs the random number server.
!
        if ( error < tolerance .or. point_max <= point_total ) then

          done = 0
          tag = tag_exit

          call MPI_Send ( dummy, 1, MPI_INTEGER, server, tag, MPI_COMM_WORLD, &
            ierr )

        end if

      end if
!
!  The Master process broadcasts to all workers whether to quit or go on.
!
      call MPI_Bcast ( done, 1, MPI_INTEGER, master, worker_comm, ierr )

      if ( done == 0 ) then
        exit
      end if

    end do

  end if

  if ( world_id == master ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'MONTE_CARLO - Master process:'
    write ( *, '(a,g14.6)' ) '  Estimate for PI = ', pi_est
    write ( *, '(a,g14.6)' ) '  Error =           ', error
! ERROR: the variable, point_total, was never defined.
    write ( *, '(a,i6)' ) '  Number of points = ', point_total
  end if
!
!  We free the worker communicator.
!
  call MPI_Comm_free ( worker_comm, ierr )
!
!  Shut down MPI.
!
  if ( world_id == master ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'MONTE_CARLO - Master process:'
    write ( *, '(a)' ) '  Normal end of execution.'
  end if

  call MPI_Finalize ( ierr )

  stop
end
municator.
!
  call MPI_Comm_free ( worker_comm, ierr )
!
!  Shut down MPI.
!
  if ( world_id == master ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'MONTE_CARLO - Master process:'
    write ( *, '(a)' ) '  Normal end of execution.'
  end if

  call MPI_Finalize ( ierr )

  stop
end