# Monte Carlo Pi with R + Fortran + MPI

I created a R package to calculate Pi with a Monte Carlo simulation.

My goal is to optimize the parallel function SnowPi as much as possible. This is the SnowPi.R code:

FSnowPi <- function(DARTS, ROUNDS, proc_num, numprocs) {
retvals <- .Fortran("snowpi", avepi = as.numeric(1),
DARTS =  as.integer(DARTS),
ROUNDS =  as.integer(ROUNDS),
proc_num = as.integer(proc_num),
numprocs = as.integer(numprocs))
return(retvals$avepi) } #'@export SnowPi <- function(DARTS, ROUNDS, numprocs){ cl <- snow::makeCluster(numprocs, type = "MPI") r<<-ROUNDS d<<-DARTS snow::clusterExport(cl, c("r", "d")) x<-snow::clusterEvalQ(cl, MyPi::FSnowPi(DARTS = d, ROUNDS = r, proc_num = mpi.comm.rank(), numprocs = mpi.comm.size()-1)) snow::stopCluster(cl) return(mean(unlist(x))) }  One thing that I think may not be optimal is the use of <<-. Alas, I don't know how to get the code to work without that (I think the problem is how to pass the correct environment to clusterExport) Finally, these are the 2 relevant Fortran subroutines: subroutine snowpi(avepi, DARTS, ROUNDS, proc_num, numprocs) bind(C, name="snowpi_") use, intrinsic :: iso_c_binding, only : c_double, c_int real(c_double), intent(out) :: avepi integer(c_int), intent(in) :: DARTS, ROUNDS, proc_num, numprocs integer :: i, n, mynpts integer, allocatable :: seed(:) double precision :: pi_est if (proc_num .eq. 1) then mynpts = ROUNDS - (numprocs-1)*(ROUNDS/numprocs) else mynpts = ROUNDS/numprocs endif ! initialize the random number generator ! we make sure the seed is different for each task call random_seed() call random_seed(size = n) allocate(seed(n)) seed = 12 + proc_num*11 call random_seed(put=seed(1:n)) deallocate(seed) avepi=0.0d0 do i = 0, mynpts-1 call dboard(darts, pi_est) ! calculate the average value of pi over all iterations avepi = ((avepi*i) + pi_est)/(i + 1) end do end subroutine snowpi  and subroutine dboard(darts, dartsscore) integer, intent(in) :: darts double precision, intent(out) :: dartsscore double precision :: x_coord, y_coord integer :: score, n score = 0 do n = 1, darts call random_number(x_coord) call random_number(y_coord) if ((x_coord**2 + y_coord**2) <= 1.0d0) then score = score + 1 end if end do dartsscore = 4.0d0*score/darts end subroutine dboard  I'm not sure if there is something that I could improve there. • I'd be curious to know how fast your Fortran code is compared to something written only in R as this one for example: stackoverflow.com/a/15186385/1201032. I bet R is not terribly slower since this implementation is all vectorized (i.e. using C under the hood.) – flodel Aug 13 '15 at 0:01 • @flodel I just used microbenchmark (10 times), Fortran is 2.1 times faster than the vectorized R code. – Ignacio Aug 13 '15 at 4:31 ## 2 Answers Here are a few pointers. I should first say that I am more familiar with R; I know next to nothing about Fortran but I could read your code and guess about a few things that will likely improve your code: both in simplicity and speed. 1) Why the need for two inputs DARTS and ROUNDS? You know that in the end you will have a total of DARTS * ROUNDS simulations: could there just be a single input num_simulations that will be split evenly among your numprocs? This should get you rid of the loop inside snowpi so it could only call dboard once. 2) You should make your slave (snowpi) as simple as possible. IMO, it should not know about things like proc_num or numprocs if you just tell it how many simulations it needs to do. Let R do the math upstream and pass it down to the slave. I don't know much about random generators but maybe the same can apply with the random seeds: if you could let R compute a seed for each slave and pass it down as an (optional) argument to the slave, that would be best. (However, if knowing the proc_num and numprocs is crucial for seed selection and can only happen on the Fortran side, then disregard my comment...) 3) This one I feel is important: do not compute averages too early. They are costly and can make you lose precision. This line in your code (avepi = ((avepi*i) + pi_est)/(i + 1)) is particularly problematic: you are doing so many unnecessary multiplications and divisions... Instead, you should just be counting the number of darts, i.e. incrementing a counter. Your slave should just be returning an integer. Let R add all these integers together and do a single division by num_simulations at the end. To summarize 1, 2, and 3, it would be great (if possible) if you could just tell your slave: here is your seed number n; now throw x darts and tell me how many ended up in the circle. 4) Now looking at your R code. Yes, the use of <<- in your code is not very elegant. That's because you are trying to use clusterEvalQ which may not be the best option from the snow package. Instead, clusterApply or clusterMap are better suited for functional programming. Your code could look like this, I let you figure out the rest: FSnowPi <- function(num_simulations, seed_num = 123456L) { retvals <- .Fortran("snowpi", num_inside = 1L, # an integer num_simulations = as.integer(num_simulations), seed_num = as.integer(seed_num)) return(retvals$num_inside)
}

