Numpy array to ctypes with FORTRAN ordering

By default, Numpy creates C-ordered arrays (because it is written in C), that is, row-major arrays. Doing a transposition with A.T creates a view of the array where the strides are reversed (ie. no copy). That being said, np.ascontiguousarray does a copy because the array is now not contiguous anymore and copies are expensive. This is why slow_conversion is slow. Note that the contiguity can be tested with yourarray.flags['F_CONTIGUOUS'] and by checking yourarray.strides. Note also that yourarray.flags and yourarray.__array_interface__ provide information about whether the array has been copied and also information about strides.

np.asfortranarray returns an array laid out in Fortran order in memory regarding the documentation. It can perform a copy if needed. In fact, np.asfortranarray(A) does a copy while np.asfortranarray(A.T) does not. You can check the C code of the function for more information about this behaviour. Since both are seen as FORTRAN-contiguous, it is better to use np.asfortranarray(A.T) that do not make any copy in this case.

Regarding ctypes, it deals with C arrays which are stored with a row-major ordering as opposed to FORTRAN array that are stored with a column-major ordering. Furthermore, C arrays do not support strides natively as opposed to FORTRAN ones. This means a row is basically a memory view of contiguous data stored in memory. Since slow_conversion(A)[1][0] == 1 is required to be true, this means the returned object should have the first item of the second row is 1 and thus that the columns are mandatory stored contiguously in memory. The thing is that the initial array is not FORTRAN-contiguous but C-contiguous and thus a transposition is required. Transpositions are pretty expensive (although the Numpy implementation is suboptimal).

Assuming you do not want to pay the overhead of a copy/transposition, the problem need to be relaxed. There are few possible options:

• Create FORTRAN ordered array directly with Numpy using for example np.zeros((n,n), dtype=np.int8, order='F'). This create a C array with transposed strides so to behave like a FORTRAN array where computations operating on columns are fast (remember that Numpy is written in C so row-major ordered array are the reference). With this, the first row in ctypes is actually a column. Note that you should be very careful when mixing C and FORTRAN ordered array for sake of performance as non-contiguous access are much slower.
• Use a strided FORTRAN array. This solution basically means that basic column-based computations will be much slower and one need to write row-based computation that are quite unusual in FORTRAN. You can extract the pointer to the C-contiguous array with A.ctypes.data_as(POINTER(c_double)), the strides with A.strides and the shape with A.shape. That being said, this solution appears not to be really portable/standard. The standard way appears to use C binding in FORTRAN. I am not very familiar with this but you can find a complete example in this answer.

A last solution is to transpose data manually using a fast transposition algorithm. An in-place transposition is faster than an out-of-place transposition (typically 2 times) but this requires a squared array and it mutates the input array that should not be used later (unless it is fine to operate on a transposed array). A fast transposition cannot be written using Numpy directly. One solution is to do in in Numba (this is a sub-optimal out-of-place solution but pretty fast already), or to do it in C or directly in FORTRAN (using a wrapper function in all cases). This should be significantly faster than what Numpy does but still far slower than a basic ctypes wrapping.

If your starting array must be C-ordered:

As you have pointed out, np.ctypeslib.as_ctypes(...) is fast.

The bottleneck of your computation is np.ascontiguousarray(A.T) - which is equivalent to np.asfortranarray(A) which is equally slow.

This leads me to believe that this can't be made faster using only numpy. I mean, since a whole dedicated function to do it already exists ( np.asfortranarray ) - we'd assume it has been optimized to have the best possible performance.

If you can start with an F-ordered array:

In this case, problem solved! You can get the performance just like that of the as_ctypes call because the array is already in the desired order in memory:

n = 10000
A = np.zeros((n, n), dtype=np.int8, order='F')
A[0, 1] = 1

def conversion(A):
    return np.ctypeslib.as_ctypes(np.ascontiguousarray(A.T))

assert conversion(A)[1][0] == 1

%%timeit
conversion(A)

# 4.76 µs ± 68.4 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

Thanks to @JérômeRichard and @StephanSchlecht for inspiring this idea. Read @JérômeRichard's answer for more insights!

There is one aspect that can be improved about it.

The operation is not only making a copy, but because it is loading and storing to main memory instead of using cache.

Usually processors will access multiple blocks multiple contiguous bytes whenever memory and use it from cache. If you run out of cache space some old block is evicted.

For sake of the argument let's say that your CPU works on blocks of 8 bytes and that the rows are contiguous. In one matrix you will access columns and the other you will access rows. When you are writing down a column you are loading multiple columns but updating only one. This overhead of one single column can be seen by copying a few columns

n = 2**14
A = np.random.randint(0, 100, (n,n), dtype=np.int8)
B = np.empty_like(A)
%%timeit
B[:1,:] = A[:1,:]
%%timeit
B[:4,:] = A[:4,:]

If you do the same on the rows you should notice something roughly linear. If you copy columns the cost of copying one column is very close to the cost of copying two columns or even 8, or 16, depending on the hardware.

I will use n=2**14 to make the things easier, but the principle applies to any dimension.

• If you have a small enough let's say 8 x 8 the entire matrix fit in the cache, so you can transpose it without accessing any cache.
• If you have are copying large continuous block of data even if you can't do the entire operation on cache you reduce the number of times a given data will be loaded from/to memory again.

Based on this what I tried is to rearrange the matrix in a matrix of smaller contiguous blocks, first I transpose the elements in a block and then the blocks in the matrix.

For the baseline

B = np.ascontiguousarray(A.T)

3.12 s ± 446 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Using 8x8 blocks

T0 = A.reshape(2048,8,2048,8)
T1 = np.ascontiguousarray(T0.transpose(0,2,3,1))
T2 = np.ascontiguousarray(T1.transpose(1,0,2,3))
T3 = np.ascontiguousarray(T2.transpose(0,2,1,3))
B = T3.reshape(A.shape)

786 ms ± 54.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

assert np.all(B == A.T) # 2.8s

It is still 200x slower than a simple copy, but it is already 4x faster than the original approach.

Allocating only two instead of 3 temporary arrays help as well

T0 = np.empty_like(A)
T1 = np.empty_like(A)
T0.reshape(2048,2048,8,8)[:] = A.reshape(2048,8,2048,8).transpose(0,2,3,1)
T1.reshape(2048,2048,8,8)[:] = T0.reshape(2048,2048,8,8).transpose(1,0,2,3)
T0.reshape(2048,8,2048,8)[:] = T1.reshape(2048,2048,8,8).transpose(0,2,1,3)
B = T0

686 ms ± 60.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)