X-Git-Url: http://git.tremily.us/?a=blobdiff_plain;f=src%2Fmonte_carlo%2Fgood.c;h=041ed50eaa9460bf48567812d17a03fc39642df5;hb=a62594e22be4b17c064b5adcff8ee0240539672f;hp=d165f0b1d666b232cfcb22eddff8861cde152d91;hpb=04d4d7e462615f8f4b7d3626ba8567b4acebe3e5;p=parallel_computing.git diff --git a/src/monte_carlo/good.c b/src/monte_carlo/good.c deleted file mode 100644 index d165f0b..0000000 --- a/src/monte_carlo/good.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