NCSA Home
Contact Us | Intranet | Search

ncsa

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 PIX
      parameter (PIX = 4*atan(1.0))

      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"

      if (myid .eq. 0) then
         print *,"Enter the number of intervals: (0 to quit)"
         read(5,*) n
      endif

      call MPI_BCAST(n, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)

c check for n > 0
      IF (N.GT.0) THEN
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 - PIX)
 97          format('  pi is approximately: ', F18.16,
     +              '  Error is: ', F18.16)
        endif

      ENDIF

      call MPI_FINALIZE(ierr)

      stop
      end