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))
SnowPi  <- function(DARTS, ROUNDS, numprocs){
  cl <- snow::makeCluster(numprocs, type = "MPI")
  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))

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)
      mynpts = ROUNDS/numprocs

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

    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 


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.

  • \$\begingroup\$ 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.) \$\endgroup\$ – flodel Aug 13 '15 at 0:01
  • \$\begingroup\$ @flodel I just used microbenchmark (10 times), Fortran is 2.1 times faster than the vectorized R code. \$\endgroup\$ – Ignacio Aug 13 '15 at 4:31

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))

SnowPi  <- function(num_simulations, num_procs) {
  # input checks
            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

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

balance_evenly <- function(total, among) {
   min_load <- floor(total / among)
   max_load <- min_load + 1L
   rest <- total - min_load * among
   load <- c(rep(max_load, rest),
             rep(min_load, among - rest))

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.

  • \$\begingroup\$ 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. \$\endgroup\$ – 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


! 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)
   seed = 12 + proc_num*11
   call random_seed(put=seed(1:n))

   ! 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

Edit: Just realized, that flodel already addressed the averaging in his answer.

If you are concerned about elimination errors in the averaging, you would need to use some tree for the summation, I think. The updating average does not really help for that. It should be also fine to put the division outside the loop, which of course is much better for the performance.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.