I'm attempting to write a cython interface to the complex version of the MUMPS solver (zmumps). I'm running into some problems as I have no previous experience with either C or cython. Following the example of the pymumps package I was able to get the real version of the code (dmumps) to work.
I believe that my problem are the pointers to the ZMUMPS_COMPLEX structures. For the
So far I have the following (lifted heavily from 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)
In my python code I'm able to set the irn and jcn using
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]
However, I have no idea how to set the values of a or rhs. Any help would be greatly appreciated.
This might help:
The following example lets you get at the C-level members of Python’s built-in “complex” object:
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
Grabbed from here.
As ZMUMPS_COMPLEX and the builtin Py_complex structs have exactly the same structure, you should be able to do the trick by creating a bridge between those two types (using typedefs and/or cast or a function that turns a Py_complex into a ZMUMPS_COMPLEX)...
I'd love to help more but I don't currently have mumps installed...
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With