将复杂的numpy阵列传递到Cython中的C

Pass complex numpy array to C++ in Cython

本文关键字：Cython 中的 复杂 numpy 阵列      更新时间：2023-10-16

我想对pyx脚本的cythonize部分进行cythize，该部分涉及使用具有复数数字的Numpy阵列。Python脚本的相关部分如下所示：

M = np.dot(N , Q)

在我的工作中，N，Q和M是具有复杂数字条目的Numpy数组。

具体来说，我想将矩阵N和Q传输到C++代码，然后在C++中执行矩阵乘法。

我知道使用指针转移到C++脚本的方法的方法，然后使用Cython，但我对如何处理具有复杂值的Numpy阵列的事物感到有些困惑。

这就是我目前试图将数组从pyx转移到C++的方式。

import numpy as np
cimport numpy as np
cdef extern from "./matmult.h" nogil:
    void mult(double* M, double* N, double* Q)
def sim():    
    cdef:
        np.ndarray[np.complex128_t,ndim=2] N = np.zeros(( 2 , 2 ), dtype=np.float64)
        np.ndarray[np.complex128_t,ndim=2] Q = np.zeros(( 2 , 2 ), dtype=np.float64)  
        np.ndarray[np.complex128_t,ndim=2] M = np.zeros(( 2 , 2 ), dtype=np.float64)          
    N = np.array([[1.1 + 2j,2.2],[3.3,4.4]])
    Q = np.array([[3.3,4.4+5j],[5.5,6.6]])  
    mult(&M[0,0], &N[0,0], &Q[0,0])
    print M

这是我的C 代码：

#include "matmult.h"
using namespace std;
int main(){}
void mult(double *M, double *N, double *Q)
{
  double P[2][2], A[2][2], B[2][2];
  for (int i=0; i<2; i++)
  {
    for (int j=0; j<2; j++)
    {
      A[i][j] = *( N + ((2*i) + j) );
      B[i][j] = *( Q + ((2*i) + j) );
      P[i][j] = 0;      
    }
  }
  for (int i=0; i<2; i++)
  {
    for (int j=0; j<2; j++)
    {
      for (int k=0; k<2; k++)
      {
         P[i][j] += A[i][k]*B[k][i];  
      }
    }
  }  
  for (int i=0; i<2; i++)
  {
    for (int j=0; j<2; j++)
    {
      *( M + ((2*i) + j) ) = P[i][j];      
    }
  }
}

当我使用Cython对此进行编译时，我会收到以下错误

mat.pyx:17:27: Cannot assign type 'double complex *' to 'double *'

我将很高兴在这里得到一些帮助。

此错误消息告诉您出了什么问题：

mat.pyx：17：27：无法将'双重复杂 *'分配给'double *'

也就是说，您有一个双重复杂的指针，从numpy(指针到复杂的128 numpy dtype(，您正在尝试使用双重指针将其传递到C 函数中。C 需要能够处理复数，因此，如果您更改double* -> std ::复合物，则可以解决问题

void mult(double *M, double *N, double *Q)

变成

#include <complex>
void mult(std::complex<double> *M, std::complex<double> *N, std::complex<double> *Q)

numpy矩阵是否足以满足您的用例？Cython可能过于杀伤。

编辑：好的，我终于得到了一些东西，与C STD :: Complex and C Double _complex类型相处有些奇怪。

cppmul.pyx：

import numpy as np
cimport numpy as np
cdef extern from "./matmult.h" nogil:
    void mult(np.complex128_t* M, np.complex128_t* N, np.complex128_t* Q)
def sim():
    cdef:
        np.ndarray[np.complex128_t,ndim=2] N = np.zeros(( 2 , 2 ), dtype=np.complex128)
        np.ndarray[np.complex128_t,ndim=2] Q = np.zeros(( 2 , 2 ), dtype=np.complex128)
        np.ndarray[np.complex128_t,ndim=2] M = np.zeros(( 2 , 2 ), dtype=np.complex128)
    N = np.array([[1.1 + 2j,2.2],[3.3,4.4]])
    Q = np.array([[3.3,4.4+5j],[5.5,6.6]])
    mult(&M[0,0], &N[0,0], &Q[0,0])
    print M

matmul.c：

#include "matmult.h"
void mult(complex_t *M, complex_t *N, complex_t *Q)
{
  complex_t P[2][2], A[2][2], B[2][2];
  for (int i=0; i<2; i++)
  {
    for (int j=0; j<2; j++)
    {
      A[i][j] = *( N + ((2*i) + j) );
      B[i][j] = *( Q + ((2*i) + j) );
      P[i][j] = 0;
    }
  }
  for (int i=0; i<2; i++)
  {
    for (int j=0; j<2; j++)
    {
      for (int k=0; k<2; k++)
      {
         P[i][j] += A[i][k]*B[k][i];
      }
    }
  }
  for (int i=0; i<2; i++)
  {
    for (int j=0; j<2; j++)
    {
      *( M + ((2*i) + j) ) = P[i][j];
    }
  }
}

matmult.h：

#include <complex.h>
typedef double _Complex complex_t;
void mult(complex_t *M, complex_t *N, complex_t *Q);

setup.py：

from distutils.core import setup
from Cython.Build import cythonize
from distutils.extension import Extension
import numpy as np
sourcefiles = ['cppmul.pyx', 'matmult.c']
extensions = [Extension("cppmul",
                        sourcefiles,
                        include_dirs=[np.get_include()],
                        extra_compile_args=['-O3']
)]
setup(
    ext_modules = cythonize(extensions)
)

运行python setup.py build_ext --inplace后，它会导入并按预期运行

import cppmul
cppmul.sim()

结果：

[[15.73 +6.6j 15.73 +6.6j]
 [43.56+16.5j 43.56+16.5j]]

尝试此

    #include "matmult.h"
 using namespace std;
int main(){}
void mult(double *M, double *N, double *Q)
{
 double P[2][2], A[2][2], B[2][2];
 for (int i=0; i<2; i++)
 {
for (int j=0; j<2; j++)
{
  A[i][j] = *( N + ((2*i) + j) );
  B[i][j] = *( Q + ((2*i) + j) );
  P[i][j] = 0;      
  }
  }
  for (int i=0; i<2; i++)
  {
    for (int j=0; j<2; j++)
    {
    for (int k=0; k<2; k++)
    {
       P[i][j] += A[i][k]*B[k][i];  
     }
 }
}  
for (int i=0; i<2; i++)
 {
  for (int j=0; j<2; j++)
  {
  *(  ((2*i) + j) )+ M = P[i][j];      
  }
 }
}
相关文章：