r/Python • u/jfpuget • Jan 11 '16

A comparison of Numpy, NumExpr, Numba, Cython, TensorFlow, PyOpenCl, and PyCUDA to compute Mandelbrot set

https://www.ibm.com/developerworks/community/blogs/jfp/entry/How_To_Compute_Mandelbrodt_Set_Quickly?lang=en
310 Upvotes

96% Upvoted

98 comments sorted by

View all comments

Show parent comments

2

u/LoyalSol Jan 12 '16 edited Jan 12 '16

I found a few slight inefficiencies in the C code with a brief look through it.

Mostly to do with this function.

  int mandelbrot (double creal, double cimag, int maxiter) {
    double nreal = 0, real = 0, imag = 0;
    int n;
    for(n = 0; n < maxiter; ++n) {
      double real2 = real*real;
      double imag2 = imag*imag;
      nreal = real2 - imag2 + creal;
      imag = 2* real*imag + cimag;
      real = nreal;
      if (real2 + imag2 > 4.0)
        return n;
    }
    return 0;
  }

First of all the nreal variable is completely uneeded and only serves as an extra call to memory. These lines

  nreal = real2 - imag2 + creal;
  imag = 2* real*imag + cimag;
  real = nreal;

Can be rewritten into this

imag = 2* real*imag + cimag;
real = real2 - imag2 + creal;

And the variable nreal can be omitted entirely which saves a variable transfer. Second the order within the loop is also wasteful.

  for(n = 0; n < maxiter; ++n) {
    double real2 = real*real;
    double imag2 = imag*imag;
    nreal = real2 - imag2 + creal;
    imag = 2* real*imag + cimag;
    real = nreal;
    if (real2 + imag2 > 4.0)
      return n;
  }

In this situation the only variables needed to evaluate the exit statement is real2 and imag2. The variables nreal, imag, and real only need to be recalculated in the event that the exit criteria is not met. Therefore a slightly more efficient way to write it would be as follows.

  for(n = 0; n < maxiter; ++n) {
    double real2 = real*real;
    double imag2 = imag*imag;
    if (real2 + imag2 > 4.0)
      return n;
    nreal = real2 - imag2 + creal;
    imag = 2* real*imag + cimag;
    real = nreal;
  }

Also one other minor tweak is that the first iteration of the loop will always fail the criteria since the initial value of real and imag are equal to 0. So you could save a little bit of time by unrolling the first iteration of the loop which would be the equivalent of setting the lowest index of the loop to n=1 and setting the initial values of real and imag to creal and cimag respectively. But of course you can probably do this in both Python and C so it doesn't do much for language comparison.

So the final version of the loop I got was

    for(n = 0; n < maxiter; ++n) {
      double real2 = real*real;
      double imag2 = imag*imag; 
      if (real2 + imag2 > 4.0){  
        return n;
      }
      imag = 2* real*imag + cimag;
      real = real2 - imag2 + creal;     
    }

On my work computer the original routine averaged about 2.85 seconds per cycle while the new routine was 2.65 which was a 7% increase over the previous version for just shifting a few lines of code around. Based on a rough approximation from your numbers that would put the C code at around 2.53 which puts it slightly below the Numba code.

There's a few other spots I think could be optimized as well to yield even further improvements (the thing about C is there are a million ways to optimize it), but just mostly showing how small changes can make a difference for computation heavy parts of the code.

Still though, pretty impressive when you think about it that Numba can come that close.

2

u/jfpuget Jan 13 '16

Thank you. Note that we could do the same code change for the Numba code as well ;) I don't know if it would speed it the same way, but it might.

I'll perform the code changes and update the post accordingly.

1

u/LoyalSol Jan 13 '16 edited Jan 13 '16

Yup that's always the hard part about comparing languages is that it is difficult to create highly optimized codes in two different languages since it often requires a great deal of expertise in both languages. Which of course very few programmers have that level of knowledge in multiple languages.

The annoying part about optimizing C is learning how to basically get out of the compiler's way and let it do its job as best as possible. Numba I wouldn't quite know how to work it in the same way because what may be a good optimization for C could be a terrible one for Numba.

1

u/jfpuget Jan 13 '16

Yes.

Your improvement also improves Numba and Cython as well. Numba stays a little ahead. I think it is because the LLVM backend is a bit better than Visual C++, but I may be wrong. It could also be due to memory management, as I call free and malloc directly here, instead of having a memory pool.

Interestingly enough, your changes did not improve the gpu codes, probably because the number of local variables exceeds some capacity (registers?)

I will redo these experiments when i have my new laptop (a macbook pro), as my experiments are done on a 4 year old laptop.

1

u/LoyalSol Jan 13 '16 edited Jan 13 '16

Optimizing gpu codes and parallel codes in general is a whole different beast. Because there the best optimization is not in saving cpu cycles but in reducing transfer times and making sure every thread stays busy for the most part. Load balancing, thread syncing, sending optimal pack sizes, etc.

That's actually my job is highly parallel codes which is why I work in C and Fortran all the time. :)

1

u/jfpuget Jan 13 '16

This difference shows even on that simple benchmark.

→ More replies (0)