Mandelbrot fractal with MPI

Question

I'm not very familiar with MPI. I've written this little piece of code to draw the Mandelbrot fractal processing each row of the image in parallel. It works but it's really slow. I've written code to perform the same task with OpenMP and with no parallelization, and both run instantly with the default parameters, whereas the MPI version runs in about 2s in my computer. Not only that, the OpenMP version used 33% less code and less memory. I suppose some of this improvements come from the fact that using threads for this task might be better (the workload distribution was different with threads too, not one line per thread), but I was wondering if I could improve my MPI code in any way.

/*
 * mandelbrot.c
 *
 * Draws the Mandelbrot fractal using MPI. Optionally write it to a csv file.
 *
 * Compiling instructions:
 * mpicc mandelbrot.c -std=c99 -o output_file
 *
 * Running instructions:
 * mpirun -np NUM_RANKS output_file [FNAME] [W] [H] [MAX_ITER] [X0] [Y0] [XN]
 * [YN]
 *
 * For every integer point in [(0, 0), (W, H)], check if its mapping to
 * [(X0, Y0), (XN, YN)] falls in the Mandelbrot set with at most MAX_ITER.
 * Write the results to a CSV file called FNAME only if the argument FNAME
 * is passed and valid (invalid e.g. on *nix: "/"). For defaults, see main().
 *
 * Note: NUM_RANKS must be >= 2
 *
 * Note: W, H and MAX_ITER are unsigned but bounded by the max value of int.
 * X0, Y0, XN and YN are floats.
 *
 */

#include <complex.h>
#include <stdio.h>
#include <stdlib.h>
#include "mpi.h"

const int MASTER = 0, TAG = 1;
int STOP = -1;

unsigned int bailout(float complex z, unsigned int max_iters) {
  /*
   * Returns i < max_iters if z is NOT in the Mandelbrot set, else max_iters.
   * Test up to max_iters. i is the number of iterations it took to determine z
   * is not in the set.
   */
  unsigned int i = 0;
  float complex zi = 0.0 + 0.0 * I;
  while (creal(zi)*creal(zi)+cimag(zi)*cimag(zi) < 4 && i++ < max_iters) {
    zi = zi*zi + z;
  }
  if (i != max_iters) {
    return i;
  } else {
    return max_iters;
  }
}

int main(int argc, char **argv) {
  /* MPI init (must come before processing argc argv) */
  int err, num_ranks, rank;
  MPI_Status status;
  err = MPI_Init(&argc, &argv);
  if (err != MPI_SUCCESS) {
    fprintf(stderr, "Couldn't start MPI.\n");
    MPI_Abort(MPI_COMM_WORLD, err);
  }
  MPI_Comm_size(MPI_COMM_WORLD, &num_ranks);
  if (num_ranks < 1) {
    fprintf(stderr, "Need at least one rank besides master.\n");
    MPI_Finalize();
    exit(EXIT_FAILURE);
  }
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  /* All ranks must have this information */
  float x0 = -2.0, y0 = -1.0, xn = 1.0, yn = 1.0;
  unsigned int w = 800, h = 600, max_iters = 20;
  if (argc >= 3) w = (unsigned)atoi(argv[2]);
  if (argc >= 4) h = (unsigned)atoi(argv[3]);
  if (argc >= 5) max_iters = (unsigned)atoi(argv[4]);
  if (argc >= 6) x0 = atof(argv[5]);
  if (argc >= 7) y0 = atof(argv[6]);
  if (argc >= 8) xn = atof(argv[7]);
  if (argc >= 9) yn = atof(argv[8]);
  float scale_x = (xn - x0) / (float)w;
  float scale_y = (yn - y0) / (float)h;
  unsigned int **img = NULL; /* Used by master */
  if (rank == MASTER) {
    /* Allocate space for the img */
    img = (unsigned int**)malloc(h*sizeof(unsigned int*));
    for (size_t i = 0; i < h; ++i) {
      img[i] = (unsigned int*)malloc(w*sizeof(unsigned int));
    }
    /* While there are rows to process */
    int i = 0; /* row number and stop signal (-1) on the same send, so int */
    while (i < h) {
      int k = 1;
      do {
        /* Assign one row to each rank */
        err = MPI_Send(&i, 1, MPI_INT, k, TAG, MPI_COMM_WORLD);
        if (err != MPI_SUCCESS) MPI_Abort(MPI_COMM_WORLD, err);
        /* Write to img the processed row from rank k */
        err = MPI_Recv(img[i], w, MPI_UNSIGNED, k, TAG, MPI_COMM_WORLD,
            &status);
        if (err != MPI_SUCCESS) MPI_Abort(MPI_COMM_WORLD, err);
      } while (++i < h && ++k < num_ranks);
    }
    /* All rows done, tell all ranks to stop receiving */
    for (int k = 1; k < num_ranks; ++k) {
      MPI_Send(&STOP, 1, MPI_INT, k, TAG, MPI_COMM_WORLD);
    }
  } else {
    /* While we don't recieve a STOP signal from MASTER */
    int i = !STOP;
    while (i != STOP) {
      /* Get the line number MASTER assgined us */
      err = MPI_Recv(&i, 1, MPI_INT, MASTER, TAG, MPI_COMM_WORLD, &status);
      if (err != MPI_SUCCESS) MPI_Abort(MPI_COMM_WORLD, err);
      if (i == STOP) break; /* Stop if it's the stop signal! */
      /* Allocate memory to process that line */
      unsigned int *row = (unsigned int*)malloc(w*sizeof(unsigned int));
      /* Fill that row with the bailout values */
      float im = y0 + i * scale_y;
      for (unsigned int j = 0; j < w; ++j) {
        float re = x0 + j * scale_x;
        float complex z = re + im * I;
        row[j] = bailout(z, max_iters);
      }
      /* Send the row for MASTER to copy */
      err = MPI_Send(row, w, MPI_UNSIGNED, MASTER, TAG, MPI_COMM_WORLD);
      /* Free it */
      free(row);
      if (err != MPI_SUCCESS) MPI_Abort(MPI_COMM_WORLD, err);
    }
  }
  /* All done, print to file if requested and free. */
  if (rank == MASTER) {
    char const *fname = NULL;
    if (argc >= 2) fname = argv[1];
    FILE *f = fopen(fname, "w");
    for (size_t i = 0; i < h; ++i) {
      if (f) {
        for (size_t j = 0; j < w; ++j) {
          fprintf(f, "%zu,", img[i][j]);
        }
        fprintf(f, "\n");
      }
      free(img[i]);
    }
    free(img);
    if (f) {
      fclose(f);
    }
  }
  MPI_Finalize();
  return EXIT_SUCCESS;
}

