program main

c*********************************************************************72
c
cc MAIN is the main program for QUAD.
c
c  Discussion:
c
c    QUAD is a sketch of a FORTRAN90 program to estimate the integral of
c    a function F(X) between A and B.
c
c    This program is only a sketch.  Your job is to finish it,
c    debug it, and get it to run.
c
c    Experiment to get a value of N that gives you a 'reasonable' accuracy.
c
c    NOW modify the program so that it will run under OpenMP.  Replace
c    CPU time by wall clock time.
c
c    Show that your program STILL gets the right answer.
c
c    Then determine the speedup when you run the OpenMP version on
c    1, 2, 4, 8 processors.
c
c  For an OpenMP version, you'll need 
c
      include 'omp_lib.h'

      integer n
      parameter ( n = 1000000 )

      double precision a
      double precision b
      double precision error
      double precision exact
      integer flops
      integer i
      double precision mflops
      double precision pi
      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...'
c
c  Step 1:
c  Load the array X with evenly spaced values between A and B.
c
      do i = 1, n
        x(i) =  ( dble ( n - i     ) * a   
     &          + dble (     i - 1 ) * b ) 
     &          / dble ( n     - 1 )
      end do
c
c  Step 2:
c
c  Use a DO loop to add to TOTAL the value of the function at each X.
c
c  Get your estimate by multiplying TOTAL by ( B - A ) and dividing by N.
c
c  Set FLOPS, the number of floating point operations
c
c  Add the OpenMP include file.
c
c  Insert an OpenMP directive before this loop of the form
c
c    c$omp parallel do private ( ... ) shared ( ... ) reduction ( + : ... )
c
c  Insert an OpenMP directive after the loop of the form
c
c    c$omp end parallel do
c
      wtime1 = omp_get_wtime ( )

      total = 0.0D+00

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

      wtime2 = omp_get_wtime ( )

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

      flops = n * 6
c
c  Step 3:
c  Print quadrature estimate, Error, CPU time, MFLOPS rate
c
      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