8 minutes

Written: 2021-10-02 04:19 +0000

# Simple Fortran Derived Types and Python

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.

### Series

Without intending it to be so, it appears that a series has been established, though not in strict order, the progression of the standards and their modern equivalents for Fortran-Python can be traced as seen in:

- NumPy Meson Fortran
**Simple Fortran Derived Types and Python**<– You are here!

## 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.

This is more surprising with each passing year ↩︎