Merge pull request #100 from JuliaParallel/eschnett/montecarlo

Clean up Monte Carlo example

Author: eschnett

examples/07-pi-impl.jl

@@ -0,0 +1,63 @@
+import MPI
+function montecarlo(mc_eval::Function, mc_monitor::Function,
+                    reps::Integer, n_returns::Integer, pooled::Integer=1,
+                    sleeptime::Real=0)
+    # set up the MPI communicator
+    MPI.Init()
+    comm = MPI.COMM_WORLD
+    rank = MPI.Comm_rank(comm)
+    commsize = MPI.Comm_size(comm)
+
+    # containers and bookkeeping
+    contrib = zeros(pooled, n_returns)
+    results = zeros(reps, n_returns)
+    pernode = div(reps, commsize - 1)
+    if mod(reps, commsize - 1) != 0
+        if rank == 0
+            println("Error (montecarlo):")
+            println("    reps=$reps, commsize=$commsize")
+            println("Choose commsize so that reps/(commsize-1) is integer")
+        end
+        MPI.Finalize()
+        exit(1)
+    end
+
+    if mod(pernode, pooled) != 0
+        pooled = div(pernode, commsize-1)
+    end
+
+    if rank > 0
+        # workers
+        @inbounds for i = 1:div(pernode, pooled)
+            # do work
+            for j = 1:pooled
+                contrib[j,:] = mc_eval()
+            end
+            MPI.Isend(contrib, 0, rank, comm)
+        end
+
+    else
+        # manager
+        sofar = 0 # results collected so far
+        while sofar < div(reps, pooled)
+            sleep(sleeptime) # flood control if job is costly
+            @inbounds for node = 1:commsize-1
+                # check for results
+                ready, _ = MPI.Iprobe(node, node, comm)
+                if ready # get them if they're ready
+                    sofar +=1
+                    if sofar*pooled <= reps
+                        MPI.Recv!(contrib, node, node, comm)
+                        results[sofar*pooled-pooled+1:sofar*pooled,:] = contrib
+                        mc_monitor(sofar*pooled, results)
+                    end
+                    if sofar == div(reps, pooled)
+                        break
+                    end
+                end
+            end
+        end
+    end
+
+    MPI.Finalize()
+end

examples/07-pi.jl

@@ -7,31 +7,26 @@ like to use, subject to  X-1 being an even divisor
 of 1e6, e.g., set X=5.
 =#
 
-include("montecarlo.jl")
+include("07-pi-impl.jl")
 
 function pi_wrapper()
-    pihat = 4.*float((norm(rand(2,1)) .< 1.))
+    4.0 * (norm(rand(2)) < 1)
 end
 
 # this function reports intermediate results during MC runs
 function pi_monitor(sofar, results)
-    # examine results every 2.5*10^5 draws
-    if mod(sofar,2.5e5)==0.
+    # Examine results every 12.5*10^5 draws
+    if mod(sofar, 10^6) == 0
         m = mean(results[1:sofar,:],1)
-        println("reps so far: ", sofar)
-        println("pihat: ", m)
-        println()
-        #if sofar == size(results,1)
-        #    writedlm("mcresults.out", results)
-        #end
+        println("reps: $sofar, pihat: $m")
     end
-end    
+end
 
 # do the monte carlo: 10^6 reps of single draws
 function main()
-    reps = Int(1e6)  # desired number of MC reps
+    reps = 10^7 # desired number of MC reps
     nreturns = 1
-    pooled = Int(50000)
+    pooled = 10^5
     montecarlo(pi_wrapper, pi_monitor, reps, nreturns, pooled)
 end