MPI:Exemplo Trapezoidal

program trapezoidal
implicit none
include 'mpif.h'

integer ::      my_rank, p,local_n, source,k
integer ::      status(MPI_STATUS_SIZE)
integer ::      ierr
real*8  ::      h,local_a,local_b,integral,total,Trap

!  Wall Time Declarations:
real*8  ::      inicio,fim, totalwalltime, traperror

real*8  ::      a   = 0.d0,  b   = 1.d0
integer ::      n   = 1048576
integer ::      dest= 0
integer ::      tag = 0
integer ::      repeticoes =1000

call MPI_Init(ierr)
call MPI_Comm_rank(MPI_COMM_WORLD, my_rank, ierr)
call MPI_Comm_size(MPI_COMM_WORLD, p, ierr)

if (my_rank == 0) then
  print *,'Processos:     p=',p
  print *,'Repeticoes:    r=',repeticoes
  print *,'Num.Trapezios: n=',n
  print *,'Intervalo: (a,b)=(',a,b,')'
endif

inicio   = MPI_WTIME(ierr)

do k=1,repeticoes
 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 == 0) then
   total = integral
   do source = 1, p-1
     call MPI_RECV(integral, 1, MPI_DOUBLE_PRECISION, source, tag, MPI_COMM_WORLD, status, ierr)
     total = total + integral
   enddo
 else
   call MPI_SEND(integral, 1, MPI_DOUBLE_PRECISION, dest, tag, MPI_COMM_WORLD, ierr)
 endif
enddo

fim = MPI_WTIME(ierr)

if (my_rank == 0) then
  traperror    = total  - 1.0d0
  totalwalltime= fim - inicio

  print *,'Estimativa=:',total
  print *,'Erro      =',traperror
  print *,'Walltime  =',totalwalltime,' s '
endif

call MPI_FINALIZE(ierr)
end program
!----------------------------------------------------
real*8 function Trap(local_a, local_b, local_n, h)
real*8  ::      local_a,local_b,h,integral,x,f
integer ::      local_n,i

integral = (f(local_a) + f(local_b))/2.0
x = local_a
do i = 1, local_n-1
   x = x + h
   integral = integral + f(x)
enddo
Trap = integral*h
end function

!----------------------------------------------------
real*8 function f(x)
real*8  ::      x

f = 5*x*x*x*x
end function