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.

Why not use f2py? Or even NumPy if that's what you're going for? — Ignacio Vazquez-Abrams, Sep 22 '15 at 10:59
Or, if you want to use ctypes instead of f2py, use iso_c_bindig internal module in fortran to expose a plain C-like interface to the routines. — Edmondo Giovannozzi, Sep 22 '15 at 11:15
Plain C-like interface is exposed by `bind(c)`, not by the module. The module just defines some types and auxiliary procedures to convert pointers. Using the module does NOT change the calling convention. — Vladimir F Героям слава, Sep 22 '15 at 12:58

score 5 · Accepted Answer · answered Sep 22 '15 at 15:10

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.

Thanks Edmondo!!! The point is that I want to avoid the modification of Fortran source code (I understandt that make the Fortran output compatible with C types is the best way to use ctypes). Is there any way to obtain this with just python tricks?? (the point is that I have tons of fortran modules, and I think that it is not affordable to modify all of them, validate them .... ) — Jfreixa, Sep 23 '15 at 08:45
As suggested by other people, you can use f2py. You may not be able to access user defined component (I haven't checked if it is possible), but there is no problem with arrays, and module variables. At least you may have to write some additional routines to access the components. I have used it many time, it is quite easy to use. Have a look and see if it can help you. — Edmondo Giovannozzi, Sep 23 '15 at 08:57

Can't get fortran function output from ctypes

1 Answers1