0

I tried to use OpenMP in my C program for creating the Mandelbrot set. I use two functions f(z) and d(z) defined in the file. When I use them inside a parallel section direct code:

dc = 5*z*z*z*z*dc + 1; 
z = z*z*z*z*z +c;

it works. When I use functions:

dc = d(z)*dc +1;
z = f(z);

It works without OpenMP (good image) but not with OpenMP (bad image).

How should it be done?

See below full code:

#include <complex.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <omp.h> //OpenM

/*
fork of 
mandelbrot-book how to write a book about the Mandelbrot set by Claude Heiland-Alle
https://code.mathr.co.uk/mandelbrot-book/blob/HEAD:/book/


gcc e.c -lm -Wall -fopenmp

./a.out >ed.ppm   // P6 = binary Portable PixMap see https://en.wikipedia.org/wiki/Netpbm#File_formats


*/

const double pi = 3.141592653589793;

double _Complex c;
double _Complex z;
double _Complex dc;
/*
 int q = 5 ;
complex double f(complex double z, int q){ return cpow(z,q) + c;}
complex double d(complex double z, int q) {return q*cpow(z, q-1); }

*/
complex double f(complex double z){ return z*z*z*z*z + c;}
complex double d(complex double z) {return 5*z*z*z*z; }




double cnorm(double _Complex z) // https://stackoverflow.com/questions/6363247/what-is-a-complex-data-type-and-an-imaginary-data-type-in-c
{
  return creal(z) * creal(z) + cimag(z) * cimag(z);
}

void hsv2rgb(double h, double s, double v, int *red, int *grn, int *blu) {
  double i, f, p, q, t, r, g, b;
  int ii;
  if (s == 0.0) { r = g = b = v; } else {
    h = 6 * (h - floor(h));
    ii = i = floor(h);
    f = h - i;
    p = v * (1 - s);
    q = v * (1 - (s * f));
    t = v * (1 - (s * (1 - f)));
    switch(ii) {
      case 0: r = v; g = t; b = p; break;
      case 1: r = q; g = v; b = p; break;
      case 2: r = p; g = v; b = t; break;
      case 3: r = p; g = q; b = v; break;
      case 4: r = t; g = p; b = v; break;
      default:r = v; g = p; b = q; break;
    }
  }
  *red = fmin(fmax(255 * r + 0.5, 0), 255);
  *grn = fmin(fmax(255 * g + 0.5, 0), 255);
  *blu = fmin(fmax(255 * b + 0.5, 0), 255);
}

