!********************************************************************
! program    : mpi_samp5.f90
! programmer : makoto nakajima
! description: compute approximate of pi using monte carlo method.
! date       : april 25, 2006
!********************************************************************
program main
	implicit none     !prohibit implicit declaration of variables
	include 'mpif.h'  !include mpi library

	!******* variables related mpi *******
	integer:: ierror !return error message from mpi subroutines
	integer:: id     !identification number of each processor
	integer:: nproc  !total number of processors

	!******* other variable declaration *******
	integer,parameter:: idmom=0   !node which works as the frontend.
	integer,parameter:: n=1000000000 !total number of random numbers used.
	integer:: n_each  !number if integers assigned to each node
	integer:: n_top   !for each node, sum(n_top:n_end) is caomputed.
	integer:: n_end
	integer:: n0        !intdex for loop
	real(8):: rnumber1  !draw of random number (1)
	real(8):: rnumber2  !draw of random number (2)
	real(8):: tot_ind   !total number of draws
	real(8):: hit_ind   !number inside the unit circle
	real(8):: tot_all   !sum of total number of draws
	real(8):: hit_all   !sum of number inside the unit circle
	real(8):: pi_approx !approximate of pi

	!******* initialization *******
	!initialization of mpi environment
	call mpi_init(ierror)

	!obtain id for each node
	call mpi_comm_rank(mpi_comm_world, id, ierror)

	!returns the number of nodes
	call mpi_comm_size(mpi_comm_world, nproc, ierror)

	!******* compute Neach *******
	!determies number of integer assigned to each node
	do n0=1,n
		if (n0*nproc>=n) exit
	end do

	!******* compute Ntop and Nend *******
	n_top=id*n0+1
	n_end=min((id+1)*n0,n) !notice n0*nproc might be above n

	!******* do random draws *******
	do n0=n_top, n_end
		call random_number(rnumber1) !draw random numbers
		call random_number(rnumber2)
		tot_ind=tot_ind+1.0_8        !increase total number of draws
		if (rnumber1**2.0_8+rnumber2**2.0_8<1.0_8) then !inside unit circle
			hit_ind=hit_ind+1.0_8      !increase number of hits
		end if
	end do

	!******* gather all sum using mpi_reduce *******
	tot_all=0.0_8     !initialize total number of draws
	call mpi_reduce(tot_ind,tot_all,1,mpi_double_precision,&
	mpi_sum,idmom,mpi_comm_world,ierror)

	hit_all=0.0_8     !initialize total number of hits
	call mpi_reduce(hit_ind,hit_all,1,mpi_double_precision,&
	mpi_sum,idmom,mpi_comm_world,ierror)

	!******* idmom comnputes the approximated pi and print out *******
	if (id==idmom) then
		pi_approx=hit_all/tot_all*4.0_8
		print *,'I am id ',id
		print *,'Obtained approximate of pi is ',pi_approx
	end if

	!******* finalize mpi environment *******
	call mpi_finalize(ierror)

end program main