integral1_openmp.c 3.03 KB
 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 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 ``````// calculate integral using master-workers model // compile: gcc -fopenmp integral1_openmp.c -lm // usage: ./a.out trial# thread# #include #include #include #include #include // w is 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 = 0.0 ; int count = atoi(argv[1]); // how many random number will generated int tnumber = atoi(argv[2]); // tnumber is total thread number, myrank is rank of each thread x_array = malloc(sizeof(double) * count); y_array = malloc(sizeof(double) * count); int ready_thread[tnumber] ; //this array holds which thread is ready to run #pragma omp parallel num_threads(tnumber) reduction ( +:sum ) { int myrank = omp_get_thread_num(); int i, j ; int number_for_each = count / (tnumber-1) ; // how many random number is used in each thread, master thread is not counted if(myrank == 0) //if thread no=0 then generate random numbers, else use generated random numbers { x = (-1.0) + 2 * ((double) rand() / RAND_MAX) ; // assign initial random number between -1 and 1 y = (-1.0) + 2 * ((double) rand() / RAND_MAX) ; // assign initial random number between -1 and 1 for(j=1 ; j < tnumber ; j++) { for(i=0 ; i < count ; i++ ) { x_array[i] = x ; y_array[i] = y ; generator(&x, &y); } ready_thread[ j ] = 1; } } else { int start_index = (myrank - 1) * number_for_each; while(ready_thread[myrank] != 1) { //hold thread until ready bit becomes 1 } for(i=0 ; i < number_for_each ; i++ ) { sum += f( x_array[start_index+i] , y_array[start_index+i] ) ; } } #pragma omp barrier } //end pragma printf("%f\n", sum / count ); return 0; }``````