使用ctypes将一个二维numpy数组传递给C

Question

什么是使用ctypes将numpy 2d数组传递给c函数的正确方法？以下是我目前的尝试（导致了段错误）：

C代码：

void test(double **in_array, int N) {
    int i, j;
    for(i = 0; i<N; i++) {
        for(j = 0; j<N; j++) {
            printf("%e \t", in_array[i][j]);
        }
        printf("\n");
    }
}

Python 代码：

from ctypes import *
import numpy.ctypeslib as npct

array_2d_double = npct.ndpointer(dtype=np.double,ndim=2, flags='CONTIGUOUS')
liblr = npct.load_library('libtest.so', './src')

liblr.test.restype = None
liblr.test.argtypes = [array_2d_double, c_int]

x = np.arange(100).reshape((10,10)).astype(np.double)
liblr.test(x, 10)

- jrsm

---

你知道在C语言中，double **和double [N][N]是不能互换的吗？ - Filipe Gonçalves

我不懂Python也不懂NumPy，但如果它是一个NxN的数组，你应该将in_array声明为double (*in_array)[N]，其中N是第二维的大小。 - Filipe Gonçalves

如果N在运行时不是固定的，那么这个程序如何工作？ - jrsm

我认为在in_array参数中使用双**不是正确的做法，所以问题是如何正确处理C语言中的2d ndpointer。 - jrsm

你看过这个问题吗？https://dev59.com/Y2025IYBdhLWcg3wr4B6?rq=1 - Filipe Gonçalves

Answer 1

这可能是一个晚回答，但我终于使它工作了。所有功劳都归于Sturla Molden，他在这个链接中提供了帮助。

关键是，要注意double **是一个np.uintp类型的数组。因此，我们有

xpp = (x.ctypes.data + np.arange(x.shape[0]) * x.strides[0]).astype(np.uintp)
doublepp = np.ctypeslib.ndpointer(dtype=np.uintp)

然后使用doublepp作为类型，将xpp传入。请查看附加的完整代码。

C 代码：

// dummy.c 
#include <stdlib.h> 

__declspec(dllexport) void foobar(const int m, const int n, const 
double **x, double **y) 
{ 
    size_t i, j; 
    for(i=0; i<m; i++) 
        for(j=0; j<n; j++) 
            y[i][j] = x[i][j]; 
} 

Python代码：

# test.py 
import numpy as np 
from numpy.ctypeslib import ndpointer 
import ctypes 

_doublepp = ndpointer(dtype=np.uintp, ndim=1, flags='C') 

_dll = ctypes.CDLL('dummy.dll') 

_foobar = _dll.foobar 
_foobar.argtypes = [ctypes.c_int, ctypes.c_int, _doublepp, _doublepp] 
_foobar.restype = None 

def foobar(x): 
    y = np.zeros_like(x) 
    xpp = (x.__array_interface__['data'][0] 
      + np.arange(x.shape[0])*x.strides[0]).astype(np.uintp) 
    ypp = (y.__array_interface__['data'][0] 
      + np.arange(y.shape[0])*y.strides[0]).astype(np.uintp) 
    m = ctypes.c_int(x.shape[0]) 
    n = ctypes.c_int(x.shape[1]) 
    _foobar(m, n, xpp, ypp) 
    return y 

if __name__ == '__main__': 
    x = np.arange(9.).reshape((3, 3)) 
    y = foobar(x) 

希望这有所帮助，

谢恩

- Yuxiang Wang

---

谢谢您的回复。不过我不太明白这里发生了什么。这如何扩展到三维及以上？ - navjotk

@navjotk 我非常推荐使用 Cython 来处理更复杂的情况，以便更轻松地进行维护 :) - Yuxiang Wang

警告：请检查 x.flags['C_CONTIGUOUS']==True。如果不是，请在执行之前使用 x = np.ascontiguousarray(x, dtype= x.dtype) 进行转换。否则，如果在传递给 foobar 之前对 x 进行了转置，则 xpp 将会出错。 - benbo

Answer 2

#include <stdio.h>

void test(double (*in_array)[3], int N){
    int i, j;

    for(i = 0; i < N; i++){
        for(j = 0; j < N; j++){
            printf("%e \t", in_array[i][j]);
        }
        printf("\n");
    }
}

int main(void)
{
    double a[][3] = {
        {1., 2., 3.},
        {4., 5., 6.},
        {7., 8., 9.},
    };

    test(a, 3);
    return 0;
}

如果您想在函数中使用double **，则必须传递指向double指针的数组（而不是2D数组）：

#include <stdio.h>

void test(double **in_array, int N){
    int i, j;

    for(i = 0; i < N; i++){
        for(j = 0; j< N; j++){
            printf("%e \t", in_array[i][j]);
        }
        printf("\n");
    }
}

int main(void)
{
    double a[][3] = {
        {1., 2., 3.},
        {4., 5., 6.},
        {7., 8., 9.},
    };
    double *p[] = {a[0], a[1], a[2]};

    test(p, 3);
    return 0;
}

另一种方法（由@eryksun建议）：传递一个单一的指针并进行一些算术运算以获取索引：

#include <stdio.h>

void test(double *in_array, int N){
    int i, j;

    for(i = 0; i < N; i++){
        for(j = 0; j< N; j++){
            printf("%e \t", in_array[i * N + j]);
        }
        printf("\n");
    }
}

