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
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
// calculate integral using client-server model
// compile: gcc -fopenmp integral2_openmp.c -lm
// usage: ./a.out trial# thread#
// note: trial# must be multiple of CHUNK_SIZE
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include <omp.h>
#define CHUNK_SIZE 1000
// 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 generated = 0 ;
int consumed = 0 ;
#pragma omp parallel num_threads(tnumber)
{
int myrank = omp_get_thread_num();
int i, j ;
int start_index;
double local_sum = 0 ;
if(myrank == 0 ) //if master thread 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(generated=0 ; generated < count ; generated++)
{
x_array[generated] = x ;
y_array[generated] = y ;
generator(&x, &y);
}
}
else
{
while( 1 )
{
# pragma omp critical
{
if( (generated - consumed) >= CHUNK_SIZE )
{
start_index = consumed ;
consumed = consumed + CHUNK_SIZE ;
for( i=0 ; i<CHUNK_SIZE ; i++)
{
local_sum = local_sum + f(x_array[start_index+i], y_array[start_index+i]) ;
}
sum = sum + local_sum ;
local_sum = 0;
}
} // end pragma omp critical
if ( (generated == consumed) && (generated > 0) )
{
break ;
}
}
}
#pragma omp barrier
} //end pragma
printf("%f\n", sum / count);
return 0;
}