This post is part of the Bridging Fortran & Python series.

Moving simple Fortran derived types to Python and back via C

Background

Object oriented programming has been part of Fortran for longer than I have been alive 1. Fortran has derived types now. They’ve been around for around for over three decades. The standards at the same time, have been supporting more and more interoperable operations. Details of these pleasant historical improvements are pretty much the most the Fortran standards committee have managed to date in the 21st century. The point of this exploration is to design interfaces to Python, mostly for the long term advancement of f2py, and just for the heck of it. Simple in the context of the title means that we consider only derived types whose members have direct C equivalents as defined in the Fortran 2003 Reid (2003) standard.

C Structures and Derived Types - I

Consider, first an example of what shall be a rather simple collection of records, the prototypical darling of most simulations, the point in Euclidean space, \(\mathbb{R}^{3}\).

1type, bind(c) :: cartesian
2   real(c_double) :: x,y,z
3end type cartesian

Now by definition, ever since the ISO_C_BINDING of the Fortran 2003 Reid (2003) standard, we know this is equivalent to:

1typedef struct {
2  double x,y,z;
3} cartesian;

To try this out, we need to have a procedure which operates on the type, due to a chronic lack of imagination, let us simply take a point and add one to each axis.

1subroutine unit_move(cartpoint) bind(c)
2  type(cartesian), intent(inout) :: cartpoint
3  print*, "Modifying the derived type now!"
4  cartpoint%x=cartpoint%x+1
5  cartpoint%y=cartpoint%y+1
6  cartpoint%z=cartpoint%z+1
7end subroutine unit_move

Now we can wrap these into a nice little module:

 1module vec
 2  use, intrinsic :: iso_c_binding
 3  implicit none
 4
 5  type, bind(c) :: cartesian
 6     real(c_double) :: x,y,z
 7  end type cartesian
 8
 9  contains
10
11  subroutine unit_move(array) bind(c)
12    type(cartesian), intent(inout) :: array
13    print*, "Modifying the derived type now!"
14    array%x=array%x+1
15    array%y=array%y+1
16    array%z=array%z+1
17  end subroutine unit_move
18
19end module vec

Call it with a simple main.

 1! main.f90
 2program main
 3  use, intrinsic :: iso_c_binding
 4  use vec
 5  implicit none
 6  type(cartesian) :: cartvec
 7  cartvec = cartesian(0.2,0.5,0.3)
 8  print*, cartvec
 9  call unit_move(cartvec)
10  print*, cartvec
11end program main

Finally we need a build system, lets say meson.

1# meson.build
2project('fcbase', 'fortran',
3  version : '0.1',
4  default_options : ['warning_level=3'])
5
6executable('fcbase',
7           'vec.f90',
8           'main.f90',
9           install : true)

Now we can finally try this out.

1meson setup testfc
2meson compile -C testfc
3./testfc/fcbase
4  0.20000000298023224       0.50000000000000000       0.30000001192092896
5 Modifying the derived type now!
6   1.2000000029802322        1.5000000000000000        1.3000000119209290

This probably surprised no one. So let’s get to something more interesting.

C and Fortran

We can start with a simple C header.

1#ifndef VECFC_H
2#define VECFC_H
3
4typedef struct {
5  double x,y,z;
6} cartesian;
7
8#endif /* VECFC_H */

Along with its associated translation unit.

 1/* vecfc.c */
 2#include<stdlib.h>
 3#include<stdio.h>
 4#include<vecfc.h>
 5
 6void* unit_move(cartesian *word);
 7
 8int main(int argc, char* argv[argc+1]) {
 9    puts("Initializing the struct");
10    cartesian a={3.0, 2.5, 6.2};
11    printf("%f %f %f",a.x,a.y,a.z);
12    puts("\nFortran function with derived type from C:");
13    unit_move(&a);
14    puts("\nReturned from Fortran");
15    printf("%f %f %f",a.x,a.y,a.z);
16    return EXIT_SUCCESS;
17}

The function declaration is straightforward, but for the fact that since we are passing something both in and out of the function, it must be a pointer, and as we have no real error propagation or feasible return type, void * makes sense too. The calling semantics are equally transparent, creating an instance of the struct and then passing it by reference.

The meson extension is fairly trivial, we only need to add the right file and declare the project language correctly.

 1project('cderived', 'c', 'fortran',
 2  version : '0.1',
 3  default_options : ['warning_level=3'])
 4
 5executable('cderived',
 6           'vec.f90',
 7           # 'vecfc.c',
 8           # 'vecfc.h',
 9           'main.f90',
10           install : true)

This seems to work.

 1meson setup testfc
 2meson compile -C testfc
 3./testfc/fcbase
 4Initializing the struct
 53.000000 2.500000 6.200000
 6Fortran function with derived type from C:
 7 Modifying the derived type now!
 8
 9Returned from Fortran
104.000000 3.500000 7.200000%

