/* -*- Mode: C++; c-basic-offset:4 ; -*- */

/*  

 *  (C) 2004 by Argonne National Laboratory.

 *      See COPYRIGHT in top-level directory.

 */

#include "mpi.h"

/* Support both new (e.g., has iostream) and old (requires iostream.h) 

   C++ compilers */

#ifdef HAVE_CXX_IOSTREAM

#include <iostream>

#ifdef HAVE_NAMESPACE_STD

// Those that do often need the std namespace; otherwise, a bare "cout"

// is likely to fail to compile

using namespace std;

#endif

#else

#include <iostream.h>

#endif

#ifdef HAVE_CXX_MATH

#include <math>

#else

#include <math.h>

#endif

double f(double);

double f(double a)

{

    return (4.0 / (1.0 + a*a));

}

int main(int argc,char **argv)

{

    int n, myid, numprocs, i;

    double PI25DT = 3.141592653589793238462643;

    double mypi, pi, h, sum, x;

    double startwtime = 0.0, endwtime;

    int  namelen;

    char processor_name[MPI_MAX_PROCESSOR_NAME];

    MPI::Init(argc,argv);

    numprocs = MPI::COMM_WORLD.Get_size();

    myid     = MPI::COMM_WORLD.Get_rank();

    MPI::Get_processor_name(processor_name,namelen);

    cout << "Process " << myid << " of " << numprocs << " is on " <<

	processor_name << endl;

    n = 10000;			/* default # of rectangles */

    if (myid == 0)

	startwtime = MPI::Wtime();

    MPI::COMM_WORLD.Bcast(&n, 1, MPI_INT, 0);

    h   = 1.0 / (double) n;

    sum = 0.0;

    /* A slightly better approach starts from large i and works back */

    for (i = myid + 1; i <= n; i += numprocs)

    {

	x = h * ((double)i - 0.5);

	sum += f(x);

    }

    mypi = h * sum;

    MPI::COMM_WORLD.Reduce(&mypi, &pi, 1, MPI_DOUBLE, MPI_SUM, 0);

    if (myid == 0) {

	endwtime = MPI::Wtime();

	cout << "pi is approximately " << pi << " Error is " <<

	    fabs(pi - PI25DT) << endl;

	cout << "wall clock time = " << endwtime-startwtime << endl;

    }

    MPI::Finalize();

    return 0;

}