pi

/* a MC example calculate value of pi*/

//const double totaltry=10000;

#include <stdlib.h>
#include <iostream>
#include <string>
#include <mpi.h>

using namespace std;

int main(int argc, char** argv) {

  double totaltry(10000);

  if (argc != 2) {
    cout << "Usage: \n" << argv[0] << " <# of trys>\n";
    return(1);
  }

  totaltry = stof(argv[1], 0);
  // Initialize the MPI environment
  MPI_Status mpistat;
  double HITSBUF[1];
  double starttime, endtime;

  MPI_Init(0, 0);

  // Get the number of processes
  int world_size;
  MPI_Comm_size(MPI_COMM_WORLD, &world_size);
  // Get the rank of the process
  int world_rank;
  MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);

  srand(time(0)+world_rank);

  for (int num_nodes = 1; num_nodes <= world_size; num_nodes <<= 1) {

    starttime = MPI_Wtime();
    MPI_Barrier(MPI_COMM_WORLD);
    double eachtry = totaltry/num_nodes;
    double localhits=0;
    double lbound = (eachtry * world_rank)/totaltry;
    double vy, vx;
    for (int i=0; i< eachtry ; i++) {
      vy = rand()*1.0/RAND_MAX;
      vx = lbound+i/totaltry;
      localhits += (vx*vx + vy*vy) < 1;
    }
    MPI_Barrier(MPI_COMM_WORLD);

    int mask = 0xFFFFFFFE;
    for (int level=num_nodes; level>1; level/=2) {
      int sourcenode = world_rank & mask;
      int targetnode = world_rank | (~mask);
      if (world_rank == sourcenode) {
        *HITSBUF = localhits;
        MPI_Send(HITSBUF, 1, MPI_DOUBLE, targetnode, 1, MPI_COMM_WORLD);
      } else {
        MPI_Recv(HITSBUF, 1, MPI_DOUBLE, sourcenode, 1, MPI_COMM_WORLD, &mpistat
);
        localhits += *HITSBUF;
      }
      mask = mask << 1 | 1;
    }
    endtime = MPI_Wtime();

    if (world_rank == num_nodes - 1)
      cout << " pi: "<< localhits*4.0/totaltry
        << " number of nodes: " << num_nodes
        << " taking " << endtime-starttime << " seconds" << endl;
  }

  // Finalize the MPI environment.
  MPI_Finalize();
}