#ifdef LINUX_CLUSTER
  #include <mpi++.h>
#else
  #include <mpi.h>
#endif
/* iostream must not have suffix .h; all others OK as is */
#include <iostream>
#include <math.h>
#include <stdio.h>
using namespace std;
float fct(float x)
{
      return cos(x);
}
/* Prototype */
float integral(float a, int i, float h, int n);
int main(int argc, char* argv[])
{
/***********************************************************************
 *                                                                     *
 * This is one of the MPI versions on the integration example          *
 * It demonstrates the use of :                                        *
 *                                                                     *
 * 1) MPI_Init                                                         *
 * 2) MPI_Comm_rank                                                    *
 * 3) MPI_Comm_size                                                    *
 * 4) MPI_Recv                                                         *
 * 5) MPI_Send                                                         *
 * 6) MPI_Finalize                                                     *
 *                                                                     *
 * Dr. Kadin Tseng                                                     *
 * Scientific Computing and Visualization                              *
 * Boston University                                                   *
 * 1998                                                                *
 *                                                                     *
 ***********************************************************************/
      int n, p;
      float h, result, a, b, pi;

      int myid, src, dest, tag;
      MPI::Status status;
      float my_result;

/* Starts MPI processes ... */

      MPI::Init(argc, argv);
      myid = MPI::COMM_WORLD.Get_rank();
      p = MPI::COMM_WORLD.Get_size();

      pi = acos(-1.0);  /* = 3.14159... */
      a = 0.;           /* lower limit of integration */
      b = pi*1./2.;     /* upper limit of integration */
      n = 500;          /* number of increment within each process */

      dest = 0;         /* define the process that computes the final result */
      tag = 123;        /* set the tag to identify this particular job */

      h = (b-a)/n/p;    /* length of increment */

      my_result = integral(a,myid,h,n);

      cout << "Process " << myid <<" has the partial result " 
         << my_result << endl;
  
      MPI::COMM_WORLD.Send(&my_result, 1, MPI::FLOAT, dest, tag);
                            /* send my_result to dest. */

      if(myid == 0) {
        result = 0.0;
        for (src=0;src<p;src++) {
          MPI::COMM_WORLD.Recv(&my_result, 1, MPI::FLOAT, src, tag, status);  /* not safe */
          result += my_result;
        }
        cout << "The result = " << result << endl;
      }
      MPI::Finalize();                       /* let MPI finish up ... */
}
float integral(float a, int i, float h, int n)
{
      int j;
      float h2, aij, integ;

      integ = 0.0;                 /* initialize integral */
      h2 = h/2.;
      for (j=0;j<n;j++) {             /* sum over all "j" integrals */
        aij = a + (i*n + j)*h;        /* lower limit of "j" integral */
        integ += fct(aij+h2)*h;
      }  
      return integ;
}