program main

  use mpi

  double precision error
  double precision f
  double precision f1
  double precision f1_part
  double precision f2
  double precision f2_part
  integer i
  integer ierr
  integer j
  integer m
  integer n
  integer np
  integer p
  integer seed
  double precision stdev
  double precision sterr
  double precision value
  double precision x(2)

  call MPI_Init ( ierr )
  call MPI_Comm_size ( MPI_COMM_WORLD, p, ierr )
  call MPI_Comm_rank ( MPI_COMM_WORLD, id, ierr )

  m = 2
  if ( id == 0 ) then
    n = 1000000
    np = n / ( p - 1 )
    n = ( p - 1 ) * np
  end if
!
!  The Broadcast command says that the master (process 0) is sending 
!  an integer "NP" to all the processes, and that each process should
!  receive this value.
!
  call MPI_Bcast ( np, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr )
!
!  At this point, each process knows the value NP.
!  The workers should work, and the master can be idle.
!
  f1_part = 0.0
  f2_part = 0.0

  if ( 0 < id ) then
!
!  FORTRAN90 is very peculiar about how to control the random number seed.
!  Just call this function instead, for convenience.
!
    seed = 12345 + id
    call srand_f90 ( seed )

    do i = 1, np
      call random_number ( harvest = x(1:m) )
      value = f ( m, x )
      f1_part = f1_part + value
      f2_part = f2_part + value * value
    end do
  end if
!
!  The master zeros out F1 and F2 so they can be used to gather the results.
!
  if ( id == 0 ) then
    f1 = 0.0
    f2 = 0.0
  end if
!
!  The partial results are gathered into F1 and F2.
!
  call MPI_Reduce ( f1_part, f1, 1, MPI_DOUBLE_PRECISION, MPI_SUM, 0, MPI_COMM_WORLD, ierr )
  call MPI_Reduce ( f2_part, f2, 1, MPI_DOUBLE_PRECISION, MPI_SUM, 0, MPI_COMM_WORLD, ierr )
!
!  Print the results.
!
  if ( id == 0 ) then
    f1 = f1 / dble ( n )
    f2 = f2 / dble ( n - 1 )
    stdev = sqrt ( f2 - f1 * f1 )
    sterr = stdev / sqrt ( dble ( n ) )
    error = abs ( f1 - 1.0D+00 )
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '      N          F           stdev   sterr      err'
    write ( *, '(a)' ) ' '
    write ( *, '(i7,2x,g14.6,2x,f10.4,2x,f10.4,2x,g8.2)' ) n, f1, stdev, sterr, error
  end if
!
!  Terminate MPI.
!
  call MPI_Finalize ( ierr )

  stop
end
function f ( m, x )

  integer m

  double precision f
  integer i
  double precision value
  double precision x(m)

  value = 1.0
  do i = 1, m 
    value = value * abs ( 4.0 * x(i) - 2.0 )
  end do

  f = value

  return
end
subroutine srand_f90 ( seed )

  implicit none

  integer ( kind = 4 ) seed
  integer ( kind = 4 ), allocatable :: seed_vector(:)
  integer ( kind = 4 ) seed_size

  call random_seed ( size = seed_size )
  allocate ( seed_vector(seed_size) )
  seed_vector(1:seed_size) = seed
  call random_seed ( put = seed_vector(1:seed_size) )
  deallocate ( seed_vector )

  return
end