Need correction to my OpenMP Mandelbrot set code

Question

I have the following Mandelbrot set code in OpenMP. My C code works just fine, and the picture that it produces is perfect. But using OpenMP, it compiles and runs correctly, but unfortunately I am not able to open the output .ppm file, simply Gimp cannot read it.

// mandopenmp.c
// to compile: gcc -fopenmp mandopenmp.c -o mandopenmp -lm
// usage: ./mandopenmp <no_of_iterations> > output.ppm

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

typedef struct {
    int r, g, b;
} rgb;


void color(rgb **m, int x, int y, int red, int green, int blue)
{
    m[x][y].r = red;
    m[x][y].g = green;
    m[x][y].b = blue;
}

void mandelbrot(int niterations, rgb **m)
{
    int w = 600, h = 400, x, y, i;
    // each iteration, it calculates: newz = oldz*oldz + p, 
    // where p is the current pixel, and oldz stars at the origin
    double pr, pi;                   // real and imaginary part of the pixel p
    double newRe, newIm, oldRe, oldIm;   // real and imaginary parts of new and old z
    double zoom = 1, moveX = -0.5, moveY = 0; // you can change these to zoom and change position

    printf("P6\n# AUTHOR: Erkan Tairi\n");
    printf("%d %d\n255\n",w,h);

    //loop through every pixel
    #pragma omp parallel for private(x,i,pr,pi,newRe,newIm,oldRe,oldIm) schedule(dynamic, 1)
    for(y = 0; y < h; y++) {
        for(x = 0; x < w; x++) {
            // calculate the initial real and imaginary part of z, 
            // based on the pixel location and zoom and position values
            pr = 1.5 * (x - w / 2) / (0.5 * zoom * w) + moveX;
                pi = (y - h / 2) / (0.5 * zoom * h) + moveY;
                newRe = newIm = oldRe = oldIm = 0; //these should start at 0,0
                // start the iteration process
                for(i = 0; i < niterations; i++) {
                        // remember value of previous iteration
                        oldRe = newRe;
                        oldIm = newIm;
                        // the actual iteration, the real and imaginary part are calculated
                        newRe = oldRe * oldRe - oldIm * oldIm + pr;
                        newIm = 2 * oldRe * oldIm + pi;
                        // if the point is outside the circle with radius 2: stop
                        if((newRe * newRe + newIm * newIm) > 4) break;
                }
                if(i == niterations)
                color(m, x, y, 0, 0, 0); // black
            else
            {
                // normalized iteration count method for proper coloring
                double z = sqrt(newRe * newRe + newIm * newIm);
                int brightness = 256. * log2(1.75 + i - log2(log2(z))) / log2((double)niterations);
                color(m, x, y, brightness, brightness, 255);
            }
            }
    }
}

int main(int argc, char *argv[])
{
    int niterations, i, j;

    if(argc != 2)
    {
        printf("Usage: %s <no_of_iterations> > output.ppm\n", argv[0]);
        exit(1);
    }

    niterations = atoi(argv[1]);

    rgb **m;
    m = malloc(600 * sizeof(rgb *));
    for(i = 0; i < 600; i++)
        m[i] = malloc(400 * sizeof(rgb));

    double begin = omp_get_wtime();
    mandelbrot(niterations, m);

    for(i = 0; i < 600; i++) {
        for(j = 0; j < 400; j++) {
            fputc((char)m[i][j].r, stdout);
            fputc((char)m[i][j].g, stdout);
            fputc((char)m[i][j].b, stdout);
        }
    }

    double end = omp_get_wtime();

    double time_spent = end - begin;
    fprintf(stderr, "Elapsed time: %.2lf seconds.\n", time_spent);

    for(i = 0; i < 600; i++)
        free(m[i]);
    free(m);

    return 0;
}

Answer 1

I don't know the internals of the Mandrelbot set, but I will give a shot based on your program workflow.

