12

Experience:
fortran for about 3 months
python - intermediate : never used the ctypes module in python before this

I was looking for a way to use the fortran code for my doctoral work in python - subsequently using the computation for visualizations on-the-fly for visualizations using matplotlib.

THIS POST helped (this tells that fortran code can be used/invoked in python using the ctypes module - and given that the fortran functions have alternate names bound to them - this much makes sense to me logically although i do not know how this works in detail. But we DO choose our battles wisely!).
Then this SO post deals with calling fortran functions from python as well.

The next logical step was to look up the documentation for the python module ctypes. This talks about how the shared library can be accessed using python at an API level.

I had all the pieces to make a minimal working example, which another answer has already done. But i wanted to look at the output mechanism and mathematical operations involving real floats. Here is the test case i made.

test.f90

function prnt(s)
    character(80):: s
    logical :: prnt
    print*, s
    prnt = .true.
end function prnt

function sin_2(r)
    real:: r,sin_2
    sin_2 = sin(r)**2
end function sin_2

make the shared object

$gfortran -shared  -g -o test.so test.f90

EDIT: for some reason my work computer needs the -fPIC option to compile

To make sure my two functions prnt and sin_2 are in there somewhere, i checked with nm: 

$ nm test.so | tail -3  
0000067f T prnt_
0000065c T sin_2_
         U sinf@@GLIBC_2.0

So far so good. My functions prnt and sin_2 have been mapped onto prnt_ and sin_2_ inside the library.

invoke fortran functions from the python interpreter

this is where all of this gets a bit soggy. Using the table in python-ctypes documentation, I typed in the following - 

>>> from ctypes import byref, cdll, c_float,c_char_p
>>> t = cdll.LoadLibrary('./test.so')
>>> c = c_char_p("Mary had a little lamb")
>>> t.prnt_('Mary had a little lamb')
 Mary had a little lambe
1
>>> t.prnt_("Mary had a little lamb")
 Mary had a little lambe
