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 &lt;- function(DARTS, ROUNDS, proc_num, numprocs) {
  retvals &lt;- .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  &lt;- function(DARTS, ROUNDS, numprocs){
  cl &lt;- snow::makeCluster(numprocs, type = "MPI")
  r&lt;&lt;-ROUNDS
  d&lt;&lt;-DARTS
  snow::clusterExport(cl, c("r", "d"))
  x&lt;-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 &lt;&lt;-. 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) &lt;= 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.
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 &lt;&lt;- 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 &lt;- function(num_simulations, seed_num = 123456L) {
  retvals &lt;- .Fortran("snowpi",
                      num_inside      = 1L,  # an integer
                      num_simulations = as.integer(num_simulations),
                      seed_num = as.integer(seed_num))
  return(retvals$num_inside)
}

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

  pi_estimate &lt;- 4 * Reduce(`+`, num_inside) / num_simulations
  return(pi_estimate)
}

balance_evenly &lt;- function(total, among) {
   min_load &lt;- floor(total / among)
   max_load &lt;- min_load + 1L
   rest &lt;- total - min_load * among
   load &lt;- c(rep(max_load, rest),
             rep(min_load, among - rest))
   return(load)
}


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.