monte_carlo.c

/* compute pi using Monte Carlo method */
#include <stdio.h>
#include <limits.h>
#include <stdlib.h>
#include <math.h>
#include "mpi.h"
#define CHUNKSIZE      1000
/* message tags */
#define REQUEST  1
#define REPLY    2

int main(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 status;

    MPI_Init(&argc, &argv);
    world = MPI_COMM_WORLD;
    MPI_Comm_size(world, &numprocs);
    MPI_Comm_rank(world, &myid);
    server = numprocs-1;	/* last proc is server */
    if (myid == 0) {
        if (argc < 2) {
            fprintf(stderr, "Usage: %s epsilon\n", argv[0] );
            MPI_Abort(MPI_COMM_WORLD, 1);
        }
        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) {	/* I am the rand server */
        do {
            MPI_Recv(&request, 1, MPI_INT, MPI_ANY_SOURCE, REQUEST,
                     world, &status);
            if (request) {
                for (i = 0; i < CHUNKSIZE; ) {
                    rands[i] = random();
                    if (rands[i] <= INT_MAX) i++;
                }
                MPI_Send(rands, CHUNKSIZE, MPI_INT,
                         status.MPI_SOURCE, REPLY, world);
            }
        }
        while(request > 0);
    }
    else {			/* I am a worker process */
        request = 1;
        done = in = out = 0;
        max  = INT_MAX;         /* max int, for normalization */
        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, MPI_STATUS_IGNORE);
            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) > 100000000);
            request = (done) ? 0 : 1;
            if (myid == 0) {
                printf( "\rpi = %23.20f", 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, <ret> to exit\n",
                totalin+totalout, totalin, totalout );
        getchar();
    }
    MPI_Finalize();
    return 0;
}