Scientific Computing & Visualization | Documentation | Tutorials | Multiprocessing by Message Passing MPI

This Site

Multiprocessing by Message Passing MPI

Example 1.4 More on Integration with MPI Collective Communications

This example demonstrates parallel integration (see Example 1 for detail) via the following MPI library routines. Click on them for details:

Example 1.4 Fortran code

Example 1.4 Fortran Code

      Program Example1_4
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_Gather
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, myid
      real h, integral_sum, a, b, integral, pi, my_int, buf(50)
      include "mpif.h"  ! brings in pre-defined MPI constants, ...
      parameter (master=0)  ! processor computing final sum

c**Starts MPI processes ...

      call MPI_Init(ierr)                               ! starts MPI
      call MPI_Comm_rank(MPI_COMM_WORLD, myid, 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
      n = 500           ! number of intervals in (b-a)/p
      h = (b-a)/n/p     ! length of increment

      my_int = integral(a,myid,h,n)

      write(*,"('Process ',i2,' has the partial sum of',f10.6)")
     &          myid,my_int

      call MPI_Gather(my_int, 1, MPI_REAL, buf, 1, MPI_REAL, master,
     &                  MPI_COMM_WORLD, ierr)

      if(myid .eq. master) then
        integral_sum = 0.0
        do i=1,p
          integral_sum = integral_sum + buf(i)
        enddo
        print *,'The Integral =',integral_sum
      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, fct, x

      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

      real function fct(x)
      implicit none
      real x
      fct = cos(x)
      end

Example 1.4 C code

#include <mpi.h>
#include <math.h>
#include <stdio.h>
float fct(float x)
{
      return cos(x);
}
/* Prototype */
float integral(float a, int i, float h, int n);
int main(int argc, char* argv[])
/*######################################################################
 #
 # This is an MPI example on parallel integration
 # It demonstrates the use of :
 #
 # * MPI_Init
 # * MPI_Comm_rank
 # * MPI_Comm_size
 # * MPI_Gather
 # * MPI_Finalize
 #
 # Dr. Kadin Tseng
 # Scientific Computing and Visualization
 # Boston University
 # 1998
 #
 #######################################################################*/
{
      int n, p, myid, i, master;
      float h, integral_sum, a, b, pi, my_int, buf[50];

/* Starts MPI processes ... */

      MPI_Init(&argc,&argv);                 /* starts MPI */
      MPI_Comm_rank(MPI_COMM_WORLD, &myid);  /* get current process id */
      MPI_Comm_size(MPI_COMM_WORLD, &p);     /* get number of processes */

      master = 0;       /* processor computing total sum */
      pi = acos(-1.0);  /* = 3.14159... */
      a = 0.;           /* lower limit of integration */
      b = pi*1./2.;     /* upper limit of integration */
      n = 500;          /* number of increment within each process */
      h = (b-a)/n/p;    /* length of increment */

      my_int = integral(a,myid,h,n);

      printf("Process %d has the partial sum of %f\n", myid,my_int);

      MPI_Gather(&my_int, 1, MPI_FLOAT, buf, 1, MPI_FLOAT, master, MPI_COMM_WORLD);

      if(myid == master) {
        integral_sum = 0.0;
        for (i=0; i<p; i++) {
          integral_sum += buf[i];
        }
        printf("The Integral =%f\n",integral_sum);
      }

      MPI_Finalize();                       /* let MPI finish up ... */
}

float integral(float a, int i, float h, int n)
{
      int j;
      float h2, aij, integ;

      integ = 0.0;                 /* initialize integral */
      h2 = h/2.;
      for (j=0;j<n;j++) {             /* sum over all "j" integrals */
        aij = a + (i*n + j)*h;        /* lower limit of "j" integral */
        integ += fct(aij+h2)*h;
      }
      return integ;
}

OIT | CCS | July 17, 2008