Probably it is because you are writing the color to your output file while in the parallel section. What this this means is that your pixels are being written as the computing process finishes, but this does not mean that the computing process of pixel X will end before the processing of pixel X+1.

This way, while writing to the file, you will end up writing first (for example) pixel X+1 and then pixel X, mixing the colors.

Try writing the output result to a matrix. You will have to change your color function, adding two parameters i and j with the coordinates of the pixel to be written.

After the whole processing finishes and every pixel is computed, then will should write the pixels of the matrix to the output file.

The code:

typedef struct {
    int r, g, b;
} rgb;

void color(rgb **m, int x, int y, int red, int green, int blue) {
    m[x][y].r = red;
    m[x][y].g = green;
    m[x][y].b = blue;
}

void mandelbrot(rgb **m, int niterations) { // note the new argument, m.
    // and your code goes on and on... until:
            if ( i == niterations )
                color(m, x, y, 0, 0, 0);
            else {
                // normalized iteration count method for proper coloring
                double z = sqrt(newRe * newRe + newIm * newIm);
                int brightness = 256. * log2(1.75 + i - log2(log2(z))) / log2((double)niterations);
                color(m, x, y, brightness, brightness, 255);
            }
        }
    }
}

int main(int argc, char *argv[]) {
    // everything ok until...

    double begin = omp_get_wtime();

    rgb **m;
    m = malloc(sizeof(rgb*) * 600);
    for ( i = 0; i < 600; i++ ) {
        m[i] = malloc(400 * sizeof(rgb));

    // finally call mandelbrot!
    mandelbrot(m, niterations);
    double end = omp_get_wtime();

    // now that you have computed your set, you just walk the array writing the output to the file.

    for ( i = 0; i < 600; i++ ) {
        free(m[i]);
    }
    free(m);

    double time_spent = end - begin;
    fprintf(stderr, "Elapsed time: %.2lf seconds.\n", time_spent);

    return 0;
}

Answer 2

Your implementation is flawed. You have declared many variables that have to be private to be shared instead. This includes pr, pi, newRe, newIm. Also oldRe and oldIm are shared by default as they are declared in a scope that is outer to the parallel region. These all should be private instead:

#pragma omp parallel for private(x,i,pr,pi,newRe,newIm,oldRe,oldIm)

Also the default scheduling for parallel for loops is often (but not necessarily always) static. This is not the optimal one for things like fractals as it takes different time to compute each line or column in the image. Therefore you should apply the schedule(dynamic,1) clause and play with the chunk size (1 in this case) until you get the best speed-up.

#pragma omp parallel for private(x,i,pr,pi,newRe,newIm,oldRe,oldIm) \
            schedule(dynamic,1)

Answer 3

If you're writing to a file sequentially (which you were doing in your original code before you edited it) then you can use the ordered pragma before you write to the file. This will make the image correct using your original code. See the following link http://bisqwit.iki.fi/story/howto/openmp/#ExampleCalculatingTheMandelbrotFractalInParallel

However, it's not the optimal solution. The optimal solutions is to write to a memory buffer first and after the mandelbrot code finishes filling the buffer then write to a file (as you do in your new code).

I have a few suggestions to speed your code up. Fuse your x and y loop (as shown in the link) and use schedule dynamic (also shown in that link) since each pixel takes different time. Lastly, use SSE/AVX to operate on two (SSE) or four (AVX) pixels at once. In total you should get a speed up of over 20x with OpenMP and SSE/AVX.

#pragma omp ordered {
    if(i == niterations)
        color(m, x, y, 0, 0, 0); // black - use original color function which writes to file
    else
    {
        // normalized iteration count method for proper coloring
        double z = sqrt(newRe * newRe + newIm * newIm);
        int brightness = 256. * log2(1.75 + i - log2(log2(z))) / log2((double)niterations);
        color(m, x, y, brightness, brightness, 255); //use original color function which writes to file
    }
}