program int_pi2
c
c For API
c
c This simple program approximates pi by computing pi = integral
c from 0 to 1 of 4/(1+x*x)dx which is approximated by sum from
c k=1 to N of 4 / ((1 + (k-1/2)**2 ).  The only input data required is N.
c
       include "mpif.h"
       integer ierr, status(MPI_STATUS_SIZE)
       real err, f, pi, sum, w
       integer i, N, info, mynum, nprocs, source, dest, type
       integer nbytes, len, dontcare, nbuf(4)
       f(x) = 4.0/(1.0+x*x)
       pi = 4.0*atan(1.0)

       dest = 0
       type = 2
       len = 4
       nbytes = 0

c All instances call the startup routine to get their instance number (mynum)
       call MPI_INIT(ierr)
       call MPI_COMM_RANK(MPI_COMM_WORLD, mynum, ierr)
       call MPI_COMM_SIZE(MPI_COMM_WORLD, nprocs, ierr)

c -------  Each new approximation to pi begins here. -------------------
c (Step 1) Get value N for a new run
5      call solicit (N,nprocs,mynum)

c Step (2): check for exit condition.
       if (N .le. 0) then
	  call exit
       endif

c Step (3): do the computation in N steps
c Parallel Version: there are "nprocs" instances participating.  Each
c instance should do 1/nprocs of the calculation.  Since we want
c i = 1..n but mynum = 0, 1, 2..., we start off with mynum+1.
 
       w = 1.0/N
       sum = 0.0
       do i = mynum+1,N,nprocs
	  sum = sum + f((i-0.5)*w)
       enddo
       sum = sum * w

c Step (4): print the results  
c (Parallel version: collect partial results and let master instance print it)
       if (mynum.eq.0) then
	  print *,'host calculated x=',sum
 	  do i = 1,nprocs-1
             nbytes = len
             call MPI_RECV(x, 1, MPI_REAL, MPI_ANY_SOURCE,
     .        type, MPI_COMM_WORLD, status, ierr)
             print *,'host got x=',x
 	     sum=sum+x
 	  enddo
          err = sum - pi
          print *, 'sum, err =', sum, err
c Other instances just send their sum and go back to 5 and wait for more input
       else
          call MPI_SEND(sum, 1, MPI_REAL, dest, type, 
     .       MPI_COMM_WORLD, ierr)
	  print *,'instance',mynum,' sent partial sum',sum,
     .            ' to instance 0'
       endif
       go to 5
       call MPI_FINALIZE(ierr)
       end



       subroutine solicit (N,nprocs,mynum)
c Get a value for N, the number of intervals in the approximation
c (Parallel versions: master instance reads in N and then
c broadcasts N to all the other instances of the program)
       include "mpif.h"
       integer ierr
       integer source, msglen,nbuf(4)
       source = 0
       msglen = 4
       if (mynum .eq. 0) then
          print *,'Enter number of approximation intervals:(0 to exit)'
          read *, N
       endif
       call MPI_BCAST(N, 1, MPI_INTEGER, source, MPI_COMM_WORLD, ierr)
       return
       end