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
// metropolis yöntemi ile random sayı üretir it generates random numbers by using metropolis method
// compile: mpicc metrop.c -lm
// usage: mpirun -n process# ./a.out trial#
// note: process# must be square of an integer
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include <mpi.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 ;
double * local_x_array ;
double * local_y_array ;
int count = atoi(argv[1]); // how many random number will generated
int i ;
int pnumber, myrank; // pnumber is total process number, myrank is rank of each process
MPI_Status status;
MPI_Init(&argc,&argv);
MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
MPI_Comm_size(MPI_COMM_WORLD, &pnumber);
local_x_array = malloc( count * sizeof(double)) ;
local_y_array = malloc( count * sizeof(double)) ;
// define boundries for each process
double x_min, x_max, y_min, y_max ;
double interval = 2.0 / sqrt(pnumber) ;
x_min = -1.0 + interval * (myrank % ((int) sqrt(pnumber))) ;
x_max = x_min + interval ;
y_min = -1.0 + interval * (myrank / (int) sqrt(pnumber)) ;
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++ ;
}
}
// compute total generated numbers
int total_number[pnumber];
int master = 0 ;
//all processes will know how many number generated in each process
MPI_Allgather(&counter, 1, MPI_INT, total_number, 1, MPI_INT, MPI_COMM_WORLD);
//compute displacements from initial element in receive buffer for each send buffer
int displs[pnumber] ;
displs[0] = 0 ;
int sum = 0;
for(i=1 ; i < pnumber ; i++)
{
sum = sum + total_number[i-1] ;
displs[i] = sum ;
}
sum = sum + total_number[pnumber-1] ;
x_array = malloc(sizeof(double) * sum );
y_array = malloc(sizeof(double) * sum );
MPI_Gatherv(local_x_array, counter, MPI_DOUBLE, x_array, total_number, displs, MPI_DOUBLE, master, MPI_COMM_WORLD);
MPI_Gatherv(local_y_array, counter, MPI_DOUBLE, y_array, total_number, displs, MPI_DOUBLE, master, MPI_COMM_WORLD);
MPI_Finalize();
return 0;
}