/* COMPILATION: mpicc integral1d_mpi.c -lm USAGE: mpirun -n processor# ./a.out trial# This program computes integral of ( 1 / (1 + x*x) ) in x=[0,1] by using importance sampling Monte Carlo */ #include #include #include #include #include // f() is the function to be integrated double f(double x) { return (1.0 / (1.0 + x*x) ) ; } // random x values are generated according to probability distribution function w() double w(double x) { return ((1.0/3.0) * (4 - 2 * x)) ; } // my_rand() function generates random numbers according to p.d.f w() // this function uses inverse of the p.d.f w() double my_rand() { double y = (double)rand() / (double)RAND_MAX ; return (2.0 - sqrt(4.0 - 3.0*y) ) ; } int main(int argc, char* argv[]) { int trials, trials_per_p; // trials is total number of random points and trials_per_p number of random points for each processor int i; double partial_sum = 0; // each processor uses its variable partial_sum to sum function value double total_sum; // at the end of the computation all processors' partial_sum gathered in total_sum double x; int pnumber, myrank; // pnumber is the number of processor in the computation and myrank is the rank of the processor MPI_Init(&argc,&argv); // Initialize the MPI execution environment MPI_Comm_rank(MPI_COMM_WORLD, &myrank); // Each processor learns its rank number MPI_Comm_size(MPI_COMM_WORLD, &pnumber); // Each processor learns how many processor participate to the calculation trials = atoi(argv[1]); // command line argument is readed for total random number trials_per_p = trials / pnumber ; // number of points for each processor is calculated srand(myrank); // It is desired every processor to generate different random numbers so // each processor seeded with different number (rank oprocessor) for(i=0 ; i < trials_per_p ; i++ ) { x = my_rand() ; partial_sum += ( f(x) / w(x) ) ; } // Assume the processor which has rank=0 is master int master_p = 0; //total sum and total random number are gathered to the master processor MPI_Reduce(&partial_sum, &total_sum, 1, MPI_DOUBLE, MPI_SUM, master_p, MPI_COMM_WORLD); MPI_Reduce(&trials, &trials_per_p, 1, MPI_INT, MPI_SUM, master_p, MPI_COMM_WORLD); //Master processor prints the result if(myrank == master_p) { double integral = total_sum / trials ; printf("computed value of the integral= %f \n", integral); printf("true value of the integral= %f \n", atan(1.0)); printf("error= %f \n", fabs( integral - atan(1.0))); } MPI_Finalize(); //Terminates MPI execution environment return 0; }