在Cython中传递结构体指针

Question

我正在尝试编写一个Cython接口，用于复杂版本的MUMPS求解器（zmumps）。由于我之前没有使用过C或Cython，所以遇到了一些问题。参考pymumps package的示例，我成功地让实数版本的代码（dmumps）运行起来了。

我认为我的问题在于指向ZMUMPS_COMPLEX结构体的指针。目前我有以下代码（大量借鉴自pymumps）：

zmumps_c.pxd:

from libc.string cimport strncpy

cdef extern from "mumps_c_types.h":

    ctypedef struct ZMUMPS_COMPLEX "ZMUMPS_COMPLEX":
        double   r
        double   i

cdef extern from "zmumps_c.h":

    ctypedef int MUMPS_INT

    ctypedef struct c_ZMUMPS_STRUC_C "ZMUMPS_STRUC_C":
        MUMPS_INT      sym, par, job
        MUMPS_INT      comm_fortran    # Fortran communicator 
        MUMPS_INT      n

        # Assembled entry
        MUMPS_INT      nz
        MUMPS_INT      *irn
        MUMPS_INT      *jcn
        ZMUMPS_COMPLEX *a

        # RHS and statistics 
        ZMUMPS_COMPLEX *rhs
        MUMPS_INT      infog[40]

    void c_zmumps_c "zmumps_c" (c_ZMUMPS_STRUC_C *)

zmumps_c.pyx

cdef class ZMUMPS_STRUC_C:
    cdef c_ZMUMPS_STRUC_C ob

    property sym:
        def __get__(self): return self.ob.sym
        def __set__(self, value): self.ob.sym = value
    property par:
        def __get__(self): return self.ob.par
        def __set__(self, value): self.ob.par = value
    property job:
        def __get__(self): return self.ob.job
        def __set__(self, value): self.ob.job = value
    property comm_fortran:
        def __get__(self): return self.ob.comm_fortran
        def __set__(self, value): self.ob.comm_fortran = value
    property n:
        def __get__(self): return self.ob.n
        def __set__(self, value): self.ob.n = value
    property nz:
        def __get__(self): return self.ob.nz
        def __set__(self, value): self.ob.nz = value
    property irn:
        def __get__(self): return <long> self.ob.irn
        def __set__(self, long value): self.ob.irn = <MUMPS_INT*> value
    property jcn:
        def __get__(self): return <long> self.ob.jcn
        def __set__(self, long value): self.ob.jcn = <MUMPS_INT*> value
    property a:
        def __get__(self): return <long> self.ob.a
        def __set__(self, long value): self.ob.a = <ZMUMPS_COMPLEX*> value
    property rhs:
        def __get__(self): return <long> self.ob.rhs
        def __set__(self, long value): self.ob.rhs = <ZMUMPS_COMPLEX*> value

    property infog:
        def __get__(self):
            cdef MUMPS_INT[:] view = self.ob.infog
            return view


def zmumps_c(ZMUMPS_STRUC_C s not None):
    c_zmumps_c(&s.ob)

在我的Python代码中，我可以使用以下方式设置irn和jcn：

import zmumps_c
import numpy as np

MUMPS_STRUC_C = staticmethod(zmumps_c.ZMUMPS_STRUC_C)
id = MUMPS_STRUC_C()
x = np.r_[1:10]
id.irn = x.__array_interface__['data'][0]

然而，我不知道如何设置a或rhs的值。非常感谢您的帮助。

- dwfm

Answer 1

有多种方法可以实现这个目标 - 这里介绍一种方法。

下面的代码定义了一个ZMUMPS_COMPLEX的包装器。然后定义了一个ZMUMPS_STRUC_C的包装器，其中__get__和__set__方法用于rhs属性，使其能够接受ZMUMPS_COMPLEX的包装器。

zmumps_c.pyx

