program main

!*****************************************************************************80
!
!! MAIN is the main program for QUAD.
!
!  Discussion:
!
!    QUAD is a FORTRAN90 program to estimate the integral of
!      f = 50.0 / ( pi * ( 2500.0 * x(i) * x(i) + 1.0 ) )
!    from 0 to 10.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    03 September 2013
!
!  Author:
!
!    John Burkardt
!
  use omp_lib

  double precision a
  double precision b
  double precision error
  double precision exact
  integer i
  integer n
  double precision pi
  double precision total
  double precision wtime
  double precision x(1000000)

  a =  0.0D+00
  b = 10.0D+00
  exact = 0.49936338107645674464D+00
  n = 1000000
  pi = 3.141592653589793D+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 from 0 to 10 is 0.49936338107645674464...'
!
!  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
!
!  Evaluate F at each node, sum and average.
!
  wtime = omp_get_wtime ( )

  total = 0.0D+00

!$omp parallel shared ( n, pi, x ) private ( i )
!$omp do 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 do
!$omp end parallel

  total = ( b - a ) * total / dble ( n )

  wtime = omp_get_wtime ( ) - wtime
!
!  Print quadrature estimate, error, time
!
  error = abs ( total - exact )

  write ( *, * ) ' '
  write ( *, * ) '  Estimate = ', total
  write ( *, * ) '  Error    = ', error
  write ( *, * ) '  Wallclock time   = ', wtime

  stop
end