int main(void)
{
    double a[][3] = {
        {1., 2., 3.},
        {4., 5., 6.},
        {7., 8., 9.},
    };

    test(a[0], 3);
    return 0;
}

- David Ranieri

---

@DavidHeffernan，我不知道ctypes代码部分，我试图解释为什么在传递2D数组时double **无法工作。 - David Ranieri

问题在于N在运行时不是固定的。 - jrsm

@jrsm：如果你真的想要使用double **方法，你可以像这样创建指针数组：c_double_p = POINTER(c_double); in_array_ptrs = (c_double_p * len(in_array))(*(r.ctypes.data_as(c_double_p) for r in in_array))。但我不建议这样做，因为它效率低下。 - Eryk Sun

Answer 3

虽然回复可能有些晚了，但我希望它能够帮助到其他遇到相同问题的人。

由于NumPy数组在内部保存为1D数组，因此可以在C中简单地重建2D形状。 这里是一个小的MWE示例：

// libtest2d.c
#include <stdlib.h> // for malloc and free
#include <stdio.h>  // for printf

// create a 2d array from the 1d one
double ** convert2d(unsigned long len1, unsigned long len2, double * arr) {
    double ** ret_arr;

    // allocate the additional memory for the additional pointers
    ret_arr = (double **)malloc(sizeof(double*)*len1);

    // set the pointers to the correct address within the array
    for (int i = 0; i < len1; i++) {
        ret_arr[i] = &arr[i*len2];
    }

    // return the 2d-array
    return ret_arr;
}

// print the 2d array
void print_2d_list(unsigned long len1,
    unsigned long len2,
    double * list) {

    // call the 1d-to-2d-conversion function
    double ** list2d = convert2d(len1, len2, list);

    // print the array just to show it works
    for (unsigned long index1 = 0; index1 < len1; index1++) {
        for (unsigned long index2 = 0; index2 < len2; index2++) {
            printf("%1.1f ", list2d[index1][index2]);
        }
        printf("\n");
    }

    // free the pointers (only)
    free(list2d);
}

和

# test2d.py

import ctypes as ct
import numpy as np

libtest2d = ct.cdll.LoadLibrary("./libtest2d.so")
libtest2d.print_2d_list.argtypes = (ct.c_ulong, ct.c_ulong,
        np.ctypeslib.ndpointer(dtype=np.float64,
            ndim=2,
            flags='C_CONTIGUOUS'
            )
        )
libtest2d.print_2d_list.restype = None

arr2d = np.meshgrid(np.linspace(0, 1, 6), np.linspace(0, 1, 11))[0]

libtest2d.print_2d_list(arr2d.shape[0], arr2d.shape[1], arr2d)

如果您使用gcc -shared -fPIC libtest2d.c -o libtest2d.so编译代码，然后运行python test2d.py，它应该会打印出数组。

我希望这个例子更或多或少是自解释的。思路是将形状也传递给C代码，然后创建一个double **指针，为其保留附加指针的空间。然后将这些指针设置为指向原始数组的正确部分。

PS：我在C方面还是初学者，如果有不应该这样做的原因，请评论告知。

- Tim

Answer 4

我这里传递了两个二维的numpy数组，并打印了一个数组的值作为参考

您可以在C++中使用并编写自己的逻辑

cpp_function.cpp

使用以下命令进行编译：g++ -shared -fPIC cpp_function.cpp -o cpp_function.so

#include <iostream>
extern "C" {
void mult_matrix(double *a1, double *a2, size_t a1_h, size_t a1_w, 
                  size_t a2_h, size_t a2_w, int size)
{
    //std::cout << "a1_h & a1_w" << a1_h << a1_w << std::endl; 
    //std::cout << "a2_h & a2_w" << a2_h << a2_w << std::endl; 
    for (size_t i = 0; i < a1_h; i++) {
        for (size_t j = 0; j < a1_w; j++) {
            printf("%f ", a1[i * a1_h + j]);
        }
        printf("\n");
    }
    printf("\n");
  }

}

Python文件 main.py

import ctypes
import numpy
from time import time

libmatmult = ctypes.CDLL("./cpp_function.so")
ND_POINTER_1 = numpy.ctypeslib.ndpointer(dtype=numpy.float64, 
                                      ndim=2,
                                      flags="C")
ND_POINTER_2 = numpy.ctypeslib.ndpointer(dtype=numpy.float64, 
                                    ndim=2,
                                    flags="C")
libmatmult.mult_matrix.argtypes = [ND_POINTER_1, ND_POINTER_2, ctypes.c_size_t, ctypes.c_size_t]
# print("-->", ctypes.c_size_t)

def mult_matrix_cpp(a,b):
    shape = a.shape[0] * a.shape[1]
    libmatmult.mult_matrix.restype = None
    libmatmult.mult_matrix(a, b, *a.shape, *b.shape , a.shape[0] * a.shape[1])

size_a = (300,300)
size_b = size_a

a = numpy.random.uniform(low=1, high=255, size=size_a)
b = numpy.random.uniform(low=1, high=255, size=size_b)

t2 = time()
out_cpp = mult_matrix_cpp(a,b)
print("cpp time taken:{:.2f} ms".format((time() - t2) * 1000))
out_cpp = numpy.array(out_cpp).reshape(size_a[0], size_a[1])

- Devil