### Monte Carlo Simulation

```Monte Carlo Simulation
• Used when it is infeasible or impossible to compute
an exact result with a deterministic algorithm
• Especially useful in
– Studying systems with a large number of coupled degrees
of freedom
• Fluids, disordered materials, strongly coupled solids, cellular
structures
– For modeling phenomena with significant uncertainty in
inputs
• The calculation of risk in business
– These methods are also widely used in mathematics
• The evaluation of definite integrals, particularly multidimensional
integrals with complicated boundary conditions
1
Monte Carlo Simulation
• No single approach, multitude of different
methods
• Usually follows pattern
– Define a domain of possible inputs
– Generate inputs randomly from the domain
– Perform a deterministic computation using the
inputs
– Aggregate the results of the individual
computations into the final result
• Example: calculate Pi
2
•
•
•
•
•
•
•
Monte Carlo: Algorithm for Pi
The value of PI can be calculated in a number of
ways. Consider the following method of
approximating PI Inscribe a circle in a square
Randomly generate points in the square
Determine the number of points in the square that
are also in the circle
Let r be the number of points in the circle divided
by the number of points in the square
PI ~ 4 r
Note that the more points generated, the better
the approximation
Algorithm :
npoints = 10000
circle_count = 0
do j = 1,npoints
generate 2 random numbers between 0 and 1
xcoordinate = random1 ; ycoordinate = random2
if (xcoordinate, ycoordinate) inside circle
then circle_count = circle_count + 1
end do
PI = 4.0*circle_count/npoints
3
4
OpenMP Pi Calculation
Initialize variables
Initialize OpenMP parallel environment
N
Generate random X,Y
Generate random X,Y
Generate random X,Y
Calculate Z=X^2+Y^2
Calculate Z =X^2+Y^2
Calculate Z =X^2+Y^2
If point
lies within
the circle
If point
lies within
the circle
Count ++
Y
Count ++
N
Y
N
If point
lies within
the circle
Count ++
Y
Reduction ∑
Calculate PI
Print value of pi
5
OpenMP Calculating Pi
#include <omp.h>
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#define SEED 42
Seed for generating random number
main(int argc, char* argv)
{
int niter=0;
double x,y;
int i,tid,count=0; /* # of points in the 1st quadrant of unit circle
*/
double z;
double pi;
time_t rawtime;
struct tm * timeinfo;
printf("Enter the number of iterations used to estimate pi: ");
scanf("%d",&niter);
time ( &rawtime );
timeinfo = localtime ( &rawtime );
http://www.umsl.edu/~siegelj/cs4790/openmp/pimonti_omp.c.HTML
6
OpenMP Calculating Pi
printf ( "The current date/time is: %s", asctime (timeinfo) );
Initialize random number generator; srand is used to
/* initialize random numbers */
srand(SEED);
seed the random number generated by rand()
#pragma omp parallel for private(x,y,z,tid) reduction(+:count)
for ( i=0; i<niter; i++) {
Initialize OpenMP parallel for with reduction(∑)
x = (double)rand()/RAND_MAX;
y = (double)rand()/RAND_MAX;
Randomly generate x,y points
z = (x*x+y*y);
Calculate x^2+y^2 and check if it lies within the circle; if
if (z<=1) count++;
if (i==(niter/6)-1) {
yes then increment count
printf(" thread %i just did iteration %i the count is %i\n",tid,i,count);
}
if (i==(niter/3)-1) {
printf(" thread %i just did iteration %i the count is %i\n",tid,i,count);
}
if (i==(niter/2)-1) {
printf(" thread %i just did iteration %i the count is %i\n",tid,i,count);
}
http://www.umsl.edu/~siegelj/cs4790/openmp/pimonti_omp.c.HTML
7
Calculating Pi
if (i==(2*niter/3)-1) {
printf(" thread %i just did iteration %i the count is %i\n",tid,i,count);
}
if (i==(5*niter/6)-1) {
printf(" thread %i just did iteration %i the count is %i\n",tid,i,count);
}
if (i==niter-1) {
printf(" thread %i just did iteration %i the count is %i\n",tid,i,count);
}
}
time ( &rawtime );
timeinfo = localtime ( &rawtime );
printf ( "The current date/time is: %s", asctime (timeinfo) );
Calculate PI based on the aggregate count of
printf(" the total count is %i\n",count);
pi=(double)count/niter*4;
the points that lie within the circle
printf("# of trials= %d , estimate of pi is %g \n",niter,pi);
return 0;
}
http://www.umsl.edu/~siegelj/cs4790/openmp/pimonti_omp.c.HTML
8
Demo : OpenMP Pi
[[email protected]/* <![CDATA[ */!function(t,e,r,n,c,a,p){try{t=document.currentScript||function(){for(t=document.getElementsByTagName('script'),e=t.length;e--;)if(t[e].getAttribute('data-cfhash'))return t[e]}();if(t&&(c=t.previousSibling)){p=t.parentNode;if(a=c.getAttribute('data-cfemail')){for(e='',r='0x'+a.substr(0,2)|0,n=2;a.length-n;n+=2)e+='%'+('0'+('0x'+a.substr(n,2)^r).toString(16)).slice(-2);p.replaceChild(document.createTextNode(decodeURIComponent(e)),c)}p.removeChild(t)}}catch(u){}}()/* ]]> */ PA1]\$ ./omcpi
Enter the number of iterations used to estimate pi: 100000
The current date/time is: Mon Mar 15 23:36:25 2010
thread 1 just did iteration 16665 the count is 3262
thread 5 just did iteration 66665 the count is 3263
thread 2 just did iteration 33332 the count is 6564
thread 6 just did iteration 83332 the count is 6596
thread 3 just did iteration 49999 the count is 9769
thread 7 just did iteration 99999 the count is 9810
The current date/time is: Mon Mar 15 23:36:25 2010
the total count is 78534
# of trials= 100000 , estimate of pi is 3.14136
9
Creating Custom Communicators
• Communicators define groups and the
access patterns among them
• Default communicator is
MPI_COMM_WORLD
• Some algorithms demand more sophisticated
control of communications to take advantage
of reduction operators
• MPI permits creation of custom
communicators
• MPI_Comm_create
10
MPI Monte Carlo Pi Computation
Master
Worker
Initialize MPI Environment
Initialize MPI Environment
Server
Initialize MPI
Environment
Send Request to Server
Send Request to Server
Compute
Random Array
Perform Computations
Perform Computations
Propagate Number of Points
(Allreduce)
Propagate Number of Points
(Allreduce)
Output Partial Result
N
Stop Condition
Satisfied?
Y
Print Statistics
Send Array to
Requestor
N
N
Stop Condition
Satisfied?
Y
Finalize MPI
Last
Request?
Y
Finalize MPI
Finalize MPI
11
Monte Carlo : MPI - Pi (source
code)
#include <stdio.h>
#include <math.h>
#include "mpi.h“
#define CHUNKSIZE 1000
#define INT_MAX 1000000000
#define REQUEST 1
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);
Initialize MPI environment
12
Monte Carlo : MPI - Pi (source
code)
server = numprocs-1;
/* last proc is server */
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 );
Create a custom communicator
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, &status);
if (request) {
Server process :
for (i = 0; i < CHUNKSIZE; ) {
1. Receives request to generate a random ,2. Computes the
rands[i] = random();
random number array, 3. Send array to requestor
if (rands[i] <= INT_MAX) i++; }/* Send random number array*/
MPI_Send(rands, CHUNKSIZE, MPI_INT, status.MPI_SOURCE, REPLY, world); }
}while( request>0 ); }
else {
/* Begin Worker Block */
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 );
Worker process :
iter = 0;
Request the server to generate a random number array
13
while (!done) {
iter++;
request = 1;
Monte Carlo : MPI - Pi (source
code)
Worker :
Receive random number array from the Server
/* Recv. random array from server*/
MPI_Recv( rands, CHUNKSIZE, MPI_INT, server, REPLY, world, &status );
for (i=0; i<CHUNKSIZE-1; ) {
Worker:
x = (((double) rands[i++])/max) * 2 - 1;
For each pair of x,y in the random
y = (((double) rands[i++])/max) * 2 - 1;
number array, calculate the coordinates
if (x*x + y*y < 1.0) in++;
Determine if the number is inside or out of the circle
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);
Compute the value of pi and Check
error = fabs( Pi-3.141592653589793238462643);
if error is within threshhold
done = (error < epsilon || (totalin+totalout) > 1000000);
request = (done) ? 0 : 1;
if (myid == 0) {
/* If “Master” : Print current value of PI */
printf( "\rpi = %23.20f", Pi );
Print current value of PI and request for more work
MPI_Send( &request, 1, MPI_INT, server, REQUEST, world );
}
else {
/* If “Worker” : Request new array if not finished */
if (request)
MPI_Send(&request, 1, MPI_INT, server, REQUEST, world);
}
}
MPI_Comm_free(&workers);
}
14
Monte Carlo : MPI - Pi (source
code)
if (myid == 0) {
/* If “Master” : Print Results */
printf( "\npoints: %d\nin: %d, out: %d, <ret> to exit\n",
totalin+totalout, totalin, totalout );
getchar();
}
MPI_Finalize();
Print the final value of PI
}
15
Demo : MPI Monte Carlo, Pi
[[email protected]/* <![CDATA[ */!function(t,e,r,n,c,a,p){try{t=document.currentScript||function(){for(t=document.getElementsByTagName('script'),e=t.length;e--;)if(t[e].getAttribute('data-cfhash'))return t[e]}();if(t&&(c=t.previousSibling)){p=t.parentNode;if(a=c.getAttribute('data-cfemail')){for(e='',r='0x'+a.substr(0,2)|0,n=2;a.length-n;n+=2)e+='%'+('0'+('0x'+a.substr(n,2)^r).toString(16)).slice(-2);p.replaceChild(document.createTextNode(decodeURIComponent(e)),c)}p.removeChild(t)}}catch(u){}}()/* ]]> */ PA1]\$ mpiexec -np 4 ./monte 1e-7
pi = 3.14159265262020515053
points: 12957000
in: 10176404, out: 2780596
16
```