c
c     compute pi by integrating f(x) = 4/(1 + x**2)     
c
c     each process:
c          - receives the # of intervals used in the approximation
c          - calculates the areas of it's rectangles
c          - synchronizes for a global summation
c     process 0 prints the result and the time it took
c
      program main

      include 'mpif.h'

      double precision PI25DT
      parameter (PI25DT = 3.141592653589793238462643d0)

      double precision  mypi, pi, h, sum, x, f, a
      integer n, myid, numprocs, i, ierr
c
c     function to integrate
c
      f(a) = 4.d0 / (1.d0 + a*a)

      call MPI_INIT(ierr)
      call MPI_COMM_RANK(MPI_COMM_WORLD, myid, ierr)
      call MPI_COMM_SIZE(MPI_COMM_WORLD, numprocs, ierr)
      print *, "Process ", myid, " of ", numprocs, " is alive"

      sizetype = 1
      sumtype = 2
      
      if (myid .eq. 0) then
           n = 100
      endif
      
      call MPI_BCAST(n, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
c
c     calculate the interval size
c
      h = 1.0d0 / n
      sum  = 0.0d0

      do 20 i = myid + 1, n, numprocs
           x = h * (dble(i) - 0.5d0)
           sum = sum + f(x)
 20   continue

      mypi = h * sum
c
c     collect all the partial sums
c
      call MPI_REDUCE(mypi, pi, 1, MPI_DOUBLE_PRECISION,
     +                MPI_SUM, 0, MPI_COMM_WORLD, ierr)
c
c     process 0 prints the result
c
      if (myid .eq. 0) then
           write(6, 97) pi, abs(pi - PI25DT)
 97        format('  pi is approximately: ', F18.16,
     +            '  Error is: ', F18.16)
      endif

      call MPI_FINALIZE(ierr)

      stop
      end
