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;
}
for (i = 0tofor (int i=0will fix your problem. Thefor (jloop is also bad style, but OpenMP privatizes that one for you.schedule(dynamic) private (c,dc,i,j,z) shared(w,h,n,r,px,r2)? It would have avoid the race condition oni(mentioned by Victor Eijkhout), but also the one ondc,zandcwhich 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.