!****************************************************************************
!   pi_mpi.f - compute pi by integrating f(x) = 4/(1 + x**2) from 0 to 1
!   This is an MPI g77 compatible program
!****************************************************************************
      PROGRAM pi_mpi

        IMPLICIT NONE

        INCLUDE 'mpif.h'

        DOUBLE PRECISION     PI25DT, mypi, pi, h, sum, x
        REAL                 t(2), t_start, t_end
        INTEGER              n, myid, numprocs, ierr, i, rc
        PARAMETER(PI25DT = 3.141592653589793238462643D0)

        CALL MPI_INIT( ierr )
        CALL MPI_COMM_RANK( MPI_COMM_WORLD, myid, ierr )
        CALL MPI_COMM_SIZE( MPI_COMM_WORLD, numprocs, ierr )
        IF ( myid==0 ) THEN
          WRITE(*,*) "MPI initialized with ",numprocs, " processors"
        END IF

        DO

           IF ( myid==0 ) then
             WRITE(*,*) "Enter the number of intervals (0 to quit)"
             READ(*,*) n
             t_start = ETIME(t)
           END IF

           CALL MPI_BCAST(n,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
           IF ( n<=0 ) EXIT

           h = 1.0d0/n
           sum  = 0.0d0
           DO i = myid+1, n, numprocs
             x = h * (DBLE(i) - 0.5d0)
             sum = sum + 4.d0/(1.d0 + x*x)
           END DO
           mypi = h * sum

           CALL MPI_REDUCE(mypi,pi,1,MPI_DOUBLE_PRECISION,MPI_SUM,0,
     +                     MPI_COMM_WORLD,ierr)

           IF ( myid==0 ) THEN
             WRITE(*,*) "Approximation for pi is ", pi
             WRITE(*,*) "The error is ", abs(pi - PI25DT)
             t_end = ETIME(t)
             WRITE(*,*) "Elapsed time is", t_end-t_start
           END IF

        END DO

        CALL MPI_FINALIZE(rc)

      END PROGRAM pi_mpi