Mandelbrot image generator 2.0

Question

I have rewritten my C Mandelbrot set image generator, including many suggestions given in my previous question (see: Mandelbrot image generator), as well as some recomendations by a friend.

The total list of additions is as follows:

• Parallelisation (omp.h)
• typedef-ing of structures
• Implementing a better file-type (P6 PPM instead of P3 PPM)
• Extensive use of command line options and arguments (getopt)
• Testing of user-defined regions of the complex plane
• Various bug fixes (possible data leak, weird image dimensions, etc.)
• Generally better coding and layout (more concise and easier to read/understand)

My question is the same as in my last post, what things have I done well here, what things have I done badly, and how can I improve in the future? I would appreciate all and any feedback on both my use of C, and my implementation of the algorithm itself.

Note: for anyone who compiles this code, I have had best results using gcc's -O3 and -fopenmp flags

//  mandelbrot.c - generates a .PPM (Portable Pixmap format, P6) file of the mandelbrot set with shading

//  Options:
//  -f [output file name] (required)
//  -h [image height in px] (required)
//  -t [max imaginary component]
//  -b [min imaginary component]
//  -v [max real component]
//  -n [min real component]

//  Return/exit codes:
//   0 - successful execution
//  -1 - argument error
//  -2 - file access error
//  -3 - image height error

#include <stdio.h>
#include <complex.h>
#include <math.h>
#include <stdlib.h>
#include <omp.h>
#include <ctype.h>
#include <unistd.h>
#define MAX_TESTS 2000

typedef struct
{
    unsigned int height;
    unsigned int width;
    double ymax, ymin;
    double xmax, xmin;
    char *file_name;
} image_meta;

typedef struct
{
    unsigned char red;
    unsigned char green;
    unsigned char blue;
} colour;

short mandelbrot_test(double complex c); //Calculates the number of iterations of the algorithm required for a given complex number to reach a magnitude >= 2
colour rgb_gen(short iterations); //Generates an RGB value based on the rate of divergence
image_meta image_meta_gen(int argc, char *argv[]); //Generates the meta data for the image by parsing the input arguments

int main(int argc, char *argv[])
{
    image_meta image;
    image = image_meta_gen(argc, argv);
    printf("Image dimensions: %dx%d\n", image.width, image.height);

    int xpx, ypx;
    double a, b, xdiff, ydiff;
    double complex num;
    FILE *file;

    xdiff = image.xmax - image.xmin;
    ydiff = image.ymax - image.ymin;

    if((file = fopen(image.file_name, "w")) != NULL)
    {
        fprintf(file, "P6 %d %d 255\n", image.width, image.height);
        colour rgb[image.width];

        #pragma omp parallel
        for(ypx = 0; ypx < image.height; ypx++)
        {
            for(xpx = 0; xpx < image.width; xpx++)
            {
                a   = image.xmin + xpx * xdiff / image.width;
                b   = image.ymax - ypx * ydiff / image.height;
                num = a + b * I;
                rgb[xpx] = rgb_gen(mandelbrot_test(num));
            }
            fwrite(rgb, sizeof(colour), image.width, file);
        }
        fclose(file);
    }
    else
    {
        fprintf(stderr, "Unable to access file!\n");
        exit(-2);
    }
    exit(0);
}

short mandelbrot_test(double complex c)
{
    double complex x = 0;
    double abs       = c * conj(c);

    if(abs * (8.0 * abs - 3.0) < 3.0 / 32.0 - creal(c)) //Quick test to see if we can bail out early by checking if the number lies within the main cardioid
    {
        return MAX_TESTS;
    }

    for(int i = 1; i < MAX_TESTS; i++)
    {
        x *= x;
        x += c;

        if(cabs(x) >= 2)
        {
            return i;
        }
    }
    return MAX_TESTS;
}

colour rgb_gen(short iterations)
{
    colour rgb;
    int brightness;

    if(iterations == MAX_TESTS)
    {
        rgb.red     = 255;
        rgb.green   = 255;
        rgb.blue    = 255;
    }
    else
    {
        brightness  = 256.0 * log2(iterations) / log2(MAX_TESTS - 1);
        rgb.red     = brightness;
        rgb.green   = brightness;
        rgb.blue    = 255;
    }
    return rgb;
}

