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.
       
      integer dim_num
      parameter ( dim_num = 4 )

      double precision f
      double precision q
      double precision q_error
      double precision q_exact
      parameter ( q_exact = 1.0D+00 )
      integer sample
      integer sample_num
      parameter ( sample_num = 1000 )
      integer seed
      double precision x(dim_num)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'MC'
      write ( *, '(a)' ) '  FORTRAN77 version.'

      seed = 123456789

      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 )

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 
     &  '   Samples Dimension      Estimate          Error'
      write ( *, '(a)' ) ' '

      write ( *, '(2x,i8,2x,i8,2x,f16.10,2x,g16.6)' ) 
     &  sample_num, dim_num, q, q_error

      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