pi.cpp

#include <math.h>
#include <mpich-devel-clang/mpi.h>
#include <iostream>

using namespace std;

int main(int argc, char* argv[]){
    int n, rank, size, i;
    // PI25DT is used for computing error
    double PI25DT = 3.141592653589793238462643;
    double mypi, pi, h, sum, x;
    // Initialize cluster
    MPI::Init(argc, argv);
    // Get cluster size ( number of nodes in the cluster )
    size = MPI::COMM_WORLD.Get_size();
    // Get the node number of currnent node.
    rank = MPI::COMM_WORLD.Get_rank();
    // NOTE: Since C doesn't have namespaces, the namings of MPI for C is different,
    // MPI::COMM_WORLD.Get_size() => MPI_Comm_size(MPI_COMM_WORLD, ...)
    while(1){
        // If this is the first node(mothership), we prompt the number of intervals user wanted.
        if (rank == 0) {
            cout << "Enter the number of intervals: (0 quits)"
                << endl;
            cin >> n;
        }

        // Broadcast the number of intervals to all other nodes
        MPI::COMM_WORLD.Bcast(&n, 1, MPI::INT, 0);
        if (n==0)
            break; // The mothership node doesn't do anything, it only waits for all other nodes finish.
        else {
            h = 1.0 / (double) n;
            sum = 0.0;
            for (i = rank + 1; i <= n; i += size) {
                x = h * ((double)i - 0.5);
                sum += (4.0 / (1.0 + x*x));
            }
            mypi = h * sum;
            // Reduce(fold) values of all nodes to a single value
            // Here we only computes the sum of all values.
            MPI::COMM_WORLD.Reduce(&mypi, &pi, 1, MPI::DOUBLE,
                    MPI::SUM, 0);
            if (rank == 0)
                cout << "pi is approximately " << pi
                    << ", Error is " << fabs(pi - PI25DT)
                    << endl;
        }
    }
    // Clean and exit
    MPI::Finalize();
    return 0;
}