image_meta image_meta_gen(int argc, char *argv[])
{
    int c;
    image_meta image;

    image.file_name = NULL;
    opterr = 0;
    image.xmax = image.ymax = image.xmin = image.ymin = -1;
    image.height = 0;

    while((c = getopt(argc, argv, "h:f:t:b:v:n:")) != -1)
    {
        switch(c)
        {
            case 'h':
                image.height = atoi(optarg);
                break;
            case 't':
                image.ymax = atof(optarg);
                break;
            case 'b':
                image.ymin = atof(optarg);
                break;
            case 'v':
                image.xmax = atof(optarg);
                break;
            case 'n':
                image.xmin = atof(optarg);
                break;
            case 'f':
                image.file_name = optarg;
                break;
            case '?':
                if(optopt == 'f' || optopt == 'h')
                {
                    fprintf(stderr, "Option -%c requires an argument.\n", optopt);
                }
                else if(isprint(optopt))
                {
                    fprintf(stderr, "Unknown option `-%c'.\n", optopt);
                }
                else
                {
                    fprintf(stderr, "Unknown option character `\\x%x'.\n", optopt);
                }
                fprintf(stderr, "Usage: ./mandelbrot -f file_name -h height [options]\nOptional: -t y-max -b y-min -v x-min -n x-max\n");
                exit(-1);
        }
    }

    if(argc < 4)
    {
        fprintf(stderr, "Error: too few args\nUsage: ./mandelbrot -f file_name -h height [options]\nOptional: -t y-max -b y-min -v x-min -n x-max\n");
        exit(-1);
    }

    if(image.height < 30)
    {
        fprintf(stderr, "Height can't be less than 30!\n");
        exit(-3);
    }

    if(image.xmax == image.xmin) 
    {
        image.xmax = 0.8;
        image.xmin = -2.0;
        printf("Using default x values...\n");
    }

    if(image.ymax == image.ymin) 
    {
        image.ymax = 1.2;
        image.ymin = -1.2;
        printf("Using default y values...\n");
    }

    if(image.xmin > image.xmax)
    {
        double temp = image.xmin;
        image.xmin  = image.xmax;
        image.xmax  = temp;
    }

    if(image.ymin > image.ymax)
    {
        double temp = image.ymin;
        image.ymin  = image.ymax;
        image.ymax  = temp;
    }

    image.width = image.height * (image.xmax - image.xmin) / (image.ymax - image.ymin);

    return image;
}

Answer 1

Current program not working correctly

