[Numpy-discussion] The NumPy Mandelbrot code 16x slower than Fortran

[Samuel John]

I'd like to add http://git.tiker.net/pyopencl.git/blob/HEAD:/examples/demo_mandelbrot.py to the discussion, since I use pyopencl  (http://mathema.tician.de/software/pyopencl) with great success in my daily scientific computing. Install with pip.

PyOpenCL does understand numpy arrays. You write a kernel (small c-program) directly into a python triple quoted strings and get a pythonic way to program GPU and core i5 and i7 CPUs with python Exception if something goes wrong. Whenever I hit a speed bottleneck that I cannot solve with pure numpy, I code a little part of the computation for GPU. The compilation is done just in time when you run the python code.

Especially for the mandelbrot this may be a _huge_ gain in speed since its embarrassingly parallel.

Samuel


On 23.01.2012, at 14:02, Robert Cimrman wrote:

> On 01/23/12 13:51, Sturla Molden wrote:
>> Den 23.01.2012 13:09, skrev Sebastian Haase:
>>> 
>>> I would think that interactive zooming would be quite nice
>>> ("illuminating")  .... and for that 13 secs would not be tolerable....
>>> Well... it's not at the top of my priority list ... ;-)
>>> 
>> 
>> Sure, that comes under the 'fast enough' issue. But even Fortran might
>> be too slow here?
>> 
>> For zooming Mandelbrot I'd use PyOpenGL and a GLSL fragment shader
>> (which would be a text string in Python):
>> 
>> madelbrot_fragment_shader = """
>> 
>> uniform sampler1D tex;
>> uniform vec2 center;
>> uniform float scale;
>> uniform int iter;
>> void main() {
>>      vec2 z, c;
>>      c.x = 1.3333 * (gl_TexCoord[0].x - 0.5) * scale - center.x;
>>      c.y = (gl_TexCoord[0].y - 0.5) * scale - center.y;
>>      int i;
>>      z = c;
>>      for(i=0; i<iter; i++) {
>>          float x = (z.x * z.x - z.y * z.y) + c.x;
>>          float y = (z.y * z.x + z.x * z.y) + c.y;
>>          if((x * x + y * y)>   4.0) break;
>>          z.x = x;
>>          z.y = y;
>>      }
>>      gl_FragColor = texture1D(tex, (i == iter ? 0.0 : float(i)) / 100.0);
>> }
>> 
>> """
>> 
>> The rest is just boiler-plate OpenGL...
>> 
>> Sources:
>> 
>> http://nuclear.mutantstargoat.com/articles/sdr_fract/
>> 
>> http://pyopengl.sourceforge.net/context/tutorials/shader_1.xhtml
> 
> Off-topic comment: Or use some algorithmic cleverness, see [1]. I recall Xaos 
> had interactive, extremely fast a fluid fractal zooming more than 10 (or 15?) 
> years ago (-> on a laughable hardware by today's standards).
> 
> r.
> 
> [1] http://wmi.math.u-szeged.hu/xaos/doku.php
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion