python - 如何在 C 中编写一个函数，它将两个矩阵作为 numpy 数组的输入

Question

我想在 C 中实现一个函数，它将两个矩阵作为 numpy 数组的输入并返回另一个矩阵。

首先我在文件 f_2_mat_in.c 中有 C 函数:

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

double *f_2_mat_in(
    const double *matrix_one, 
    size_t n, 
    size_t m,
    const double *matrix_two,
    size_t u, 
    size_t v
    )
{
  double *matrix_out = malloc(sizeof *matrix_out * n * m);

  for (size_t i = 0; i < n; i++){
    for (size_t j = 0; j < m; j++){
      printf("%f ", matrix_one[i * m + j]);
      matrix_out[i * m + j] = matrix_one[i * m + j];
    printf("--\n");
    }
  }

  for (size_t i = 0; i < u; i++){
    for (size_t j = 0; j < v; j++){
      printf("%f ", matrix_two[u * v + j]);
    printf("--\n");
    }
  }

  printf("\n");
  return matrix_out;
}

然后我有了使用 f_2_mat_in 的 main.py python 脚本:c

from ctypes import c_void_p, c_double, c_size_t, c_int, cdll, POINTER
from numpy.ctypeslib import ndpointer
import pdb
import numpy as np

c_double_p1 = POINTER(c_double)
c_double_p2 = POINTER(c_double)

n = 5 
m = 4
kernel_size = 3
u = 3 
v = 6
matrix_one = np.random.randn(n, m).astype(c_double)
matrix_one[0][1] = np.nan
lib = cdll.LoadLibrary("f_2_mat_in.so")
f_2_mat_in = lib.f_2_mat_in
f_2_mat_in.argtypes = c_int, c_double_p1, c_size_t, c_size_t, c_double_p2, c_size_t, c_size_t
f_2_mat_in.restype = ndpointer(
        dtype=c_double,
        shape=(n,m),
        flags='C')
f_2_mat_in.restype = ndpointer(
        dtype=c_double,
        shape=(u,v),
        flags='C')

matrix_out = f_2_mat_in(
        c_int(kernel_size),
        matrix_one.ctypes.data_as(c_double_p1),
        c_size_t(n),
        c_size_t(m),
        matrix_one.ctypes.data_as(c_double_p2),
        c_size_t(u),
        c_size_t(v))

print("matrix_one:", matrix_one)
print("matrix_two:", matrix_two)
print("matrix_out:", matrix_out)
print("matrix_one.shape:", matrix_one.shape)
print("matrix_two.shape:", matrix_two.shape)
print("matrix_out.shape:", matrix_out.shape)
print("in == out", matrix_one == matrix_out)
pdb.set_trace()

最后我有了编译 C 脚本并运行 pyhon 脚本的 bash 脚本:

a=f_2_mat_in
cc -fPIC -shared -o $a.so $a.c && python main.py

我收到以下错误(在 MacOS 11.6.1 中):

Python(53609,0x111c84e00) malloc: can't allocate region
:*** mach_vm_map(size=5626830804037632, flags: 100) failed (error code=3)
Python(53609,0x111c84e00) malloc: *** set a breakpoint in malloc_error_break to debug

我做错了什么？

Answer 1

检查 [SO]: C function called from Python via ctypes returns incorrect value (@CristiFati's answer)对于使用 CTypes 时的常见错误(关于argtypes、restype)。< br/>但是自从使用 NumPy 数组后，可以根据 [NumPy]: C-Types Foreign Function Interface (numpy.ctypeslib) 进一步改进。 .

原代码有问题:

• C:

  • 索引在第 2^nd 循环内越界(U未定义 Behavior)

  • 矩阵未被释放(内存泄漏)

• Python:

  • argtypes 和 restype 未正确设置

    • 函数 header 不匹配 C

  • 2^nd 矩阵未定义

  • 还有很多

  所以我决定从头开始。

我假设(因为我没有运行原始代码)malloc 错误是因为它试图取消引用 Python st 参数em>' 传递了一个 int，而 C 将它解释为一个 double* - 这是(再次)UB。

dll00.c:

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>

#if defined(_WIN32)
#  define DLL00_EXPORT_API __declspec(dllexport)
#else
#  define DLL00_EXPORT_API
#endif

#define ELEMENT_TYPE double

#if defined(__cplusplus)
extern "C" {
#endif

DLL00_EXPORT_API ELEMENT_TYPE* matrixFunc(const ELEMENT_TYPE *pMat0, size_t rows0, size_t cols0, const ELEMENT_TYPE *pMat1, size_t rows1, size_t cols1);
DLL00_EXPORT_API void deallocArray(ELEMENT_TYPE *pMat);

#if defined(__cplusplus)
}
#endif


ELEMENT_TYPE* matrixFunc(
    const ELEMENT_TYPE *pMat0,
    size_t rows0,
    size_t cols0,
    const ELEMENT_TYPE *pMat1,
    size_t rows1,
    size_t cols1)
{
    size_t matSize = sizeof(ELEMENT_TYPE) * cols0 * rows0;
    ELEMENT_TYPE *pRet = (ELEMENT_TYPE*)(malloc(matSize));

    printf("\nC: - 1st matrix (%zu, %zu):\n", rows0, cols0);
    memcpy(pRet, pMat0, matSize);  // Faster than setting each element individually
    for (size_t i = 0; i < rows0; ++i) {
        for (size_t j = 0; j < cols0; ++j) {
            printf("%3.3f ", pMat0[i * cols0 + j]);
            //pRet[i * cols0 + j] = pMat0[i * cols0 + j];
        }
        printf("\n");
    }

    printf("\nC: - 2nd matrix (%zu, %zu):\n", rows1, cols1);
    for (size_t i = 0; i < rows1; ++i) {
        for (size_t j = 0; j < cols1; ++j) {
            printf("%3.3f ", pMat1[i * cols1 + j]);
        }
        printf("\n");
    }

    printf("\n");
    return pRet;
}


void deallocArray(ELEMENT_TYPE *pMat)
{
    if (pMat)
        free(pMat);
}

code00.py:

#!/usr/bin/env python

import ctypes as cts
import sys

import numpy as np


DLL_NAME = "./dll00.{:s}".format("dll" if sys.platform[:3].lower() == "win" else "so")


def np_mat_type(rows, cols, element_type=float):
    return np.ctypeslib.ndpointer(dtype=element_type, shape=(rows, cols), flags="C_CONTIGUOUS")


def main(*argv):
    np.set_printoptions(precision=3)
    rows0 = 5
    cols0 = 4
    kernel_size = 3
    rows1 = 3
    cols1 = 6
    mat0 = np.random.randn(rows0, cols0).astype(cts.c_double)
    mat1 = np.random.randn(rows1, cols1).astype(cts.c_double)

    dll = cts.CDLL(DLL_NAME)
    matrix_func = dll.matrixFunc
    matrix_func.argtypes = (
        np_mat_type(rows0, cols0), cts.c_size_t, cts.c_size_t,
        np_mat_type(rows1, cols1), cts.c_size_t, cts.c_size_t)
    matrix_func.restype = np_mat_type(rows0, cols0)

    dealloc_array = dll.deallocArray
    dealloc_array.argtypes = (np_mat_type(rows0, cols0),)
    dealloc_array.restype = None

    print("mat0:")
    print(mat0)
    print("\nmat1:")
    print(mat1)
    mat_res = matrix_func(mat0, rows0, cols0, mat1, rows1, cols1)
    print("result:")
    print(mat_res)
    print("\nequality:", np.all(mat_res == mat0))

    dealloc_array(mat_res)


if __name__ == "__main__":
    print("Python {:s} {:03d}bit on {:s}\n".format(" ".join(elem.strip() for elem in sys.version.split("\n")),
                                                   64 if sys.maxsize > 0x100000000 else 32, sys.platform))
    rc = main(*sys.argv[1:])
    print("\nDone.\n")
    sys.exit(rc)

输出:

• Nix(Linux):

  (qaic-env) [cfati@cfati-5510-0:/mnt/e/Work/Dev/StackOverflow/q073806188]> ~/sopr.sh
  ### Set shorter prompt to better fit when pasted in StackOverflow (or other) pages ###
  
  [064bit prompt]> ls
  code00.py  dll00.c
  [064bit prompt]> gcc -fPIC -shared -o dll00.so dll00.c
  [064bit prompt]> ls
  code00.py  dll00.c  dll00.so
  [064bit prompt]>
  [064bit prompt]> python ./code00.py
  Python 3.8.10 (default, Jun 22 2022, 20:18:18) [GCC 9.4.0] 064bit on linux
  
  mat0:
  [[-0.83  -0.892 -1.137  0.333]
   [-0.608 -1.403 -0.74   0.652]
   [-0.299 -1.987 -1.009 -0.35 ]
   [-0.383  1.173  0.772 -0.616]
   [-1.063 -0.33   0.124 -0.639]]
  
  mat1:
  [[ 0.017  0.207  0.051  0.269  0.131  0.282]
   [ 0.911  0.686  0.491 -1.846 -2.196  0.287]
   [ 1.897 -0.174  2.99  -0.104 -0.376 -0.404]]
  
  C: - 1st matrix (5, 4):
  -0.830 -0.892 -1.137 0.333
  -0.608 -1.403 -0.740 0.652
  -0.299 -1.987 -1.009 -0.350
  -0.383 1.173 0.772 -0.616
  -1.063 -0.330 0.124 -0.639
  
  C: - 2nd matrix (3, 6):
  0.017 0.207 0.051 0.269 0.131 0.282
  0.911 0.686 0.491 -1.846 -2.196 0.287
  1.897 -0.174 2.990 -0.104 -0.376 -0.404
  
  result:
  [[-0.83  -0.892 -1.137  0.333]
   [-0.608 -1.403 -0.74   0.652]
   [-0.299 -1.987 -1.009 -0.35 ]
   [-0.383  1.173  0.772 -0.616]
   [-1.063 -0.33   0.124 -0.639]]
  
  equality: True
  
  Done.
  

• 获胜:

  [cfati@CFATI-5510-0:e:\Work\Dev\StackOverflow\q073806188]> sopr.bat
  ### Set shorter prompt to better fit when pasted in StackOverflow (or other) pages ###
  
  [prompt]> "c:\Install\pc032\Microsoft\VisualStudioCommunity\2019\VC\Auxiliary\Build\vcvarsall.bat" x64 > nul
  
  [prompt]> dir /b
  code00.py
  dll00.c
  dll00.so
  
  [prompt]> cl /nologo /MD /DDLL dll00.c  /link /NOLOGO /DLL /OUT:dll00.dll
  dll00.c
     Creating library dll00.lib and object dll00.exp
  
  [prompt]> dir /b
  code00.py
  dll00.c
  dll00.dll
  dll00.exp
  dll00.lib
  dll00.obj
  dll00.so
  
  [prompt]>
  [prompt]> "e:\Work\Dev\VEnvs\py_pc064_03.09_test0\Scripts\python.exe" ./code00.py
  Python 3.9.9 (tags/v3.9.9:ccb0e6a, Nov 15 2021, 18:08:50) [MSC v.1929 64 bit (AMD64)] 064bit on win32
  
  mat0:
  [[-0.883  1.441 -1.563  1.452]
   [-1.751 -0.625 -0.782 -0.326]
   [-0.257 -0.838 -0.356 -0.187]
   [ 0.503  0.981  0.041  0.383]
   [-0.637 -1.098  0.122  1.086]]
  
  mat1:
  [[ 1.139  1.059  0.621  1.128  0.872 -0.163]
   [-0.159  0.15   0.163  1.379  0.486  0.349]
   [ 0.855  0.447 -1.12   0.257 -0.22   0.099]]
  
  C: - 1st matrix (5, 4):
  -0.883 1.441 -1.563 1.452
  -1.751 -0.625 -0.782 -0.326
  -0.257 -0.838 -0.356 -0.187
  0.503 0.981 0.041 0.383
  -0.637 -1.098 0.122 1.086
  
  C: - 2nd matrix (3, 6):
  1.139 1.059 0.621 1.128 0.872 -0.163
  -0.159 0.150 0.163 1.379 0.486 0.349
  0.855 0.447 -1.120 0.257 -0.220 0.099
  
  result:
  [[-0.883  1.441 -1.563  1.452]
   [-1.751 -0.625 -0.782 -0.326]
   [-0.257 -0.838 -0.356 -0.187]
   [ 0.503  0.981  0.041  0.383]
   [-0.637 -1.098  0.122  1.086]]
  
  equality: True
  
  Done.