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

using namespace std;

# include "mpi.h"

int main ( int argc, char *argv[] );
double f ( int dim_num, double x[] );
void r8vec_uniform_01 ( int n, int *seed, double r[] );

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

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

//****************************************************************************80
//
//  Purpose:
//
//    MAIN is the main program for MC.
//
//  Discussion:
//
//    Things you can vary include:
//
//    * the function to be integrated (see the list below)
//    * the number of samples used, SAMPLE_NUM.
//    * the spatial dimension D_NUM (try 1, 2, 4, 8)
//
//    The function to be integrated is
//
//      F(X) = Product ( 1 <= DIM <= DIM_NUM ) ( 2 * abs ( 2 * X(DIM) - 1 ) ) 
//      Exact integral is 1.
//
{
  int dim_num = 4;
  int id;
  int ierr;
  int master = 0;
  int p;
  double q;
  double q_error;
  double q_exact = 1.0;
  double q_total;
  int sample;
  int sample_num = 1000;
  int sample_total;
  int seed;
  double wtime;
  double wtime1;
  double wtime2;
  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 ( );
//
//  Say hello.
//
  if ( id == master ) 
  {
    printf ( "\n" );
    printf ( "MC\n" );
    printf ( "  C++ version\n" );
    printf ( "  An MPI program to estimate an integral\n" );
    printf ( "  The number of processes is %d\n", p );
    wtime1 = MPI::Wtime ( );
  }

  seed = 123456789 + id;

  x = new double[dim_num];

  q = 0.0;

  for ( sample = 1; sample <= sample_num; sample++ )
  {
    r8vec_uniform_01 ( dim_num, &seed, x );

    q = q + f ( dim_num, x );
  }

  q = q / ( double ) ( sample_num );
  q_error = fabs ( q - q_exact );

  if ( id == master ) 
  {
    printf ( "\n" );
    printf ( "   Process  Samples Dimension      Estimate          Error\n" );
    printf ( "\n" );
  }

  printf ( "  %8d  %8d  %8d  %16f  %16e\n", id, sample_num, dim_num, q, q_error );
//
//  Have each process send results to process MASTER for reduction to
//  final result.
//
  MPI::COMM_WORLD.Reduce ( &q, &q_total, 1, MPI::DOUBLE, MPI::SUM, master );

  if ( id == master )
  {
    wtime2 = MPI::Wtime ( );

    q_total = q_total / ( double ) ( p );
    q_error = fabs ( q_total - q_exact );
    sample_total = p * sample_num;
    printf ( "\n" );
    printf ( "     Total  %8d  %8d  %16f  %16e\n", 
      sample_total, dim_num, q_total, q_error );
    printf ( "\n" );
    printf ( "  Elapsed wall clock seconds = %f\n", wtime );
  }

  delete [] x;

  MPI::Finalize ( );

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

double f ( int dim_num, double x[] )

//****************************************************************************80
//
//  Purpose:
//
//    F evaluates the function F(X) which we are integrating.
//
//  Modified:
//
//    04 September 2008
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the number of spatial dimensions.
//
//    Input, double X[DIM_NUM], the point at which to evaluate F.
//
//    Output, double F, the value of F(X).
//
{
  int dim;
  double value;

  value = 1.0;
  for ( dim = 0; dim < dim_num; dim++ )
  {
    value = value * 2.0 * fabs ( 2.0 * x[dim] - 1.0 );
  }

  return value;
}
//****************************************************************************80

void r8vec_uniform_01 ( int n, int *seed, double r[] )

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_UNIFORM_01 returns a unit pseudorandom vector.
//
//  Parameters:
//
//    Input, int N, the number of entries in the vector.
//
//    Input/output, int *SEED, a seed for the random number generator.
//
//    Output, double R[N], the vector of pseudorandom values.
//
{
  int i;
  int i4_huge = 2147483647;
  int k;

  if ( *seed == 0 )
  {
    cerr << "\n";
    cerr << "R8VEC_UNIFORM_01 - Fatal error!\n";
    cerr << "  Input value of SEED = 0.\n";
    exit ( 1 );
  }

  for ( i = 0; i < n; i++ )
  {
    k = *seed / 127773;

    *seed = 16807 * ( *seed - k * 127773 ) - k * 2836;

    if ( *seed < 0 )
    {
      *seed = *seed + i4_huge;
    }

    r[i] = ( double ) ( *seed ) * 4.656612875E-10;
  }

  return;
}