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

int main ( int argc, char *argv[] );
double cpu_time ( void );

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

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

/******************************************************************************/
/*
  Purpose:

    MAIN is the main program for QUAD.

  Discussion:

    QUAD is a sketch of a C program to estimate the integral of
    a function F(X) between A and B.

    This program is only a sketch.  Your job is to finish it,
    debug it, and get it to run.

    Experiment to get a value of N that gives you a 'reasonable' accuracy.

    NOW modify the program so that it will run under OpenMP.  Replace
    CPU time by wall clock time.

    Show that your program STILL gets the right answer.

    Then determine the speedup when you run the OpenMP version on
    1, 2, 4, 8 processors.
*/
{
# define N 1000000

  double a = 0.0;
  double b = 10.0;
  double error;
  double exact = 0.49936338107645674464;
  int flops;
  int i;
  double mflops;
  double pi = 3.141592653589793;
  double total;
  double wtime;
  double wtime1;
  double wtime2;
  double x[N];

  printf ( "\n" );
  printf ( "QUAD:\n" );
  printf ( "  C version\n" );
  printf ( "\n" );
  printf ( "  Estimate the integral of f(x) from A to B.\n" );
  printf ( "  f(x) = 50 / ( pi * ( 2500 * x * x + 1 ) ).\n" );
  printf ( "  A = %f\n", a );
  printf ( "  B = %f\n", b );
  printf ( "  Exact integral is 0.49936338107645674464...\n" );
/*
  Step 1:
  Load the array X with evenly spaced values between A and B.
*/
  for ( i = 0; i < N; i++ )
  {
    x[i] =  ( ( double ) ( N - i - 1 ) * a   
            + ( double ) (     i     ) * b ) 
            / ( double ) ( N     - 1 );
  }
/*
  Step 2:

  Use a FOR loop to add to TOTAL the value of the function at each X.

  Set FLOPS, the number of floating point operations, by counting the
  number of floating operations in the statement that updates TOTAL,
  and multiplying by N (number of times you did the loop).

  add the OpenMP include file,

  insert an OpenMP directive before the loop of the form

    # pragma omp parallel for private ( ... ) shared ( ... ) reduction ( + : ... )
*/
  wtime1 = omp_get_wtime ( );
  total = 0.0;
/*
  Note that "N" is not a variable!  It's a symbolic constant.
*/

# pragma omp parallel for private ( i ) shared ( pi, x ) reduction ( + : total )
  for ( i = 0; i < N; i++ )
  {
    total = total + 50 / pi / ( 2500.0 * x[i] * x[i] + 1.0 );
  }

  flops = 6 * N;
  wtime2 = omp_get_wtime ( );
  total = ( b - a ) * total / ( double ) N;
/*
  Step 3:
  Print quadrature estimate, Error, W time, MFLOPS rate
*/
  error = fabs ( total - exact );
  wtime = wtime2 - wtime1;
  mflops = ( double ) ( flops ) / 1000000.0 / wtime;

  printf ( "\n" );
  printf ( "  Estimate = %f\n", total );
  printf ( "  Error    = %e\n", error );
  printf ( "  W time   = %f\n", wtime );
  printf ( "  FLOPS    = %d\n", flops );
  printf ( "  MFLOPS   = %f\n", mflops );

  return 0;
# undef N
}
/*******************************************************************************/

double cpu_time ( void )

/*******************************************************************************/
/*
  Purpose:
 
    CPU_TIME reports the total CPU time for a program.

  Modified:

    27 September 2005

  Author:

    John Burkardt

  Parameters:

    Output, double CPU_TIME, the current total elapsed CPU time in second.
*/
{
  double value;

  value = ( double ) clock ( ) / ( double ) CLOCKS_PER_SEC;

  return value;
}