1
>>> t.prnt_(c)
 Mary had a little lambe[�  .prnt_(c)

1

I suppose that the 1 printed at the end of each output is python's way of letting me know that the boolean output to t.prnt_ is .true..
I am a bit worried about how the situation gets worse with t.prnt_ when i switch to the proper datatype for the string. The literals print okay, only with a e at the end. Is that the EOL character?

Then there is the t.sin_2_ function. I have decided to use it to compute sin(4.56)**2. here is how that went -

>>> f = c_float(4.56)
>>> t.sin_2_(4.56)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
 ctypes.ArgumentError: argument 1: <type 'exceptions.TypeError'>: Don't know how to convert parameter 1
>>> t.sin_2_(f)
Segmentation fault (core dumped)

Where did i go wrong here? I tried to explain how I approached the problem so that it would be apparent if I made an obvious gaffe somewhere.
The multitude of links to other SO posts is to help other people who ask the same question as I do now.

Community
  • 1
  • 1
Debanjan Basu
  • 834
  • 1
  • 8
  • 29

3 Answers3

12

In Fortran, arguments are passed by reference. Fortran character arrays aren't null terminated; the length is passed by value as an implicit long int argument. Also, Python's float type is a double, so you may want to use Fortran real(8) for consistency.

test.f90:

function prnt(s)          ! byref(s), byval(length) [long int, implicit]
    character(len=*):: s  ! variable length input
    logical :: prnt
    write(*, "(A)") s     ! formatted, to remove initial space
    prnt = .true.
end function prnt

function sin_2(r)         ! byref(r)
    real:: r, sin_2       ! float; use real(8) for double
    sin_2 = sin(r)**2
end function sin_2

Remember to set ctypes argtypes for functions, and restype where appropriate. In this case sin_2 takes a float pointer and returns a float.

ctypes example:

>>> from ctypes import *
>>> test = CDLL('./test.so')

>>> test.prnt_.argtypes = [c_char_p, c_long]
>>> test.sin_2_.argtypes = [POINTER(c_float)]
>>> test.sin_2_.restype = c_float

>>> s = 'Mary had a little lamb'
>>> test.prnt_(s, len(s))
Mary had a little lamb
1

>>> x = c_float(4.56)
>>> test.sin_2_(byref(x))
0.9769567847251892
Eryk Sun
  • 33,190
  • 5
  • 92
  • 111
  • a tiny little follow up question - how is a `-shared` object different from the generic object file? this is the only part of the workflow i dont have a logic for. – Debanjan Basu Apr 11 '13 at 10:05
  • 2
    Compile with the -v option to see what `gfortran` is doing. It calls the compiler `f951` to create a .s assembler file. Next it calls the [assembler `as`](http://sourceware.org/binutils/docs-2.23.1/as) to create a .o object file. Finally it calls [`collect2`](http://gcc.gnu.org/onlinedocs/gccint/Collect2.html), which calls the [linker `ld`](http://sourceware.org/binutils/docs-2.23.1/ld) to create an ELF .so shared object. Look at the output of `readelf -d test.so` to see the linked libs such as libgfortran. – Eryk Sun Apr 11 '13 at 13:17
  • 1
    Compiling with the -S option will output the .s assembler file. Compiling with -c will output an ELF .o RELocatable object file. The .so shared object has the needed dynamic linking headers for the runtime linker/loader `ld.so` to load it into the process image at runtime (or using `dlopen` in Python to import an extension module, or with ctypes). If you add a `program` section to test.f90, you can remove the -shared option to instead compile an ELF executable. Use `readelf` to compare the headers and symbol tables for all 3 ELF file types. – Eryk Sun Apr 11 '13 at 13:31
  • you have been very patient with my questions and very helpful. I have some added questions to ask you for understanding your answer that i couldn't find using google. Here they are - maybe you would consider editing your answer to explain the following. 1. when you define `function prnt(s)` - there is no len(s) argument there - optional or not! 2. Why did you use a reference to the argument for `sin_2` but not for `prnt`? Is there a rule about this? – Debanjan Basu Apr 12 '13 at 09:52
  • The length argument is implicit, not optional. It's a behind the scenes implementation detail of Fortran. – Eryk Sun Apr 12 '13 at 14:59
  • 1
    Argument `s` is also by reference. CPython uses C null-terminated strings (a pointer to a null-terminated array of `char`). So the `c_char_p` constructor (I set it in `argtypes`) can simply reuse the string object's internal pointer, which is why you never use a `c_char_p` for a function that mutates the string. Instead use `create_string_buffer` to create a mutable `c_char` array. – Eryk Sun Apr 12 '13 at 15:02
  • 1
    Since I set `sin_2_.argtypes = [POINTER(c_float)]`, I could have used `sin_2_(c_float(4.56))`. ctypes is smart enough to create the pointer implicitly, but in this case I think it's clearer in the code to use `byref` explicitly. Notice also that the `c_float` return type is automatically converted to a Python `float`, which has double precision. It would be better to use double precision in Fortran, too, to keep things consistent. Then use a ctypes `c_double`. Alternatively use NumPy on the Python side, which has a `float32` dtype. – Eryk Sun Apr 12 '13 at 15:10
5

I use f2py pretty frequently, its quite simple. If your code complies to Fortran 90 declarations (e.g. double precision, intent(in) :: myvar), f2py will wrap your code automatically to C and is compiled into a shared object that is directly callable through Python. You can import your f2py created module as you would any other python module. The C wrapper is transparent and handles the type interfacing between Python and Fortran, you don't have to mess with c_types at all. I would definitely recommend this route.

John Warner
  • 51
  • 1
  • 1
    while you are technically correct, f2py should be easier to handle. But at the same time, python-ctypes makes more sense to me - although i have to deal with the learning curve now - i can assign datatypes to input and output more explicitly this way - and explicit is good for me if i want to code in python!! :) this is why i went for ctypes rather than f2py. it is a personal preference - and i hope your answer benefits those who come looking for answers to the same question and have the choice of a different way. NOTE: i am not an experienced programmer at all!! – Debanjan Basu Apr 11 '13 at 08:41
  • I would agree with you, you always choose the method that makes sense to you. If you're worried about specifying data types, these are dictated by the Fortran code and are enforced by the C wrapper f2py generates (i.e. you get a type error if your input/output variables don't match what is expected by the fortran code). – John Warner Apr 11 '13 at 14:45
3

I would just like to point out that one should always use the iso_c_binding module and the bind attribute when interfacing Python and Fortran using ctypes. That way you can be sure that your variable types match. And the compiler is free to mangle names however it likes when building the shared library so you might get undefined references when switching compiler vendor or even just using a different version of your compiler and trying to call t.sin_2_ as per the example. See also here.

So you should do something like:

module test
  use iso_c_binding
  implicit none

contains

  function sin_2(r) bind(c, name='c_sin_2')
    real(c_float) :: sin_2
    real(c_float), intent(in), value :: r

    sin_2 = sin(r)**2
  end function sin_2
end module test

And then

from ctypes import *

test = CDLL('./test.so')

test.c_sin_2.restype = c_float

print(test.c_sin_2(c_float(3)))

I'm not quite sure, though, why I still have to set the restype, maybe someone with more experience can tell us that.

Here you can find some additional info.

user35915
  • 1,222
  • 3
  • 17
  • 23