无法从ctypes获取fortran函数输出

Question

我正在尝试使用ctypes从Python调用Fortran函数.我试图从子程序和函数(两者都具有相同的功能)获得结果,但我无法从函数中获得预期的输出,而子程序运行良好.问题是我有很多库使用Fortran函数而不是子例程.Fortran函数和ctypes有什么问题吗？

一段Fortran代码:

MODULE Vector
! Public types
 TYPE VectorType
    PRIVATE
    DOUBLE PRECISION, DIMENSION(3):: components = 0.0d0
 END TYPE VectorType
!---------------------------------------------------------------------    
 CONTAINS 
!---------------------------------------------------------------------
 SUBROUTINE newVect(this,vectorIn)
 TYPE (VectorType),       INTENT(OUT):: this
 DOUBLE PRECISION, DIMENSION(3), INTENT(IN)::vectorIn

   this%components = (/vectorIn(1), vectorIn(2), vectorIn(3)/)

 END SUBROUTINE newVect
!---------------------------------------------------------------------
 SUBROUTINE subVect(this,vectorOut)

 TYPE(VectorType), INTENT (OUT):: vectorOut
 TYPE(VectorType), INTENT (IN) :: this

   vectorOut%components = this%components

 END SUBROUTINE subVect
!----------------------------------------------------------------------
 TYPE(VectorType) FUNCTION getVect(this) RESULT(vectorOut)

 TYPE(VectorType), INTENT (IN) :: this

   vectorOut%components = this%components

  END FUNCTION getVect
!--------------------------------------------------------------------
END MODULE Vector

我正在使用的Python代码是:

import ctypes
import numpy as np

class _VectorType(ctypes.Structure):
    _fields_ = [('components',  ctypes.c_double*3)]


lib_gen_ctypes = '/local/scratch/jfreixa/lib/lib_ctypes_vector.so'

try_ctypes = ctypes.CDLL(lib_gen_ctypes,ctypes.RTLD_GLOBAL)

class vector(object):

    _ctypes_newVect = try_ctypes['Vector.newVect_']
    _ctypes_subVect = try_ctypes['Vector._subVect_']
    _ctypes_getVect = try_ctypes['Vector.getVect_']

    vector_pointer = ctypes.POINTER(_VectorType)

    _ctypes_getVect.argtypes = [vector_pointer,]
    _ctypes_getVect.restype  = _VectorType


    def __init__(self,*args):
        self._vector = _VectorType()
        self._newVect(*args)

    def _newVect(self,vectIn):
        pdb.set_trace()
        c_vect = (ctypes.c_double*3)(*vectIn)
        self._ctypes_newVect(self._vector,c_vect)

    def subVect(self):
        pdb.set_trace()
        c_vect =  _VectorType()
        self._ctypes_subVect(ctypes.byref(self._vector),ctypes.byref(c_vect))
        print c_vect.components[:]
        return np.array(c_vect.components[:])

    def getVect(self):
        pdb.set_trace()
        c_vect = self._ctypes_getVect(ctypes.byref(self._vector))
        vect = self._ctypes_getVect(self.vector_pointer.from_address(ctypes.addressof(c_vect)))
        print vect.components[:]
        return np.array(vect.components[:])

对于这个功能,我尝试过很多东西,但是我从未获得过正确的结果.要运行该程序,请尝试:

import pyctp.vector
newVect = pyctp.vector.vector((1.0,2.0,3.0))
newVect.subVect()
newVect.getVect()

子程序调用返回预期向量,而函数调用返回空向量或满垃圾的向量.

Answer 1

首先,您应该将bind(C)属性放在您希望从Python中可见的所有过程和类型中.所有Fortran类型都应该来自iso_c_binding,例如使用real(c_double)而不是double precision你将确保它是与C可互操作的类型.

MODULE Vector
use iso_c_binding
! Public types
TYPE,bind(C) :: VectorType
    real(c_double), DIMENSION(3):: components = 0.0d0
END TYPE VectorType

CONTAINS

!---------------------------------------------------------------------
SUBROUTINE newVect(this,vectorIn) bind(c,name="newVect")
   TYPE (VectorType),       INTENT(OUT):: this
   real(c_double), DIMENSION(3), INTENT(IN)::vectorIn

   this%components = (/vectorIn(1), vectorIn(2), vectorIn(3)/)

END SUBROUTINE newVect
!---------------------------------------------------------------------
SUBROUTINE subVect(this,vectorOut) bind(c,name="subVect")
   TYPE(VectorType), INTENT (OUT):: vectorOut
   TYPE(VectorType), INTENT (IN) :: this

   vectorOut%components = this%components

END SUBROUTINE subVect
!----------------------------------------------------------------------
TYPE(VectorType) FUNCTION getVect(this) RESULT(vectorOut) bind(c,name="getVect")

TYPE(VectorType), INTENT (IN) :: this

    vectorOut%components = this%components

END FUNCTION getVect
!--------------------------------------------------------------------
END MODULE Vector

然后在python中传递所有参数作为参考:

import ctypes
import numpy as np

class _VectorType(ctypes.Structure):
    _fields_ = [('components',  ctypes.c_double*3)]

lib_gen_ctypes = 'lib_ctypes_vector.so'
try_ctypes = ctypes.CDLL(lib_gen_ctypes,ctypes.RTLD_GLOBAL)

class vector(object):

    _ctypes_newVect = try_ctypes['newVect']
    _ctypes_subVect = try_ctypes['subVect']
    _ctypes_getVect = try_ctypes['getVect']

    vector_pointer = ctypes.POINTER(_VectorType)

    _ctypes_getVect.argtypes = [vector_pointer,]
    _ctypes_getVect.restype  = _VectorType

    def __init__(self,*args):
        self._vector = _VectorType()
        self._newVect(*args)

    def _newVect(self,vectIn):
        c_vect = (ctypes.c_double*3)(*vectIn)
        self._ctypes_newVect(ctypes.byref(self._vector),ctypes.byref(c_vect))

    def subVect(self):
        c_vect =  _VectorType()
        self._ctypes_subVect(ctypes.byref(self._vector),ctypes.byref(c_vect))
        return np.array(c_vect.components[:])

    def getVect(self):
        c_vect = self._ctypes_getVect(ctypes.byref(self._vector))
        return np.array(vect.components[:])

getVect的结果已经是VectorType,因此您可以直接访问其组件.