Passing numpy arrays in Cython to a C function that requires dynamically allocated arrays

Question

I have some C code that has the following declaration:

int myfunc(int m, int n, const double **a, double **b, double *c);

So a is a constant 2D array, b is a 2D array, and c is a 1D array, all dynamically allocated. b and c do not need to be anything specifically before they are passed to myfunc, and should be understood as output information. For the purposes of this question, I'm not allowed to change the declaration of myfunc.

Question 1: How do I convert a given numpy array a_np into an array a with the format required by this C function, so that I can call this C function in Cython with a?

Question 2: Are the declarations for b and c below correct, or do they need to be in some other format for the C function to understand them as a 2D and 1D array (respectively)?

My attempt:

myfile.pxd

cdef extern from "myfile.h":
    int myfunc(int p, int q, const double **a, double **b, double *c)

mytest.pyx

cimport cython
cimport myfile
import numpy as np
cimport numpy as np

p = 3
q = 4
cdef:
    double** a = np.random.random([p,q])
    double** b
    double* c

myfile.myfunc(p, q, a, b, c)

Then in iPython I run

import pyximport; pyximport.install()
import mytest

The line with the definition of a gives me the error message Cannot convert Python object to 'double **'. I don't get any error messages regarding b or c, but since I'm unable to run the C function at this time, I'm not sure the declarations of b and c are written correctly (that is, in a way that will enable the C function to output a 2D and a 1D array, respectively).

Other attempts: I've also tried following the solution here, but this doesn't work with the double-asterisk type of arrays I have in the myfunc declaration. The solution here does not apply to my task because I can't change the declaration of myfunc.

Answer 1

Create a helper array in cython

To get a double** from a numpy array, you can create a helper-array of pointers in your *.pyx file. Further more, you have to make sure that the numpy array has the correct memory layout. (It might involve creating a copy)

Fortran order

If your C-function expects fortran order (all x-coordinates in one list, all y coordinates in another list, all z-coordinates in a third list, if your array a corresponds to a list of points in 3D space)

N,M = a.shape
# Make sure the array a has the correct memory layout (here F-order)
cdef np.ndarray[double, ndim=2, mode="fortran"] a_cython =
                         np.asarray(a, dtype = float, order="F")
#Create our helper array
cdef double** point_to_a = <double **>malloc(M * sizeof(double*))
if not point_to_a: raise MemoryError
try:
    #Fillup the array with pointers
    for i in range(M): 
        point_to_a[i] = &a_cython[0, i]
    # Call the C function that expects a double**
    myfunc(... &point_to_a[0], ...)
finally:
    free(point_to_a)

C-order

If your C-function expects C-order ([x1,y1,z1] is the first list, [x2,y2,z2] the second list for a list of 3D points):

N,M = a.shape
# Make sure the array a has the correct memory layout (here C-order)
cdef np.ndarray[double, ndim=2, mode="c"] a_cython =
                         np.asarray(a, dtype = float, order="C")
#Create our helper array
cdef double** point_to_a = <double **>malloc(N * sizeof(double*))
if not point_to_a: raise MemoryError
try:
    for i in range(N): 
        point_to_a[i] = &a_cython[i, 0]
    # Call the C function that expects a double**
    myfunc(... &point_to_a[0], ...)
finally:
    free(point_to_a)

Answer 2

Reply 1: You can pass NumPy array via Cython to C using the location of the start of the array (see code below).

Reply 2: Your declarations seem correct but I don't use this approach of explicit memory management. You can use NumPy to declare cdef-ed arrays.

Use

cdef double[:,::1] a = np.random.random([p, q])
cdef double[:,::1] b = np.empty([p, q])
cdef double[::1] b = np.empty(q)

Then pass &a[0], the location of the start of the array, to your C function. The ::1 is to ensure contiguousness.

A good reference for this is Jake Vanderplas' blog: https://jakevdp.github.io/blog/2012/08/08/memoryview-benchmarks/

Finally, typically one creates functions in Cython and calls them in Python, so your Python code would be:

import pyximport; pyximport.install()
import mytest
mytest.mywrappedfunc()

where mywrappedfunc is a Python (def and not cdef) function defined in the module that can do the array declaration show above.