Program Example1_2
c#######################################################################
c#
c# This is a SPMD MPI example on parallel integration
c# It demonstrates the use of :
c#
c# * <A HREF="http://www-unix.mcs.anl.gov/mpi/www/www3/MPI_Init.html">MPI_Init</A>
c# * <A HREF="http://www-unix.mcs.anl.gov/mpi/www/www3/MPI_Comm_rank.html">MPI_Comm_rank</A>
c# * <A HREF="http://www-unix.mcs.anl.gov/mpi/www/www3/MPI_Comm_size.html">MPI_Comm_size</A>
c# * <A HREF="http://www-unix.mcs.anl.gov/mpi/www/www3/MPI_Recv.html">MPI_Recv</A>
c# * <A HREF="http://www-unix.mcs.anl.gov/mpi/www/www3/MPI_Send.html">MPI_Send</A>
c# * <A HREF="http://www-unix.mcs.anl.gov/mpi/www/www3/MPI_Finalize.html">MPI_Finalize</A>
c#
c#######################################################################
      implicit none
      integer n, p, i, j, k, ierr, master
      real h, a, b, integral, pi

      include "mpif.h"  !! This brings in pre-defined MPI constants, ...
      integer myid, source, dest, tag, status(MPI_STATUS_SIZE)
      real my_result, Total_result, result

      data master/0/

c**Starts MPI processes ...

      call MPI_Init(ierr)                             !! starts MPI
      call MPI_Comm_rank(MPI_COMM_WORLD, myid, ierr)  !! get current proc id
      call MPI_Comm_size(MPI_COMM_WORLD, p, ierr)     !! get number of procs

      pi = acos(-1.0)   !!  = 3.14159...
      a = 0.0           !! lower limit of integration
      b = pi/2.         !! upper limit of integration
      n = 500           !! number of increment within each process
      dest = master     !! define the process that computes the final result
      tag = 123         !! set the tag to identify this particular job
      h = (b-a)/n/p     !! length of increment

      my_result = integral(a,myid,h,n)
      write(*,*)'myid=',myid,',  my_result=',my_result

      if(myid .eq. master) then               ! the following is serial
        result = my_result
        do k=1,p-1
          call MPI_Recv(my_result, 1, MPI_REAL, 
     &       MPI_ANY_SOURCE, tag,    !! more efficient, less prone to deadlock
     &       MPI_COMM_WORLD, status, ierr) !! root receives my_result from proc
          result = result + my_result
        enddo
      else
        call MPI_Send(my_result, 1, MPI_REAL, dest, tag,
     &     MPI_COMM_WORLD, ierr)     !! send my_result to intended dest.
      endif

c**results from all procs have been collected and summed ...

      if(myid .eq. 0) then
        write(*,*)'Final Result =',result
      endif
      call MPI_Finalize(ierr)                       !! let MPI finish up ...

      stop
      end
      real function integral(a,i,h,n)
      implicit none
      integer n, i, j
      real h, h2, aij, a
      real fct, x

      fct(x) = cos(x)               !! kernel of the integral
      integral = 0.0                !! initialize integral
      h2 = h/2.
      do j=0,n-1                    !! sum over all "j" integrals
        aij = a + (i*n +j)*h        !! lower limit of "j" integral
        integral = integral + fct(aij+h2)*h
      enddo

      return
      end