Each process allocates memory for the row it is processing, while rank 0 (master) allocates memory for the whole image, which is redundant. I thought about only allocating memory in rank 0, but I read that processes don't share memory and indeed I had some segfaults when I tried to pass a pointer from one process to the other. Oh and I know everything is a little bit clustered inside main.

I've made a script to draw the csv file to a .png.

Example image:

Answer 1

I can see two reasons why your program might be slow.

The first one is that for each row you allocate memory and then free it. If you move the allocation outside of the while loop and reuse the memory then this will get rid of a whole bunch of unnecessary memory allocation operations.

The second one is probably more important: MPI stands for Message Passing Interface. This means that communication between nodes is done by passing messages around. The easiest to imagine an MPI program is to assume that each node is a computer and that the program is running on all these computers in parallel and that the programs can exchange messages over a network.

Hence one goal for a performant MPI program is to minimize the number of message exchanges while maximising the work each individual node is performing.

What you have done is to split your problem in very small chunks which you are then distributing to the workers. What you should do instead: If you have N workers then split up your program into N parts and send each part to a worker. This is essentially the minimal number of messages to exchange if you want all N workers to participate in the solving of the problem.

So rather than sending one row at a time send height / N rows to each worker, let them compute all of them and gather the result back.

Answer 2

At a glance, the MPI code specifically doesn't look too bad. The processes start and stop as needed, and you don't appear to accidentally have a process freeing the wrong memory. As I'm not familiar with this fractal, I'll just comment on the MPI and some general best practices.

• You have STOP as a global variable. Try to keep it local so that it won't be prone to bugs.

  In addition, rename it to stop since all-caps are already used for the constants.

• There could be more spacing in main() to allow more readability. All the code and comments are crammed together, almost making it look like a wall of text. For instance, have every section preceded by a comment have its own section.

  You can start by removing needless comments such as this:

  /* Allocate space for the img */
  img = (unsigned int**)malloc(h*sizeof(unsigned int*));
  

  This reveals nothing already visible in the code, so no comment is needed.

  From there, it may be necessary to separate more of this functionality into additional functions. Ideally, main() should do as little as possible, but it looks like it's doing too much.

• You have separate exit cases that call both MPI_Abort() and MPI_Finalize(). You may just need the latter as the former is for extreme cases and doesn't exit gracefully.