From: W. Trevor King Date: Tue, 14 Sep 2010 19:24:19 +0000 (-0400) Subject: Merged src/monte_carlo/transfer/ into src/monte_carlo/. X-Git-Url: http://git.tremily.us/?a=commitdiff_plain;h=4f45f333e9590117eea2d8ac9f6c2295188227f8;p=parallel_computing.git Merged src/monte_carlo/transfer/ into src/monte_carlo/. --- diff --git a/src/monte_carlo/transfer/README b/src/monte_carlo/README similarity index 100% rename from src/monte_carlo/transfer/README rename to src/monte_carlo/README diff --git a/src/monte_carlo/transfer/compile_all b/src/monte_carlo/compile_all similarity index 100% rename from src/monte_carlo/transfer/compile_all rename to src/monte_carlo/compile_all diff --git a/src/monte_carlo/transfer/Metropolis_CS.c b/src/monte_carlo/transfer/Metropolis_CS.c deleted file mode 100644 index ef122fe..0000000 --- a/src/monte_carlo/transfer/Metropolis_CS.c +++ /dev/null @@ -1,347 +0,0 @@ - -/* Metropolis_CS.c */ -/* Demonstration of integration by random numbers */ -/* Client-Server */ -/* */ -/* Metropolis Algorithm */ -/* Gaussian Random Numbers distribution */ -/* */ - -/* Syntax: */ -/* mpiexec -np 5 ./Metropolis_CS -p > data_file */ - - /* Michel Vallieres */ - -#include -#include -#include -#include - -#define NMAX 1000 /* trial size of randomarrays x & y */ - - -/* function to integrate */ -double func( double x, double y ) -{ - return M_PI * ( x*x + y*y ); -} - - - -/* probability profile of the random numbers */ -double weight( double x, double y ) -{ - return exp( - ( x*x + y*y ) ); -} - - - -/* Metropolis algorithm to generate Random */ -/* Points in the plane with probability profile */ -/* weight(x,y) */ -void Metropolis ( double x[], double y[], double (*weight)(double,double), - double delta, int N, int *accepted_count ) -{ - double xt, yt, xx, yy, eta, ratio; - int local_count, count; - - count = 0; /* accepted count */ - - for( local_count=0 ; local_count < N-1 ; local_count++ ) - { - xx = (double)rand()/(double)RAND_MAX; /* random # in 0,1 */ - yy = (double)rand()/(double)RAND_MAX; - xt = x[local_count] + delta * ( 2*xx - 1.0 ); /* in range -1 to 1 */ - yt = y[local_count] + delta * ( 2*yy - 1.0 ); /* trial step */ - /* probability ratio */ - x[local_count+1] = x[local_count]; /* new step = old step */ - y[local_count+1] = y[local_count]; - ratio = weight( xt, yt ) / weight( x[local_count] , y[local_count] ); - eta = (double)rand()/(double)RAND_MAX; /* random # in 0-1 */ - if ( ratio > eta ) /* Metropolis selection */ - { - count++; - x[local_count+1] = xt; /* selected step */ - y[local_count+1] = yt; - } - - } - *accepted_count = *accepted_count + count; -} - - - - -/* code executed by node 0 */ -/* control of calculation & final Monte Carlo sums */ -void node_zero_code( int numprocs ) -{ - double f, sumf, sumf2, integral, sigma, store[2]; - int N_average, N_input, leftover, N_node, N_total, count, first; - int node, kill; - double x0, y0; - MPI_Status recv_status; - - x0 = 0.5; /* initial position in plane*/ - y0 = 0.6; - /* transmit to server node */ - MPI_Ssend( &x0, 1, MPI_DOUBLE, numprocs-1, 2123, MPI_COMM_WORLD ); - MPI_Ssend( &y0, 1, MPI_DOUBLE, numprocs-1, 3123, MPI_COMM_WORLD ); - - f = func( x0, y0 ); /* Monte-Carlo sums */ - sumf = f; - sumf2 = f*f; - - N_total = 1; /* Monte-Carlo counters */ - first = 1; - - /* accuracy improvement */ - fprintf( stderr, "\n Enter # of random numbers to use ( 0 to quit ): " ); - while ( scanf( "%d", &N_input ) != EOF && N_input > 0 ) - { - /* adjust N_input to # nodes */ - N_average = N_input/(numprocs-2); - leftover = N_input % (numprocs-2); - /* add to total */ - N_total = N_total + N_input - first; - /* send to slave nodes */ - for ( node=1 ; node ~0.5 */ - fprintf( stderr, "\n Acceptance ratio: %d %d %f\n\n", - N_total, count, (double)count/(double)N_total ); - - /* end all */ - MPI_Finalize(); - exit(0); - -} - - - - -/* worker code: requests & receives random numbers and */ -/* calculates MC partial sums */ -void worker_nodes_code( int numprocs ) -{ - static double *x, *y; /* static: permanent storage */ - static int current_x_y_size; - double f, sumf, sumf2, store[2]; - int N, ir; - MPI_Status recv_status; - /* current size of x & y */ - /* arrays - reserve memory */ - current_x_y_size = NMAX; - x = (double *)malloc( current_x_y_size * sizeof(double) ); - y = (double *)malloc( current_x_y_size * sizeof(double) ); - - /* wait for # of MC steps */ - N = 1; - while ( N > 0 ) - { - /* receive N from node 0 */ - MPI_Recv( &N, 1, MPI_INT, 0, 51, MPI_COMM_WORLD, &recv_status ); - - /* suicide point */ - if ( N < 0 ) - { - MPI_Finalize(); - exit(0); - } - /* request MC steps */ - MPI_Ssend( &N, 1, MPI_INT, numprocs-1, 123, MPI_COMM_WORLD ); - - /* space for x & y arrays */ - if ( N > current_x_y_size || - current_x_y_size == 0 ) /* enough ? */ - { /* reserve more */ - free ( x ); - free ( y ); - current_x_y_size = N + 100; - x = (double *)malloc( current_x_y_size * sizeof(double) ); - y = (double *)malloc( current_x_y_size * sizeof(double) ); - } - - /* receive random #s */ - MPI_Recv( x, N, MPI_DOUBLE, numprocs-1, 201, - MPI_COMM_WORLD, &recv_status ); - MPI_Recv( y, N, MPI_DOUBLE, numprocs-1, 202, - MPI_COMM_WORLD, &recv_status ); - - /* Monte-Carlo integral */ - sumf = 0.0; - sumf2 = 0.0; - for ( ir=1 ; ir 0 ) /* loop over requests */ - /* from nodes */ - { - /* receive request */ - MPI_Recv( &N, 1, MPI_INT, MPI_ANY_SOURCE, 123, - MPI_COMM_WORLD, &recv_status ); - - /* in case of early suicide */ - if ( N <= 0 ) break; - - /* space for x & y arrays */ - if ( N > current_x_y_size || - current_x_y_size == 0 ) /* enough ? */ - { /* reserve more */ - free ( x ); - free ( y ); - current_x_y_size = N + 100; - x = (double *)malloc( current_x_y_size * sizeof(double) ); - y = (double *)malloc( current_x_y_size * sizeof(double) ); - } - - /* Metropolis Algorithm */ - x[0] = x00; - y[0] = y00; - Metropolis ( x, y, weight, delta, N , &accepted_count ); - - /* plot info */ - if ( plot_control == 1 ) - for ( ir=1 ; ir 0 && myid < numprocs-1 ) - worker_nodes_code( numprocs ); - /***** SERVER ****/ - else - { - /* plot control */ - plot_control = 0; /* only server needs it */ - if ( argc == 2 ) - if ( strcmp( argv[1], "-p") == 0 ) - plot_control = 1; - /* call the server code */ - server_code( plot_control ); - } - -} - diff --git a/src/monte_carlo/transfer/Metropolis_Gaussian_Client_Server.c b/src/monte_carlo/transfer/Metropolis_Gaussian_Client_Server.c deleted file mode 100644 index d165f0b..0000000 --- a/src/monte_carlo/transfer/Metropolis_Gaussian_Client_Server.c +++ /dev/null @@ -1,280 +0,0 @@ - -/* Metropolis_Gaussian_Client_Server.c */ -/* Demonstration of integration by random numbers */ -/* Metropolis Algorithm */ -/* Gaussian Random Numbers distribution */ -/* */ - -/* Syntax: */ -/* Metropolis_Gaussian_Client_Server -p > data_file */ -/* */ - - /* Michel Vallieres */ - -#include -#include -#include -#include - -#define NMAX 1000 - -/* function to integrate */ -double func( double x, double y ) -{ - return M_PI * ( x*x + y*y ); -} - - - -/* probability profile of the random numbers */ -double weight( double x, double y ) -{ - return exp( - ( x*x + y*y ) ); -} - - -/* Metropolis algorithm to generate Random */ -/* Points in the plane with probability profile */ -/* weight(x,y) */ -void Metropolis ( double x[], double y[], double (*weight)(double,double), - double delta, int N, int *accepted_count ) -{ - double xt, yt, xx, yy, eta, ratio; - int local_count, count; - - count = 0; /* accepted count */ - - for( local_count=0 ; local_count < N-1 ; local_count++ ) - { - xx = (double)rand()/(double)RAND_MAX; /* random # 0,1 */ - yy = (double)rand()/(double)RAND_MAX; - xt = x[local_count] + delta * ( 2*xx - 1.0 ); /* in range -1 to 1 */ - yt = y[local_count] + delta * ( 2*yy - 1.0 ); /* trial step */ - /* probability ratio */ - x[local_count+1] = x[local_count]; /* new step = old step */ - y[local_count+1] = y[local_count]; - ratio = weight( xt, yt ) / weight( x[local_count] , y[local_count] ); - eta = (double)rand()/(double)RAND_MAX; /* random # in 0-1 */ - if ( ratio > eta ) /* Metropolis selection */ - { - count++; - x[local_count+1] = xt; /* selected step */ - y[local_count+1] = yt; - } - - // fprintf(stderr,"%d %f %f %f %f %d\n", local_count,xt,yt,ratio,eta,count); - - } - *accepted_count = *accepted_count + count; -} - - -/* demo program to illustrate Monte-Carlo integration using */ -/* Gaussian random numbers in the plane */ -int main( int argc, char *argv[] ) -{ - int ir, accepted_count, N, N_total, N_input, - N_node, N_average, leftover, first; - int current_x_y_size; - double *x, *y, delta, x0, y0; - int plot_control; - double f, sumf, sumf2, integral, sigma, store[2]; - MPI_Status recv_status; - int myid, numprocs, node, to; - - /* join the MPI virtual machine */ - MPI_Init(&argc, &argv); - MPI_Comm_rank(MPI_COMM_WORLD, &myid); - MPI_Comm_size(MPI_COMM_WORLD, &numprocs); - - if ( myid != 0 ) - { - /* current size of x & y */ - /* arrays - reserve memory */ - current_x_y_size = NMAX; - x = (double *)malloc( current_x_y_size * sizeof(double) ); - y = (double *)malloc( current_x_y_size * sizeof(double) ); - } - - x0 = 0.5; /* initial plane position */ - y0 = 0.6; - - accepted_count = 1; /* counter - selected numbers */ - - - /**** MASTER NODE ****/ - if ( myid == 0 ) - { - f = func( x0, y0 ); /* Monte-Carlo sums */ - sumf = f; - sumf2 = f*f; - - first = 1; - N_total = 1; - /* accuracy improvement */ - fprintf( stderr, "\n Enter # of random numbers to use ( 0 to quit ): " ); - while ( scanf( "%d", &N_input ) != EOF && N_input > 0 ) - { - /* adjust N_input to # nodes */ - N_average = N_input/(numprocs-2); - leftover = N_input % (numprocs-2); - /* add to total */ - N_total = N_total + N_input - first; - /* send to slave nodes */ - for ( node=1 ; node ~0.5 */ - fprintf( stderr, "\n Acceptance ratio: %d %d %f\n\n", - N_total, accepted_count, (double)accepted_count/(double)N_total ); - - /* end all */ - MPI_Finalize(); - exit(0); - - } - - /**** SLAVE NODES ****/ - else if ( myid > 0 && myid < numprocs-1 ) - { - /* wait for # of random numbers */ - N = 1; - while ( N > 0 ) - { - /* receive N from node 0 */ - MPI_Recv( &N, 1, MPI_INT, 0, 51, MPI_COMM_WORLD, &recv_status ); - - /* suicide point */ - if ( N < 0 ) - { - MPI_Finalize(); - exit(0); - } - /* request random # from server */ - MPI_Ssend( &N, 1, MPI_INT, numprocs-1, 123, MPI_COMM_WORLD ); - - /* space for x & y arrays */ - if ( N > current_x_y_size || - current_x_y_size == 0 ) /* enough ? */ - { /* reserve more */ - free ( x ); - free ( y ); - current_x_y_size = N + 100; - x = (double *)malloc( current_x_y_size * sizeof(double) ); - y = (double *)malloc( current_x_y_size * sizeof(double) ); - } - - /* receive random #s */ - MPI_Recv( x, N, MPI_DOUBLE, numprocs-1, 201, MPI_COMM_WORLD, &recv_status ); - MPI_Recv( y, N, MPI_DOUBLE, numprocs-1, 202, MPI_COMM_WORLD, &recv_status ); - - /* Monte-Carlo integral */ - sumf = 0.0; - sumf2 = 0.0; - for ( ir=1 ; ir 0 ) /* loop over requests from nodes */ - { - /* receive request */ - MPI_Recv( &N, 1, MPI_INT, MPI_ANY_SOURCE, 123, - MPI_COMM_WORLD, &recv_status ); - - /* in case of early suicide */ - if ( N <= 0 ) break; - - /* space for x & y arrays */ - if ( N > current_x_y_size || - current_x_y_size == 0 ) /* enough ? */ - { /* reserve more */ - free ( x ); - free ( y ); - current_x_y_size = N + 100; - x = (double *)malloc( current_x_y_size * sizeof(double) ); - y = (double *)malloc( current_x_y_size * sizeof(double) ); - } - - /* Metropolis Algorithm */ - x[0] = x0; - y[0] = y0; - Metropolis ( x, y, weight, delta, N , &accepted_count ); - - /* plot info */ - if ( plot_control == 1 ) - for ( ir=1 ; ir data_file */ - - /* Michel Vallieres */ - -#include -#include -#include - -#define NMAX 1000 - - -/* function to integrate */ -double func( double x, double y ) -{ - return M_PI * ( x*x + y*y ); -} - - - -/* probability profile of the random numbers */ -double weight( double x, double y ) -{ - return exp( - ( x*x + y*y ) ); -} - - -/* Metropolis algorithm to generate Random */ -/* Points in the plane with probability profile */ -/* weight(x,y) */ -void Metropolis ( double x[], double y[], double (*weight)(double,double), - double delta, int N, int *accepted_count ) -{ - double xt, yt, xx, yy, eta, ratio; - int local_count, count; - - count = 0; /* accepted count */ - - for( local_count=0 ; local_count < N-1 ; local_count++ ) - { - xx = (double)rand()/(double)RAND_MAX; /* random # in 0,1 */ - yy = (double)rand()/(double)RAND_MAX; - xt = x[local_count] + delta * ( 2*xx - 1.0 ); /* in range -1 to 1 */ - yt = y[local_count] + delta * ( 2*yy - 1.0 ); /* trial step */ - /* probability ratio */ - x[local_count+1] = x[local_count]; /* new step = old step */ - y[local_count+1] = y[local_count]; - ratio = weight( xt, yt ) / weight( x[local_count] , y[local_count] ); - eta = (double)rand()/(double)RAND_MAX; /* random # in 0-1 */ - if ( ratio > eta ) /* Metropolis selection */ - { - count++; - x[local_count+1] = xt; /* selected step */ - y[local_count+1] = yt; - } - - // fprintf(stderr,"<-> %d %f %f %f %f %d\n", local_count,xt,yt,ratio,eta,count); - - } - *accepted_count = *accepted_count + count; -} - - -/* demo program to illustrate Monte-Carlo integration using */ -/* Gaussian random numbers in the plane */ -int main( int argc, char *argv[] ) -{ - int ir, accepted_count, N_total, N_input, first; - double delta, x00, y00; - int plot_control, current_x_y_size; - double f, sumf, sumf2, integral, sigma; - double *x, *y; - /* plot control */ - plot_control = 0; - if ( argc == 2 ) - if ( strcmp( argv[1], "-p") == 0 ) - plot_control = 1; - - delta = 0.2; /* Metropolis delta (size of jump) */ - - N_total = 1; /* total # of random numbers */ - - accepted_count = 1; /* accepted steps counter */ - - first = 1; /* first bunch of MC steps */ - - /* space for the MC steps */ - current_x_y_size = NMAX + 100; - x = (double *)malloc( current_x_y_size * sizeof(double) ); - y = (double *)malloc( current_x_y_size * sizeof(double) ); - - x00 = 0.5; /* initial position in plane */ - y00 = 0.6; - - f = func( x00, y00 ); /* Monte-Carlo sums */ - sumf = f; - sumf2 = f*f; - /* accuracy improvement */ - fprintf( stderr, "\n Enter # of random number to use ( 0 to quit ): " ); - while ( scanf( "%d", &N_input ) != EOF && N_input > 0 ) - { - N_total = N_total + N_input - first; /* add to total */ - - if ( first == 0 ) N_input++; /* adjust count on MC steps */ - first = 0; /* on subsequent calls */ - - if ( N_input > current_x_y_size ) /* enough space? */ - { /* if not -- reserve more */ - free ( x ); - free ( y ); - current_x_y_size = N_input + 100; - x = (double *)malloc( current_x_y_size * sizeof(double) ); - y = (double *)malloc( current_x_y_size * sizeof(double) ); - } - - x[0] = x00; /* initial step */ - y[0] = y00; - /* Metropolis Algorithm */ - Metropolis ( x, y, weight, delta, N_input, &accepted_count ); - - /* plot info */ - if ( plot_control == 1 ) - for ( ir=0 ; ir ~0.5 */ - fprintf( stderr, "\n Acceptance ratio: %d %d %f\n\n", - N_total, accepted_count, (double)accepted_count/(double)N_total ); - -} diff --git a/src/monte_carlo/transfer/Metropolis_Gaussian_Parallel.c b/src/monte_carlo/transfer/Metropolis_Gaussian_Parallel.c deleted file mode 100644 index 5dce7c6..0000000 --- a/src/monte_carlo/transfer/Metropolis_Gaussian_Parallel.c +++ /dev/null @@ -1,279 +0,0 @@ - -/* Metropolis_Gaussian_Parallel.c */ -/* Demonstration of integration by random numbers */ -/* Metropolis Algorithm */ -/* Gaussian Random Numbers distribution */ -/* Parallel version */ -/* */ - -/* Syntax: */ -/* mpiexec -np N Metropolis_Gaussian_Parallel -p > data_file */ - - /* Michel Vallieres */ - -#include -#include -#include -#include - -#define NMAX 1000 /* initial array dimensions */ - /* will be modified if need be */ - -double x00, y00; /* initial position of random walk */ - - -/* function to integrate */ -double func( double x, double y ) -{ - return M_PI * ( x*x + y*y ); -} - - - -/* probability profile of the random numbers */ -double weight( double x, double y ) -{ - return exp( - ( x*x + y*y ) ); -} - - - -/* Metropolis algorithm to generate Random */ -/* Points in the plane with probability profile */ -/* weight(x,y) */ -void Metropolis ( double x[], double y[], double (*weight)(double,double), - double delta, int N, int *accepted_count ) -{ - double xt, yt, xx, yy, eta, ratio; - int local_count, count; - - count = 0; /* accepted count */ - - for( local_count=0 ; local_count < N-1 ; local_count++ ) - { - xx = (double)rand()/(double)RAND_MAX; /* random # in 0,1 */ - yy = (double)rand()/(double)RAND_MAX; - xt = x[local_count] + delta * ( 2*xx - 1.0 ); /* in range -1 to 1 */ - yt = y[local_count] + delta * ( 2*yy - 1.0 ); /* trial step */ - /* probability ratio */ - x[local_count+1] = x[local_count]; /* new step = old step */ - y[local_count+1] = y[local_count]; - ratio = weight( xt, yt ) / weight( x[local_count] , y[local_count] ); - eta = (double)rand()/(double)RAND_MAX; /* random # in 0-1 */ - if ( ratio > eta ) /* Metropolis selection */ - { - count++; - x[local_count+1] = xt; /* selected step */ - y[local_count+1] = yt; - } - - // fprintf(stderr,"%d %f %f %f %f %d %d\n", - // local_count,xt,yt,ratio,eta,N,count); - - } - *accepted_count = *accepted_count + count; -} - - - -/* demo program to illustrate Monte-Carlo integration using */ -/* Gaussian random numbers in the plane */ -/* Parallel version */ -int main( int argc, char *argv[] ) -{ - int round, ir, ir_min, accepted_count, N, N_total, N_input; - int plot_control, first, ir_first; - double *x, *y, *x_local, *y_local; - double delta; - double f, sumf, sumf2, local_sumf, local_sumf2, - total_sumf, total_sumf2, integral, sigma; - int node, myid, numprocs; - int current_x_y_size, N_average, leftover, number, nmax; - MPI_Status recv_status; - - /* join the MPI virtual machine */ - MPI_Init(&argc, &argv); - MPI_Comm_rank(MPI_COMM_WORLD, &myid); - MPI_Comm_size(MPI_COMM_WORLD, &numprocs); - - /* adjust NMAX for # nodes */ - nmax = NMAX - ( NMAX % ( numprocs-1) ); - /* current size of x & y */ - /* arrays - reserve memory */ - current_x_y_size = nmax; - x = (double *)malloc( current_x_y_size * sizeof(double) ); - y = (double *)malloc( current_x_y_size * sizeof(double) ); - x_local = (double *)malloc( current_x_y_size * sizeof(double) ); - y_local = (double *)malloc( current_x_y_size * sizeof(double) ); - - x00 = 0.5; /* initial position in plane */ - y00 = 0.6; - - /* MASTER NODE */ - if ( myid == 0 ) - { - /* plot control */ - plot_control = 0; - if ( argc == 2 ) - if ( strcmp( argv[1], "-p") == 0 ) - plot_control = 1; - - delta = 0.2; /* Metropolis delta - size of jump */ - - N_total = 1; /* total # of random numbers */ - accepted_count = 1; /* counter for selected numbers */ - first = 1; - - f = func(x00,y00); /* Monte-Carlo sums */ - sumf = f; - sumf2 = f*f ; - - /* accuracy improvement */ - fprintf( stderr, "\n Enter # of random number to use ( 0 to quit ): " ); - while ( scanf( "%d", &N_input ) != EOF && N_input > 0 ) - { - - /* find # per processor */ - N_average = N_input / ( numprocs -1 ); - leftover = N_input % ( numprocs -1 ); - /* add to total */ - N_total = N_total + N_input - first; - /* generate & send random */ - /* numbers to each node */ - for ( node=1 ; node current_x_y_size ) - { - free( x ); - free( y ); - current_x_y_size = number; - x = (double *)malloc( (current_x_y_size+100) * sizeof(double) ); - y = (double *)malloc( (current_x_y_size+100) * sizeof(double) ); - } - /* Metropolis Algorithm */ - x[0] = x00; - y[0] = y00; - Metropolis ( x, y, weight, delta, number , &accepted_count ); - - /* send x & y to nodes */ - /* send x & y to nodes */ - MPI_Ssend( x, number, MPI_DOUBLE, node, 123, MPI_COMM_WORLD ); - MPI_Ssend( y, number, MPI_DOUBLE, node, 124, MPI_COMM_WORLD ); - - /* plot info */ - if ( plot_control == 1 ) - { - ir_first = 0; - if ( first == 0 ) ir_first = 1; - for ( ir=ir_first ; ir nodes to suicide */ - number = -1; - for ( node=1 ; node ~0.5 */ - fprintf( stderr, "\n Acceptance ratio: %d %d %f\n\n", - N_total, accepted_count, (double)accepted_count/(double)N_total ); - - /* end main code */ - free( x ); - free( y ); - MPI_Finalize(); - exit(0); - - } /* if my_node */ - /* WORKER NODES */ - else - { - /* loop over all bunches of */ - /* random numbers */ - for ( ; ; ) - { - /* receive signal or # of numbers */ - MPI_Recv( &number, 1, MPI_INT, 0, 121, - MPI_COMM_WORLD, &recv_status); - - /* end code - suicide */ - if ( number <0 ) - { - free( x_local ); - free( y_local ); - MPI_Finalize(); - exit(0); - } - else if ( number > 0 ) /* perform MC integral */ - { - /* more space if needed */ - if ( number > current_x_y_size ) - { - free( x_local ); - free( y_local ); - current_x_y_size = number; - x_local = (double *)malloc( current_x_y_size * sizeof(double) ); - y_local = (double *)malloc( current_x_y_size * sizeof(double) ); - } - - MPI_Recv( x_local, number, MPI_DOUBLE, 0, 123, - MPI_COMM_WORLD, &recv_status); - MPI_Recv( y_local, number, MPI_DOUBLE, 0, 124, - MPI_COMM_WORLD, &recv_status); - /* Monte-Carlo integral */ - local_sumf = 0.0; - local_sumf2 = 0.0; - for ( ir=1 ; ir data_file */ - - /* Michel Vallieres */ - -#include -#include -#include - -#define NMAX 1000 - - -/* function to integrate */ -double func( double x, double y ) -{ - return M_PI * ( x*x + y*y ); -} - - - -/* probability profile of the random numbers */ -double weight( double x, double y ) -{ - return exp( - ( x*x + y*y ) ); -} - - -/* Metropolis algorithm to generate Random */ -/* Points in the plane with probability profile */ -/* weight(x,y) */ -void Metropolis ( double x[], double y[], double (*weight)(double,double), - double delta, int N, int *accepted_count ) -{ - double xt, yt, xx, yy, eta, ratio; - int local_count, count; - - count = 0; /* accepted count */ - - for( local_count=0 ; local_count < N-1 ; local_count++ ) - { - xx = (double)rand()/(double)RAND_MAX; /* random # in 0,1 */ - yy = (double)rand()/(double)RAND_MAX; - xt = x[local_count] + delta * ( 2*xx - 1.0 ); /* in range -1 to 1 */ - yt = y[local_count] + delta * ( 2*yy - 1.0 ); /* trial step */ - /* probability ratio */ - x[local_count+1] = x[local_count]; /* new step = old step */ - y[local_count+1] = y[local_count]; - ratio = weight( xt, yt ) / weight( x[local_count] , y[local_count] ); - eta = (double)rand()/(double)RAND_MAX; /* random # in 0-1 */ - if ( ratio > eta ) /* Metropolis selection */ - { - count++; - x[local_count+1] = xt; /* selected step */ - y[local_count+1] = yt; - } - - // fprintf(stderr,"<-> %d %f %f %f %f %d\n", local_count,xt,yt,ratio,eta,count); - - } - *accepted_count = *accepted_count + count; -} - - -/* demo program to illustrate Monte-Carlo integration using */ -/* Gaussian random numbers in the plane */ -int main( int argc, char *argv[] ) -{ - int ir, accepted_count, N_total, N_input, first; - double delta, x00, y00; - int plot_control, current_x_y_size; - double f, sumf, sumf2, integral, sigma; - double *x, *y; - /* plot control */ - plot_control = 0; - if ( argc == 2 ) - if ( strcmp( argv[1], "-p") == 0 ) - plot_control = 1; - - delta = 0.2; /* Metropolis delta (size of jump) */ - - N_total = 1; /* total # of random numbers */ - - accepted_count = 1; /* accepted steps counter */ - - first = 1; /* first bunch of MC steps */ - - /* space for the MC steps */ - current_x_y_size = NMAX + 100; - x = (double *)malloc( current_x_y_size * sizeof(double) ); - y = (double *)malloc( current_x_y_size * sizeof(double) ); - - x00 = 0.5; /* initial position in plane */ - y00 = 0.6; - - f = func( x00, y00 ); /* Monte-Carlo sums */ - sumf = f; - sumf2 = f*f; - /* accuracy improvement */ - fprintf( stderr, "\n Enter # of random number to use ( 0 to quit ): " ); - while ( scanf( "%d", &N_input ) != EOF && N_input > 0 ) - { - N_total = N_total + N_input - first; /* add to total */ - - if ( first == 0 ) N_input++; /* adjust count on MC steps */ - first = 0; /* on subsequent calls */ - - if ( N_input > current_x_y_size ) /* enough space? */ - { /* if not -- reserve more */ - free ( x ); - free ( y ); - current_x_y_size = N_input + 100; - x = (double *)malloc( current_x_y_size * sizeof(double) ); - y = (double *)malloc( current_x_y_size * sizeof(double) ); - } - - x[0] = x00; /* initial step */ - y[0] = y00; - /* Metropolis Algorithm */ - Metropolis ( x, y, weight, delta, N_input, &accepted_count ); - - /* plot info */ - if ( plot_control == 1 ) - for ( ir=0 ; ir ~0.5 */ - fprintf( stderr, "\n Acceptance ratio: %d %d %f\n\n", - N_total, accepted_count, (double)accepted_count/(double)N_total ); - -}