Pi Computation

• $\pi$ by numerically evaluating the integral
  $$\int_0^1\frac{1}{1+x^2}dx=\frac{\pi}{4}$$

• Midpoint Rule for $\int_a^b f(x)dx \approx (b-a)f(x_m)$

  Figure 4.4: Midpoint Rule.

  
• Midpoint Rule becomes
  $$\int_0^1\frac{1}{1+x^2}dx \approx \sum_{i=1}^n \frac{1}{1+\left( \frac{i-0.5}{n}\right)^2}$$

Sequential Code:

    1 #include <stdio.h>
    2 #include <math.h>
    3 int main(int argc, char* argv[])
    4 {
    5     int done = 0, n, i;
    6     double PI25DT = 3.141592653589793238462643;
    7     double mypi, h, sum, x;
    8     while (!done)
    9     {
   10         printf("Enter the number of intervals: (0 quits) ");
   11         scanf("%d",&n);
   12 		if (n == 0) break; /* Quit when "0" entered*/
   13         /* Integral limits are from 0 to 1 */            
   14         h   = (1.0-0.0)/(double)n; /* Step length*/
   15         sum = 0.0;           /* Initialize sum variable */
   16         /* loop over interval for integration*/
   17         for (i = 1; i <= n; i += 1) 
   18         {
   19             x = h * ((double)i - 0.5); /* Middle point at step */
   20             sum += 4.0 / (1.0 + x*x);  /* Sum up at each step */
   21 //("i=%d x=%f sum=%f \n",i,x,sum); /* print intermediate steps */
   22         }
   23         mypi = h * sum; /* Obtain resulting pi number */
   24         printf("pi is approximately %.16f, Error is %.16f\n",mypi, \\
   25         fabs(mypi - PI25DT));
   26     }
   27 }

Figure 4.5: Sequential Code Output.

• Parallel Code:
  • The master process reads number of intervals from standard input, this number is then sent to the processes.
  • Having received the number of intervals, each process evaluates the total area of n/size rectangles under the curve.
  • The contributions to the total area under the curve are collected from participating processes by the master process, which then adds them up, and prints the result on standard output.

    1 #include <stdio.h>
    2 #include <math.h>
    3 #include "mpi.h"
    4 
    5 int main(int argc, char* argv[])
    6 {
    7   int done = 0, n, i;
    8   double PI25DT = 3.141592653589793238462643;
    9   double mypi, h, sum, x;
   10   int size, rank, me;
   11   int tag=11;         
   12   MPI_Status  status; 
   13   double mysum;       
   14   double pi;          
   15     
   16   MPI_Init(&argc, &argv);              /* Initialize MPI */          
   17   MPI_Comm_size(MPI_COMM_WORLD, &size);/* Get number of processes */
   18   MPI_Comm_rank(MPI_COMM_WORLD, &rank);/* Get own identifier */      
   19     
   20   while (!done)
   21     {
   22       if (rank == 0) { /* Process 0 does this */ 
   23 	printf("Enter the number of intervals: (0 quits) ");
   24 	scanf("%d",&n);
   25 	/* Send a message containing number of intervals to all other processes */
   26 	for (i=1; i<size; i++) {
   27 	  MPI_Send(&n, 1, MPI_INT, i, tag, MPI_COMM_WORLD); /* Blocking send */ 
   28 	}
   29 	if (n == 0) break; /* Quit when "0" entered */
   30 	/* Computing local pi number for rank 0 process*/
   31 	/* Integral limits are from 0 to 1 */            
   32 	h   = (1.0-0.0)/(double)n; /* Step length*/       
   33 	mysum = 0.0; /* Initialize sum variable */        
   34 	for (i = rank+1; i <= n; i += size) /* Loop over interval for integration */ 
   35 	  {
   36 	    x = h * ((double)i - 0.5);  /* Middle point at step */
   37 	    mysum += 4.0 / (1.0 + x*x); /* Sum up at each step */
   38 	    //printf("i=%d x=%f sum=%f \n",i,x,sum); /* Intermediate steps */
   39 	  }
   40 	mypi = h * mysum; /* Obtain local resulting pi number */ 
   41 	/* Receive a message containing local resulting pi number from all other processes */
   42 	for (i=1; i<size; i++) {
   43 	  MPI_Recv (&pi, 1, MPI_DOUBLE, i, tag, MPI_COMM_WORLD, &status); /* Blocking recieve */
   44 	  printf("Process 0 : Received local resulting pi number: %.16f from process %d \n",pi,i);
   45 	  mypi=mypi+pi; /* Reduce all local values to mypi variable */
   46 	}
   47 	printf("pi is approximately %.16f, Error is %.16f\n",mypi, fabs(mypi - PI25DT));  
   48       }
   49       else /* Other processes do this */
   50         { 
   51 	  MPI_Recv (&n, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status); /* Blocking recieve */
   52 	  printf("Process %d : Received number of intervals as %d from process 0 \n",rank, n); 
   53 	  if (n == 0) break; /* Quit when "0" entered*/
   54 	  /* Computing local pi number for other processes*/ 
   55 	  /* Integral limits are from 0 to 1 */             
   56 	  h   = (1.0-0.0)/(double)n; /* Step length*/       
   57 	  mysum = 0.0; /* Initialize sum variable */        
   58 	  for (i = rank+1; i <= n; i += size) /* Loop over interval for integration */ 
   59             {
   60 	      x = h * ((double)i - 0.5); /* Middle point at step */
   61 	      mysum += 4.0 / (1.0 + x*x);  /* Sum up at each step */
   62 	      //printf("i=%d x=%f sum=%f \n",i,x,sum); /* Intermediate steps */
   63             }
   64 	  mypi = h * mysum; /* Obtain local resulting pi number */ 
   65 	  /* Send a message containing local resulting pi number to master processes */
   66 	  MPI_Send(&mypi, 1, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD); /* Blocking send */
   67         }
   68     }
   69   MPI_Finalize();   
   70 }

Figure 4.6: Parallel Code Output.