!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!									!
!  pi_f90.f
!
!  31 March, 98 							!
!  Amy C. Sun, AHPCC/UNM						!
!									!	
!  Modified by John Grondalski AHPCC/UNM	6/21/99			!
!									!	
!  Compute pi as 4.*Integral{0->1}[1/(1+x^2)] on multiple processors	!
!  using MPI.  See "Using MPI", W. Gropp, E. Lusk, and A. Skjellum,	!
!  MIT Press, p23, (1994).						!
!  Modified from Paul Alsing's original F77 version                     !
!									!
!  This code can be used interactively or with LoadLeveler with input	!
!  file 'pi.input'. The code will ask user for the # of intervals N 	!	
!  used to compute pi.  						!
!									!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

program pi_f90

	implicit none
	include "mpif.h"

	double precision  pi_part, pi, h, sum, x, f, a, pi25dt
	double precision  starttime, endtime, runtime, runtime_max
	integer N, rank, nprocrs, i, answer, ierr
	parameter (PI25DT = 3.141592653589793238462643d0)

!  function to integrate

	f(a) = 4.0 / (1.0 + a*a)

!  Initialize MPI

	call MPI_INIT(ierr)
	call MPI_COMM_RANK(MPI_COMM_WORLD, rank, ierr)
	call MPI_COMM_SIZE(MPI_COMM_WORLD, nprocrs, ierr)

!  Message to User

       if(rank == 0) then
	   print*,'This program computes pi as'
	   print*,'4.*Integral(0->1)[1/(1+x^2)]'
 	   print*,'in parallel on ',nprocrs,' processors.'
      endif

!  Loop over User inputs and pi computations
         
	do

!  User inputs

  	   if(rank  == 0) then
	      print*,'Enter the number of intervals N: '
	      read(5,*) N
	      print*,'pi will be computed with ',N,' intervals on ', & 
                    nprocrs,' processors'
	      print*,''
	   endif

!  Procr 0 = P(0) gives n to all other procrs

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

!  P(rank) computes pi on intervals rank+1 -> N modulo nprocrs.
!  Note: rank={0,1,2,....,nprocrs-1}
!  Note: we will time how long each procr takes with MPI_Wtime() 

	   starttime=MPI_Wtime()

	   sum=0.0
	   h=1.0/N

	   do i = rank+1, N, nprocrs
	      x=h*(i-0.5)
	      sum = sum + f(x)
	   enddo

	   pi_part = sum*h

	   endtime=MPI_Wtime()
	   runtime=endtime-starttime

!  Reduce final answer of pi onto procr P(0)

	   call MPI_REDUCE(pi_part,pi,1,MPI_DOUBLE_PRECISION, &
     		      MPI_SUM,0,MPI_COMM_WORLD,ierr)
	   call MPI_REDUCE(runtime,runtime_max,1,MPI_DOUBLE_PRECISION, &
     		      MPI_MAX,0,MPI_COMM_WORLD,ierr)

	   if (rank == 0) then
	      print*,'For ',N,' intervals, pi = ',pi
	      print*,'with error = ',abs(pi-PI25DT)
	      print*,'computed in ',runtime_max,' secs.'
	      print*,' ' 
              print*,'Do you wish to try another run? (y=1;n=0)'
              read(5,*) answer
              print*,'answer = ',answer
	   endif

!  Program does not continue until all procrs have reached this barrier
!  (procrs P(>0) are waiting for P(0) to get here)
           
           call MPI_BCAST(answer,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)

	   if (answer == 0) then
		exit
	   endif
	
	enddo

	print*,'Thank you for computing pi from processor ',rank

	   call MPI_FINALIZE(ierr)


end program pi_f90
