program main

c*********************************************************************72
c
cc MAIN is the main program for QUAD.
c
c  Discussion:
c
c    QUADRATURE estimates an integral using quadrature.
c
c    We break up the interval [A,B] into N subintervals, evaluate
c    F(X) at the midpoint of each subinterval, and divide the
c    sum of these values by N to get an estimate for the integral.
c
c    If we have M processes available because we are using MPI, then
c    we can ask processes 0, 1, 2, ... M-1 to handle the subintervals
c    in the following order:
c
c          0      1       2            M-1  <-- Process numbers begin at 0
c     ------ ------  ------  -----  ------
c          1      2       3    ...       M
c        M+1    M+2     M+3    ...     2*M
c      2*M+1    2*M+2 2*M+3    ...     3*M
c                              
c    and so on up to subinterval N.  The partial sums collected by 
c    each process are then sent to the master process to be added 
c    together to get the estimated integral.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    22 August 2008
c
c  Author:
c
c    John Burkardt
c
      implicit none

      include 'mpif.h'

      double precision a
      parameter ( a = 0.0D+00 )
      character * ( 80 ) arg
      double precision b
      parameter ( b = 1.0D+00 )
      double precision f
      integer i
      integer iargc
      integer id
      integer ierr
      integer n
      integer n_part
      integer p
      double precision q
      double precision q_diff
      double precision q_exact
      parameter ( q_exact = 2.0D+00 )
      double precision q_part
      double precision wtime
      double precision x
c
c  Establish the MPI environment.
c
      call MPI_Init ( ierr )
c
c  Get this process's ID.
c
      call MPI_Comm_rank ( MPI_COMM_WORLD, id, ierr )
c
c  Find out how many processes are available.
c
      call MPI_Comm_size ( MPI_COMM_WORLD, p, ierr )

      n = 1000000

      wtime = MPI_Wtime ( )
c
c  The master process must broadcast the value of N to all processes.
c
      call MPI_Bcast ( n, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr )
c
c  Every process, including the master, now adds up its terms.
c
c  There are N evaluation points total.
c  Each evaluation point is in the center of an interval of width 1/N.
c
      q_part = 0.0D+00
      n_part = 0

      do i = id+1, n, p

        x = ( dble ( 2 * n - 2 + i + 1 ) * a
     &      + dble (         2 * i - 1 ) * b ) 
     &      / dble ( 2 * n     )

        n_part = n_part + 1
        q_part = q_part + f ( x )

      end do
c
c  All the partial sums are collected and summed by the master process.
c
      call MPI_Reduce ( q_part, q, 1, MPI_DOUBLE_PRECISION, MPI_SUM, 
     &  0, MPI_COMM_WORLD, ierr )
c
c  Every process scales the total and prints the answer.
c
      q = q * ( b - a ) / dble ( n )
      q_diff = abs ( q - q_exact )
      wtime = MPI_Wtime ( ) - wtime

      write ( *, '(a)' ) ' '
      write ( *, '(a,g24.16)' ) '  Integral estimate  ', q
      write ( *, '(a,f14.6)' ) '  Elapsed time (seconds) = ', wtime
c
c  Finish up.
c
      call MPI_Finalize ( ierr )

      stop
      end
      function f ( x )

c*********************************************************************72
c
cc F is the function we are integrating.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    22 August 2008
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    Input, double precision X, the argument of the function.
c
c    Output, double precision F, the value of the function.
c
      implicit none

      double precision f
      double precision x

      if ( x .le. 0.0D+00 ) then
        f = 0.0D+00
      else
        f = 1.0D+00 / sqrt ( x )
      end if

      return
      end
      subroutine s_to_i4 ( s, s_len, ival )

c*********************************************************************72
c
cc S_TO_I4 reads an I4 from a string.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    28 April 2008
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    Input, character * ( * ) S, a string to be examined.
c
c    Input, integer S_LEN, the length of S.
c
c    Output, integer IVAL, the integer value read from the string.
c    If the string is blank, then IVAL will be returned 0.
c
      implicit none

      integer s_len

      character c
      integer i
      integer ierror
      integer isgn
      integer istate
      integer ival
      integer length
      character * ( s_len ) s

      ierror = 0
      istate = 0
      isgn = 1
      ival = 0

      do i = 1, s_len

        c = s(i:i)
c
c  Haven't read anything.
c
        if ( istate .eq. 0 ) then

          if ( c .eq. ' ' ) then

          else if ( c .eq. '-' ) then
            istate = 1
            isgn = -1
          else if ( c .eq. '+' ) then
            istate = 1
            isgn = + 1
          else if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then
            istate = 2
            ival = ichar ( c ) - ichar ( '0' )
          else
            ierror = 1
            return
          end if
c
c  Have read the sign, expecting digits.
c
        else if ( istate .eq. 1 ) then

          if ( c .eq. ' ' ) then

          else if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then
            istate = 2
            ival = ichar ( c ) - ichar ( '0' )
          else
            ierror = 1
            return
          end if
c
c  Have read at least one digit, expecting more.
c
        else if ( istate .eq. 2 ) then

          if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then
            ival = 10 * ival + ichar ( c ) - ichar ( '0' )
          else
            ival = isgn * ival
            length = i - 1
            return
          end if

        end if

      end do
c
c  If we read all the characters in the string, see if we're OK.
c
      if ( istate .eq. 2 ) then
        ival = isgn * ival
        length = s_len
      else
        ierror = 1
        length = 0
      end if

      return
      end