Multiprocessing by Message Passing MPI

Example 1.3 Integration with MPI Collective Communications

Two new routines are introduced in this example. First, we wish to share data that is available on a specific process with all participating processes. Rather than setting up a send/recv pair of point-to-point communications between the root process and another process, this can best be achieved by a broadcast routine. In addition, after the local sums have been computed by all participating processes, they are to be summed by the root process to obtain the final sum. This again can be best acheived by a reduction routine. These type of communications that involve multiple processes are called collective communications. As they are expected to perform very specific "group" communication and are implemented by the vendors, they are more efficient as well as more convenient for the application programmers.

Examples of reduction operations are: MPI_SUM (summing), MPI_MAX and MPI_MIN ( find maximum and minimum values) and MPI_MAXLOC and MPI_MINLOC (find locations giving maximum and minimum values).

      Program Example1_3
c#######################################################################
c#
c# This is an MPI example on parallel integration
c# It demonstrates the use of :
c#
c# * MPI_Init
c# * MPI_Comm_rank
c# * MPI_Comm_size
c# * MPI_Bcast
c# * MPI_Reduce
c# * MPI_SUM
c# * MPI_Finalize
c#
c# Dr. Kadin Tseng
c# Scientific Computing and Visualization
c# Boston University
c# 1998
c#
c#######################################################################
      implicit none
      integer n, p, i, j, ierr, master
      real h, result, a, b, integral, pi

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

      data master/0/


c**Starts MPI processes ...

      call MPI_Init(ierr)                               !! starts MPI
      call MPI_Comm_rank(MPI_COMM_WORLD, Iam, ierr)    !! get current process id
      call MPI_Comm_size(MPI_COMM_WORLD, p, ierr)       !! get number of processes

      pi = acos(-1.0)   !!  = 3.14159...
      a = 0.0           !! lower limit of integration
      b = pi*1./2.      !! upper limit of integration

      dest = 0          !! define the process that computes the final result
      tag = 123         !! set the tag to identify this particular job

      if(Iam .eq. master) then
        print *,'The requested number of processors =',p
        print *,'enter number of increments within each process'
        read(*,*)n
      endif

c**Broadcast "n" to all processes

      call MPI_Bcast(n, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)

      h = (b-a)/n/p     !! length of increment

      my_result = integral(a,Iam,h,n)

      write(*,"('Process ',i2,' has the partial result of',f10.6)")
     &          Iam,my_result

      call MPI_Reduce(my_result, result, 1, MPI_REAL, MPI_SUM, dest,
     &                  MPI_COMM_WORLD, ierr)

      if(Iam .eq. master) then
        print *,'The 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