ctypes 中的复数

Question

这可能有点愚蠢，但我正在尝试使用 ctypes 来调用接收复数向量作为参数的函数。但在 ctypes 中没有 c_complex 类。有人知道如何解决这个问题吗？

编辑：我指的是python的ctypes，以防还有其他......

Answer 1

我不太喜欢上面的 @Biggsy 答案，因为它需要用 c/fortran 编写的附加包装器。对于简单的情况来说，这并不重要，但如果你想使用外部库（如 lapack），这将需要你为你想要调用的每个函数/库制作一些代理 dll，这可能不会很有趣，而且很多事情可能会出错（如链接、可重用性等）。我的方法仅依赖于 python 并借助 ctypes 和 numpy。希望这会有所帮助。

在下面的示例中，我展示了如何在 python、numpy（两者本质上都是 C）和纯 C 或 fortran（我使用 fortran，但对于 C 而言，python 端完全相同）之间传递复数/数组。

首先让我们创建一个虚拟的 Fortran (f90) 库，例如 f.f90，带有 C 接口：

! print complex number
subroutine print_c(c) bind(c)
  use iso_c_binding
  implicit none
  complex(c_double_complex), intent(in) :: c
  print "(A, 2F20.16)", "@print_c, c=", c
end subroutine print_c

! multiply real numbers
real(c_double) function mtp_rr(r1,r2) bind(c)
  use iso_c_binding
  implicit none
  real(c_double), intent(in) :: r1,r2
  print "(A)", "@mpt_rr" ! make sure its fortran
  mtp_rr = r1 * r2
end function mtp_rr

! multiply complex numbers
complex(c_double_complex) function mtp_cc(c1,c2) bind(c)
  use iso_c_binding
  implicit none
  complex(c_double_complex), intent(in) :: c1,c2
  print "(A)", "@mpt_cc" ! make sure its fortran
  mtp_cc = c1 * c2
  return
end function mtp_cc

! print array of complex numbers
subroutine print_carr(n, carr) bind(c)
  use iso_c_binding
  implicit none
  integer(c_long), intent(in) :: n
  integer(c_long) :: i
  complex(c_double_complex), intent(in) :: carr(n)
  print "(A)", "@print_carr"
  do i=1,n
    print "(I5,2F20.16)", i, carr(i)
  end do
end subroutine print_carr

该库可以像往常一样简单地编译（在本例中编译为 libf.so）。请注意，每个子例程都包含自己的“use”语句，如果在模块级别声明“use iso_c_binding”，则可以避免这种情况。（我不知道为什么 gfortran 可以在不全局使用 iso_c_binding 的情况下理解函数类型，但它可以工作并且对我来说很好。）

gfortran -共享 -fPIC f.f90 -o libf.so

然后我创建了一个 python 脚本，例如 p.py，其中包含以下内容：

#!/usr/bin/env python3
import ctypes as ct
import numpy as np
## ctypes 1.1.0
## numpy 1.19.5
# print("ctypes", ct.__version__)
# print("numpy", np.__version__)
from numpy.ctypeslib import ndpointer
## first we prepare some datatypes
c_double = ct.c_double
class c_double_complex(ct.Structure): 
    """complex is a c structure
    https://docs.python.org/3/library/ctypes.html#module-ctypes suggests
    to use ctypes.Structure to pass structures (and, therefore, complex)
    """
    _fields_ = [("real", ct.c_double),("imag", ct.c_double)]
    @property
    def value(self):
        return self.real+1j*self.imag # fields declared above
c_double_complex_p = ct.POINTER(c_double_complex) # pointer to our complex
## pointers
c_long_p = ct.POINTER(ct.c_long) # pointer to long (python `int`)
c_double_p = ct.POINTER(ct.c_double) # similar to ctypes.c_char_p, i guess?
# ndpointers work well with unknown dimension and shape
c_double_complex_ndp = ndpointer(np.complex128, ndim=None)
### ct.c_double_complex
### > AttributeError: module 'ctypes' has no attribute 'c_double_complex'

## then we prepare some dummy functions to simplify argument passing
b_long = lambda i: ct.byref(ct.c_long(i))
b_double = lambda d: ct.byref(ct.c_double(d))
b_double_complex = lambda c: ct.byref(c_double_complex(c.real, c.imag))

## now we load the library
libf = ct.CDLL("libf.so")

## and define IO types of its functions/routines
print_c = libf.print_c
print_c.argtypes = [c_double_complex_p] # fortran accepes only references
print_c.restype = None # fortran subroutine (c void)

mtp_rr = libf.mtp_rr
mtp_rr.argtypes = [c_double_p, c_double_p] 
mtp_rr.restype = c_double # ctypes.c_double

mtp_cc = libf.mtp_cc
mtp_cc.argtypes = [c_double_complex_p, c_double_complex_p] 
mtp_cc.restype = c_double_complex # custom c_double_complex

print_carr = libf.print_carr
print_carr.argtypes = [c_long_p, c_double_complex_ndp]
print_carr.restype = None

## finally we can prepare some data and test what we got
print("test print_c")
c = 5.99j+3.1234567890123456789
print_c(b_double_complex(c)) # directly call c/fortran function
print(c)

print("\ntest mtp_rr")
d1 = 2.2
d2 = 3.3
res_d = mtp_rr(b_double(d1),b_double(d2))
print("d1*d2 =", res_d)

