Newer
Older
atuncer @raptor
committed
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 <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include <omp.h>
// 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;
}