Can't get fortran function output from ctypes

Question

I'm trying to call a Fortran function from Python using ctypes. I have tried to obtain the result from a subroutine and from a function (both with the same functionality), but I can't obtain the expected output from the function whereas the subroutine works well. The problem is that I have a lot of libraries with Fortran functions instead of subroutines. Is there any problem with Fortran functions and ctypes?

Piece of Fortran code:

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

The Python code I'm using is:

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[:])

For the function I have tried a lot of things but I have never obtained the correct result. To run the piece of program try with:

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

The subroutine call returns the expected vector while the function call returns a null vector or a vector full of garbage.

Answer 1

First of all you should put the bind(C) attributes to all the procedures and types that you want to be visible from Python. All the Fortran type should be taken from iso_c_binding, for example using real(c_double) instead of double precision you will be sure that it is of the type interoperable with 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

Then in python pass all the arguments as reference:

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[:])

the results from getVect is already of VectorType, so you can access its components directly.