Cython中返回C malloc数组指针

Question

怎么才能把一个malloc数组指针（或者numpy数组指针）从cython高效地返回到python3呢？

只要我不返回数组指针，cython代码运行得非常好。

我想要的是：

def double complex* randn_zig(int n):
  ...
  r = malloc(n*n*sizeof(double complex))
  ...
  return r

c11（gcc 11）的等效代码是：

double complex* randn_zig(int n){

    r = malloc(n*n*sizeof(double complex))

    return r
}

我尝试过

<double complex*> randn_zig(int n):

还有 randn_zig(<double complex*> r, int n):

以及其他组合，但到现在为止都没有成功。如果我能找到一种方法把指向大约10^6到10^10的双精度复数数组的指针返回，那么c和cython的代码版本速度是NumPy/pylab randn版本的5倍。

Answer 1

对于使用C-11标准的gcc 5及以上版本（gcc -std=gnu11 ...），多维数组的malloc和calloc语法发生了很大的变化。

现在，如果你想在main()函数中创建一个大小为n=1024的二维双精度复数calloc数组r[n][n]，代码如下：

long n = 1024;
complex double (*r)[n] = calloc(n, sizeof *r);

下面是一个使用指向这个calloc数组r[n][n]的指针的高斯随机数生成器randn_box_muller()的例子：

inline static void randn_box_muller(long n, complex double r[][n])
{
    long i, j; 
    register double x, y;

    for(i = 0; i < n; i++){
        for(j = 0; j < n; j++){  
            x = 2.*M_PI*dsfmt_genrand_close_open(&dsfmt);
            y = sqrt(-2.*log(dsfmt_genrand_close_open(&dsfmt)));
            r[i][j] = (cos(x) + I*sin(x))*y;
        }
     }
     return;
}

这种相对较新的calloc分配语法有点奇怪。不过它在处理一维、二维甚至n维的calloc和malloc数组时都能很好地工作。希望它也能和Python3一起使用。我希望很快就能进行测试。

Answer 2

我觉得最好的方法是把一个已经在Python中通过NumPy创建的数组的指针传递给Cython。否则，你就得把用malloc创建的数组的内容复制到另一个数组里，就像这个简单的例子所展示的那样：

import numpy as np
cimport numpy as np

from libc.stdlib cimport malloc, free

def main():
  cdef int i, n=40
  cdef double complex *r
  cdef np.ndarray[np.complex128_t, ndim=1] a
  a = np.zeros(n*n, dtype=np.complex128)
  r = <double complex *>malloc(n*n*sizeof(double complex))
  for i in range(n*n):
      r[i] = 1.
  for i in range(n*n):
      a[i] = r[i]
  free(r)
  return a

Answer 3

除了上面提到的两个选项：PyArray_SimpleNewFromData 和直接返回类型化的内存视图（不处理内存），还有一个进一步的选择，就是使用这个 cython.view.array 类。

这个类比较底层，可以用来包装已经存在的内存。它有一个属性叫 callback_free_data，你可以在这里设置一个函数，当这个对象被销毁时会调用这个函数，从而释放内存（这里的示例代码是从文档中复制过来的）：

cdef view.array my_array = view.array(..., mode="fortran", allocate_buffer=False)
my_array.data = <char *> my_data_pointer

# define a function that can deallocate the data (if needed)
my_array.callback_free_data = free

这个类支持缓冲协议，这样你就可以对它进行索引，使用类型化的内存视图，或者用 np.asarray 将它包装成一个 Numpy 数组（不需要复制）。后面这个功能可能比 PyArray_SimpleNewFromData 更容易使用。

Answer 4

Numpy C API

你的问题跟这个帖子很相似。

你可以使用下面的函数把一个C语言的指针传给Numpy数组。这样，当Numpy数组不再使用时，它会自动释放内存。如果你想手动释放这个指针，就不要设置NPY_OWNDATA这个标志。

import numpy as np
cimport numpy as np

cdef pointer_to_numpy_array_complex128(void * ptr, np.npy_intp size):
    '''Convert c pointer to numpy array.
    The memory will be freed as soon as the ndarray is deallocated.
    '''
    cdef extern from "numpy/arrayobject.h":
        void PyArray_ENABLEFLAGS(np.ndarray arr, int flags)
    cdef np.ndarray[np.complex128, ndim=1] arr = \
            np.PyArray_SimpleNewFromData(1, &size, np.NPY_COMPLEX128, ptr)
    PyArray_ENABLEFLAGS(arr, np.NPY_OWNDATA)
    return arr

供参考：

• PyArray_SimpleNewFromData
• Numpy 数据类型 API

Cython 类型内存视图

当然，你也可以使用Cython 的内存视图。

import numpy as np
cimport numpy as np

cdef np.complex128_t[:,:] view = <np.complex128_t[:n,:n]> c_pointer
numpy_arr = np.asarray(view)

上面的代码会把C指针转成一个Numpy数组。不过，这样不会自动释放内存，你必须自己去释放内存，否则会导致内存泄漏！