Computing $\pi$

Computing PI using Sequential Code
An OpenMP version of a threaded program to compute PI number using random numbers.
Computing PI using OpenMP directives - Random Numbers. Follow the steps below for executing OpenMP code;

gcc -o code21 code21.c -fopenmp 
./code21

    1 /**************************************************
    2  An OpenMP version of a threaded program to compute PI.
    3  **************************************************/
    4 #pragma omp parallel default(none) private(rand_no_x, rand_no_y, num_threads, sample_points_per_thread) shared(npoints) reduction(+: sum) num_threads(8)
    5     {
    6         num_threads = omp_get_num_threads();
    7         sample_points_per_thread  = npoints / num_threads;
    8 
    9         sum = 0;
   10         for (int i = 0; i < sample_points_per_thread; i++) {
   11             rand_no_x = (double)rand() / (double)RAND_MAX;
   12             rand_no_y = (double)rand() / (double)RAND_MAX;
   13             if (((rand_no_x - 0.5) * (rand_no_x - 0.5) + (rand_no_y - 0.5) * (rand_no_y - 0.5)) < 0.25) {
   14                 sum = sum + 1;
   15             }
   16         }
   17     }

• Note that this program is much easier to write in terms of specifying creation and termination of threads compared to the corresponding POSIX threaded program.
• The omp_get_num_threads() function returns the number of threads in the parallel region
• The omp_get_thread_num() function returns the integer id of each thread (recall that the master thread has an id 0).
• The parallel directive specifies that all variables except npoints, the total number of random points in two dimensions across all threads, are local.
• Furthermore, the directive specifies that there are eight threads, and the value of sum after all threads complete execution is the sum of local values at each thread.
• A for loop generates the required number of random points (in two dimensions) and determines how many of them are within the prescribed circle of unit diameter.

An OpenMP version of a threaded program t o compute PI number by numerical integration without reduction clause.
Computing PI using OpenMP directives - Numerical Integration. Follow the steps below;

export OMP_NUM_THREADS=4
gcc -o code22 code22.c -fopenmp 
./code22

    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 
   13         if (n == 0) break; /* Quit when "0" entered*/
   14         /* Integral limits are from 0 to 1 */            
   15         h   = (1.0-0.0)/(double)n; /* Step length*/
   16         sum = 0.0;           /* Initialize sum variable */
   17         for (i = 1; i <= n; i += 1) /* loop over interval for integration*/
   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             //printf("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, fabs(mypi - PI25DT));
   25     }
   26 }

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

Access to Sequential Code & OpenMP Parallel Code

Table 5.2: Sequential and OpenMP outputs for computing $\pi$ number.

An OpenMP version of a threaded program to compute PI number by numerical integration with reduction clause.
Computing PI using OpenMP directives - Reduction. Follow the steps below for executing OpenMP code;

export OMP_NUM_THREADS=8
gcc -o code23 code23.c -fopenmp 
./code23

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