int main()
{
  int aa = 4;
  int w = 800 * aa;
  int h = 800 * aa;
  int n = 1024;
  double r = 2;
  double px = r / (h/2);
  double r2 = 25 * 25;
  unsigned char *img = malloc(3 * w * h);
  int i,j;
  #pragma omp parallel for //schedule(dynamic) private (c,dc,i,j,z) shared(w,h,n,r,px,r2)
  
  for ( j = 0; j < h; ++j)
  {
    double y = (h/2 - (j + 0.5)) / (h/2) * r;
    for (i = 0; i < w; ++i)
    {
      double x =  (i + 0.5 - w/2) / (h/2) * r; // for q=2 add -0.5
      c = x + I * y;
      //double _Complex 
      dc = 0; // first derivative of zn with respect to c
      //double _Complex z = 0;
      z = 0;
      int k;
      for (k = 0; k < n; ++k)
      { 
      
        //complex double temp = z*z*z*z; // optimisation ?
        
        // works for openmp
        //dc = 5*z*z*z*z*dc + 1; 
        //z = z*z*z*z*z +c;
        
        
        // not works for openmp
        dc = d(z)*dc +1;
        z = f(z);
        
        if (cnorm(z) > r2)
          break;
      }
      
      // color
      double hue = 0, sat = 0, val = 1; // interior color = white
      
      if (k < n) 
      { // exterior and boundary color
        double _Complex de = 2 * z * log(cabs(z)) / dc;
        hue = fmod(1 + carg(de) / (2 * pi), 1); // ? slope of de
        sat = 0.25;
        val = tanh(cabs(de) / px / aa);
      }
      
      // hsv to rgb conversion
      int red, grn, blu;
      hsv2rgb(hue, sat, val, &red, &grn, &blu);
      // save rgb color to array
      img[3*(j * w + i)+0] = red;
      img[3*(j * w + i)+1] = grn;
      img[3*(j * w + i)+2] = blu;
    }
  }
  
  //
  printf("P6\n%d %d\n255\n", w, h);
  fwrite(img, 3 * w * h, 1, stdout);
  free(img);
  
  
  return 0;
}
2
  • 2
    You should always define loop variables in the loop header. I suspect that changing for (i = 0 to for (int i=0 will fix your problem. The for (j loop is also bad style, but OpenMP privatizes that one for you. Commented Sep 19, 2023 at 16:30
  • 2
    Why did you comment the schedule(dynamic) private (c,dc,i,j,z) shared(w,h,n,r,px,r2)? It would have avoid the race condition on i (mentioned by Victor Eijkhout), but also the one on dc, z and c which are actually shared global variables... You should avoid global variables like the plague in parallel computing, especially if you write them. Global variables are a serious concern in software engineering. Commented Sep 19, 2023 at 16:42

1 Answer 1

3

Here's a version of your code that works and produces the image correctly:

#include <complex.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <omp.h> //OpenM

/*
fork of 
mandelbrot-book how to write a book about the Mandelbrot set by Claude Heiland-Alle
https://code.mathr.co.uk/mandelbrot-book/blob/HEAD:/book/


gcc e.c -lm -Wall -fopenmp

./a.out >ed.ppm   // P6 = binary Portable PixMap see https://en.wikipedia.org/wiki/Netpbm#File_formats


*/

const double pi = 3.141592653589793;

/*
 int q = 5 ;
complex double f(complex double z, int q){ return cpow(z,q) + c;}
complex double d(complex double z, int q) {return q*cpow(z, q-1); }

*/
complex double f(complex double z, complex double c){ return z*z*z*z*z + c;}
complex double d(complex double z) {return 5*z*z*z*z; }




double cnorm(double _Complex z) // https://stackoverflow.com/questions/6363247/what-is-a-complex-data-type-and-an-imaginary-data-type-in-c
{
  return creal(z) * creal(z) + cimag(z) * cimag(z);
}

void hsv2rgb(double h, double s, double v, int *red, int *grn, int *blu) {
  double i, f, p, q, t, r, g, b;
  int ii;
  if (s == 0.0) { r = g = b = v; } else {
    h = 6 * (h - floor(h));
    ii = i = floor(h);
    f = h - i;
    p = v * (1 - s);
    q = v * (1 - (s * f));
    t = v * (1 - (s * (1 - f)));
    switch(ii) {
      case 0: r = v; g = t; b = p; break;
      case 1: r = q; g = v; b = p; break;
      case 2: r = p; g = v; b = t; break;
      case 3: r = p; g = q; b = v; break;
      case 4: r = t; g = p; b = v; break;
      default:r = v; g = p; b = q; break;
    }
  }
  *red = fmin(fmax(255 * r + 0.5, 0), 255);
  *grn = fmin(fmax(255 * g + 0.5, 0), 255);
  *blu = fmin(fmax(255 * b + 0.5, 0), 255);
}

int main()
{
  const int aa = 4;
  const int w = 800 * aa;
  const int h = 800 * aa;
  const int n = 1024;
  const double r = 2;
  const double px = r / (h/2);
  const double r2 = 25 * 25;
  unsigned char *img = malloc(3 * w * h);

  #pragma omp parallel for schedule(dynamic)
  for (int j = 0; j < h; ++j)
  {
    double _Complex c;
    double _Complex z;
    double _Complex dc;
    double y = (h/2 - (j + 0.5)) / (h/2) * r;
    for (int i = 0; i < w; ++i)
    {
      double x =  (i + 0.5 - w/2) / (h/2) * r; // for q=2 add -0.5
      c = x + I * y;
      //double _Complex 
      dc = 0; // first derivative of zn with respect to c
      //double _Complex z = 0;
      z = 0;
      int k;
      for (k = 0; k < n; ++k)
      { 
      
        //complex double temp = z*z*z*z; // optimisation ?
        
        // works for openmp
        //dc = 5*z*z*z*z*dc + 1; 
        //z = z*z*z*z*z +c;
        
        
        // not works for openmp
        dc = d(z)*dc +1;
        z = f(z,c);
        
        if (cnorm(z) > r2)
          break;
      }
      
      // color
      double hue = 0, sat = 0, val = 1; // interior color = white
      
      if (k < n) 
      { // exterior and boundary color
        double _Complex de = 2 * z * log(cabs(z)) / dc;
        hue = fmod(1 + carg(de) / (2 * pi), 1); // ? slope of de
        sat = 0.25;
        val = tanh(cabs(de) / px / aa);
      }
      
      // hsv to rgb conversion
      int red, grn, blu;
      hsv2rgb(hue, sat, val, &red, &grn, &blu);
      // save rgb color to array
      img[3*(j * w + i)+0] = red;
      img[3*(j * w + i)+1] = grn;
      img[3*(j * w + i)+2] = blu;
    }
  }
  
  //
  printf("P6\n%d %d\n255\n", w, h);
  fwrite(img, 3 * w * h, 1, stdout);
  free(img);
  
  
  return 0;
}

The important changes is that I declared i, j, c, z, dc withing the parallel region and passed c to the function f. Why this works now, and didn't before, is a bit of a mystery, but I assume that, because c, z, dc were global, they could not be made private.

Sign up to request clarification or add additional context in comments.

Comments

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.