program main

c*********************************************************************72
c
cc MAIN is the main program for MC.
c
c  Discussion:
c
c    Things you can vary include:
c
c    * the function to be integrated (see the list below)
c    * the number of samples used, SAMPLE_NUM.
c    * the spatial dimension D_NUM (try 1, 2, 4)
c
c    The function to be integrated is:
c
c      F = Product ( 1 <= DIM <= DIM_NUM ) ( 2 * abs ( 2 * X(DIM) - 1 ) ) 
c      Exact integral is 1.
c
c  You will need the "INCLUDE" statement here in order to access MPI functions
c  and variables.
c
      include 'mpif.h'

      integer dim_num
      parameter ( dim_num = 4 )

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

      if ( id .eq. master ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'MC:'
        write ( *, '(a)' ) '  FORTRAN77 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
c
c  Each process must use a different seed.
c
      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 / dble ( sample_num )
      q_error = abs ( q - q_exact )

      if ( id .eq. 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
c
c  Have each process sent results to process MASTER for reduction
c  to final result.
c
      call MPI_Reduce ( q, q_total, 1, MPI_DOUBLE_PRECISION, MPI_SUM, 
     &  master, MPI_COMM_WORLD, ierr )
c
c  "Clean up" the result.
c
      if ( id .eq. 0 ) then
        q_total = q_total / dble ( p )
        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 )

c*********************************************************************72
c
cc F evaluates the function F(X) which we are integrating.
c
c  Modified:
c
c    04 September 2008
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    Input, integer DIM_NUM, the number of spatial dimensions.
c
c    Input, double precision X[DIM_NUM], the point at which to evaluate F.
c
c    Output, double precision F, the value of F(X).
c
      integer dim_num

      integer dim
      double precision f
      double precision x(dim_num)
      double precision value

      value = 1.0D+00

      do dim = 1, dim_num

        value = value * 2.0D+00 * abs ( 2.0D+00 * x(dim) - 1.0D+00 )

      end do

      f = value

      return
      end
      subroutine r8vec_uniform_01 ( n, seed, r )

c*********************************************************************72
c
cc R8VEC_UNIFORM_01 returns a unit pseudorandom vector.
c
c  Parameters:
c
c    Input, integer N, the number of entries in the vector.
c
c    Input/output, integer SEED, the "seed" value, which 
c    should NOT be 0.  On output, SEED has been updated.
c
c    Output, double precision R(N), the vector of pseudorandom values.
c
      implicit none

      integer n

      integer i
      integer i4_huge
      parameter ( i4_huge = 2147483647 )
      integer k
      integer seed
      double precision r(n)

      if ( seed .eq. 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'R8VEC_UNIFORM_01 - Fatal errorc'
        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 .lt. 0 ) then
          seed = seed + i4_huge
        end if

        r(i) = dble ( seed ) * 4.656612875D-10

      end do

      return
      end