# Rmpi code to find prime numbers

# for illustration purposes, not intended to be optimal

# returns vector of all primes in 2..n; the vector "divisors" is used as
# a basis for a Sieve of Erathosthenes operation; must have n <=
# (max(divisors)^2) and n even

# the argument "msgsize" controls the chunk size in communication from
# the manager to the first worker, node 1

primepipe <- function(n,divisors,msgsize) {
   mpi.bcast.Robj2slave(dowork)
   mpi.bcast.Robj2slave(dosieve)
   # start workers; note nonblocking call
   mpi.bcast.cmd(dowork,n,divisors,msgsize)
   # remove the evens right away
   odds <- seq(from=3,to=n,by=2)
   nodd <- length(odds)
   # send odds to node 1, in chunks 
   startmsg <- seq(from=1,to=nodd,by=msgsize)
   for (s in startmsg) {
      rng <- s:min(s+msgsize-1,nodd)
      mpi.send.Robj(odds[rng],tag=0,dest=1)
   }
   # send end-data sentinel
   mpi.send.Robj(NA,tag=0,dest=1)
   # receive results from last node
   lastnode <- mpi.comm.size()-1
   # return te result; don't forget the 2
   c(2,mpi.recv.Robj(tag=0,source=lastnode))
}

# worker code
dowork <- function(n,divisors,msgsize) {
   # which are my divisors?
   me <- mpi.comm.rank()
   lastnode <- mpi.comm.size()-1
   ld <- length(divisors)
   tmp <- floor(ld / lastnode)
   mystart <- (me-1) * tmp + 1
   myend <- mystart + tmp - 1
   if (me == lastnode) myend <- ld
   mydivs <- divisors[mystart:myend]
   if (me == lastnode) out <- NULL
   repeat {
      msg <- mpi.recv.Robj(tag=0,source=me-1)
      if (me < lastnode) {
         if (!is.na(msg[1])) {
            sieveout <- dosieve(msg,mydivs)
            mpi.send.Robj(sieveout,tag=0,dest=me+1)
         } else {  # no more coming in, so send sentinel
            mpi.send.Robj(NA,tag=0,dest=me+1)
            return()
         }
      } else {
         if (!is.na(msg[1])) {
            sieveout <- dosieve(msg,mydivs)
            out <- c(out,sieveout)
         } else {  # no more coming in, so send results to manager
            mpi.send.Robj(out,tag=0,dest=0)
            return()
         }
      }
   }
}

# check divisibility of the current chunk x
dosieve <- function(x,divs) {
   for (d in divs) {
      x <- x[x %% d != 0 | x == d]
   }
   x
}

# serial prime finder; can be used to generate divisor list of primepipe
serprime <- function(n) {
   nums <- 1:n
   x <- rep(1,n)
   maxdiv <- ceiling(sqrt(n))
   for (d in 2:maxdiv) {
      if (x[d])
         x[x !=0 & nums > d & nums %% d == 0] <- 0
   }
   nums[x != 0 & nums >= 2]
}