program main !*****************************************************************************80 ! !! MAIN is the main program for QUAD. ! ! Discussion: ! ! QUAD estimates an integral using quadrature. ! ! We break up the interval [A,B] into N subintervals, evaluate ! F(X) at the midpoint of each subinterval, and divide 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. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 22 August 2008 ! ! Author: ! ! John Burkardt ! use mpi implicit none real ( kind = 8 ), parameter :: a = 0.0D+00 character ( len = 80 ) arg real ( kind = 8 ), parameter :: b = 1.0D+00 real ( kind = 8 ) f integer ( kind = 4 ) i integer ( kind = 4 ) iargc integer ( kind = 4 ) id integer ( kind = 4 ) ierr integer ( kind = 4 ) n integer ( kind = 4 ) n_part integer ( kind = 4 ) p real ( kind = 8 ) q real ( kind = 8 ) q_diff real ( kind = 8 ), parameter :: q_exact = 2.0D+00 real ( kind = 8 ) q_part real ( kind = 8 ) wtime real ( kind = 8 ) x ! ! Initialize MPI. ! call MPI_Init ( ierr ) ! ! Get this process's ID. ! call MPI_Comm_rank ( MPI_COMM_WORLD, id, ierr ) ! ! Find out how many processes are available. ! call MPI_Comm_size ( MPI_COMM_WORLD, p, ierr ) n = 1000000 wtime = MPI_Wtime ( ) ! ! The master process broadcasts the value of N to all other processes. ! call MPI_Bcast ( n, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr ) ! ! Every process, including the master, now adds up its terms. ! ! There are N evaluation points total. ! Each evaluation point is in the center of an interval of width 1/N. ! q_part = 0.0D+00 n_part = 0 do i = id + 1, n, p x = ( real ( 2 * n - 2 * i + 1, kind = 8 ) * a & + real ( 2 * i - 1, kind = 8 ) * b ) & / real ( 2 * n, kind = 8 ) n_part = n_part + 1 q_part = q_part + f ( x ) end do ! ! Each process sends its value of Q_PART to the master process, to ! be summed in Q. ! call MPI_Reduce ( q_part, q, 1, MPI_DOUBLE_PRECISION, MPI_SUM, 0, & MPI_COMM_WORLD, ierr ) ! ! Every process scales the sum and prints the answer. ! q = q * ( b - a ) / real ( n, kind = 8 ) q_diff = abs ( q_exact - q ) wtime = MPI_Wtime ( ) - wtime write ( *, '(a)' ) ' ' write ( *, '(a,g24.16)' ) ' Integral estimate ', q write ( *, '(a,f14.6)' ) ' Elapsed wall clock time = ', wtime ! ! Finish up. ! call MPI_Finalize ( ierr ) stop end function f ( x ) !*****************************************************************************80 ! !! F is the function we are integrating. ! ! Modified: ! ! 22 August 2008 ! ! 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 if ( x <= 0.0D+00 ) then f = 0.0D+00 else f = 1.0D+00 / sqrt ( x ) end if return end subroutine s_to_i4 ( s, value ) !*****************************************************************************80 ! !! S_TO_I4 reads an integer value from a string. ! ! Discussion: ! ! Instead of ICHAR, we now use the IACHAR function, which ! guarantees the ASCII collating sequence. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 12 January 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, a string to be examined. ! ! Output, integer ( kind = 4 ) VALUE, the integer value read from the string. ! If the string is blank, then VALUE will be returned 0. ! ! Output, integer ( kind = 4 ) IERROR, an error flag. ! 0, no error. ! 1, an error occurred. ! ! Output, integer ( kind = 4 ) LENGTH, the number of characters ! of S used to make the integer. ! implicit none character c integer ( kind = 4 ) i integer ( kind = 4 ) ierror integer ( kind = 4 ) isgn integer ( kind = 4 ) length character ( len = * ) s integer ( kind = 4 ) state character :: TAB = achar ( 9 ) integer ( kind = 4 ) value value = 0 ierror = 0 length = 0 state = 0 isgn = 1 do i = 1, len_trim ( s ) c = s(i:i) ! ! STATE = 0, haven't read anything. ! if ( state == 0 ) then if ( c == ' ' .or. c == TAB ) then else if ( c == '-' ) then state = 1 isgn = -1 else if ( c == '+' ) then state = 1 isgn = +1 else if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then state = 2 value = iachar ( c ) - iachar ( '0' ) else ierror = 1 return end if ! ! STATE = 1, have read the sign, expecting digits or spaces. ! else if ( state == 1 ) then if ( c == ' ' .or. c == TAB ) then else if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then state = 2 value = iachar ( c ) - iachar ( '0' ) else ierror = 1 return end if ! ! STATE = 2, have read at least one digit, expecting more. ! else if ( state == 2 ) then if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then value = 10 * value + iachar ( c ) - iachar ( '0' ) else value = isgn * value ierror = 0 length = i - 1 return end if end if end do ! ! If we read all the characters in the string, see if we're OK. ! if ( state == 2 ) then value = isgn * value ierror = 0 length = len_trim ( s ) else value = 0 ierror = 1 length = 0 end if return end