# include <cstdlib>
# include <iostream>
# include <iomanip>
# include <cmath>

using namespace std;

# include "mpi.h"

int main ( int argc, char *argv[] );
double f ( double x );

//****************************************************************************80

int main ( int argc, char *argv[] )

//****************************************************************************80
//
//  Purpose:
//
//    MAIN is the main program for QUAD.
//
//  Discussion:
//
//    QUADRATURE estimates an integral using quadrature.
//
//    We break up the interval [A,B] into N subintervals, evaluate
//    F(X) at the midpoint of each subinterval, and divide the
//    sum of these values by N to get an estimate for the integral.
//
//    If we have M processes available because we are using MPI, then
//    we can ask processes 0, 1, 2, ... M-1 to handle the subintervals
//    in the following order:
//
//          0      1       2            M-1  <-- Process numbers begin at 0
//     ------ ------  ------  -----  ------
//          1      2       3    ...       M
//        M+1    M+2     M+3    ...     2*M
//      2*M+1    2*M+2 2*M+3    ...     3*M
//                              
//    and so on up to subinterval N.  The partial sums collected by 
//    each process are then sent to the master process to be added 
//    together to get the estimated integral.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    23 September 2008
//
//  Author:
//
//    John Burkardt
//
{
  double a = 0.0;
  double b = 1.0;
  int i;
  int id;
  int ierr;
  int n;
  int n_part;
  int p;
  double q;
  double q_diff;
  double q_exact = 2.0;
  double q_part;
  double wtime;
  double x;
//
//  Initialize MPI.
//
  MPI::Init ( argc, argv );
//
//  Get the number of processes.
//
  p = MPI::COMM_WORLD.Get_size (  );
//
//  Determine this processes's rank.
//
  id = MPI::COMM_WORLD.Get_rank ( );

  n = 1000000;
//
//  Record the starting time.
//
  wtime = MPI::Wtime ( );
//
//  Broadcast the number of intervals to the other processes.
//
  MPI::COMM_WORLD.Bcast ( &n, 1, MPI::INT, 0 );
//
//  Integrate F(X) over a subinterval determined by your process ID.
//
  q_part = 0.0;
  n_part = 0;

  for ( i = id + 1; i <= n; i = i + p ) 
  {
    x = ( ( double ) ( 2 * n - 2 * i + 1 ) * a
        + ( double ) (         2 * i - 1 ) * b )
        / ( double ) ( 2 * n     );

    n_part = n_part + 1;
    q_part = q_part + f ( x );
  }
//
//  Each process sends its value Q_PART to the MASTER process, to be added to
//  the global result Q.
//
  MPI::COMM_WORLD.Reduce ( &q_part, &q, 1, MPI::DOUBLE, MPI::SUM, 0 );
//
//  Every process scales the sum and prints the result.
//
  q = q * ( b - a ) / ( double ) n;
  q_diff = fabs ( q - q_exact );
  wtime = MPI::Wtime ( ) - wtime;

  cout << "\n";  
  cout << "  The estimated integral " << q << "\n";
  cout << "  Wall clock elapsed seconds = " << wtime << "\n";
//
//  Shut down MPI.
//
  MPI::Finalize ( );

  return 0;
}
//****************************************************************************80

double f ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    F evaluates the function F(X) which we are integrating.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    22 August 2008
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double X, the argument of the function.
//
//    Output, double F, the value of the function.
//
{
  double value;

  if ( x <= 0.0 )
  {
    value = 0.0;
  }
  else
  {
    value = 1.0 / sqrt ( x );
  }

  return value;
}