!
! Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
! Use is subject to license terms.
!
! Sun, Sun Microsystems, the Sun logo, Sun HPC ClusterTools, Sun PFS,
! Sun C++, Sun MPI, Prism, Sun Prism, and all Sun-based trademarks and logos,
! are trademarks or registered trademarks of Sun Microsystems, Inc. in
! the United States and other countries.
!

!
! Estimate pi via Monte-Carlo method.
!   		
! Each process sums how many of samplesize random points generated 
! in the square (-1,-1),(-1,1),(1,1),(1,-1) fall in the circle of 
! radius 1 and center (0,0), and then estimates pi from the formula
! pi = (4 * sum) / samplesize.
! The final estimate of pi is calculated at rank 0 as the average of 
! all the estimates.
!
	program monte

	include 'mpif.h'

	double precision drand
	external drand

	double precision x, y, pi, pisum
        integer*4 ierr, rank, np
	integer*4 incircle, samplesize

	parameter(samplesize=2000000)

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

!	seed random number generator
	x = drand(2 + 11*rank)

	incircle = 0
	do i = 1, samplesize
	   x = drand(0)*2.0d0 - 1.0d0           ! generate a random point
	   y = drand(0)*2.0d0 - 1.0d0

	   if ((x*x + y*y) .lt. 1.0d0) then
	      incircle = incircle+1	        ! point is in the circle
	   endif
	end do

	pi = 4.0d0 * DBLE(incircle) / DBLE(samplesize)

!	sum estimates at rank 0
	call MPI_REDUCE(pi, pisum, 1, MPI_DOUBLE_PRECISION, MPI_SUM, 
     &		0 , MPI_COMM_WORLD, ierr)

	if (rank .eq. 0) then
!	   final estimate is the average
	   pi = pisum / DBLE(np)
	   print '(A,I4,A,F8.6,A)','Monte-Carlo estimate of pi by ',np,
     &		 ' processes is ',pi,'.'
	endif

	call MPI_FINALIZE(ierr)
	end