SnowPi  <- function(num_simulations, num_procs) {
# input checks
stopifnot(is.integer(num_simulations),
is.integer(num_procs),
num_simulations > 0L,
num_procs > 0L)
# initialize cluster
cl <- snow::makeCluster(num_procs, type = "MPI")
# balance the number of simulations per node
node_num_simu <- balance_evenly(num_simulations, num_procs)
# compute random seeds
node_seeds <- sample.int(1e6, num_procs)   # please review...
# run in parallel
num_inside <- snow::clusterMap(cl, MyPi::FSnowPi, node_num_simu, node_seeds)
# close the cluster
snow::stopCluster(cl)

pi_estimate <- 4 * Reduce(+, num_inside) / num_simulations
return(pi_estimate)
}

balance_evenly <- function(total, among) {
rest <- total - min_load * among
}


I hope some of it is helpful. Please let me know. And thanks for sharing your Fortran code; I know at least one other R package that uses it: fields whose rdist() function for computing euclidean distances is extremely fast. I want to give a try one day.

• A reason to use darts and rounds could be to allow an adjustment to cache sizes and allow for vectorized operations. I have no idea wether this would pay off in this case, but I try to illustrate this in my answer. Might be worth a try, optimal settings for darts/rounds would probably be system dependent. – haraldkl Sep 4 '15 at 20:55

I have no clue about R, but have some comments on the Fortran and the MPI. You are distributing your work such, that one process takes the rest from an otherwise equal splitting, this might result in one process doing significantly different work than the others. A better load balancing with a difference of at most 1 between any two processes could be achieved by something like this:

mynpts = rounds / numprocs
if (proc_num < mod(rounds, numprocs)) mynpts = myntps + 1


Here:

! initialize the random number generator
! we make sure the seed is different for each task
call random_seed()
call random_seed(size = n)


the first call to random_seed is unnecessary.

When aiming for performance, try to avoid calls to small routines in loops. I would collapse your two routines and end up with the following:

!> Compute Pi via the Monte-Carlo method.
!! Do you really need to mimic the name-mangling here, or would bind(C, name="snowpi")
!! work as well?
subroutine snowpi(avepi, DARTS, ROUNDS, proc_num, numprocs) bind(C, name="snowpi_")
use, intrinsic :: iso_c_binding, only : c_double, c_int

!> Number of darts to use per round.
!! This should be number sufficiently small to fit into L1-Cache,
!! but sufficiently large to allow vectorized processing.
integer(kind=c_int), intent(in)  :: DARTS

!> Total number of rounds all processes together will need to compute.
!! Each process will compute a fraction of rounds/numprocs out of these.
integer(kind=c_int), intent(in)  :: ROUNDS

!> Rank of the calling MPI process.
integer(kind=c_int), intent(in)  :: proc_num

!> Total number of processes to work on this.
integer(kind=c_int), intent(in)  :: numprocs

!> Resulting approximation of Pi.
real(kind=c_double), intent(out) :: avepi

! Local variables
integer              :: i, n, mynpts
real(kind=c_double)  :: pi_est
integer, allocatable :: seed(:)
real(kind=c_double)  :: coord(darts, 2)
real(kind=c_double)  :: dart_weight

mynpts = rounds/numprocs
if (proc_num < modulo(rounds, numprocs)) mynpts = mynpts + 1

! initialize the random number generator
! we make sure the seed is different for each task
call random_seed(size = n)
allocate(seed(n))
seed = 12 + proc_num*11
call random_seed(put=seed(1:n))
deallocate(seed)

! Avoid division inside the loop.
dart_weight = 4.0_c_double / real(darts, kind=c_double)

! Actual computation of Pi
avepi = 0.0_c_double
do i = 0, mynpts-1
call random_number(coord)
! No idea, wether count would actually be faster than the
! self implemented do loop, but it is more concise.
score = count(coord(:,1)**2 + coord(:,2)**2 <= 1.0_c_double)
pi_est = dart_weight * score
avepi = avepi + pi_est
end do
avepi = avepi/real(mynpts, kind=c_double)

end subroutine snowpi