MPI: Calculation of Pi Using the Monte Carlo Method

The “Monte Carlo Method” is a stochastic method of solving a problem.

I use this method in calculating the value of PI for my next MPI code. But first, I will explain the idea behind the method and how I implemented it in MPI and the performance of my code on my Hackintosh.

If a circle of radius is inscribed inside a square with side length , then the area of the circle will be and the area of the square will be . So the ratio of the area of the circle to the area of the square will be .

So we pick points at random inside the square. The random selection of points makes the method stochastic. We checkt to see of the point is in fact inside the circle, this is done by checking the following inequality

We keep track of the number of points that reside inside the circle, . PI is them approximated as

The Monte Carlo method is often used to approximate solutions to problems in physics and mathematics. There is a draw back using this method however, it is very slow. So I thought to parallelize this method using MPI. It is quite simple. It only took me 5min to code and debug (I’m not a great programmer either)

The MPI code is in C as alway:

#include <stdio.h>
#include <stdlib.h>
#include <mpi.h>
#include <math.h>
#include <time.h>

int main(int argc, char* argv[]){
int i,id, np,N;
double x, y,double_N,eTime,sTime,pTime;
int lhit;
MPI_Init(&argc, &argv);
MPI_Comm_rank(MPI_COMM_WORLD, &id);
MPI_Comm_size(MPI_COMM_WORLD, &np);

if( argc !=2){
if (id==0){
fprintf(stderr,"Incorrect NARGIN\n");
fflush(stderr);
}
MPI_Abort(MPI_COMM_WORLD,1);
}

sscanf(argv[1], "%lf", &double_N);
N = lround(double_N);
MPI_Barrier(MPI_COMM_WORLD);
sTime = MPI_Wtime();
lhit = 0;
srand((unsigned)(time(0)));
int lN = N/np;

for(i = 0; i<lN;i++){
x = ((double)rand())/((double)RAND_MAX);
y = ((double)rand())/((double)RAND_MAX);
if (((x*x) + (y*y)) <= 1) lhit++;
}

int hit=0;
MPI_Allreduce(&lhit,&hit,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
double est;
est = (hit*4)/((double)N);
MPI_Barrier(MPI_COMM_WORLD);

eTime = MPI_Wtime();
pTime = fabs(eTime - sTime);

if (id == 0) {
printf("Number of Points Used:		%d\n",N);
printf("Estimate of Pi:			%24.16f\n",est);
printf("Elapsed Wall time:		%5.3e\n",pTime);
}

MPI_Finalize();
return 0;
}

The runtime result of the above code on my quad-core I7
It is obvious fromt the graph that this method is much faster running on multiple processors. The second graph is the same code, only using GSL random number generator. It appears to be more efficient, and it has a better approximation to PI. One thing troubles me though. The I7 is hyperthreads, so I can use up to 8-cores, but the result is not “stable”
regardless, my machine still rocks

About these ads