c  trap.f -- Parallel Trapezoidal Rule, first version
*  Fall97:   Timing, real*8, new f  added to Pacheco's PPMI source
c 
c  Input: None.
c  Output:  Estimate of the integral from a to b of f(x) 
c     using the trapezoidal rule and n trapezoids.
c 
c  Algorithm:
c     1.  Each process calculates "its" interval of 
c         integration.
c     2.  Each process estimates the integral of f(x)
c         over its interval using the trapezoidal rule.
c     3a. Each process != 0 sends its integral to 0.
c     3b. Process 0 sums the calculations received from
c         the individual processes and prints the result.
c 
c  Note:  f(x), a, b, and n are all hardwired.
c
c  See Chap. 4, pp. 56 & ff. in PPMPI.
c
	program trapezoidal
c
	include 'mpif.h'
c
	integer   my_rank
	integer   p
	real*8      a
	real*8      b
	integer   n
	real*8      h
	real*8      local_a
	real*8      local_b
	integer   local_n
	real*8      integral
	real*8      total 
	integer   source
	integer   dest
	integer   tag
	integer   status(MPI_STATUS_SIZE)
        integer   ierr
c
        real*8      Trap
c
*  Wall Time Declarations:
        real*8       start, finish, totalwalltime, traperror
c
	data a, b, n, dest, tag /0.0, 1.0, 1024, 0, 0/
	
	call MPI_INIT(ierr)
	call MPI_COMM_RANK(MPI_COMM_WORLD, my_rank, ierr)
	call MPI_COMM_SIZE(MPI_COMM_WORLD, p, ierr)

* Fall 97:  Wall Timer Added:
        start = MPI_WTIME(ierr)
	
	h = (b-a)/n
	local_n = n/p
	
	local_a = a + my_rank*local_n*h
	local_b = local_a + local_n*h
	integral = Trap(local_a, local_b, local_n, h)
	
	if (my_rank .EQ. 0) then
	    total = integral
	    do 100 source = 1, p-1
	        call MPI_RECV(integral, 1, MPI_DOUBLE_PRECISION, source,
     +                        tag, MPI_COMM_WORLD, status, ierr)
	        total = total + integral
 100        continue
	else
	    call MPI_SEND(integral, 1, MPI_DOUBLE_PRECISION, dest, 
     +          tag, MPI_COMM_WORLD, ierr)
        endif

* Fall 97:  Wall Timer Added:
        finish = MPI_WTIME(ierr)
	
	if (my_rank .EQ. 0) then
            write(6,200) n
 200        format(' ','With n = ',I4,' trapezoids, our estimate')
            traperror = total - 1.0e0
            write(6,300) a, b, total, traperror
 300        format(' ','of the integral on (',f6.2,',',f6.2,
     +             ') = ',f11.5,'; Error =', e12.4)
            totalwalltime=finish - start
            write(6,400) totalwalltime, p
 400        format(' ','Integration Wall Time =',f9.6,' Seconds on'
     &    i4,' Processors')
        endif

	call MPI_FINALIZE(ierr) 
        end 
c	
c	
	real*8 function Trap(local_a, local_b, local_n, h)
	real*8     local_a
	real*8     local_b
	integer  local_n
	real*8     h 
c
	real*8     integral
	real*8     x 
	real*8     i 
c
	real*8     f
c
	integral = (f(local_a) + f(local_b))/2.0 
	x = local_a 
	do 100 i = 1, local_n-1
	   x = x + h 
	   integral = integral + f(x) 
 100    continue
	Trap = integral*h 
        return
        end
c	
c	
	real*8 function f(x)
	real*8 x
	real*8 return_val 
*  Fall 97:  f changes for more computational effort:
        return_val = 5*x*x*x*x
        f = return_val
	return 
	end