print("\ntest mtp_cc")
c1 = 2j+3
c2 = 3j
res = mtp_cc(b_double_complex(c1), b_double_complex(c2))
print(res.value)

print("\ntest print_carr")
n = 10
arr = np.arange(n) + 10j * np.arange(10)
print("arr.shape =",arr.shape)
print_carr(b_long(n), arr)
arr = arr.reshape((5,2)) # still contiguous chunk of memory
print("arr.shape =",arr.shape)
print_carr(b_long(n), arr) # same result
print("done")

一切正常。我不知道为什么 ctypes 中还没有实现这样的东西。

还与此处的其他建议相关：您可以通过创建新的 ctypes 类型来传递复数

__c_double_complex_p = ct.POINTER(2*ct.c_double)
dll.some_function.argtypes = [__c_double_complex_p, ...]

但你不能用这种方式检索结果！要通过类设置正确的 restype only 方法对我有用。

Answer 2

正如 OP 在对 @Armin Rigo 的回答的评论中指出的那样，执行此操作的正确方法是使用包装函数。另外，正如我在原始问题的评论中指出的那样，这对于 C++ 和 Fortran 也是可能的。然而，让它在所有三种语言中工作的方法并不一定是显而易见的。因此，这个答案为每种语言提供了一个工作示例。

假设您有一个采用标量复数参数的 C/C++/Fortran 过程。然后，您需要编写一个包装程序，该程序采用两个浮点数/双精度数（实部和虚部），将它们组合成一个复数，然后调用原始过程。显然，这可以扩展到复数数组，但现在让我们用单个复数来保持简单。

C

例如，假设您有一个 C 函数来打印格式化的复数。所以我们有my_complex.c：

#include <stdio.h>
#include <complex.h>

void print_complex(double complex z)
{
    printf("%.1f + %.1fi\n", creal(z), cimag(z));
}

然后我们必须添加一个包装函数（在同一文件的末尾），如下所示：

void print_complex_wrapper(double z_real, double z_imag)
{
    double complex z = z_real + z_imag * I;
    print_complex(z);
}

以通常的方式将其编译到库中：

gcc -shared -fPIC -o my_complex_c.so my_complex.c

然后从 Python 调用包装函数：

from ctypes import CDLL, c_double
c = CDLL('./my_complex_c.so')
c.print_complex_wrapper.argtypes = [c_double, c_double]
z = complex(1.0 + 1j * 2.0)
c.print_complex_wrapper(c_double(z.real), c_double(z.imag))

C++

在 C++ 中同样的事情有点复杂，因为extern "C"需要定义一个接口以避免名称修改，并且我们需要在 Python 中处理该类（均按照此 SO Q&A）。

所以，现在我们有了my_complex.cpp（我已经添加了包装函数）：

#include <stdio.h>
#include <complex>

class ComplexPrinter
{
    public:
        void printComplex(std::complex<double> z)
        {   
            printf("%.1f + %.1fi\n", real(z), imag(z));
        }

        void printComplexWrapper(double z_real, double z_imag)
        {   
            std::complex<double> z(z_real, z_imag);
            printComplex(z);
        }   
};

我们还需要添加一个extern "C"接口（在同一文件的末尾），如下所示：

extern "C" 
{
    ComplexPrinter* ComplexPrinter_new()
    {   
        return new ComplexPrinter();
    }   
    void ComplexPrinter_printComplexWrapper(ComplexPrinter* printer, double z_real, double z_imag)
    {   
        printer->printComplexWrapper(z_real, z_imag);
    }   
}

按照通常的方式编译成库：

g++ -shared -fPIC -o my_complex_cpp.so my_complex.cpp

并从 Python 调用包装器：

from ctypes import CDLL, c_double
cpp = CDLL('./my_complex_cpp.so')
cpp.ComplexPrinter_printComplexWrapper.argtypes = [c_double, c_double]


class ComplexPrinter:

    def __init__(self):
        self.obj = cpp.ComplexPrinter_new()

    def printComplex(self, z):
        cpp.ComplexPrinter_printComplexWrapper(c_double(z.real), c_double(z.imag))


printer = ComplexPrinter()
z = complex(1.0 + 1j * 2.0)
printer.printComplex(z)

福特兰语言

最后，在 Fortran 中，我们有原始子例程my_complex.f90：

subroutine print_complex(z)
  use iso_c_binding, only: c_double_complex
  implicit none
  complex(c_double_complex), intent(in) :: z
  character(len=16) :: my_format = "(f4.1,a3,f4.1,a)"
  print my_format, real(z), " + ", aimag(z), "i" 
end subroutine print_complex

我们向其中添加包装函数（在同一文件的末尾）：

subroutine print_complex_wrapper(z_real, z_imag) bind(c, name="print_complex_wrapper")
  use iso_c_binding, only: c_double, c_double_complex
  implicit none
  real(c_double), intent(in) :: z_real, z_imag
  complex(c_double_complex) :: z
  z = cmplx(z_real, z_imag)
  call print_complex(z)
end subroutine print_complex_wrapper

然后按照通常的方式编译成库：

gfortran -shared -fPIC -o my_complex_f90.so my_complex.f90

并从 Python 调用（注意 POINTER 的使用）：

from ctypes import CDLL, c_double, POINTER
f90 = CDLL('./my_complex_f90.so')
f90.print_complex_wrapper.argtypes = [POINTER(c_double), POINTER(c_double)]
z = complex(1.0 + 1j * 2.0)
f90.print_complex_wrapper(c_double(z.real), c_double(z.imag))