I ran your program with and without OpenMP to see how much of a difference it made (I commented out the #pragma to test the slow case). I noticed the following:

1. The two programs outputted different files. Not only different contents but different sizes as well.
2. On a 2000 height image, the slow version ran in 100 seconds. The parallel version ran in 7.5 seconds. That raised a red flag because my machine only has 4 cores so there shouldn't have been a 13x speedup.
3. When I compared the images in an image viewer, at first it seemed that both images were ok. But when I zoomed into the image created by the parallel program, it showed that many of the lines were repeated (i.e. 3-4 lines in a row were the same line).

What is the current program doing?

Let's examine the parallel part:

    #pragma omp parallel
    for(ypx = 0; ypx < image.height; ypx++)
    {
        for(xpx = 0; xpx < image.width; xpx++)
        {
            a   = image.xmin + xpx * xdiff / image.width;
            b   = image.ymax - ypx * ydiff / image.height;
            num = a + b * I;
            rgb[xpx] = rgb_gen(mandelbrot_test(num));
        }
        fwrite(rgb, sizeof(colour), image.width, file);
    }

What this part is doing is it is spawning N threads (4 in my case), and it is running the same code in all of the threads using the same shared variables. Any variable declared before the #pragma is considered a shared variable, so in this case, xpx and ypx are being shared.

Therefore, what is happening is that each of the 4 threads is computing the line at ypx=0. In the xpx loop, the 4 threads all work together to complete that loop 4x as fast. Then all 4 threads are writing the same rgb buffer for ypx=0 out to the file (4 times total), and then each advancing ypx by 1 (+4 total). In other words, they wrote the first line out to the file 4 times in a row and didn't compute lines 1-3 at all.

Since the x and y directions are each being sped up by 4x, this results in a 16x speedup (close to the 13x speedup I observed). But of course, only one out of every 4 lines was actually being generated.

Other findings

Writing to file in parallel

This line:

        fwrite(rgb, sizeof(colour), image.width, file);

shouldn't be included in the parallelized region. Since the multiple threads are all generating data at the same time, calling fwrite() could result in out of order writes to the file. I suggest creating one master buffer for the whole file and then only writing the file out at the end of the parallel execution. Theoretically, you could leave that fwrite() in the parallel part if you synchronized it such that the lines were generated and written out in order. But that would slow down execution of your program.

Scope of variables

One thing that I noticed when I was trying to get your program to work is that all of your variables were declared up top. This actually hindered things when I was fixing the parallel part because all of the variables became shared by default. I had to move the variables to the correct scope to get things to work right (it could also be done using some #pragma options).

Corrected program

Here is the corrected main() function. It ran in 26 seconds, so it got the expected 4x speedup.

int main(int argc, char *argv[])
{
    image_meta image = image_meta_gen(argc, argv);
    FILE *file;

    printf("Image dimensions: %dx%d\n", image.width, image.height);
    if ((file = fopen(image.file_name, "w")) != NULL) {
        colour *rgb = calloc(image.width * image.height, sizeof(colour));
        if (rgb == NULL)
            exit(-1);

        fprintf(file, "P6 %d %d 255\n", image.width, image.height);

        #pragma omp parallel
        {
            double xdiff = image.xmax - image.xmin;
            double ydiff = image.ymax - image.ymin;

            #pragma omp for schedule(dynamic)
            for(int ypx = 0; ypx < image.height; ypx++) {
                for(int xpx = 0; xpx < image.width; xpx++) {
                    double a = image.xmin + xpx * xdiff / image.width;
                    double b = image.ymax - ypx * ydiff / image.height;
                    double complex num = a + b * I;
                    rgb[ypx*image.width + xpx] = rgb_gen(mandelbrot_test(num));
                }
            }
        }
        fwrite(rgb, sizeof(colour), image.width*image.height, file);
        fclose(file);
    } else {
        fprintf(stderr, "Unable to access file!\n");
        exit(-2);
    }
    exit(0);
}

Explanation of program

1. I changed the one line rgb buffer into a whole image buffer. The buffer is written to the file outside of the parallel block.
2. I moved all of the variables into the scope of the parallel block so that none of the variables would be shared. Only image is shared but it is never written to, so it's fine.
3. The #pragma omp for divides the next for loop (the ypx loop) among the threads.
4. The schedule(dynamic) part is an optimization that helps all the threads distribute the work. At first when it wasn't dynamic, the loop was split evenly among the threads by number of lines, but some threads had more work to do than others (because some parts of the mandlebrot image are easier to compute than others). Some threads finished much earlier than others and became idle, while others struggled through the hard parts. The whole thing took 35+ seconds instead of 26 seconds. By using schedule(dynamic), each loop iteration is scheduled one at a time, so threads do not become idle until there is no more work to be done.

Answer 2

This looks pretty good! I've never used OpenMP, so I can't comment on that, and I have to admit this is the first time I've seen the use of complex numbers in C. Very cool!

I don't have a lot to add to the discussion except for this:

If the output image is n x m pixels, then you're currently generating the coordinates (a and b) n*m times. You could calculate them n + m times by calculating all the as outside the loop, and all the bs outside the loop like this:

        for(xpx = 0; xpx < image.width; xpx++)
        {
            as[xpx]   = image.xmin + xpx * xdiff / image.width;
        }

        for(ypy = 0; ypy < image.height; ypy++)
        {
            bs[ypy] = image.ymin + ypy * ydiff / image.height;
        }

Then, in the inner loop, calculate num like this:

num = as[xpx] + I * bs[ypy];

It's a small win, but it's something.

You could also do something similar for the coloring. You could generate a table of MAX_TESTS RGB values and just return the appropriate one in rgb_gen().

If you really want to make things even more parallel, you also have the option of using SIMD instructions for your target CPU. On Intel, that would be MMX/SSE/AVX instructions. On ARM I think it would be NEON. And if you want to go even further you could write OpenCL, CUDA, or OpenGL Computation shaders and do the calculations on the GPU.