# include <stdlib.h>
# include <stdio.h>
# include <math.h>

# 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[] );

/******************************************************************************/

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

/******************************************************************************/
/*
  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.
*/
  ierr = MPI_Init ( &argc, &argv );
/*
  Get the number of processes.
*/
  ierr = MPI_Comm_size ( MPI_COMM_WORLD, &p );
/*
  Determine this processes's rank.
*/
  ierr = MPI_Comm_rank ( MPI_COMM_WORLD, &id );
/*
  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 = ( double * ) malloc ( dim_num * sizeof ( double ) );

  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.
*/
  ierr = MPI_Reduce ( &q, &q_total, 1, MPI_DOUBLE, MPI_SUM, master, 
    MPI_COMM_WORLD );

  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 );
  }

  free ( x );

  ierr = MPI_Finalize ( );

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

double f ( int dim_num, double x[] )

/******************************************************************************/
/*
  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;
}
/******************************************************************************/

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

/******************************************************************************/
/*
  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 )
  {
    printf ( "\n" );
    printf ( "R8VEC_UNIFORM_01 - Fatal error!\n" );
    printf ( "  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;
}