How to pass multidimensional Numpy arrays to C function using CFFI?

I am trying to speed up my Python program by implementing a function in C++ and embedding it in my code using CFFI. The function takes two 3x3 arrays and computes a distance. The Python code is the following:

import cffi
import numpy as np

ffi = cffi.FFI()
ffi.cdef("""
    extern double dist(const double s[3][3], const double t[3][3]);
""")
lib = ffi.dlopen("./dist.so")

S = np.array([[-1.63538,  0.379116, -1.16372],[-1.63538, 0.378137, -1.16366 ],[-1.63193, 0.379116, -1.16366]], dtype=np.float32)
T = np.array([[-1.6467834, 0.3749715, -1.1484985],[-1.6623441, 0.37410975, -1.1647063 ],[-1.6602284, 0.37400728, -1.1496595 ]], dtype=np.float32)
Sp = ffi.cast("double(*) [3]", S.ctypes.data)
Tp = ffi.cast("double(*) [3]", T.ctypes.data)

dd = lib.dist(Sp,Tp);

This solution doesn't work as intended. Indeed the arguments printed by the C function are:

 Sp=[[0.000002,  -0.270760,  -0.020458]
    [0.000002,  0.000000,  0.000000]
    [0.000000,  0.000000,  0.000000]]
 Tp=[[0.000002,  -0.324688,  -0.020588]
    [0.000002,  0.000000,  0.000000]
    [0.000000,  0.000000,  -nan]]

I have also tried the following to initialize the pointers:

Sp = ffi.new("double *[3]")
for i in range(3):
    Sp[i] = ffi.cast("double *", S[i].ctypes.data)
Tp = ffi.new("double *[3]")
for i in range(3):
    Tp[i] = ffi.cast("double *", T[i].ctypes.data)
dd = lib.dist(Sp,Tp);

But this solution rises an error in dist(Sp,Tp):

TypeError: initializer for ctype 'double(*)[3]' must be a pointer to same type, not cdata 'double *[3]'

Do you have any idea on how to make it work? Thanks.

Solution

The types double[3][3] and double *[3] are not equivalent. The former is an 2D array of double 3x3, stored as 9 contiguous double. The latter is an array of 3 double pointers, which is not how 2D static arrays are implemented in C or C++.

Both numpy arrays and C++ static arrays are represented in memory as a contiguous block of elements, it's just the double *[3] type in the middle that's throwing a wrench in the works. What you'd want is to use double[3][3] proper or double[3]* (a pointer to a row of three doubles). Note that if you use the latter you may need to change the function prototype to take double [][3].