ctypeslib | Generate python ctypes classes from C headers | Compiler library

First of all, a scoped C array allocated on the stack (like in somefunction) must never be returned by a function. The space of the stack will be reused by other function like the one of CPython for example. The returned array must be allocated on the heap instead.

Moreover, writing a function working with Numpy arrays using ctypes is pretty cumbersome. As you found out, you need to pass the full shape in parameter. But the thing is you also need to pass the strides for each dimension and for each input arrays in parameter of the function since they may not be contiguous in memory (for example np.transpose change this). That being said, we can assume that the input array is contiguous for sake of performance and sanity. This can be enforced with np.ascontiguousarray. The pointer of the views a and b can be extracted using numpy.ctypeslib.as_ctypes, but hopefully ctype can do that automatically. Furthermore, the returned array is currently a C pointer and not a Numpy array. Thus, you need to create a Numpy array with the right shape and strides from it numpy.ctypeslib.as_array. Because the resulting shape is not known from the caller, you need to retrieve it from the callee function using several integer pointers (one per dimension). In the end, this results in a pretty-big ugly highly-bug-prone code (which will often silently crash if anything goes wrong not to mention the possible memory leaks if you do not pay attention to that). You can use Cython to do most of this work for you.

Assuming you do not want to use Cython or you cannot, here is an example code with ctypes:

Your C function takes int * arguments but your np.array is using 64-bit integers. In C, int is 4 bytes (32-bit). Try np.int32.

Numba currently uses LLVM-Lite to compile the code efficiently to a binary (after the Python code has been translated to an LLVM intermediate representation). The code is optimized like en C++ code would be using Clang with the flags -O3 and -march=native. This last parameter is very important as is enable LLVM to use wider SIMD instructions on relatively-recent x86-64 processors: AVX and AVX2 (possible AVX512 for very recent Intel processors). Otherwise, by default Clang and GCC use only the SSE/SSE2 instructions (because of backward compatibility).

Another difference come from the comparison between GCC and the LLVM code from Numba. Clang/LLVM tends to aggressively unroll the loops while GCC often don't. This has a significant performance impact on the resulting program. In fact, you can see that the generated assembly code from Clang:

With Clang (128 items per loops):

So after many failures and countless tries, I came up with this, which seems to work:

Listing [Python.Docs]: ctypes - A foreign function library for Python.

According to [NumPy]: numpy.ctypeslib.ndpointer(dtype=None, ndim=None, shape=None, flags=None) (emphasis is mine):

An ndpointer instance is used to describe an ndarray in restypes and argtypes specifications.

argtypes and restype only have meaning in the context of a function (call), and that is when ndpointer should be used (exemplified in the URL).

As a workaround, you can have an ndpointer structure member, but you have to initialize it in the same manner as the ctypes.POINTER one (in MyStruct_Working). However, bear in mind that this will come with some limitations (you won't have indexing) afterwards.
My suggestion is to stick with ctypes.POINTER(ctypes.c_double).

I finally found the way:

Computing transpositions efficiently is hard. This primitive is not compute-bound but memory-bound. This is especially true for big matrices stored in the RAM (and not CPU caches).

Numpy use a view-based approach which is great when only a slice of the array is needed and quite slow the computation is done eagerly (eg. when a copy is performed). The way Numpy is implemented results in a lot of cache misses strongly decreasing performance when a copy is performed in this case.

I found this function virtually always run in single-threaded, whatever how large the transposed matrix/array is.

This is clearly dependant of the Numpy implementation used. AFAIK, some optimized packages like the one of Intel are more efficient and take more often advantage of multithreading.

I think a parallelized transpose function should be a resolution. The problem is how to implement it.

Yes and no. It may not be necessary faster to use more threads. At least not much more, and not on all platforms. The algorithm used is far more important than using multiple threads.

On modern desktop x86-64 processors, each core can be bounded by its cache hierarchy. But this limit is quite high. Thus, two cores are often enough to nearly saturate the RAM throughput. For example, on my (4-core) machine, a sequential copy can reach 20.4 GiB/s (Numpy succeed to reach this limit), while my (practical) memory throughput is close to 35 GiB/s. Copying A takes 72 ms while the naive Numpy transposition A takes 700 ms. Even using all my cores, a parallel implementation would not be faster than 175 ms while the optimal theoretical time is 42 ms. Actually, a naive parallel implementation would be much slower than 175 ms because of caches-misses and the saturation of my L3 cache.

Naive transposition implementations do not write/read data contiguously. The memory access pattern is strided and most cache-lines are wasted. Because of this, data are read/written multiple time from memory on huge matrices (typically 8 times on current x86-64 platforms using double-precision). Tile-based transposition algorithm is an efficient way to prevent this issue. It also strongly reduces cache misses. Ideally, hierarchical tiles should be used or the cache-oblivious Z-tiling pattern although this is hard to implement.

Here is a Numba-based implementation based on the previous informations:

Listing [Python.Docs]: ctypes - A foreign function library for Python and [NumPy]: C-Types Foreign Function Interface (numpy.ctypeslib). The latter states about the following 2 items (emphasis is mine):

1. as_array

   Create a numpy array from a ctypes array or POINTER.

2. ndpointer

   An ndpointer instance is used to describe an ndarray in restypes and argtypes specifications. This approach is more flexible than using, for example, POINTER(c_double), since several restrictions can be specified, which are verified upon calling the ctypes function.

imageBuffer_ptr is not a "ctypes array or POINTER" (literally: it's not a ctypes.POINTER instance, but ctypes.c_void_p one (there's a subtle difference between the 2)).

The C layer (which in our case lies in the middle) doesn't know about NumPy arrays and so on, so from its PoV, that argument is simply an unsigned short* (doesn't have information about shape, ...). So, on it's other side (in the Python callback), you have to reconstruct the array from the plain C pointer. That can be achieved in 2 ways, like I did in the example below:

• callback_ct: Use good-old C way: modify the callback's signature to use ct.POINTER(ct.c_ushort). That is the most straightforward, but the drawback is that it loses all the extra checks when called, meaning that it takes any pointer not just a ndarray one

• callback_np: Manually recreate the array from the argument (using a series of conversions). It's more advanced (some might even consider it a bit hackish)

code00.py:

I can't compile your code as is, but the problem is that np.shape returns (rows,columns) or the equivalent (height,width), not (width,height):