Cython and Fortran

So far so good. Let’s get to something slightly more interesting. Although in the long run cython isn’t really the approach numpy is going towards at the moment, it still makes sense to give it a shot as a first approximation.

 1cdef struct cartesian:
 2    double x,y,z
 3
 4cdef extern:
 5    void* unit_move(cartesian *word)
 6
 7def unitmv(dpt={'x':0,'y':0,'z':0}):
 8    cdef cartesian A=dpt
 9    unit_move(&A)
10    return(A)

The most interesting aspect is that, having defined our C structure and the external function, we need something to offer to Python. Since cython structures can be mapped from simple python dictionaries, we take an input of one, construct the structure on the fly in the function and then use our Fortran subroutine. The return type can also be coerced into a list, and that’s what can be printed or manipulated further. This isn’t the most efficient way to work with this, but it will do for now.

The meson update follows from the explanation in the previous post, namely, that we need to include our Python sources.

 1project('pointeuc', 'c', 'fortran',
 2  version : '0.1',
 3  default_options : ['warning_level=3'])
 4
 5add_languages('cython')
 6
 7py_mod = import('python')
 8py3 = py_mod.find_installation()
 9py3_dep = py3.dependency()
10
11incdir_numpy = run_command(py3,
12  ['-c', 'import os; os.chdir(".."); import numpy; print(numpy.get_include())'],
13  check : true
14).stdout().strip()
15
16inc_np = include_directories(incdir_numpy)
17
18py3.extension_module('veccyf',
19                     'veccyf.pyx',
20                     'vec.f90',
21                     include_directories: inc_np,
22                     dependencies : py3_dep)

The NumPy dependency in this situation is really optional, but makes for a more interoperable structure. Now we’re ready to try this out in an interactive manner.

1meson setup testcpyth
2meson compile -C testcpyth
3cd testcpyth
4python -c "import veccyf; A=veccyf.unitmv({'x':3.3,'y':5.2,'z':2.6}); print(A)"
5 Modifying the derived type now!
6{'x': 4.3, 'y': 6.2, 'z': 3.6}

Taking this a bit slower, it is still as neat as ever.

1import veccyf
2point = {'x': 4.3,
3         'y': 2.6,
4         'z': 6.4}
5opoint = veccyf.unitmv(point)
6print(opoint)

Aside from the output string, this works as expected, and is pretty neat, all things considered. This might feed into an f2py workflow, but then, it adds a layer of indirection which seems overly circuitous. The design could have also been far improved by using a bind(c, name=uniquname) as well.

Note that a much more rational approach would have been to generate a class in Python.

Prognosis for F2PY

The essential information used from the fortran code for the cython types (and C structures) was essentially:

  • The types of each variable for the derived type
    • double for every one of them in this case
  • The names of each variable
    • x, y, and z here

Apart from this, for functions which operate on aforementioned types,we require:

  • The external declaration (due to the iso_c_binding)
  • A mapping from structure to a python object, typically a dictionary

It so happens that much of the information needed for generation of the relevant C code is already present in f2py, and can be inspected via the crackfortran module.

 1# python -m numpy.f2py.crackfortran vec.f90 -show
 2Reading fortran codes...
 3      Reading file 'vec.f90' (format:free)
 4{'before': '', 'this': 'use', 'after': ', intrinsic :: iso_c_binding '}
 5Line #2 in vec.f90:"  use, intrinsic :: iso_c_binding "
 6      analyzeline: Could not crack the use statement.
 7Line #5 in vec.f90:"  type, bind(c) :: cartesian "
 8      analyzeline: No name/args pattern found for line.
 9Post-processing...
10      Block: vec
11              Block: unknown_type
12              Block: unit_move
13Post-processing (stage 2)...
14[{'args': [],
15  'block': 'module',
16  'body': [{'args': [],
17            'block': 'type',
18            'body': [],
19            'entry': {},
20            'externals': [],
21            'from': 'vec.f90:vec',
22            'interfaced': [],
23            'name': 'unknown_type',
24            'parent_block': <Recursion on dict with id=4500228096>,
25            'sortvars': ['x', 'y', 'z'],
26            'varnames': ['x', 'y', 'z'],
27            'vars': {'x': {'kindselector': {'kind': 'c_double'},
28                           'typespec': 'real'},
29                     'y': {'kindselector': {'kind': 'c_double'},
30                           'typespec': 'real'},
31                     'z': {'kindselector': {'kind': 'c_double'},
32                           'typespec': 'real'}}},
33           {'args': ['array'],
34            'block': 'subroutine',
35            'body': [],
36            'entry': {},
37            'externals': [],
38            'from': 'vec.f90:vec',
39            'interfaced': [],
40            'name': 'unit_move',
41            'parent_block': <Recursion on dict with id=4500228096>,
42            'saved_interface': '\n'
43                               '      subroutine unit_move(array) \n'
44                               '          type(cartesian), intent(inout) :: '
45                               'array\n'
46                               '      end subroutine unit_move',
47            'sortvars': ['array'],
48            'vars': {'array': {'attrspec': [],
49                               'intent': ['inout'],
50                               'typename': 'cartesian',
51                               'typespec': 'type'}}}],
52  'entry': {},
53  'externals': [],
54  'from': 'vec.f90',
55  'implicit': None,
56  'interfaced': [],
57  'name': 'vec',
58  'sortvars': [],
59  'vars': {}}]

