/
calc_pi_reduce.cpp
42 lines (39 loc) · 1.11 KB
/
calc_pi_reduce.cpp
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
#include <cstdio>
#include <random>
#include <algorithm>
#include <cmath>
#include <mpi.h>
const int TRIAL = 100000;
double calc_pi(const int seed) {
std::mt19937 mt(seed);
std::uniform_real_distribution<double> ud(0.0, 1.0);
int n = 0;
for (int i = 0; i < TRIAL; i++) {
double x = ud(mt);
double y = ud(mt);
if (x * x + y * y < 1.0) n++;
}
return 4.0 * static_cast<double>(n) / static_cast<double>(TRIAL);
}
int main(int argc, char **argv) {
MPI_Init(&argc, &argv);
int rank;
int procs;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_size(MPI_COMM_WORLD, &procs);
double pi = calc_pi(rank);
double pi2 = pi * pi;
double pi_sum = 0.0;
double pi2_sum = 0.0;
printf("%f\n", pi);
MPI_Allreduce(&pi, &pi_sum, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(&pi2, &pi2_sum, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
double pi_ave = pi_sum / procs;
double pi_var = pi2_sum / (procs - 1) - pi_sum * pi_sum / procs / (procs - 1);
double pi_stdev = sqrt(pi_var);
MPI_Barrier(MPI_COMM_WORLD);
if (rank == 0) {
printf("pi = %f +- %f\n", pi_ave, pi_stdev);
}
MPI_Finalize();
}