c**********************************************************************
c   pi.f - compute pi by integrating f(x) = 4/(1 + x**2)     
c     
c   Each node: 
c    1) receives the number of rectangles used in the approximation.
c    2) calculates the areas of it's rectangles.
c    3) Synchronizes for a global summation.
c   Node 0 prints the result.
c
c  Variables:
c
c    pi  the calculated result
c    n   number of points of integration.  
c    x           midpoint of each rectangle's interval
c    f           function to integrate
c    sum,pi      area of rectangles
c    tmp         temporary scratch space for global summation
c    i           do loop index
c****************************************************************************
c     examples/int_coll

      include 'mpif.h'

      double precision  PI25DT
      parameter        (PI25DT = 3.141592653589793238462643d0)

      double precision  mypi,pi,h,sum,x,f,a
      integer n,my_rank,nproc,i,ierr
c                                 function to integrate
      f(a)=4.d0/(1.d0 + a*a)
      call init(nproc,my_rank)
      n=10
      do while(n.gt.0) 
         if(my_rank.eq.0) then
            write(*,*) 'Enter the number of intervals: (0 quits)'
            read(*,*) n
         endif
         call MPI_Bcast(n,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
         if(n.gt.0) then
            h=1.0d0/dble(n)
            sum=0.0d0
            do i=my_rank+1,n,nproc
               x=h*(dble(i)-0.5d0)
               sum=sum+f(x)
            enddo
            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            node 0 prints the answer.
            if (my_rank .eq. 0) then
               write(6, 100) pi, abs(pi - PI25DT)
 100           format('  pi is approximately: ', F18.16,
     *              '  Error is: ', F18.16)
            endif
         endif
      enddo
      call MPI_FINALIZE(ierr)
      stop
      end
c
      subroutine init(nproc,my_rank)
      include "mpif.h"
      integer my_rank,nproc
      character*30 nodename
      integer nchar
      call MPI_init(ierr)
      call MPI_Comm_rank(MPI_COMM_WORLD,my_rank,ierr)
      call MPI_Comm_size(MPI_COMM_WORLD,nproc,ierr)
      call MPI_Get_processor_name(nodename,nchar,ierr)
      write(*,*) 'my_rank=',my_rank,' processor=',nodename
      return
      end