Why do I just change the atomic clause to a reduction clause and the result is an error

angelavor 0
Viewed 25

What is the difference between a reduction clause and an atomic clause? Why do I just change the atomic clause to a reduction clause and the result is an error. Here are the source codes.

version using atomic clause

/*
**  PROGRAM: Mandelbrot area
**
**  PURPOSE: Program to compute the area of a  Mandelbrot set.
**           Correct answer should be around 1.510659.
**           WARNING: this program may contain errors
**
**  USAGE:   Program runs without input ... just run the executable
**            
**  HISTORY: Written:  (Mark Bull, August 2011).
**           Changed "comples" to "d_comples" to avoid collsion with 
**           math.h complex type (Tim Mattson, September 2011)
*/

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

# define NPOINTS 1000
# define MAXITER 1000


struct d_complex{
   double r;
   double i;
};
void testpoint(struct d_complex);

int numoutside = 0;

int main(){
   int i;
   double area, error;


//   Loop over grid of points in the complex plane which contains the Mandelbrot set,
//   testing each point to see whether it is inside or outside the set.

#pragma omp parallel for default(shared) 
   for (i=0; i<NPOINTS; i++) {
     struct d_complex c;
     double eps = 1.0e-5;
     for (int j=0; j<NPOINTS; j++) {
       c.r = -2.0+2.5*(double)(i)/(double)(NPOINTS)+eps;
       c.i = 1.125*(double)(j)/(double)(NPOINTS)+eps;
       testpoint(c);
     }
   }

// Calculate area of set and error estimate and output the results
printf("numoutside:%d\n",numoutside);
   
area=2.0*2.5*1.125*(double)(NPOINTS*NPOINTS-numoutside)/(double)(NPOINTS*NPOINTS);
   error=area/(double)NPOINTS;

   printf("Area of Mandlebrot set = %12.8f +/- %12.8f\n",area,error);
   printf("Correct answer should be around 1.510659\n");

}

void testpoint(struct d_complex c){

// Does the iteration z=z*z+c, until |z| > 2 when point is known to be outside set
// If loop count reaches MAXITER, point is considered to be inside the set

       struct d_complex z;
       int iter;
       double temp;

       z=c;
       for (iter=0; iter<MAXITER; iter++){
         temp = (z.r*z.r)-(z.i*z.i)+c.r;
         z.i = z.r*z.i*2+c.i;
         z.r = temp;
         if ((z.r*z.r+z.i*z.i)>4.0) {
           #pragma omp atomic
           numoutside++;
           break;
         }
       }

}

Because numoutside is a shared variable that increments in the parallel region, I added atomic operations to ensure correct results. But, when I try to change atomic to reduction clause, the result is error. And the error result of each execution is different. I know there must be thread race, but I don't know why?

version using reduction clause

/*
**  PROGRAM: Mandelbrot area
**
**  PURPOSE: Program to compute the area of a  Mandelbrot set.
**           Correct answer should be around 1.510659.
**           WARNING: this program may contain errors
**
**  USAGE:   Program runs without input ... just run the executable
**            
**  HISTORY: Written:  (Mark Bull, August 2011).
**           Changed "comples" to "d_comples" to avoid collsion with 
**           math.h complex type (Tim Mattson, September 2011)
*/

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

# define NPOINTS 1000
# define MAXITER 1000


struct d_complex{
   double r;
   double i;
};
void testpoint(struct d_complex);

int numoutside = 0;

int main(){
   int i;
   double area, error;


//   Loop over grid of points in the complex plane which contains the Mandelbrot set,
//   testing each point to see whether it is inside or outside the set.

#pragma omp parallel for default(shared) reduction(+:numoutside)
   for (i=0; i<NPOINTS; i++) {
     struct d_complex c;
     double eps = 1.0e-5;
     for (int j=0; j<NPOINTS; j++) {
       c.r = -2.0+2.5*(double)(i)/(double)(NPOINTS)+eps;
       c.i = 1.125*(double)(j)/(double)(NPOINTS)+eps;
       testpoint(c);
     }
   }

// Calculate area of set and error estimate and output the results
printf("numoutside:%d\n",numoutside);
   
area=2.0*2.5*1.125*(double)(NPOINTS*NPOINTS-numoutside)/(double)(NPOINTS*NPOINTS);
   error=area/(double)NPOINTS;

   printf("Area of Mandlebrot set = %12.8f +/- %12.8f\n",area,error);
   printf("Correct answer should be around 1.510659\n");

}

void testpoint(struct d_complex c){

// Does the iteration z=z*z+c, until |z| > 2 when point is known to be outside set
// If loop count reaches MAXITER, point is considered to be inside the set

       struct d_complex z;
       int iter;
       double temp;

       z=c;
       for (iter=0; iter<MAXITER; iter++){
         temp = (z.r*z.r)-(z.i*z.i)+c.r;
         z.i = z.r*z.i*2+c.i;
         z.r = temp;
         if ((z.r*z.r+z.i*z.i)>4.0) {
           numoutside++;
           break;
         }
       }

}
openmp c
1 Answers

There is an error since testpoint is compiled separately of main and numoutside accessed in testpoint is still the global variable and so it is left unprotected. To fix the problem, you need to pass a numoutside pointer to testpoint. Indeed, under the hood, reduction(+:numoutside) causes the compiler rewrite the code of the OpenMP parallel section so numoutside is typically a thread-local variable then aggregated to the global variable at the end of the section.

Here is a modified version of the code (untested):

// [...] (same)


   #pragma omp parallel for default(shared) reduction(+:numoutside)
   for (i=0; i<NPOINTS; i++) {
     struct d_complex c;
     double eps = 1.0e-5;
     for (int j=0; j<NPOINTS; j++) {
       c.r = -2.0+2.5*(double)(i)/(double)(NPOINTS)+eps;
       c.i = 1.125*(double)(j)/(double)(NPOINTS)+eps;
       testpoint(c, &numoutside);
     }
   }


// [...] (same)


void testpoint(struct d_complex c, int* numoutside){

       // Does the iteration z=z*z+c, until |z| > 2 when point is known to be outside set
       // If loop count reaches MAXITER, point is considered to be inside the set

       struct d_complex z;
       int iter;
       double temp;

       z=c;
       for (iter=0; iter<MAXITER; iter++){
         temp = (z.r*z.r)-(z.i*z.i)+c.r;
         z.i = z.r*z.i*2+c.i;
         z.r = temp;
         if ((z.r*z.r+z.i*z.i)>4.0) {
           (*numoutside)++;
           break;
         }
       }
}

A better approach may be to just return a local version of numoutside in the function testpoint (initialized to 0). You can then just use an instruction like numoutside += testpoint(c); in the main loop.

Note that global variable are generally considered as bad, especially in parallel program and even more when they are mutated (bug prone, cause performance issues and many software engineering problems).

Jérôme Richard0
Related