Computing PI using Parallel Code

The first step is to identify a chunk of code that is doing the big part of the computation...

In this elementary example it is apparent that most of the work is done by the for-loop. Therefore, the for-loop is the segment of code we should divide up for simultaneous work by separate processors.

The lower and upper limits, and/or the setpsize, specified for the for-loop can be made to be functions of each processor's ID number, myrank , and the total number of processors, numProcs .

A nice way of partitioning the data is to rewite the for-loop limits and stride as follows:

  Serial code:   for ( interval = 0; interval < numberOfIntervals; interval++ )
      ==>
  Parallel code: for ( interval = myrank; interval < numberOfIntervals; interval += numProcs )       

You may want to experiment with these limits and stride to verify that indeed all the intervals are distributed among the processors as evenly as possible and that no two processors are computing the same interval. This kind of clever use of myrank and numProcs is characteristic of data decomposition in parallel programming.

Anyway, before we go into more details, here's the parallel version of our program to compute PI, written in C:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <mpi.h>

const double pi = 3.1415926535897932385;

int main( int argc, char *argv[] ) {
  /*  This is a parallelization of serialPi.c.                             */
  /*  The for-loop is partitioned among the processors.                    */
  /*  Function MPI_Reduce sums the processors' results for processor 0.    */
  /*  Processor 0 outputs.                                                 */
  /*                                                                       */
  /*  Remember, it is convenient to reference the typical processor in the */
  /*  first person, using "I" and "me" and "my." In other words,           */
  /*  read through this program as if you were one of the processors.      */

  /*  Arguments required for executing argv[0]:      */
  /*  0                         -> serialPi          */
  /*  1                         -> numberOfIntervals */

  /*  My ID number and the total number of processors: */
  int myrank, numProcs;

  /*  Variables for the computation. */
  int numberOfIntervals, interval;
  double intervalWidth, intervalMidPoint, totalArea, myArea = 0.0;

  /*  For completeness I'll declare an MPI status variable, */
  /*  although I don't plan to use it. */
  MPI_Status status;

  /*  Initialize the MPI API: */
  MPI_Init(&argc, &argv);

  /*  Request my ID number: */
  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);

  /*  I'll also ask how many other processors are out there: */
  MPI_Comm_size(MPI_COMM_WORLD, &numProcs);

  /*  Okay. The preparations have been made. */

  /*  The number of intervals is a command line argument. */
  numberOfIntervals = atoi( argv[1] );

  /*  Compute the interval width. */
  intervalWidth = 1.0 / numberOfIntervals;

  /*  Now I'll compute my area. */
  for ( interval = myrank; interval < numberOfIntervals; interval += numProcs )    {
      intervalMidPoint = (interval + 0.5) * intervalWidth;
      myArea += 4.0 / ( 1.0 + intervalMidPoint*intervalMidPoint );
  }

  /*  Okay. Now I'll submit my results to the MPI API. */
  /*  If I'm processor 0, I'll receive the total area. */
  /*  Otherwise, I'm done. */
  MPI_Reduce( &myArea,
              &totalArea,
	      1,
	      MPI_DOUBLE,
	      MPI_SUM,
   	      0,
 	      MPI_COMM_WORLD );

  /*  If I'm processor 0, I need to normalize the total area and output. */
  if (myrank == 0)    {
      totalArea *= intervalWidth;

      printf( "The computed value of the integral is %.15f\n", totalArea );
      printf( "The exact value of the integral is    %.15f\n", pi );
  }

  /*  Close the MIP API: */
  MPI_Finalize();

  return 0;
}