#include <cassert>
#include <cstdio>
#include <mpi.h>
#include <hpc/aux/slices.h>
using namespace std;
using namespace hpc::aux;
constexpr double a = 0;
constexpr double b = 1;
constexpr unsigned int N = 100; /* number of intervals */
// numerical integration according to the Simpson rule
// for f over the interval [a,b] using n subintervals
template <typename F>
double simpson(F f, double a, double b, int n) {
assert(n > 0 && a <= b);
double value = f(a)/2 + f(b)/2;
double xleft;
double x = a;
for (int i = 1; i < n; ++i) {
xleft = x; x = a + i * (b - a) / n;
value += f(x) + 2 * f((xleft + x)/2);
}
value += 2 * f((x + b)/2);
value *= (b - a) / n / 3;
return value;
}
char* cmdname;
void usage() {
cerr << "Usage: " << cmdname << " [# intervals] " << endl;
MPI_Abort(MPI_COMM_WORLD, 1);
}
int main(int argc, char** argv) {
MPI_Init(&argc, &argv);
int nof_processes; MPI_Comm_size(MPI_COMM_WORLD, &nof_processes);
int rank; MPI_Comm_rank(MPI_COMM_WORLD, &rank);
Slices<decltype(N)> slices(nof_processes, N);
auto first_interval = slices.offset(rank);
auto nofintervals = slices.size(rank);
auto next_interval = first_interval + nofintervals;
double x0 = a + first_interval * (b - a) / N;
double x1 = a + next_interval * (b - a) / N;
double value = simpson([](double x) -> double {
return 4 / (1 + x*x);
}, x0, x1, nofintervals);
if (rank) {
MPI_Send(&value, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
} else {
double sum = value;
for (int i = 0; i + 1 < nof_processes; ++i) {
MPI_Status status;
double value;
MPI_Recv(&value, 1, MPI_DOUBLE, MPI_ANY_SOURCE,
0, MPI_COMM_WORLD, &status);
sum += value;
}
printf("%.15lf\n", sum);
}
MPI_Finalize();
}