Have pointer to data. Need Numpy array in Fortran

2020-06-22 03:05发布

问题:

This is a very common use-case for me. I have a C function that returns me a pointer to doubles:

 //myheader.h
 double *mycfuntion(...)

I know the dimensions of the data that get returned. I also know that the data are Fortran-ordered. I want to write a Cython "shim" to get the data into Python as a Numpy array:

#myshim.pyx
import numpy
cimport numpy as cnumpy

cnumpy.import_array()

cdef extern from "myheader.h" :
     double *mycfunction(...)

def mypyfunc(...) :
     cdef double *data = mycfunction(...)
          **MAGIC**
     return outarray

MAGIC ideas

(A) cdef cnumpy.ndarray[ cnumpy.double_t, mode='fortran', ...] outarray
This would be the handiest way of doing things. There's something critical that I'm missing here, though, concerning how to turn the pointer data into a buffer that I can pass to the cnumpy.ndarray constructor. I've tried:

cdef cnumpy.ndarray[ cnumpy.double_t, mode='fortran', ...] outarray
cdef bytes databuffer = <char *>data
outarray = numpy.ndarray(buffer=databuffer, dtype=numpy.double, ...)

This approach consistently fails with TypeError: buffer is too small for requested array

(B) The Numpy C-API
I've used cnumpy.PyArray_SimpleNewFromData(...) plenty from Cython. It works just fine. The problem is that it doesn't support a flags argument so I can't tell it to construct a Fortran array. The alternative that I'd used in pure C implementations is PyArray_NewFromDescr(...). It accepts flags. This approach is long-winded and painful, and means getting some symbols from numpy via an extern block that aren't already imported. There has got to be a better way.


I have been Googling my face off on this problem but nothing obvious has popped up. Maybe I'm an idiot. Or just sleep-depping. Cheers!

回答1:

This is a little hackish, but you can reverse the order of your dimensions then return the transpose. This should get your striding correct, but not copy the data.

cimport numpy as cnp
cnp.import_array()

cdef double a1[10]
for i in range(10):
    a1[i] = i

# C order
a2 = cnp.PyArray_SimpleNewFromData(2, [2, 5], cnp.NPY_FLOAT64, a1)
# fortran order
a3 = cnp.PyArray_SimpleNewFromData(2, [5, 2], cnp.NPY_FLOAT64, a1).T
# a2 and a3 point to the same data
a2[0,0] = 10.
print a2
print a3