program main

!*****************************************************************************80
!
!! MAIN is the main program for MC.
!
!  Discussion:
!
!    Things you can vary include:
!
!    * the function to be integrated (see the list below)
!    * the number of samples used, SAMPLE_NUM.
!    * the spatial dimension D_NUM (try 1, 2, 4)
!
!    The function to be integrated is:
!
!      F = Product ( 1 <= DIM <= DIM_NUM ) ( 2 * abs ( 2 * X(DIM) - 1 ) ) 
!      Exact integral is 1.
!
  use mpi

  integer ( kind = 4 ), parameter :: dim_num = 4

  real    ( kind = 8 ) f
  integer ( kind = 4 ) id
  integer ( kind = 4 ) ierr
  integer ( kind = 4 ), parameter :: master = 0
  integer ( kind = 4 ) p
  real    ( kind = 8 ) q
  real    ( kind = 8 ) q_error
  real    ( kind = 8 ), parameter :: q_exact = 1.0D+00
  real    ( kind = 8 ) q_total
  integer ( kind = 4 ) sample
  integer ( kind = 4 ), parameter :: sample_num = 1000
  integer ( kind = 4 ) sample_total
  integer ( kind = 4 ) seed
  real    ( kind = 8 ) wtime
  real    ( kind = 8 ) wtime1
  real    ( kind = 8 ) wtime2
  real    ( kind = 8 ) x(dim_num)
!
!  Establish the MPI environment.
!
  call MPI_Init ( ierr )
!
!  Get this process's ID.
!
  call MPI_Comm_rank ( MPI_COMM_WORLD, id, ierr )
!
!  Find out how many processes are available.
!
  call MPI_Comm_size ( MPI_COMM_WORLD, p, ierr )

  if ( id == master ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'MC:'
    write ( *, '(a)' ) '  FORTRAN90 version'
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  An MPI program to estimate an integral.'
    write ( *, '(a,i8)' ) '  The number of processes is ', p
    wtime1 = MPI_Wtime ( )
  end if
!
!  Each process must use a different seed.
!
  seed = 123456789 + id

  q = 0.0D+00

  do sample = 1, sample_num

    call r8vec_uniform_01 ( dim_num, seed, x )

    q = q + f ( dim_num, x )

  end do

  q = q / real ( sample_num, kind = 8 )
  q_error = abs ( q - q_exact )

  if ( id == master ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '     Process  Samples Dimension      Estimate          Error'
    write ( *, '(a)' ) ' '
  end if

  write ( *, '(2x,i8,2x,i8,2x,i8,2x,f16.10,2x,g16.6)' ) &
    id, sample_num, dim_num, q, q_error
!
!  Have each process sent results to process MASTER for reduction
!  to final result.
!
  call MPI_Reduce ( q, q_total, 1, MPI_DOUBLE_PRECISION, MPI_SUM, &
    master, MPI_COMM_WORLD, ierr )
!
!  "Clean up" the result.
!
  if ( id == 0 ) then

    q_total = q_total / real ( p, kind = 8 )
    q_error = abs ( q_total - q_exact )
    sample_total = p * sample_num
    write ( *, '(a)' ) ' '
    write ( *, '(2x,a8,2x,i8,2x,i8,2x,f16.10,2x,g16.6)' ) &
      '   Total', sample_total, dim_num, q_total, q_error

    wtime2 = MPI_Wtime ( )
    wtime = wtime2 - wtime1

    write ( *, '(a)' ) ' '
    write ( *, '(a,f14.6)' ) '  Elapsed wall clock seconds = ', wtime

  end if

  call MPI_Finalize ( ierr )

  stop
end
function f ( dim_num, x )

!*****************************************************************************80
!
!! F evaluates the function F(X) which we are integrating.
!
!  Modified:
!
!    04 September 2008
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the number of spatial dimensions.
!
!    Input, real ( kind = 8 ) X(DIM_NUM), the point at which to evaluate F.
!
!    Output, real ( kind = 8 ) F, the value of F(X).
!
  integer ( kind = 4 ) dim_num

  real ( kind = 8 ) f
  real ( kind = 8 ) x(dim_num)

  f = product ( 2.0D+00 * abs ( 2.0D+00 * x(1:dim_num) - 1.0D+00 ) )

  return
end
subroutine r8vec_uniform_01 ( n, seed, r )

!*****************************************************************************80
!
!! R8VEC_UNIFORM_01 returns a unit pseudorandom vector.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the number of entries in the vector.
!
!    Input/output, integer ( kind = 4 ) SEED, the "seed" value, which 
!    should NOT be 0.  On output, SEED has been updated.
!
!    Output, real ( kind = 8 ) R(N), the vector of pseudorandom values.
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) i
  integer ( kind = 4 ), parameter :: i4_huge = 2147483647
  integer ( kind = 4 ) k
  integer ( kind = 4 ) seed
  real    ( kind = 8 ) r(n)

  if ( seed == 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'R8VEC_UNIFORM_01 - Fatal error!'
    write ( *, '(a)' ) '  Input value of SEED = 0.'
    stop
  end if

  do i = 1, n

    k = seed / 127773

    seed = 16807 * ( seed - k * 127773 ) - k * 2836

    if ( seed < 0 ) then
      seed = seed + i4_huge
    end if

    r(i) = real ( seed, kind = 8 ) * 4.656612875D-10

  end do

  return
end