# 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] }