The Code

/*
 * %Id: randompie.c,v 1.1 2003/10/12 20:22:50 gustav Exp %
 *
 * %Log: randompie.c,v %
 * Revision 1.1  2003/10/12 20:22:50  gustav
 * Initial revision
 *
 *
 */

#include <stdio.h>
#include <stdlib.h> /* function random is defined here */
#include <limits.h> /* INT_MAX is defined here */
#include <mpi.h>

#define CHUNKSIZE 1000
#define REQUEST      1
#define REPLY        2
#define TOTAL     5000000

main(argc, argv)
     int argc;
     char *argv[];
{
  int iter;
  int in, out, i, iters, max, ix, iy, ranks[1], done, temp;
  double x, y, Pi, error, epsilon;
  int numprocs, myid, server, totalin, totalout, workerid;
  int rands[CHUNKSIZE], request;
  MPI_Comm world, workers;
  MPI_Group world_group, worker_group;
  MPI_Status stat;

  MPI_Init(&argc, &argv);
  world = MPI_COMM_WORLD;
  MPI_Comm_size(world, &numprocs);
  MPI_Comm_rank(world, &myid);
  server = numprocs - 1;
  if (myid == 0)
    sscanf( argv[1], "%lf", &epsilon);
  MPI_Bcast(&epsilon, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Comm_group( world, &world_group);
  ranks[0] = server;
  MPI_Group_excl(world_group, 1, ranks, &worker_group);
  MPI_Comm_create(world, worker_group, &workers);
  MPI_Group_free(&worker_group);

  if(myid == server) {
    do {
      MPI_Recv(&request, 1, MPI_INT, MPI_ANY_SOURCE, REQUEST, world, &stat);
      if (request) {
        for (i = 0; i < CHUNKSIZE; i++)
          rands[i] = random();
        MPI_Send(rands, CHUNKSIZE, MPI_INT, stat.MPI_SOURCE, REPLY, world);
      }
    }
    while (request> 0);
  }
  else {
    request = 1;
    done = in = out = 0;
    max = INT_MAX;
    MPI_Send(&request, 1, MPI_INT, server, REQUEST, world);
    MPI_Comm_rank(workers, &workerid);
    iter = 0;
    while(!done) {
      iter++;
      request = 1;
      MPI_Recv(rands, CHUNKSIZE, MPI_INT, server, REPLY, world, &stat);
      for (i=0; i < CHUNKSIZE; ) {
        x = (((double) rands[i++])/max) * 2 - 1;
        y = (((double) rands[i++])/max) * 2 - 1;
        if (x*x + y*y < 1.0)
          in++;
        else
          out++;
      }
      MPI_Allreduce(&in, &totalin, 1, MPI_INT, MPI_SUM, workers);
      MPI_Allreduce(&out, &totalout, 1, MPI_INT, MPI_SUM, workers);
      Pi = (4.0*totalin)/(totalin + totalout);
      error = fabs( Pi-3.141592653589793238462643);
      done = ((error < epsilon) || ((totalin+totalout) > TOTAL));
      request = (done) ? 0 : 1;
      if (myid == 0) {
        printf( "\rpi = %23.20lf", Pi);
        MPI_Send(&request, 1, MPI_INT, server, REQUEST, world);
      }
      else {
        if (request)
          MPI_Send(&request, 1, MPI_INT, server, REQUEST, world);
      }
    }
    MPI_Comm_free(&workers);
  }
  if (myid == 0) 
    printf( "\npoints: %d\nin: %d, out: %d\n", 
            totalin+totalout, totalin, totalout);
  MPI_Finalize();
  exit (0);
}

Here's how the code is compiled and installed:

gustav@bh1 $ pwd
/N/B/gustav/src/I590/randompie
gustav@bh1 $ make       
co  RCS/Makefile,v Makefile
RCS/Makefile,v  -->  Makefile
revision 1.1
done
co  RCS/randompie.c,v randompie.c
RCS/randompie.c,v  -->  randompie.c
revision 1.1
done
mpicc -o randompie randompie.c
gustav@bh1 $ make install
install randompie /N/B/gustav/bin
gustav@bh1 $

And here's how the code is run. Observe that the MPI program takes one argument from the command line, namely the $\epsilon$. It is simply typed on the command line following the name of the MPI program, randompie.

gustav@bh1 $ mpdtrace | wc
     32      32     159
gustav@bh1 $ mpiexec -n 16 randompie 0.0001 
pi =  3.14164761904761924427
points: 420000
in: 329873, out: 90127
gustav@bh1 $