Code Example

program.cxx

# include 
# include 
# include 
# include 

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 QUADRATURE.
//
//  Discussion:
//
//    QUADRATURE estimates an integral using quadrature.
//
//    The integral of F(X) = 4 / ( 1 + X * X ) from 0 to 1 is PI.
//
//    We break up the interval [0,1] into N subintervals, evaluate
//    F(X) at the midpoint of each subinterval, and multiply 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.
//
//  Modified:
//
//    15 October 2007
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    William Gropp, Ewing Lusk, Anthony Skjellum,
//    Using MPI: Portable Parallel Programming with the
//    Message-Passing Interface,
//    Second Edition,
//    MIT Press, 1999,
//    ISBN: 0262571323.
//
{
  double h;
  int i;
  int ierr;
  int master = 0;
  double my_part;
  int n;
  double pi;
  double pi_diff;
  double pi_exact = 3.141592653589793238462643;
  int process_id;
  int process_num;
  double wtime_diff;
  double wtime_end;
  double wtime_start;
  double x;
//
//  Initialize MPI.
//
  MPI::Init ( argc, argv );
//
//  Get the number of processes.
//
  process_num = MPI::COMM_WORLD.Get_size (  );
//
//  Determine this processes's rank.
//
  process_id = MPI::COMM_WORLD.Get_rank ( );
//
//  Say hello.
//
  if ( process_id == master )
  {
    cout << "\n";
    cout << "QUADRATURE - Master process:\n";
    cout << "  C++ version\n";
    cout << "\n";
    cout << "  Estimate the value of PI by approximating an integral.\n";
    cout << "  The integral is approximated by a sum,\n";
    cout << "  whose calculation is divided among a number of processes.\n";
    cout << "\n";
    cout << "  The number of processes available is " << process_num << "\n";
  }
//
//  Record the starting time.
//
  if ( process_id == master )
  {
    wtime_start = MPI::Wtime ( );
  }
//
//  Normally, we will assume, the number of intervals N is read in
//  by process 0.  We'll just simulate this by assigning a value to
//  N, but only on process 0.
//
  if ( process_id == master )
  {
    n = 1000;
    cout << "\n";
    cout << "QUADRATURE - Master process:\n";
    cout << "  Number of intervals used is " << n << "\n";
  }
//
//  Broadcast the number of intervals to the other processes.
//
  MPI::COMM_WORLD.Bcast ( &n, 1, MPI::INT, master );
//
//  Integrate F(X) over a subinterval determined by your process ID.
//
  h = 1.0 / ( double ) n;
  my_part = 0.0;
  for ( i = process_id + 1; i <= n; i = i + process_num )
  {
    x = h * ( ( double ) i - 0.5 );
    my_part = my_part + f ( x );
  }

  cout << "\n";
  cout << "QUADRATURE - Process " << process_id << "\n";
  cout << "  My contribution is " << my_part * h << ".\n";
//
//  Each process sends its value MY_PART to the MASTER process, to be added to
//  the global result PI.
//
  MPI::COMM_WORLD.Reduce ( &my_part, &pi, 1, MPI::DOUBLE, MPI::SUM, master );
//
//  The master process scales the sum and prints the result.
//
  if ( process_id == master )
  {
    pi = pi * h;

    cout << "\n";
    cout << "QUADRATURE - Master process:\n";
    cout << "  The estimate for PI is " << pi << "\n";
    cout << "  The exact value is     " << pi_exact << "\n";
    pi_diff = fabs ( pi - pi_exact );
    cout << "  The error is           " << pi_diff << "\n";

    wtime_end = MPI_Wtime ( );
    wtime_diff = wtime_end - wtime_start;

    cout << "\n";
    cout << "  Wall clock elapsed seconds = " << wtime_diff << "\n";
  }
//
//  Shut down MPI.
//
  MPI::Finalize ( );

  if ( process_id == master )
  {
    cout << "\n";
    cout << "QUADRATURE - Master process:\n";
    cout << "  Normal end of execution.\n";
  }

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

double f ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    F evaluates the function F(X) which we are integrating.
//
//  Discussion:
//
//    Integral ( 0 <= X <= 1 ) 4/(1+X*X) dX = PI
//
//  Modified:
//
//    19 March 2006
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double X, the argument of the function.
//
//    Output, double F, the value of the function.
//
{
  double value;

  value = 4.0 / ( 1.0 + x * x );

  return value;
}