The information is clearly present, which brings us to the next phase.

Python-C and Fortran

The most natural way to feed into f2py would be to directly write out an equivalent Python-C API which effectively mimics the extension module of the last section.

The crucial definition essentially takes a dictionary, unrolls into a list, creates the C structure, then passes it through Fortran and then back.

 1static PyObject* eucli_unitinc(PyObject* self, PyObject* args) {
 2    cartesian a;
 3    PyObject* dict;
 4    PyArg_ParseTuple(args, "O", &dict);
 5    PyObject* vals = PyDict_Values(dict);
 6    a.x = PyFloat_AsDouble(PyList_GetItem(vals,0));
 7    a.y = PyFloat_AsDouble(PyList_GetItem(vals,1));
 8    a.z = PyFloat_AsDouble(PyList_GetItem(vals,2));
 9    // Call Fortran on it
10    unit_move(&a);
11    //
12    PyObject* ret = Py_BuildValue("{s:f,s:f,s:f}",
13                         "x", a.x,
14                         "y", a.y,
15                         "z", a.z);
16    return ret;
17}

This also works well enough with minor modifications to meson.

1python -c "import eucli; A=eucli.unitinc({'x':3.3,'y':5.2,'z':2.6}); print(A)"
2 Modifying the derived type now!
3{'x': 4.3, 'y': 6.2, 'z': 3.6}

Boilerplate

The rest of the hand written wrapper is boilerplate, but the full file is reproduced here:

 1/* Make "s#" use Py_ssize_t rather than int. */
 2#ifndef PY_SSIZE_T_CLEAN
 3#define PY_SSIZE_T_CLEAN
 4#include <stdio.h>
 5#endif /* PY_SSIZE_T_CLEAN */
 6
 7#include "Python.h"
 8
 9// Fortran-C stuff
10
11typedef struct {
12  double x,y,z;
13} cartesian;
14
15void* unit_move(cartesian *word);
16
17// Python-C stuff
18
19static PyObject* eucli_unitinc(PyObject* self, PyObject* args) {
20    cartesian a;
21    PyObject* dict;
22    PyArg_ParseTuple(args, "O", &dict);
23    PyObject* vals = PyDict_Values(dict);
24    a.x = PyFloat_AsDouble(PyList_GetItem(vals,0));
25    a.y = PyFloat_AsDouble(PyList_GetItem(vals,1));
26    a.z = PyFloat_AsDouble(PyList_GetItem(vals,2));
27    // Call Fortran on it
28    unit_move(&a);
29    //
30    PyObject* ret = Py_BuildValue("{s:f,s:f,s:f}",
31                         "x", a.x,
32                         "y", a.y,
33                         "z", a.z);
34    return ret;
35}
36
37static PyMethodDef EucliMethods[] = {
38    {"unitinc", (PyCFunction)eucli_unitinc, METH_VARARGS,
39     "Add 1 to the axes of a point"},
40    {NULL, NULL, 0, NULL}        /* Sentinel */
41};
42
43static struct PyModuleDef eucli_module = {
44   PyModuleDef_HEAD_INIT,
45   "eucli",
46   "blah",
47   -1,
48   EucliMethods
49};
50
51PyMODINIT_FUNC PyInit_eucli(void)
52{
53    PyObject *m;
54    m = PyModule_Create(&eucli_module);
55    if (!m) {
56        return NULL;
57    }
58
59    return m;
60}

Conclusions

Slowly but surely, compilers have caught up with to the standards, and soon, with some support, numpy shall too. A subsequent post shall cover more exotic derived types, typically (until 2018) treated by non-portable opaque pointer tricks. There are still some rather awkward design decisions in this case; including things like considering the mapping and their interoperability with NumPy arrays, but by and large this would depend on the functions which use these. In terms of derived types of simple structures, without allocatable arrays, the extensions appear to be straightforward.

References

Reid, John. 2003. “The New Features of Fortran 2003,” 38. https://wg5-fortran.org/N1551-N1600/N1579.pdf.

  1. This is more surprising with each passing year ↩︎

Series info

Bridging Fortran & Python series

  1. NumPy, Meson and f2py
  2. Simple Fortran Derived Types and Python <-- You are here!
  3. Exploring ISO_C_BINDING and type-bound procedures
  4. Fortran OOP and Python
  5. Types from Fortran to Python via Opaque Pointers