// calculate integral using master-workers model // compile: mpicc integral1_mpi.c -lm // usage: mpirun -n process# ./a.out trial# #include #include #include #include #include // w is gaussian weight function double w(double x, double y) { return exp( - ( x*x + y*y ) ); } void generator(double * x, double * y) { double xt, yt, ratio, tmp ; double delta = 0.5 ; tmp = (double)rand() / (double)RAND_MAX ; // update x by adding a multiples of delta , (2 * tmp - 1) creates a number between -1.0 and 1.0 xt = (*x) + delta * (2 * tmp - 1) ; tmp = (double)rand() / (double)RAND_MAX ; // update y by adding a multiples of delta , (2 * tmp - 1) creates a number between -1.0 and 1.0 yt = (*y) + delta * (2 * tmp - 1) ; // compare updated x,y values with old x,y values, accept or reject updated values as new values according to weight function ratio = w(xt, yt ) / w(*x, *y) ; tmp = (double)rand() / (double)RAND_MAX ; if(ratio > tmp ) { *x = xt ; *y = yt ; } } //f is the function to be integrated double f(double x, double y) { return M_PI * ( x*x + y*y ) ; } int main(int argc, char* argv[]) { double x, y ; double * x_array ; double * y_array ; double sum ; double total_sum ; double integral ; int trials, trials_per_p; int i, j ; int pnumber, myrank; MPI_Status status; MPI_Init(&argc,&argv); MPI_Comm_rank(MPI_COMM_WORLD, &myrank); MPI_Comm_size(MPI_COMM_WORLD, &pnumber); x = 0.0 ; y = 0.0 ; trials = atoi(argv[1]); trials_per_p = trials / (pnumber - 1 ); x_array = malloc( trials_per_p * sizeof(double)) ; y_array = malloc( trials_per_p * sizeof(double)) ; if(myrank == 0) //master process generates random numbers and sends them to workers { for(j=1 ; j < pnumber ; j++) { for(i=0 ; i < trials_per_p ; i++) { generator(&x, &y); x_array[i] = x ; y_array[i] = y ; } MPI_Send(x_array, trials_per_p , MPI_DOUBLE, j, 1234 , MPI_COMM_WORLD) ; MPI_Send(y_array, trials_per_p , MPI_DOUBLE, j, 5678 , MPI_COMM_WORLD) ; } } if(myrank != 0) // worker processes receive random numbers from master process and calculate sum { MPI_Recv(x_array, trials_per_p, MPI_DOUBLE, 0, 1234, MPI_COMM_WORLD, &status) ; MPI_Recv(y_array, trials_per_p, MPI_DOUBLE, 0, 5678, MPI_COMM_WORLD, &status) ; sum = 0.0 ; for(i=0 ; i < trials_per_p ; i++) { sum += f(x_array[i], y_array[i]) ; } } MPI_Reduce(&sum, &total_sum, 1, MPI_DOUBLE , MPI_SUM, 0, MPI_COMM_WORLD) ; // all calculated sums are accumulated in master process if(myrank == 0) { integral = total_sum / ( (pnumber -1) * trials_per_p ) ; printf("%f %f\n", integral , M_PI); } MPI_Finalize(); return 0; }