cdef class ZMUMPS_COMPLEX:
    '''A wrapper for the ZMUMPS_COMPLEX struct'''
    cdef c_ZMUMPS_COMPLEX c_zm

    def __init__(self, double real, double imag=0):
        self.c_zm.r = real
        self.c_zm.i = imag

    property r:
        def __get__(self):
            return self.c_zm.r
        def __set__(self, value):
            self.c_zm.r = value

    property i:
        def __get__(self):
            return self.c_zm.i
        def __set__(self, value):
            self.c_zm.i = value

    def __repr__(self):
        return '({real}{imag:+}j)'.format(real=self.c_zm.r, imag=self.c_zm.i)

cdef class ZMUMPS_STRUC_C:
    '''A wrapper for the ZMUMPS_STRUC_C struct'''
    cdef c_ZMUMPS_STRUC_C ob
    cdef object _rhs

    property rhs:
        def __get__(self):
            return self._rhs
        def __set__(self, ZMUMPS_COMPLEX c):
            self._rhs = c
            self.ob.rhs = &c.c_zm

    def test(self):
        return (self.ob.rhs[0].r, self.ob.rhs[0].i,)

def main():
    z = ZMUMPS_STRUC_C()
    c = ZMUMPS_COMPLEX(-3.5, 2.0)
    z.rhs = c
    print z.rhs
    print z.test()
    c.r = 42.0
    print z.rhs
    z.rhs.i = -5.0
    print z.rhs

main() 函数演示了这两个对象的行为。输出应该是这样的:

(-3.5+2.0j)
(-3.5, 2.0)
(42.0+2.0j)
(42.0-5.0j)

我没有安装这个库，所以我使用下面的虚拟定义进行测试：

zmumps_c.pxd

cdef struct c_ZMUMPS_COMPLEX "ZMUMPS_COMPLEX":
    double r
    double i

cdef struct c_ZMUMPS_STRUC_C "ZMUMPS_STRUC_C":
    c_ZMUMPS_COMPLEX *rhs

setup.py

from distutils.core import setup
from Cython.Build import cythonize

setup(
    ext_modules = cythonize("example.pyx")
)

- Snorfalorpagus

---

谢谢您的回复，这非常有帮助。当设置实部和虚部时，它似乎仍然会遇到问题。您的示例在实部和虚部为标量时有效，但在它们是向量（numpy数组）时会崩溃。例如，运行： ' import numpy as np r = np.r_[1:5] i = np.r_[1:5] c = ZMUMPS_COMPLEX(r, i) ' 会导致崩溃。再次感谢您的帮助！ - dwfm

@dwfm 这是可以预料的。ZMUMPS_COMPLEX 结构体期望实部和虚部为双精度浮点数，而不是数组！是否有另一个结构体期望向量？ - Snorfalorpagus

@dwfm 我想我误解了。rhs 应该是一个长度为 n 或 nz 的 ZMUMPS_COMPLEX 数组吗？ - Snorfalorpagus

rhs 应该是长度为 n 的 ZMUMPS_COMPLEX 数组。a 应该是长度为 nz 的 ZMUMPS_COMPLEX 数组。 - dwfm

Answer 2

这可能会有所帮助：

以下示例允许您访问Python内置“complex”对象的C级成员：

cdef extern from "complexobject.h":

    struct Py_complex:
        double real
        double imag

    ctypedef class __builtin__.complex [object PyComplexObject]:
        cdef Py_complex cval

# A function which uses the above type
def spam(complex c):
    print "Real:", c.cval.real
    print "Imag:", c.cval.imag

从这里获取。

由于ZMUMPS_COMPLEX和内置的Py_complex 结构具有完全相同的结构，您应该能够通过创建这两种类型之间的桥梁来完成操作（使用typedef和/或强制类型转换或将Py_complex转换为ZMUMPS_COMPLEX的函数）...

我很想提供更多帮助，但目前我没有安装mumps...

- Gauthier Boaglio