integral1d_OMP.c 2.21 KB
 atuncer @raptor committed Jan 15, 2016 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 ``````/* COMPILATION: gcc -fopenmp integral1d_OMP.c -lm USAGE: ./a.out thread# 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 i; double sum = 0; double x; int thread_count = atoi(argv[1]); // number of threads participating computation is read from the command lile arguments int trials = atoi(argv[2]); // trials is total number of random points, it is read from the command line arguments // multiple threads are started after this compiler directive, num_threads() function indicates number of threads to be started // reduction(+:sum) function indicates end of the parallel region total_hits variables of each thread are gathered in master thread // private(x) indicates variable x is private to each thread instead of being shared by threads #pragma omp parallel num_threads(thread_count) reduction(+:sum) private(x) { // It is desired every thread to generate different random numbers so // each thread are seeded with different number (thread id) srand(omp_get_thread_num()); //this compiler directive divides the loop iterations between the spawned threads #pragma omp for for(i=0 ; i < trials ; i++ ) { x = my_rand() ; sum += ( f(x) / w(x) ) ; } }// all threads are joined after this block double integral = 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))); return 0; } ``````