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
| double MonteCarloSequential(
std::vector<double> a, double sd,
double (*pfunction)(std::vector<double>), bool (*parea)(std::vector<double>),
unsigned int dim, int p_count) {
std::mt19937 gen;
gen.seed(static_cast<unsigned int>(time(0)));
std::uniform_real_distribution<> urd(0, 1);
std::vector<double> p(dim * p_count);
for (unsigned int i = 0; i < p.size(); i++) {
p[i] = a[i % dim] + urd(gen) * sd;
}
int n = 0;
double val = 0;
for (int i = 0; i < p_count; ++i) {
std::vector<double> cpoint;
for (unsigned int j = 0; j < dim; ++j) {
cpoint.push_back(p[dim * i + j]);
}
if (parea(cpoint)) {
n++;
val += pfunction(cpoint);
}
}
val = val / n;
double area = std::pow(sd, dim) * n / p_count;
return area * val;
}
double MonteCarloParallel(std::vector<double> a, double sd,
double (*pfunction)(std::vector<double>),
bool (*parea)(std::vector<double>),
unsigned int dim, int p_count) {
if (a.size() != dim) throw(1);
int size, rank;
MPI_Comm_size(MPI_COMM_WORLD, &size);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
int delta = p_count / size;
std::vector<double> p(dim * p_count);
if (rank == 0) {
std::mt19937 gen;
gen.seed(static_cast<unsigned int>(time(NULL)));
std::uniform_real_distribution<> urd(0, 1);
for (unsigned int i = 0; i < p.size(); i++) {
p[i] = a[i % dim] + urd(gen) * sd;
}
}
std::vector<double> loc_p(delta * dim);
MPI_Scatter(&p[0], delta * dim, MPI_DOUBLE, &loc_p[0],
delta * dim, MPI_DOUBLE, 0, MPI_COMM_WORLD);
int loc_n = 0;
int n = 0;
double loc_val = 0;
double val = 0;
for (int i = 0; i < delta; ++i) {
std::vector<double> cpoint;
for (unsigned int j = 0; j < dim; ++j) {
cpoint.push_back(loc_p[dim * i + j]);
}
if (parea(cpoint)) {
loc_n++;
loc_val += pfunction(cpoint);
}
}
MPI_Reduce(&loc_n, &n, 1, MPI_INT, MPI_SUM, 0,
MPI_COMM_WORLD);
MPI_Reduce(&loc_val, &val, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
val = val / n;
double area = std::pow(sd, dim) * n / p_count;
return area * val;
} |
(
