c     /***************************************************************/
c     /* pi.f - compute pi by integrating f(x) = 4/(1 + x**2)        */
c     /* MPI VERSION                                                 */
c     /*                                                             */
c     /* Variables:                                                  */
c     /*                                                             */
c     /* pi   the calculated result                                  */
c     /* mypi the local sum                                          */
c     /* n    number of points of integration.                       */
c     /* x    midpoint of each rectangle's interval                  */
c     /* f    function to integrate                                  */
c     /* i    do loop index                                          */
c     /* numprocs total number of processors available               */
c     /* myid processor id (0<=myid<=numprocs-1)                     */
c     /***************************************************************/
      program mainpar
      include 'mpif.h'
      double precision  PI25DT
      parameter        (PI25DT = 3.141592653589793238462643d0)
      double precision  mypi,pi,h,sum,x,f,a,t0,t1,t2
      integer n,myid,numprocs,i,ierr
c     /* function to integrate */
      f(a) = 4.d0 / (1.d0 + a*a)
c     /* initialize MPI */
      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'
c     /* input of the number of intervals */
 10   if (myid.eq.0) then
         write(6,98)
 98      format('Enter the number of intervals: (0 quits)')
         read(5,99) n
 99      format(i10)
      endif  
c     /* broadcast the value for 'n' */
      call MPI_BCAST(n,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
c     /* check for quit signal */
      if ( n .le. 0 ) goto 30
c     /* start timing */
      t0=MPI_WTIME()
c     /* calculate the interval size */
      h = 1.0d0/n
c     /* perform the PARTIAL rectangular summation */
      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     /* collect all the partial sums */
      call MPI_REDUCE(mypi,pi,1,MPI_DOUBLE_PRECISION,MPI_SUM,0,
     $     MPI_COMM_WORLD,ierr)
c     /* stop timing & compute global maximum elapsed time*/
      t1=MPI_WTIME()
      call MPI_REDUCE(t1-t0,t2,1,MPI_DOUBLE_PRECISION,MPI_MAX,0,
     $     MPI_COMM_WORLD,ierr)
      if(myid.eq.0)write(*,*)'elapsed time: ',t2
c     /* node 0 prints result to screen */
      if (myid .eq. 0) then
         write(6, 97) pi, abs(pi - PI25DT)
 97      format('  pi is approximately: ', F18.16,
     +        '  Error is: ', E20.15)
      endif
c     
      goto 10
c     /* finalize MPI */
 30   call MPI_FINALIZE(ierr)
      stop
      end



