next up previous index
Next: The Discussion Up: Manipulating Communicators Previous: Manipulating Communicators

The Code

Here is the code:

#include <stdio.h>
#include <math.h>
#include "mpi.h"
#define CHUNKSIZE 1000
#define REQUEST      1
#define REPLY        2

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) > 1000000));
      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();
}
And here is how to compile and run this code:
gustav@sp20:../MPI 17:39:37 !550 $ mpcc -o pirand pirand.c
gustav@sp20:../MPI 17:39:50 !551 $ cat pirand.ll
# @ job_type = parallel
# @ environment = COPY_ALL; MP_EUILIB=ip; MP_INFOLEVEL=2
# @ requirements = (Adapter == "hps_ip") 
# @ min_processors = 6
# @ max_processors = 10
# @ class = test
# @ output = pirand.out
# @ error = pirand.err
# @ executable = /usr/bin/poe
# @ arguments = pirand
# @ queue
gustav@sp20:../MPI 17:39:57 !552 $ llsubmit pirand.ll
submit: The job "sp20.135" has been submitted.
gustav@sp20:../MPI 17:40:14 !553 $ cat pirand.out
pi =  3.13959541604384639868
points: 1003500
in: 787646, out: 215854
gustav@sp20:../MPI 17:41:30 !554 $
You can also invoke program pirand with an option that specifies the value of the $\varepsilon$ parameter, which is used to assess the accuracy of the computation. Let us try to do this interactively a few times:
gustav@sp20:../MPI 17:41:30 !554 $ poe pirand 0.001 -procs 8
pi =  3.14158730158730170601
points: 31500
in: 24740, out: 6760
gustav@sp20:../MPI 17:43:32 !555 $ poe pirand 0.00001 -procs 8
pi =  3.14158730158730170601
points: 31500
in: 24740, out: 6760
gustav@sp20:../MPI 17:43:46 !556 $ poe pirand 0.000001 -procs 8
pi =  3.13962437562437557403
points: 1001000
in: 785691, out: 215309
gustav@sp20:../MPI 17:44:13 !557 $ poe pirand 0.01 -procs 8
pi =  3.15085714285714280081
points: 3500
in: 2757, out: 743
gustav@sp20:../MPI 17:44:28 !558 $
Observe that asking for too good a resolution does not pay. The program runs up to the preset limit of one million points and returns a result that is less accurate than the result returned with only 31,500 points, which is accurate to within 0.00001. This tells us something about the random number generator, which, clearly, isn't random enough for long sequences.



Zdzislaw Meglicki
2001-02-26