Transcript 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 Master Thread N Generate random X,Y Calculate Z=X^2+Y^2 If point lies within the circle Y Count ++ Generate random X,Y Calculate Z =X^2+Y^2 If point lies within the circle Count ++ Y N N WorkerThreads Generate random X,Y Calculate Z =X^2+Y^2 If point lies within the circle Y Count ++ N Reduction ∑ Calculate PI Print value of pi 5
OpenMP Calculating Pi
#include
#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 numbers */
srand(SEED);
Initialize random number generator; srand is used to seed the random number generated by rand()
#pragma omp parallel for private(x,y,z,tid) reduction(+:count)
for ( i=0; i x = (double)rand()/RAND_MAX; y = (double)rand()/RAND_MAX; Initialize OpenMP parallel for with reduction(∑) Randomly generate x,y points z = (x*x+y*y); if (z<=1) count++; Calculate x^2+y^2 and check if it lies within the circle; if if (i==(niter/6)-1) { yes then increment count tid = omp_get_thread_num(); printf(" thread %i just did iteration %i the count is %i\n",tid,i,count); } if (i==(niter/3)-1) { tid = omp_get_thread_num(); printf(" thread %i just did iteration %i the count is %i\n",tid,i,count); } if (i==(niter/2)-1) { tid = omp_get_thread_num(); 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 if (i==(2*niter/3)-1) { tid = omp_get_thread_num(); } printf(" thread %i just did iteration %i the count is %i\n",tid,i,count); } if (i==(5*niter/6)-1) { tid = omp_get_thread_num(); printf(" thread %i just did iteration %i the count is %i\n",tid,i,count); } if (i==niter-1) { tid = omp_get_thread_num(); 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) ); printf(" the total count is %i\n",count); pi=(double)count/niter*4; Calculate PI based on the aggregate count of 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 [LSU760000@n00 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 10 N Server Master Worker Initialize MPI Environment Initialize MPI Environment Initialize MPI Environment Broadcast Error Bound Receive Error Bound Send Request to Server Receive Request Send Request to Server Receive Random Array Receive Random Array Perform Computations Compute Random Array Perform Computations Propagate Number of Points (Allreduce) Send Array to Requestor Propagate Number of Points (Allreduce) Output Partial Result N Stop Condition Satisfied? Y Print Statistics Finalize MPI N Stop Condition Satisfied? Y Finalize MPI Last Request? Y Finalize MPI 11 #include { #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); Initialize MPI environment 12 server = numprocs-1; /* last proc is server */ if (myid == 0) sscanf( argv[1], "%lf", &epsilon ); Broadcast Error Bounds: 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 { if (request) { for (i = 0; i < CHUNKSIZE; ) { rands[i] = random(); }while( request>0 ); } Create a custom communicator MPI_Recv(&request, 1, MPI_INT, MPI_ANY_SOURCE, REQUEST, world, &status); Server process : 1. Receives request to generate a random ,2. Computes the 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); } 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 ); iter = 0; Worker process : Request the server to generate a random number array 13 } while (!done) { iter++; request = 1; /* Recv. random array from server*/ Worker : Receive random number array from the Server MPI_Recv( rands, CHUNKSIZE, MPI_INT, server, REPLY, world, &status ); for (i=0; i Worker: For each pair of x,y in the random number array, calculate the coordinates Determine if the number is inside or out of the circle } 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); Compute the value of pi and Check 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 } if (myid == 0) { /* If “Master” : Print Results */ printf( "\npoints: %d\nin: %d, out: %d, getchar(); } MPI_Finalize(); Print the final value of PI 15 [LSU760000@n00 PA1]$ mpiexec -np 4 ./monte 1e-7 pi = 3.14159265262020515053 points: 12957000 in: 10176404, out: 2780596 16Calculating Pi
Demo : OpenMP Pi
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
MPI Monte Carlo Pi Computation
Monte Carlo : MPI - Pi (source
code)
Monte Carlo : MPI - Pi (source code)
Monte Carlo : MPI - Pi (source
code)
Monte Carlo : MPI - Pi (source code)
Demo : MPI Monte Carlo, Pi