How to pass a 2d array returned in a C function to python using ctypes

Question

I have a C function returns a 2d array (actually double**), I want to pass this double** to python using ctypes. There are questions about this but the array is modified in parameter list rather than using the return of the function. Here is my code (The C function and the corresponding python wrapper)

double** flux_function_2d(double *u, double gamma){

  double **F;
  double p, H;

  F = (double **)malloc(4*sizeof(double *));
  for (int i = 0; i < 4; i++)
    F[i] = (double *)malloc(2*sizeof(double));

  p = (gamma - 1.0) * (u[3] - 0.5 * (u[1]*u[1] + u[2]*u[2]) / u[0]);
  H = u[3] / u[0] + p / u[0];

  F[0][0] = u[1]; F[1][0] = u[1] * u[1] / u[0] + p;
  F[2][0] = u[1] * u[2] / u[0]; F[3][0] = u[1] * H;

  F[0][1] = u[2]; F[1][1] = u[2] * u[1] / u[0];
  F[2][1] = u[2] * u[2] / u[0] + p; F[3][1] = u[2] * H;

  return F;
}


def c_flux_function_2d(u, gamma):
    flux = np.ctypeslib.load_library("libs/flux.so", ".")
    flux.flux_function_2d.argtypes = [POINTER(c_double), c_double]
    flux.flux_function_2d.restype = POINTER(POINTER(c_double))
    F = flux.flux_function_2d(u.ctypes.data_as(POINTER(c_double)),
        c_double(gamma))
    F_arr = np.ctypeslib.as_array(F, shape=(4, 2))
    return F_arr

The error happens at F_arr = np.ctypeslib.as_array(F, shape=(4, 2)) which means that numpy cannot parsing the pointer of a pointer in ctypes.

Thanks in advance!

Answer 1

If you actually had a 2D array, you could call numpy.ctypeslib.as_array, but you don't.

You have a pointer to the first element of an array of pointers to the first elements of other arrays. This is a completely different thing from a 2D array, that just happens to be indexable with the same syntax in C.

The memory layout of what you have is completely incompatible with NumPy. NumPy arrays have a single memory buffer with shape and stride information telling how far you can go in each dimension and how many bytes you have to go to take a step in each dimension.

For a C-contiguous array (the most common case), the memory layout of an MxN array's buffer is the same as a C whatevertype[M][N]. In contrast, your data structure's subarrays are scattered in arbitrary memory locations, with no such thing as a stride in the first dimension.

You will need to get your data into a layout NumPy can work with. There are various ways to do this, including just allocating and working with a contiguous buffer in the C function, or creating a numpy.empty([4, 2]) and copying the data into that with a nested loop.