Monthly Archives: November 2013

Making PyOpenCL handle NumPy arrays as images

PyOpenCL Image objects take a shape tuple that gives (width, height, depth), but NumPy arrays specify shape in the order (rows, columns, …) a.k.a. (height, width, …) where the ellipsis indicates higher dimensions.  What’s important is that the width and height dimensions have been swapped.  The PyOpenCL documentation suggests creating the NumPy arrays in the Fortran or column-major order, instead of the default row major order.  While that is fine for single channel images, RGB and RGBA images still get messed up.  The solution turned out to be quite easy: swap the first and second dimension.  I assume an OpenCL context (ctx) has already been created and you have an image shape tuple specified as image_shape=(rows, columns, …).  This line creates the OpenCL 2D image and stores it in input_image_cl: Continue reading

My computer can’t add (part 2)

In My computer can’t add (part 1) introduced Kahan summation and showed how it helps to improve the accuracy of summing up a large number of floating point values.  Kahan’s algorithm is fine for this task, but what happens when one tries combining the result of two sums?  What if you need to multiply the compensated sum by some value?  This led to a number of searches on Google and online publishers that turned up relatively little.  What I did eventually find led me to the double-double precision and double-single precision libraries written in Fortran.  Currently, they reside at  I am going to focus on the addition operators again, but the Fortran code shows how to handle multiplication, transcendental functions and pretty much everything else you could want. Continue reading