program main

!*****************************************************************************80
!
!! MAIN is the main program for QUAD.
!
!  Discussion:
!
!    QUAD is a sketch of a FORTRAN90 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.
!
!
!  For an OpenMP version, you'll need 
!
  use omp_lib

  integer n
  parameter ( n = 100000 )

  double precision a
  double precision b
  double precision error
  double precision exact
  integer flops
  integer i
  double precision mflops
  double precision, parameter :: pi = 3.141592653589793D+00
  double precision total
  double precision wtime
  double precision wtime1
  double precision wtime2
  double precision x(n)

  a =  0.0
  b = 10.0
  exact = 0.49936338107645674464D+00

  write ( *, * ) ' '
  write ( *, * ) 'QUAD:'
  write ( *, * ) '  FORTRAN90 version'
  write ( *, * ) ' '
  write ( *, * ) '  Estimate the integral of f(x) from A to B.'
  write ( *, * ) '  f(x) = 50 / ( pi * ( 2500 * x * x + 1 ) ).'
  write ( *, * ) '  A = ', a
  write ( *, * ) '  B = ', b
  write ( *, * ) '  Exact integral is 0.49936338107645674464...'
!
!  Step 1:
!  Load the array X with evenly spaced values between A and B.
!
  do i = 1, n
    x(i) =  ( dble ( n - i     ) * a   &
            + dble (     i - 1 ) * b ) &
            / dble ( n     - 1 )
  end do
!
!  Step 2:
!
!  Read the wall clock timer.
!  Use a do loop to add to TOTAL the value of the function at each X.
!  Read the wall clock timer.
!  Get your estimate by multiplying TOTAL by ( B - A ) and dividing by N.
!  Set FLOPS, the number of floating point operations.
!
!  AFTER you get your code running sequentially, 
!    add the OpenMP include file,
!    insert an OpenMP directive before this loop of the form
!      !$omp parallel do private ( ... ) shared ( ... ) reduction ( + : ... )
!    switch from CPU time measurements to wall clock time, by 
!      ctime1 = omp_get_wtime ( )
!
  wtime1 = omp_get_wtime ( )

  total = 0.0D+00

!$omp parallel do shared ( n, pi, x ) private ( i ) reduction ( + : total )
  do i = 1, n
    total = total + 50.0D+00 / ( pi * ( 2500.0D+00 * x(i) * x(i) + 1.0D+00 ) )
  end do
!$omp end parallel do

  wtime2 = omp_get_wtime ( )

  total = ( b - a ) * total / dble ( n )
!
!  I count 6 floating point operations in the statement inside the loop.
!
  flops = n * 6
!
!  Step 3:
!  Print quadrature estimate, Error, W time, MFLOPS rate
!
  error = abs ( total - exact )
  wtime = wtime2 - wtime1
  mflops = dble ( flops ) / 1000000.0 / wtime

  write ( *, * ) ' '
  write ( *, * ) '  Estimate = ', total
  write ( *, * ) '  Error    = ', error
  write ( *, * ) '  W time   = ', wtime
  write ( *, * ) '  FLOPS    = ', flops
  write ( *, * ) '  MFLOPS   = ', mflops

  stop
end