Newer
Older
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
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
// metropolis yöntemi ile random sayı üretir it generates random numbers by using metropolis method
// compile: gcc -lm -fopenmp random_openmp.c
// usage: ./a.out trial# thread#
// note: process# must be square of an integer
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include <omp.h>
// w is gaussian weight function
double w(double x, double y)
{
double sigma ;
// sigma is standard deviation
sigma = 0.05 ;
return exp(-(x *x + y*y) / (2 * sigma * sigma) ) ;
}
void generator(double * x, double * y)
{
double xt, yt, ratio, tmp ;
double delta = 0.05 ;
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 ;
}
}
int main(int argc, char* argv[])
{
double x, y ;
double * x_array ;
double * y_array ;
int count = atoi(argv[1]); // how many random number will generated
int tnumber; // tnumber is total thread number, myrank is rank of each thread
tnumber = atoi(argv[2]);
int total_number[tnumber]; // this array holds how many number generated in each thread, and this is shared among all threads
#pragma omp parallel num_threads(tnumber)
{
double * local_x_array ;
double * local_y_array ;
int myrank = omp_get_thread_num();
int i ;
local_x_array = malloc( count * sizeof(double)) ;
local_y_array = malloc( count * sizeof(double)) ;
// define boundries for each thread
double x_min, x_max, y_min, y_max ;
double interval = 2.0 / sqrt(tnumber) ;
x_min = -1.0 + interval * (myrank % ((int) sqrt(tnumber))) ;
x_max = x_min + interval ;
y_min = -1.0 + interval * (myrank / (int) sqrt(tnumber)) ;
y_max = y_min + interval ;
int counter = 0 ; // counter holds number of random points in local_x_array and local_y_array
srand(myrank) ; // seed each process a different number to generate different random number
while(counter != 1)
{
x = x_min + interval * ((double) rand() / RAND_MAX) ; //assign initial random number between 0-1 to x
y = y_min + interval * ((double) rand() / RAND_MAX) ; //assign initial random number between 0-1 to y
if( (x > x_min) && (x < x_max) && (y > y_min) && (y < y_max) )
{
local_x_array[0] = x ;
local_y_array[0] = y ;
counter++ ;
}
}
for(i=1 ; i < count ; i++ )
{
generator(&x, &y);
if( (x > x_min) && (x < x_max) && (y > y_min) && (y < y_max) )
{
local_x_array[counter] = x ;
local_y_array[counter] = y ;
counter++ ;
}
}
total_number[myrank] = counter ;
// to wait until total_number array is populated place a barrier
#pragma omp barrier
if (myrank == 0)
{
int sum = 0 ; // sum holds sum of generated numbers by each thead
for(i=0 ; i < tnumber ; i++)
{
sum = sum + total_number[i] ;
}
printf("sum is %d\n", sum);
x_array = malloc(sizeof(double) * sum);
y_array = malloc(sizeof(double) * sum);
}
// to prevent threads to write before thread 0 create x_array, y_array place another barrier
#pragma omp barrier
// now each array populates global x_array, y_array by using private array local_x_array, local_y_array
// firstly, each thread computes displacement from initial element of global array
int displs = 0 ;
for(i=0 ; i < myrank ; i++)
{
displs = displs + total_number[i];
}
for(i=0 ; i < counter ; i++ )
{
x_array[displs + i] = local_x_array[i] ;
y_array[displs + i] = local_y_array[i] ;
}
} //end pragma
return 0;
}