program main

!*****************************************************************************80
!
!! MAIN is the main program for QUADRATURE.
!
!  Discussion:
!
!    QUADRATURE estimates an integral using quadrature.
!
!    The integral of F(X) = 4 / ( 1 + X * X ) from 0 to 1 is PI.
!
!    We break up the interval [0,1] into N subintervals, evaluate
!    F(X) at the midpoint of each subinterval, and multiply the
!    sum of these values by N to get an estimate for the integral.
!
!    If we have M processes available because we are using MPI, then
!    we can ask processes 0, 1, 2, ... M-1 to handle the subintervals
!    in the following order:
!
!          0      1       2            M-1  <-- Process numbers begin at 0
!     ------ ------  ------  -----  ------
!          1      2       3    ...       M
!        M+1    M+2     M+3    ...     2*M
!      2*M+1    2*M+2 2*M+3    ...     3*M
!                              
!    and so on up to subinterval N.  The partial sums collected by 
!    each process are then sent to the master process to be added 
!    together to get the estimated integral.
!
!  Modified:
!
!    15 October 2007
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    William Gropp, Ewing Lusk, Anthony Skjellum,
!    Using MPI: Portable Parallel Programming with the
!    Message-Passing Interface,
!    Second Edition,
!    MIT Press, 1999,
!    ISBN: 0262571323.
!
!    Snir, Otto, Huss-Lederman, Walker, Dongarra,
!    MPI - The Complete Reference,
!    Volume 1, The MPI Core,
!    second edition,
!    MIT Press, 1998.
!

!
!  Fortran77 include file:
!
  include 'mpif.h'
!
!  Fortran90 module:
!
!  use mpi
!
! implicit none
!
  real    ( kind = 8 ) f
  real    ( kind = 8 ) h
  integer ( kind = 4 ) i
  integer ( kind = 4 ) ierr
  integer ( kind = 4 ), parameter :: master = 0
  integer ( kind = 4 ) my_id
  real    ( kind = 8 ) my_part
  integer ( kind = 4 ) n
  integer ( kind = 4 ) num_procs
  real    ( kind = 8 ) pi
  real    ( kind = 8 ) pi_diff
  real    ( kind = 8 ), parameter :: pi_exact = 3.141592653589793238462643D+00
  real    ( kind = 8 ) sum2
  real    ( kind = 8 ) wtime_diff
  real    ( kind = 8 ) wtime_end
  real    ( kind = 8 ) wtime_start
  real    ( kind = 8 ) x
!
!  Initialize MPI.
!
  call MPI_Init ( ierr )
!
!  Get this process's ID.
!
  call MPI_Comm_rank ( MPI_COMM_WORLD, my_id, ierr )
!
!  Find out how many processes are available.
!
  call MPI_Comm_size ( MPI_COMM_WORLD, num_procs, ierr )

  if ( my_id == master ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'QUADRATURE - Master process:'
    write ( *, '(a)' ) '  FORTRAN90 version'
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  An MPI example program to compute PI.'
    write ( *, '(a,i8)' ) '  The number of processes is ', num_procs

    start_time = MPI_Wtime ( )

  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8,a)' ) 'Process ', my_id, ' is active.'
!
!  Assume that the master process just got the value of N from the user.
!  Here, we'll use an assignment statement.
!
  if ( my_id == master ) then
    n = 1000
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'QUADRATURE - Master process:'
    write ( *, '(a,i8,a)' ) 'Number of intervals being used is ', n
  end if
!
!  The master process broadcasts the value of N to all processes.
!
  call MPI_Bcast ( n, 1, MPI_INTEGER, master, MPI_COMM_WORLD, ierr )
!
!  Process MY_ID now adds up its terms.
!
  h = 1.0D+00 / real ( n, kind = 8 )
  my_part = 0.0D+00
  do i = my_id + 1, n, num_procs
    x = h * ( real ( i, kind = 8 ) - 0.5D+00 )
    my_part = my_part + f ( x )
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) 'QUADRATURE - Process ', my_id
  write ( *, '(a,g24.16)' ) '  My contribution is ', my_part * h
!
!  Each process sends its value of MY_PART to the master process, to
!  be summed in PI.
!
  call MPI_Reduce ( my_part, pi, 1, MPI_DOUBLE_PRECISION, MPI_SUM, master, &
    MPI_COMM_WORLD, ierr )
!
!  The master process scales the sum by H and prints the answer.
!
  if ( my_id == master ) then

    pi = h * pi

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'QUADRATURE - Master process:'
    write ( *, '(a,g24.16)' ) '  Estimate for PI is ', pi
    write ( *, '(a,g24.16)' ) '  Exact value is     ', pi_exact
    pi_diff = abs ( pi_exact - pi )
    write ( *, '(a,g24.16)' ) '  Error is           ', pi_diff

    wtime_end = MPI_Wtime ( )
    wtime_diff = wtime_end - wtime_start

    write ( *, '(a)' ) ' '
    write ( *, '(a,f14.6)' ) '  Elapsed wall clock seconds = ', wtime_diff

  end if
!
!  Finish up.
!
  call MPI_Finalize ( ierr )

  if ( my_id == master ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'QUADRATURE - Master process:'
    write ( *, '(a)' ) '  Normal end of execution.'
  end if

  stop
end
function f ( x )

!*****************************************************************************80
!
!! F is the function we are integrating.
!
!  Discussion:
!
!    Integral ( 0 <= X <= 1 ) 4/(1+X*X) dX = PI
!
!  Modified:
!
!    10 February 2000
!
!  Parameters:
!
!    Input, real ( kind = 8 ) X, the argument of the function.
!
!    Output, real ( kind = 8 ) F, the value of the function.
!
  implicit none

  real ( kind = 8 ) f
  real ( kind = 8 ) x

  f = 4.0D+00 / ( 1.0D+00 + x * x )

  return
end