Calling Fortran from C, C++ and Python

This content was generated with the assistance of AI and should be verified against the official documentation before being relied on in production.

The reference material consulted while preparing these pages is listed in the section bibliography.

This page is written against Fortran 2023 (ISO/IEC 1539-1:2023, J3/23-007r1) and was checked on 2026-10-19. Every Fortran source on it was compiled warning-free with gfortran -std=f2023 -Wall -Wextra -pedantic -fimplicit-none -fcheck=all (GNU Fortran 16.2) and also with LLVM Flang 23. The C and C++ clients were built with GCC 16.2 (gcc -std=c17 and g++ -std=c++23), the Python examples ran on Python 3.13 with NumPy 2.5.4, Meson 1.12.1 and cffi 2.1.1 in a virtual environment, and the CMake example used CMake 3.31.6 with Ninja. The output shown is real output. Examples that fail on purpose, to show a linker or runtime diagnostic, are labelled. The Julia and R snippets near the end are taken from the projects' own manuals and were not run here. The language side of the boundary, the bind(c) rules and the iso_c_binding module, is covered in Interoperability with C; assumed-shape, allocatable and pointer arguments seen from C are in C Descriptors and ISO_Fortran_binding.

This page is about the other direction: a program or script written in another language that wants the numerical kernel of heatflow, the running example. The kernel is the explicit FTCS step in heatflow_kernels, exposed to the outside world through the heatflow_capi module as one function, heatflow_step. Everything below calls that function, or a small variant of it, from C, from C++, from Python through f2py, and from Python through ctypes and cffi. The same function is also the entry point for MATLAB in MATLAB: MEX, Engine and MAT-File APIs. How the libraries are built is in Building with Make, CMake and Meson, and the compiler flags used on the command lines are explained in Compilers and Compiler Options.

Choosing a bridge

Fortran has no foreign-function interface of its own that other languages target; every bridge is built on the same fact: a compiled Fortran procedure is an object-code symbol with a calling convention, and a bind(c) procedure (18.10) has the symbol name and the C calling convention that the C compiler would give an equivalent C function. What differs between the bridges is who writes the glue.

Bridge Fortran side needs Glue written by Typical use

C or C++ client

bind(c) interface and a header

you, once, as a header with extern "C"

applications and libraries in C or C++ that reuse a kernel

f2py

plain Fortran, intent on every dummy; bind(c) not needed

f2py generates a CPython extension from the source or a signature file

NumPy code that wants a compiled kernel with array arguments

ctypes or cffi

bind(c) procedures in a shared library

you, as a few lines of Python declaring the C prototype

loading an existing shared library without compiling an extension

Julia ccall

bind(c) procedures in a shared library

you, at each call site

Julia programs

R .Fortran or .Call

bind(c) or legacy external procedures

a small R wrapper, or C glue for .Call

R packages with compiled code

MATLAB MEX

a gateway mexFunction

you, in Fortran

MATLAB calling Fortran, see the MATLAB page

A bind(c) library with a C header is the lowest common denominator: every other row, except f2py and MEX, can consume it, and f2py does not need it. When a project must serve several languages, write the bind(c) layer and the header first and put the others on top.

flowchart TD A["Fortran kernel: heatflow_kernels"] --> B["bind(c) wrapper: heatflow_capi"] B --> C["C header heatflow.h"] C --> D["C client"] C --> E["C++ client with extern C"] B --> F["shared library libheatflow.so"] F --> G["Python ctypes or cffi"] F --> H["Julia ccall or R .Fortran"] A --> I["f2py wrapper: kernel.f90"] I --> J["Python extension module"]

The C API of heatflow

The C-visible part of heatflow is one module and one header. The module is the one used in Modules and Submodules; it uses the FTCS sweep from heatflow_kernels and exports a single bind(c) function:

!> C API: heatflow_step advances an nx-by-ny field by one FTCS step.
module heatflow_capi
  use, intrinsic :: iso_c_binding, only: c_int, c_double
  use heatflow_kernels, only: ftcs_sweep
  implicit none
  private
  public :: heatflow_step

contains

  function heatflow_step(nx, ny, u, unew, alpha, dt, dx, dy, residual) &
      bind(c, name='heatflow_step') result(status)
    integer(c_int), value, intent(in) :: nx, ny
    real(c_double), intent(in) :: u(nx, ny)
    real(c_double), intent(inout) :: unew(nx, ny)
    real(c_double), value, intent(in) :: alpha, dt, dx, dy
    real(c_double), intent(out) :: residual
    integer(c_int) :: status
    if (nx < 3 .or. ny < 3) then
      residual = 0.0_c_double
      status = 1_c_int
      return
    end if
    call ftcs_sweep(u, unew, alpha, dt, dx, dy, residual)
    status = 0_c_int
  end function heatflow_step
end module heatflow_capi

Source: J3/23-007r1: 18.3.7 Interoperability of procedures and procedure interfaces and 18.10.1 Definition and reference of interoperable procedures, GNU Fortran manual: Interoperability with C

Four decisions in this function matter to every caller.

  • Sizes and scalars are passed by value. nx, ny, alpha, dt, dx and dy have the value attribute, so the C prototype takes int and double, not pointers. Without value a bind(c) dummy is passed by address, as in every other Fortran call (18.3.7).

  • The arrays are explicit-shape, not assumed-shape. u(nx, ny) is interoperable with a plain C pointer to the first element, because an explicit-shape array has no descriptor (18.3.6). Assumed-shape u(:, :) would need a C descriptor, and a C compiler cannot build one without ISO_Fortran_binding.h. The assumed-shape routine ftcs_sweep is called inside, after the dummy arguments have been given their shape.

  • The result is a status code. A Fortran procedure cannot throw into a C caller. Failures travel as an integer; heatflow_step returns 1 for a grid smaller than 3 by 3. stop and error stop inside a library routine would end the caller’s process, so they stay out of the C-facing layer.

  • bind(c, name='heatflow_step') fixes the symbol. The name is the binding label (18.10.2). Without name=, the label is the lowercase procedure name; with it, the label is exactly the string, whatever the compiler would otherwise generate. Without bind(c) the symbol is whatever the compiler chooses; see the section on legacy symbols below.

The header is the contract the C and C++ programmers see:

#ifndef HEATFLOW_H
#define HEATFLOW_H
#ifdef __cplusplus
extern "C" {
#endif
/* Advance the nx-by-ny column-major field u by one FTCS step into unew.
   Returns 0 on success, 1 if the grid is smaller than 3x3. */
int heatflow_step(int nx, int ny, const double *u, double *unew,
                  double alpha, double dt, double dx, double dy,
                  double *residual);
#ifdef __cplusplus
}
#endif
#endif

The correspondence is mechanical: integer(c_int) is int, real(c_double) is double, intent(in) becomes const, intent(out) and intent(inout) become non-const pointers, and an array dummy without value becomes a pointer to its first element. Keep the comment about column-major order in the header; the compiler cannot check it.

Source: J3/23-007r1: 18.3.7 Interoperability of procedures and procedure interfaces and 18.10.1 Definition and reference of interoperable procedures, GNU Fortran manual: Interoperability with C

A C client

The client below fills a 21 by 21 field with a single hot point in the middle, advances it one step and prints the result. The index expression (NX / 2) + NX * (NY / 2) is the column-major address of u(11, 11), zero-based.

#include <stdio.h>
#include <stdlib.h>
#include "heatflow.h"

int main(void) {
    enum { NX = 21, NY = 21 };
    double u[NX * NY] = {0}, unew[NX * NY] = {0}, residual = 0.0;
    u[(NX / 2) + NX * (NY / 2)] = 100.0;
    unew[(NX / 2) + NX * (NY / 2)] = 100.0;
    int status = heatflow_step(NX, NY, u, unew, 1.0, 5.0e-5, 0.05, 0.05, &residual);
    printf("status=%d residual=%.4f center=%.4f\n", status, residual, unew[(NX / 2) + NX * (NY / 2)]);
    return status;
}

unew is initialised with the same hot point on purpose: ftcs_sweep writes only the interior points, so the boundary values of unew come from the caller. This is the same convention the Fortran solvers follow, and it is the first thing a C caller forgets.

Build it in three steps: compile the Fortran into an archive, compile the client, and link. The first attempt below is the instructive one.

gfortran -std=f2023 -Wall -Wextra -pedantic -fimplicit-none -fcheck=all -c \
  src/heatflow_kinds.f90 src/heatflow_kernels.f90 src/heatflow_capi.f90
ar rcs libheatflow_capi.a heatflow_kinds.o heatflow_kernels.o heatflow_capi.o
gcc -std=c17 -Wall -Wextra -pedantic -Iinclude c_client.c -L. -lheatflow_capi -o c_client

This link fails, and the diagnostic is the real one (the first lines; the others repeat the same symbol):

/usr/bin/ld: ./libheatflow_capi.a(heatflow_capi.o): in function `heatflow_step':
heatflow_capi.f90:(.text+0x50): undefined reference to `_gfortran_runtime_error_at'
/usr/bin/ld: heatflow_capi.f90:(.text+0x17c): undefined reference to `_gfortran_runtime_error_at'
collect2: error: ld returned 1 exit status

The C client uses no Fortran feature, but the Fortran objects were compiled with -fcheck=all, which inserts calls to the GNU Fortran runtime library for the bounds and pointer checks. The linker is gcc, which knows nothing about libgfortran. Adding it fixes the link:

gcc -std=c17 -Wall -Wextra -pedantic -Iinclude c_client.c -L. -lheatflow_capi -lgfortran -o c_client
./c_client
status=0 residual=8.0000 center=92.0000

Recompiling the Fortran with -O2 and without -fcheck=all removes every reference to _gfortran_* from heatflow_step: nm -g heatflow_capi.o | grep -c ' U _gfortran' printed 0, and the client linked without -lgfortran and printed the same line. A library that works without the runtime today may need it tomorrow, when someone adds a print, an allocate, a stop or a runtime check, so a build that links Fortran objects with a C driver should always name the runtime explicitly.

Source: GNU Fortran manual: Influencing the linking step, Influencing runtime behavior

Linking: order, drivers and runtime libraries

Linking Fortran with C or C++ goes wrong in a few repeatable ways.

Order on the command line

A static archive is searched once, left to right, and only for symbols that are undefined at that point. The client object must therefore come before the archive:

gcc -std=c17 -Iinclude -L. -lheatflow_capi c_client.c -o bad
/usr/bin/ld: /tmp/ccjy7LQP.o: in function `main':
c_client.c:(.text+0x9c): undefined reference to `heatflow_step'
collect2: error: ld returned 1 exit status

Move -lheatflow_capi after c_client.c, and the runtime library after the archive that needs it, as in the working command above. A shared library given with -l is less sensitive on some systems, but the portable rule is the same: objects, then libraries, with each library after the objects that call it.

Which runtime does a Fortran object need?

The runtime library is a property of the compiler that produced the object, not of the language.

Compiler Runtime to name when the driver is not the Fortran compiler Observed here

GNU Fortran 16.2

-lgfortran (and -lm)

the driver adds -lgfortran -lm -lgcc_s -lgcc -latomic_asneeded -lm -lc after the objects (gfortran -v)

LLVM Flang 23

-lflang_rt.runtime in the directory /usr/lib/llvm-23/lib/clang/23/lib/aarch64-unknown-linux-gnu on this machine

flang -v printed -lflang_rt.runtime -lm -lgcc in the link line; the directory depends on the installation

A routine that does I/O makes the point clearly. This one prints a residual, exactly as a Fortran kernel might log one:

module hf_report
  use, intrinsic :: iso_c_binding, only: c_int, c_double
  implicit none
  private
  public :: hf_print_residual
contains
  subroutine hf_print_residual(step, residual) bind(c, name='hf_print_residual')
    integer(c_int), value, intent(in) :: step
    real(c_double), value, intent(in) :: residual
    print '(a, i0, a, es10.3)', 'step ', step, ' residual ', residual
  end subroutine hf_print_residual
end module hf_report

Source: GNU Fortran manual: Influencing the linking step, Flang documentation: Flang driver

The C side is two lines:

void hf_print_residual(int step, double residual);

int main(void) {
    hf_print_residual(7, 2.514e-3);
    return 0;
}

With gcc as the linker and no runtime, the six I/O entry points the print statement uses are all unresolved:

/usr/bin/ld: hf_report.o: in function `hf_print_residual':
hf_report.f90:(.text+0x50): undefined reference to `_gfortran_st_write'
/usr/bin/ld: hf_report.f90:(.text+0x68): undefined reference to `_gfortran_transfer_character_write'
/usr/bin/ld: hf_report.f90:(.text+0x78): undefined reference to `_gfortran_transfer_integer_write'
/usr/bin/ld: hf_report.f90:(.text+0x90): undefined reference to `_gfortran_transfer_character_write'
/usr/bin/ld: hf_report.f90:(.text+0xa0): undefined reference to `_gfortran_transfer_real_write'
/usr/bin/ld: hf_report.f90:(.text+0xa8): undefined reference to `_gfortran_st_write_done'
collect2: error: ld returned 1 exit status

Three ways to link it, all of which printed the same line on GNU Fortran 16.2:

gcc hf_report_main.c hf_report.o -lgfortran -o report      # the C driver plus the runtime
gfortran hf_report_main.c hf_report.o -o report            # the Fortran driver: it adds the runtime itself
./report
step 7 residual  2.514E-03

With Flang 23 the same object, compiled by flang -c, references _FortranAio... symbols instead of _gfortran_*, so an object compiled by flang and linked by gcc with no runtime fails the same way:

/usr/bin/ld: hf_report.o: in function `hf_print_residual':
FIRModule:(.text+0x38): undefined reference to `_FortranAioBeginExternalFormattedOutput'
/usr/bin/ld: FIRModule:(.text+0x50): undefined reference to `_FortranAioOutputAscii'

Source: GNU Fortran manual: Influencing the linking step, Flang documentation: Flang driver

and either the Flang driver (flang main.o hf_report.o -o report, where main.o was compiled by gcc -c) or gcc with -L/usr/lib/llvm-23/lib/clang/23/lib/aarch64-unknown-linux-gnu -lflang_rt.runtime -lm printed the same step 7 residual 2.514E-03. Passing a .c file to flang itself failed on this machine with error: unknown integrated tool '-cc1'; the Flang driver links, but it does not compile C.

Two rules follow. Link with the Fortran compiler’s driver whenever the program contains any Fortran, because it knows the right runtime, the right order and any further libraries it needs. And never mix runtimes: a program that uses objects from gfortran and from flang needs both runtimes, and their I/O units and exit handling do not know about each other. A C++ main program linked with g++ needs -lgfortran added by hand when the Fortran objects use the runtime, by the same rule as for gcc; linking C++ objects with gfortran needs the C++ runtime (-lstdc++). The C++ example below linked against an archive that needed no runtime, and these two link lines were not run.

Shared libraries and the loader

For Python, Julia, R or any dynamic loader, build a position-independent shared library:

gfortran -std=f2023 -Wall -Wextra -pedantic -fimplicit-none -O2 -fPIC -shared \
  src/heatflow_kinds.f90 src/heatflow_kernels.f90 src/heatflow_capi.f90 -o libheatflow.so
nm -D libheatflow.so | grep ' T '
ldd libheatflow.so | grep -o 'libgfortran[^ ]*'
0000000000000640 T __heatflow_kernels_MOD_ftcs_sweep
00000000000007c0 T __heatflow_kernels_MOD_jacobi_sweep
0000000000000920 T heatflow_step
libgfortran.so.5
libgfortran.so.5

Three things in this output are worth reading. The library depends on libgfortran.so.5 at run time (ldd matched it on two lines, the dependency and its resolved path), so the machine that loads it needs the GNU Fortran runtime installed, or the library has to ship with it. Only heatflow_step has a stable, language-neutral name; the other two symbols carry the compiler’s module-procedure mangling and are not API. And a static archive (ar rcs) has none of these problems for a C program but cannot be loaded by Python.

Source: GNU Fortran manual: Influencing the linking step, Flang documentation: Flang driver

C++: extern "C", layouts and mdspan

A C++ compiler mangles function names to encode the parameter types, so a declaration of heatflow_step without extern "C" asks the linker for a symbol that no Fortran compiler produced. The header above guards its declaration with #ifdef __cplusplus, which is what makes it usable from both languages. Without the guard, a C++ client fails like this (the header with the extern "C" lines removed):

/usr/bin/ld: /tmp/cc9FlmU5.o: in function `main':
noguard.cpp:(.text+0x78): undefined reference to `heatflow_step(int, int, double const*, double*, double, double, double, double, double*)'
collect2: error: ld returned 1 exit status

and nm on the two objects shows the difference in one line each: U heatflow_step with the guard, and U _Z13heatflow_stepiiPKdPdddddS1_ without it.

C++ can also say what the array layout is. Fortran’s u(i, j) has the first index running fastest, which is the layout that C++23 calls std::layout_left in <mdspan>. A client that wants to write u[i, j] for the Fortran element u(i + 1, j + 1) can wrap the buffer in an mdspan without copying:

#include <cstdio>
#include <mdspan>
#include <vector>

#include "heatflow.h"

int main() {
    constexpr int nx = 21;
    constexpr int ny = 11;
    std::vector<double> u(static_cast<std::size_t>(nx) * ny, 0.0);
    std::vector<double> unew(u.size(), 0.0);

    // Fortran u(i, j) is column-major: the first index runs fastest.
    // layout_left is the C++23 name for exactly that layout.
    using field_view = std::mdspan<double, std::dextents<int, 2>, std::layout_left>;
    field_view uf(u.data(), nx, ny);
    field_view unewf(unew.data(), nx, ny);
    uf[nx / 2, ny / 2] = 100.0;       // u(11, 6) in Fortran terms
    unewf[nx / 2, ny / 2] = 100.0;

    double residual = 0.0;
    const int status = heatflow_step(nx, ny, u.data(), unew.data(),
                                     1.0, 5.0e-5, 0.05, 0.05, &residual);

    std::printf("status=%d residual=%.4f center=%.4f\n", status, residual,
                unewf[nx / 2, ny / 2]);
    // A row-major (layout_right) view of the same memory swaps both extents and indices.
    std::mdspan<double, std::dextents<int, 2>, std::layout_right> row(unew.data(), ny, nx);
    std::printf("row-major view: row[%d, %d] = %.4f\n", ny / 2, nx / 2, row[ny / 2, nx / 2]);
    // The layouts differ when the field is not square.
    std::printf("stride of i: %d, stride of j: %d\n", unewf.stride(0), unewf.stride(1));
    return status;
}

The grid here is 21 by 11 so that the two index orders cannot be confused. Built with g++ -std=c++23 -Wall -Wextra -pedantic -Iinclude 02_client.cpp -L. -lheatflow_capi -o cpp_client (against the archive compiled with -O2, so no -lgfortran was needed) it printed:

status=0 residual=8.0000 center=92.0000
row-major view: row[5, 10] = 92.0000
stride of i: 1, stride of j: 21

The strides are the whole story. layout_left has stride 1 for the first index and nx for the second, which is the Fortran rule. The layout_right view over the same memory addresses the element as row[j, i] with extents (ny, nx): the dimensions and the indices are both swapped, and nothing is copied. The same program, with the archive rebuilt by flang -c, printed the same three lines. std::mdspan is not a Fortran descriptor and carries no lower bounds, so the 1-based indices of Fortran are the 0-based indices of C++ plus one.

Source: cppreference: std::mdspan, cppreference: language linkage

Arrays across the boundary: column-major and row-major

A Fortran array is stored in column-major order: the elements of a(:, 1) are contiguous, then a(:, 2), and so on (array element order, 9.5.3.3). A C array of arrays is stored row-major. The memory is the same for a pointer to the first element; only the reading of it differs. The standard states the correspondence for explicit-shape arrays (18.3.6, Note 3): the Fortran array integer(c_int) :: a(18, 3:7, *) is interoperable with the C array int b[][5][18]. The dimensions are reversed, and lower bounds do not matter, only extents. In general a Fortran a(n, m) is a C a[m][n], and a(i, j) is a[j - 1][i - 1].

This Fortran function declares a matrix with the extents it is told and prints it by rows:

module layout_demo
  use, intrinsic :: iso_c_binding, only: c_int, c_double
  implicit none
  private
  public :: show_matrix
contains
  !> Print a(i, j) for an m-by-n explicit-shape dummy and return its sum.
  function show_matrix(a, m, n) bind(c, name='show_matrix') result(total)
    integer(c_int), value, intent(in) :: m, n
    real(c_double), intent(in) :: a(m, n)
    real(c_double) :: total
    integer :: i
    do i = 1, m
      print '(a, i0, a, *(f5.1))', 'a(', i, ', :) =', a(i, :)
    end do
    total = sum(a)
  end function show_matrix
end module layout_demo

Source: J3/23-007r1: 18.3.6 Interoperability of array variables and 9.5.3.3 Array element order, GNU Fortran manual: Interoperability with C

and this C program passes a 2 by 3 row-major array twice, once with the extents swapped as the rule says, and once with the C extents copied unchanged:

#include <stdio.h>

double show_matrix(const double *a, int m, int n);

int main(void) {
    /* C: 2 rows of 3 columns, row-major. */
    double c[2][3] = {{1.0, 2.0, 3.0}, {4.0, 5.0, 6.0}};

    /* Handing it over as a Fortran 3-by-2 matrix: the dimensions swap,
       and a(i, j) is c[j - 1][i - 1]. */
    double total = show_matrix(&c[0][0], 3, 2);
    printf("total=%.1f\n", total);

    /* Handing it over as a Fortran 2-by-3 matrix is a different matrix, not a transpose of the intent. */
    total = show_matrix(&c[0][0], 2, 3);
    printf("total=%.1f\n", total);
    return 0;
}
a(1, :) =  1.0  4.0
a(2, :) =  2.0  5.0
a(3, :) =  3.0  6.0
total=21.0
a(1, :) =  1.0  3.0  5.0
a(2, :) =  2.0  4.0  6.0
total=21.0

The first call shows the rule at work: row a(1, :) of the Fortran matrix is the first column of the C array, 1.0 4.0. The second call, with the extents copied, reads the same six numbers as a different matrix, and the sum, 21, is identical in both, which is why a sum or a norm cannot tell a wrong layout from a right one. A test of an interface across this boundary needs an asymmetric input, a non-square shape and an element-by-element comparison.

There are three honest ways to deal with the mismatch:

  • Keep the Fortran layout. Allocate the field as nx * ny doubles and address it as u[i + nx * j], as the C client above does, or through layout_left. No copy, and the Fortran code stays idiomatic.

  • Swap the extents and indices. A C double a[ny][nx] is the Fortran a(nx, ny), element for element, with no copy. This is the right choice for C code that must keep row-major arrays and is willing to write a[j][i] for a(i, j).

  • Transpose at the boundary. Copy once on entry and once on exit when the calling code is row-major throughout and the kernel is called rarely. For a time-stepping loop the copies cost more than the stencil.

Python has the same three choices, and NumPy names them order="F" and order="C". They are the subject of the f2py section below.

Source: J3/23-007r1: 18.3.6 Interoperability of array variables and 9.5.3.3 Array element order, GNU Fortran manual: Interoperability with C

Legacy symbols and CMake FortranCInterface

Code that predates bind(c) is called through whatever name the compiler chose, and the compilers disagree. GNU Fortran and Flang lowercase the name and append one underscore for an external procedure (scale_vec becomes scale_vec_), but a module procedure is named very differently by each: __heatflow_legacy_MOD_halve for GNU Fortran and _QMheatflow_legacyPhalve for Flang, as the nm output further down shows. Hand-written C that spells those names is tied to one compiler version. The CMake module FortranCInterface detects the convention of the configured compilers by building a probe project, and FortranCInterface_HEADER writes a header with macros that expand to the right names.

The Fortran file has one external subroutine and one module procedure, both written in the style found in older code:

! A pre-bind(c) style external routine: its C-visible name is compiler dependent.
subroutine scale_vec(n, x, factor)
  implicit none
  integer, intent(in) :: n
  double precision, intent(inout) :: x(n)
  double precision, intent(in) :: factor
  x = factor*x
end subroutine scale_vec

module heatflow_legacy
  implicit none
contains
  subroutine halve(n, x)
    integer, intent(in) :: n
    double precision, intent(inout) :: x(n)
    x = 0.5d0*x
  end subroutine halve
end module heatflow_legacy

Source: CMake: FortranCInterface, J3/23-007r1: 18.10.2 Binding labels for procedures

The CMakeLists.txt asks FortranCInterface for a header with the two symbols and verifies that the C and Fortran compilers can call each other:

cmake_minimum_required(VERSION 3.20)
project(legacy_demo LANGUAGES C Fortran)

include(FortranCInterface)
FortranCInterface_HEADER(fc_mangle.h MACRO_NAMESPACE "FC_"
  SYMBOLS scale_vec heatflow_legacy:halve)
message(STATUS "Fortran symbol mangling: global='${FortranCInterface_GLOBAL_SUFFIX}' module prefix='${FortranCInterface_MODULE_PREFIX}'")

add_library(legacy STATIC scale_vec.f90)
add_executable(legacy_main main.c)
target_include_directories(legacy_main PRIVATE ${CMAKE_CURRENT_BINARY_DIR})
target_link_libraries(legacy_main PRIVATE legacy)

# Check that C and Fortran in the chosen compilers can call each other.
FortranCInterface_VERIFY()

Source: CMake: FortranCInterface, J3/23-007r1: 18.10.2 Binding labels for procedures

The C source includes the generated header. scale_vec and heatflow_legacy_halve are macros that the header defines; every Fortran argument is a pointer, because without value Fortran passes by reference:

#include <stdio.h>
#include "fc_mangle.h"

/* Fortran passes every argument by reference. scale_vec and heatflow_legacy_halve
   are macros from the generated header that expand to the mangled names. */
void scale_vec(const int *n, double *x, const double *factor);
void heatflow_legacy_halve(const int *n, double *x);

int main(void) {
    int n = 3;
    double x[3] = {1.0, 2.0, 3.0};
    double factor = 4.0;
    scale_vec(&n, x, &factor);
    heatflow_legacy_halve(&n, x);
    printf("%.1f %.1f %.1f\n", x[0], x[1], x[2]);
    return 0;
}

Configured with GNU Fortran 16.2, the relevant lines of the configure output and the program’s output were:

-- Detecting Fortran/C Interface
-- Detecting Fortran/C Interface - Found GLOBAL and MODULE mangling
-- Fortran symbol mangling: global='_' module prefix='__'
-- Verifying Fortran/C Compiler Compatibility
-- Verifying Fortran/C Compiler Compatibility - Success
[6/6] Linking C executable legacy_main
2.0 4.0 6.0

Source: CMake: FortranCInterface, J3/23-007r1: 18.10.2 Binding labels for procedures

The result is (1, 2, 3) * 4 / 2. The generated fc_mangle.h contains, besides comments, these definitions:

#define FC_GLOBAL(name,NAME) name##_
#define FC_GLOBAL_(name,NAME) name##_
#define FC_MODULE(mod_name,name, mod_NAME,NAME) __##mod_name##_MOD_##name
#define FC_MODULE_(mod_name,name, mod_NAME,NAME) __##mod_name##_MOD_##name
#define scale_vec FC_GLOBAL_(scale_vec, SCALE_VEC)
#define heatflow_legacy_halve FC_MODULE(heatflow_legacy,halve, HEATFLOW_LEGACY,HALVE)

and nm on the archive showed what the compiler really emitted, T scale_vec_ and T __heatflow_legacy_MOD_halve. The CMake link line for the C executable was cc main.c.o -o legacy_main liblegacy.a -lgfortran -lm -lgcc_s: CMake knew that the static library was Fortran and added the runtime to a C link, which is the part that hand-written makefiles forget.

Configured again with FC=flang-23, the same project reported module prefix='_QM' and nm showed T _QMheatflow_legacyPhalve and T scale_vec_. The C source did not change and the program printed the same 2.0 4.0 6.0.

FortranCInterface is the right tool for code that cannot be changed. For code that can, a bind(c, name=...) wrapper in one new module, as heatflow_capi does, makes the macros unnecessary, and it is portable to every compiler that implements the standard. The binding label is what the standard guarantees; the underscore convention is only folklore (18.10.2). A related legacy convention is that many compilers pass a hidden length argument at the end of the argument list for every character(len=*) dummy; the Julia manual warns about it, and a bind(c) interface avoids it.

Source: CMake: FortranCInterface, J3/23-007r1: 18.10.2 Binding labels for procedures

Python with f2py

f2py is the Fortran-to-Python wrapper generator of NumPy. It reads the Fortran source (or a signature file), works out the Python-visible argument list from the declarations, generates a C extension module that converts NumPy arrays to Fortran arrays and back, and builds it. No bind(c) is needed: f2py calls the Fortran procedure through the compiler’s own convention, which it also knows how to discover.

NumPy 2.5 and the Meson backend

NumPy’s documentation states that the old numpy.distutils backend of f2py was removed in NumPy 2.5.0 and that Meson is the only backend. Python 3.12 had already removed distutils from the standard library (PEP 632). In practice this changes four things.

  • python -m numpy.f2py -c ... still builds a module, but through Meson and Ninja. Both have to be installed, as they were in the virtual environment used here (meson 1.12.1 and ninja from PyPI). The compilers are chosen with the environment variables FC and CC, which the migration page of the f2py guide lists together with FFLAGS, CFLAGS and LDFLAGS.

  • Libraries are linked with -l and -L on the command line, and system dependencies with --dep lapack, which maps to a Meson dependency('lapack').

  • --build-dir keeps the generated Meson project, which can be committed or used as a Meson subproject. This is the migration route for projects that used setup.py and numpy.distutils.core.setup: either call python -m numpy.f2py -c from the build script, or write a meson.build as below.

  • --backend is still accepted but documented as a legacy option, since only Meson remains.

The compiler used in the examples below is whatever FC selects, GNU Fortran 16.2 here.

A first wrapper, and the kind pitfall

The FTCS step as a plain subroutine, with explicit sizes so that f2py can read the shapes:

subroutine ftcs_step(u, nx, ny, alpha, dt, dx, dy, unew, residual)
  use iso_fortran_env, only: dp => real64
  implicit none
  integer, intent(in) :: nx, ny
  real(dp), intent(in) :: u(nx, ny), alpha, dt, dx, dy
  real(dp), intent(out) :: unew(nx, ny), residual
  integer :: i, j
  unew = u
  residual = 0.0_dp
  do j = 2, ny - 1
    do i = 2, nx - 1
      unew(i, j) = u(i, j) + alpha*dt*((u(i-1, j) - 2*u(i, j) + u(i+1, j))/dx**2 &
                                     + (u(i, j-1) - 2*u(i, j) + u(i, j+1))/dy**2)
      residual = max(residual, abs(unew(i, j) - u(i, j)))
    end do
  end do
end subroutine ftcs_step

and a Python script that calls it on a 21 by 21 field:

import numpy as np
import kernel

n = 21
u = np.zeros((n, n), order="F")
u[n // 2, n // 2] = 100.0
unew, residual = kernel.ftcs_step(u, 1.0, 5e-5, 0.05, 0.05)
print(f"residual={residual:.4f} center={unew[n // 2, n // 2]:.4f}")

The first build is the obvious one:

python -m numpy.f2py -c kernel.f90 -m kernel
python demo.py

It builds without a warning and then the script dies:

free(): corrupted unsorted chunks
Aborted

Source: NumPy: Using F2PY

This is the single most important f2py pitfall for modern code. f2py reads the declaration real(dp) and has to turn the kind dp into a C type, but dp is an alias for real64 imported from iso_fortran_env, which f2py cannot evaluate. It falls back to the default real, four bytes, while Fortran writes eight-byte values into the buffer: the wrapper allocates too little memory, and the process aborts or, worse, silently returns garbage. There is no diagnostic at build time. The fix is a f2cmap file named .f2py_f2cmap in the build directory (or passed with --f2cmap) that tells f2py what dp means:

{'real': {'dp': 'double'}}

With that file next to kernel.f90, the same two commands work:

residual=8.0000 center=92.0000

Source: NumPy: Using F2PY

The same value as the C client printed, from the same kernel. The alternatives to a map file are to write the kinds that f2py already knows, real(c_double) and integer(c_int) from iso_c_binding, as the second example below does, or to declare real(kind=8) for code that is only ever built with compilers that use byte sizes. A kinds module with a parameter is invisible to f2py unless the map says otherwise.

What f2py understood about the routine is visible in the docstring of the generated function:

unew,residual = ftcs_step(u,alpha,dt,dx,dy,[nx,ny])

Wrapper for ``ftcs_step``.

Parameters
----------
u : input rank-2 array('d') with bounds (nx,ny)
alpha : input float
dt : input float
dx : input float
dy : input float

Other Parameters
----------------
nx : input int, optional
    Default: shape(u, 0)
ny : input int, optional
    Default: shape(u, 1)

Returns
-------
unew : rank-2 array('d') with bounds (nx,ny)
residual : float

Four rules visible in that signature apply to every wrapped routine. intent(out) arguments are not parameters; they become return values, in declaration order. The dimension arguments nx and ny appear in the dummy-argument list of the Fortran routine but are optional in Python, with defaults derived from the array shape, because the array declaration u(nx, ny) ties them together. Scalars are passed by value from Python and by reference to Fortran without the caller noticing. And an intent(in) array is converted with whatever copy the layout and the type require.

Signature files and f2py directives

The wrapper’s interface can be written down, inspected and edited as a signature file, a .pyf file in a Fortran-90-like syntax. It is generated by running f2py with -h. The second example has three routines that treat the array argument in three ways:

! f2py example: three ways to hand an array to Fortran.
subroutine scale_inplace(u, nx, ny, factor)
  use iso_c_binding, only: c_double, c_int
  implicit none
  integer(c_int), intent(in) :: nx, ny
  real(c_double), intent(inout) :: u(nx, ny)
  real(c_double), intent(in) :: factor
  u = factor*u
end subroutine scale_inplace

subroutine scale_returned(u, nx, ny, factor)
  use iso_c_binding, only: c_double, c_int
  implicit none
  integer(c_int), intent(in) :: nx, ny
  real(c_double), intent(inout) :: u(nx, ny)
  !f2py intent(in,out) :: u
  real(c_double), intent(in) :: factor
  u = factor*u
end subroutine scale_returned

subroutine row_sums(u, nx, ny, s, total)
  use iso_c_binding, only: c_double, c_int
  implicit none
  integer(c_int), intent(in) :: nx, ny
  real(c_double), intent(in) :: u(nx, ny)
  real(c_double), intent(out) :: s(nx), total
  s = sum(u, dim=2)
  total = sum(s)
end subroutine row_sums

!f2py intent(in,out) :: u is a comment to every Fortran compiler, and an instruction to f2py. It lets the Fortran side say intent(inout) (the truth for a Fortran caller) while the Python side sees an argument that is passed in and returned. The comment directive is the usual way to adjust the Python interface without editing the generated signature file. Generating the signature file:

python -m numpy.f2py fieldops.f90 -m fieldops -h fieldops.pyf --overwrite-signature

gave (the generated header and trailer comments are omitted; the list is otherwise the file):

python module fieldops ! in
    interface  ! in :fieldops
        subroutine scale_inplace(u,nx,ny,factor) ! in :fieldops:fieldops.f90
            use iso_c_binding, only: c_double,c_int
            real(kind=c_double) dimension(nx,ny),intent(inout) :: u
            integer(kind=c_int), optional,intent(in),check(shape(u, 0) == nx),depend(u) :: nx=shape(u, 0)
            integer(kind=c_int), optional,intent(in),check(shape(u, 1) == ny),depend(u) :: ny=shape(u, 1)
            real(kind=c_double) intent(in) :: factor
        end subroutine scale_inplace
        subroutine scale_returned(u,nx,ny,factor) ! in :fieldops:fieldops.f90
            use iso_c_binding, only: c_double,c_int
            real(kind=c_double) dimension(nx,ny),intent(in,out) :: u
            integer(kind=c_int), optional,intent(in),check(shape(u, 0) == nx),depend(u) :: nx=shape(u, 0)
            integer(kind=c_int), optional,intent(in),check(shape(u, 1) == ny),depend(u) :: ny=shape(u, 1)
            real(kind=c_double) intent(in) :: factor
        end subroutine scale_returned
        subroutine row_sums(u,nx,ny,s,total) ! in :fieldops:fieldops.f90
            use iso_c_binding, only: c_double,c_int
            real(kind=c_double) dimension(nx,ny),intent(in) :: u
            integer(kind=c_int), optional,intent(in),check(shape(u, 0) == nx),depend(u) :: nx=shape(u, 0)
            integer(kind=c_int), optional,intent(in),check(shape(u, 1) == ny),depend(u) :: ny=shape(u, 1)
            real(kind=c_double) dimension(nx),intent(out),depend(nx) :: s
            real(kind=c_double) intent(out) :: total
        end subroutine row_sums
    end interface
end python module fieldops

The attributes in the signature are the vocabulary of f2py:

Attribute Meaning in Python

intent(in)

input only; the default for dummies with no intent

intent(out)

not a parameter; returned as part of the result tuple

intent(in,out)

input that is also returned; the returned object is a new Fortran-ordered array, not necessarily the argument

intent(inout)

the Fortran code modifies the caller’s array in place; the argument must already be a NumPy array of the right type, size and layout

intent(inplace)

like inout, but a copy with write-back is made when the array does not qualify; documented as preferable to inout

intent(hide)

removed from the argument list; used for values that are computed or defaulted

intent(c)

treat the array as C-contiguous instead of Fortran-contiguous; returned multidimensional arrays are then C-contiguous

optional, check(...), depend(...)

a default value, a validity test and the dependence that orders argument evaluation

A signature file is edited by hand when the defaults are wrong, and then built together with the source: python -m numpy.f2py -c fieldops.pyf fieldops.f90. If the .pyf is on the command line, f2py ignores -m.

Array order, copies and in-place updates

The Python script exercises each behaviour. a is Fortran-ordered, c is C-ordered, i is an integer array:

import numpy as np
import fieldops

a = np.arange(6.0).reshape(2, 3, order="F")           # Fortran-ordered, float64
print("flags F/C:", a.flags.f_contiguous, a.flags.c_contiguous)
fieldops.scale_inplace(a, 2.0)                        # nx, ny deduced from the shape
print("inplace F-ordered:", a.tolist())

c = np.arange(6.0).reshape(2, 3)                      # C-ordered
try:
    fieldops.scale_inplace(c, 2.0)
except Exception as e:
    print(type(e).__name__, str(e).splitlines()[0])
print("C-ordered unchanged:", c.tolist())

r = fieldops.scale_returned(c, 2.0)                   # intent(in,out): copy in, array out
print("returned:", r.tolist(), "same object:", r is c, "F-contiguous result:", r.flags.f_contiguous)
print("c after:", c.tolist())

i = np.arange(6).reshape(2, 3, order="F")             # int64 data: converted to float64 (a copy)
s, total = fieldops.row_sums(i)
print("row sums:", s.tolist(), "total:", total)
print(fieldops.row_sums.__doc__.splitlines()[0])
print(fieldops.scale_returned.__doc__.splitlines()[0])
flags F/C: True False
inplace F-ordered: [[0.0, 4.0, 8.0], [2.0, 6.0, 10.0]]
ValueError failed to initialize intent(inout) array -- input not fortran contiguous
C-ordered unchanged: [[0.0, 1.0, 2.0], [3.0, 4.0, 5.0]]
returned: [[0.0, 2.0, 4.0], [6.0, 8.0, 10.0]] same object: False F-contiguous result: True
c after: [[0.0, 1.0, 2.0], [3.0, 4.0, 5.0]]
row sums: [6.0, 9.0] total: 15.0
s,total = row_sums(u,[nx,ny])
u = scale_returned(u,factor,[nx,ny])

Source: NumPy: Using F2PY bindings in Python

Read it line by line.

  • intent(inout) modifies the caller’s memory and is strict. The Fortran-ordered a is changed in place. The C-ordered c is refused with a ValueError, because an in-place update of a C-ordered array would need a copy and the changes would be lost; c is unchanged.

  • intent(in,out) always copies when the layout is wrong, and returns a new array. r is c is False and c still holds its original values, so code that expects the argument to be updated has to use the return value.

  • The result is Fortran-ordered, because intent(c) was not given. A C-ordered u can come back as an array that indexes the same way but is laid out differently; downstream code that assumes C order should call np.ascontiguousarray or ask for intent(c).

  • Type conversion is a copy too. The integer array i is converted to float64 before the call. The computation is correct, and invisible in the output.

f2py can say how many copies it made. The macro F2PY_REPORT_ON_ARRAY_COPY makes the generated wrapper print a line to stderr for every array larger than the given size that it copies. On the Meson backend used here the macro was applied through CFLAGS (passing -DF2PY_REPORT_ON_ARRAY_COPY=1 on the f2py command line had no effect in this test):

CFLAGS=-DF2PY_REPORT_ON_ARRAY_COPY=1 python -m numpy.f2py -c fieldops.f90 -m fieldops
python demo2.py 2>&1 | grep copied
copied an array: size=6, elsize=8
copied an array: size=6, elsize=8

The two lines are the intent(in,out) copy of c and the integer-to-double conversion of i; the in-place call on a made none. For a time-stepping kernel called thousands of times, an unreported copy of a large field per step is the usual reason a wrapped kernel is slower than the Fortran program, so allocate the arrays once as np.empty(shape, order="F") and pass the same ones every time.

Building with Meson directly

python -m numpy.f2py -c is enough for one file. A project with several sources, or one that wants the same build description for CI, can use a meson.build that runs f2py on the sources and compiles the generated C with the Fortran:

project('heatflow_kernel', 'c', 'fortran', version: '0.1.0',
  meson_version: '>= 1.1.0', default_options: ['warning_level=1', 'buildtype=release'])

py = import('python').find_installation(pure: false)
py_dep = py.dependency()

incdir_numpy = run_command(py, ['-c', 'import numpy; print(numpy.get_include())'],
  check: true).stdout().strip()
incdir_f2py = run_command(py, ['-c', 'import numpy.f2py; print(numpy.f2py.get_include())'],
  check: true).stdout().strip()
inc_np = include_directories(incdir_numpy, incdir_f2py)

kernel_source = custom_target('kernelmodule.c',
  input: ['kernel.f90'],
  output: ['kernelmodule.c'],
  command: [py, '-m', 'numpy.f2py', '@INPUT@', '-m', 'kernel', '--lower', '--f2cmap', meson.current_source_dir() / '.f2py_f2cmap'])

py.extension_module('kernel',
  ['kernel.f90', kernel_source, incdir_f2py / 'fortranobject.c'],
  include_directories: [inc_np],
  dependencies: [py_dep],
  install: false)

The custom target runs f2py to produce kernelmodule.c; the extension module is then the Fortran source, the generated C and NumPy’s fortranobject.c. Note that the map file is passed explicitly with --f2cmap, since the build runs in another directory. Configuring and building:

meson setup mbuild
meson compile -C mbuild
PYTHONPATH=mbuild python demo.py
Run-time dependency python found: YES 3.13
Build targets in project: 2
[7/7] Linking target kernel.cpython-313-aarch64-linux-gnu.so
residual=8.0000 center=92.0000

The file name of the extension depends on the Python version and the platform, as with any CPython extension. The Meson build ran with Meson 1.12.1; the project requires Meson 1.1.0 or newer.

Source: NumPy: F2PY user guide and reference manual, Using F2PY, F2PY signature file, Migrating to Meson, Using via meson

ctypes and cffi against a bind(c) library

When the library already exists, or must be usable from more than one language, a Python wrapper that loads the shared library at run time is simpler than building an extension. The library is libheatflow.so from the linking section; the function is heatflow_step from the header. There is no signature to generate: the Python code states the C prototype, and a wrong statement is a crash, not an exception.

ctypes

import ctypes
import sys

import numpy as np

lib = ctypes.CDLL(sys.argv[1])
f64_2d = np.ctypeslib.ndpointer(dtype=np.float64, ndim=2, flags="F_CONTIGUOUS")

lib.heatflow_step.restype = ctypes.c_int
lib.heatflow_step.argtypes = [
    ctypes.c_int, ctypes.c_int,                  # nx, ny (value)
    f64_2d, f64_2d,                              # u, unew
    ctypes.c_double, ctypes.c_double,            # alpha, dt (value)
    ctypes.c_double, ctypes.c_double,            # dx, dy (value)
    ctypes.POINTER(ctypes.c_double),             # residual (by reference)
]


def step(u, alpha, dt, dx, dy):
    nx, ny = u.shape
    unew = u.copy(order="F")
    residual = ctypes.c_double(0.0)
    status = lib.heatflow_step(nx, ny, u, unew, alpha, dt, dx, dy, ctypes.byref(residual))
    return status, unew, residual.value


n = 21
u = np.zeros((n, n), order="F")
u[n // 2, n // 2] = 100.0
status, unew, residual = step(u, 1.0, 5e-5, 0.05, 0.05)
print(f"status={status} residual={residual:.4f} center={unew[n // 2, n // 2]:.4f}")

try:
    step(np.ascontiguousarray(u), 1.0, 5e-5, 0.05, 0.05)
except ctypes.ArgumentError as e:
    print("C-ordered array rejected:", e)

status, _, _ = step(np.zeros((2, 2), order="F"), 1.0, 5e-5, 0.05, 0.05)
print("2x2 grid status:", status)
status=0 residual=8.0000 center=92.0000
C-ordered array rejected: argument 3: TypeError: array must have flags ['F_CONTIGUOUS']
2x2 grid status: 1

Source: Python documentation: ctypes

The ndpointer with flags="F_CONTIGUOUS" is the part that does the work: it makes the layout a precondition that is checked on every call, instead of a silent transposition. The value arguments are plain c_int and c_double, and residual, the intent(out) scalar, is passed with byref. Without restype, ctypes assumes the result is an int, which is right here only by luck; without argtypes it guesses the type of every argument from the Python object and nothing checks the call. The status is returned to Python as a number, and it is Python’s job to turn 1 into an exception.

cffi

cffi in its ABI mode takes the prototype as a string, which can be pasted from the header:

import sys

import numpy as np
from cffi import FFI

ffi = FFI()
ffi.cdef("""
int heatflow_step(int nx, int ny, const double *u, double *unew,
                  double alpha, double dt, double dx, double dy,
                  double *residual);
""")
lib = ffi.dlopen(sys.argv[1])

n = 21
u = np.zeros((n, n), order="F")
u[n // 2, n // 2] = 100.0
unew = u.copy(order="F")
residual = ffi.new("double *")
# from_buffer needs a C-contiguous buffer: the transpose of an F-ordered array is C-contiguous
# and shares the same memory, so Fortran sees u(i, j) unchanged.
status = lib.heatflow_step(n, n, ffi.from_buffer("double[]", u.T), ffi.from_buffer("double[]", unew.T),
                           1.0, 5e-5, 0.05, 0.05, residual)
print(f"status={status} residual={residual[0]:.4f} center={unew[n // 2, n // 2]:.4f}")
status=0 residual=8.0000 center=92.0000

The first attempt, passing u directly to ffi.from_buffer, failed with ValueError: ndarray is not C-contiguous, because cffi accepts only C-contiguous buffers. Passing u.T, the transpose, which for an F-ordered array is a C-contiguous view of the same memory, solves it without a copy; Fortran sees the column-major array it expects. Both scripts ran against a libheatflow.so built by gfortran and, with the same Python code, against one built with flang -fPIC -shared (for ctypes).

The trade-off with f2py is clear. f2py generates checked conversions (shapes, defaults, types, copies) and returns tuples; it needs a compiler and Meson at install time and rebuilds for every Python version. ctypes and cffi need no build step on the Python side, work with every Python version, and can load a library built by someone else, but each call site states the prototype, and nothing verifies it. For a library published to other users, a C header plus cffi or ctypes, or a Python package that ships f2py wheels per platform, are the two usual options. How Python packages and their wheels are built and distributed is in Python Virtual Environments and Packaging.

Source: Python documentation: ctypes, cffi documentation, NumPy: numpy.ctypeslib

Other languages in one line each

These bridges all consume the same bind(c) shared library. They were not run here, and the lines below are community references, not verified examples.

Language Mechanism Notes from the project’s own manual

Julia

ccall((:heatflow_step, "libheatflow"), Cint, (Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}), nx, ny, u, unew, alpha, dt, dx, dy, residual)

the Julia manual says Fortran compilers mangle names by default and that Fortran inputs must be passed as pointers; a bind(c) function with value arguments avoids both. Julia arrays are column-major, so no reordering is needed

R

.Fortran("name", ...) passes everything by reference to a legacy external routine; .Call goes through C glue code

the R manual on writing extensions says that integer and double precision correspond to int and double, and recommends c_int and c_double from iso_c_binding to say so explicitly

A bind(c) function with value arguments does not fit the .Fortran interface, which passes every argument by reference; wrap it in a subroutine without value for R, or write the .Call glue.

Source: Julia manual: Calling C and Fortran Code, Writing R Extensions

The Python and Fortran Rosetta Stone

The Fortran-lang site keeps a side-by-side page of Python and NumPy idioms next to their Fortran equivalents. It is the quickest way to find the array syntax a kernel needs when the prototype was written in NumPy. It is a mapping of syntax, not a calling convention, and it complements f2py rather than replacing it. The comparisons with C are in Fortran for Programmers of Other Languages. For the C and C++ languages themselves, this site has C: Numbers and Math and C++: Numbers and Math, and the Python side is in Python Numbers. The Rust equivalent of this page, calling C from Rust and exposing Rust to C, is in Rust: FFI and C Interop. The NumPy-based data science material that sits on top of f2py kernels is planned as a separate section (#244).

Source: Fortran Programming Language: Python Fortran Rosetta Stone

Where the books are outdated

Area Outdated in the books, current approach

Calling convention

older material explains foreign calls with %VAL, %REF and compiler-specific upper or lower case and underscore rules; today bind(c) with value and name= specifies the symbol and the passing mode, and FortranCInterface handles the code that cannot be changed

Python builds

f2py with numpy.distutils and setup.py; NumPy 2.5 removed that backend and Meson is the only one, with FC and CC selecting the compilers

Fortran 77 style wrappers

wrappers that pass a hidden length for every character(len=*) argument and fixed-size arrays by address; explicit-shape arrays with explicit sizes still work, and assumed-shape arguments are now reachable from C through ISO_Fortran_binding.h

Arrays

older advice on avoiding the row-major versus column-major copy by transposing by hand; C++23 std::mdspan with layout_left and NumPy order="F" make the layout part of the type or the flag

Compilers

older text names compilers that have since changed: ifort was replaced by ifx and the LLVM Fortran driver is now flang

Standard version provenance and compiler support

Feature Introduced Compiler support checked on 2026-10-19

bind(c) and the iso_c_binding module

Fortran 2003 (GNU Fortran manual)

GNU Fortran 16.2 and LLVM Flang 23 build and run all examples

Interoperable procedures with assumed-shape, allocatable, pointer, optional and assumed-rank dummies, and ISO_Fortran_binding.h

Fortran 2018, from ISO/IEC TS 29113:2012 (J3/23-007r1, C.1)

see C Descriptors and ISO_Fortran_binding

f_c_string, c_f_strpointer

Fortran 2023

covered with strings in Interoperability with C

f2py with the Meson backend only

NumPy 2.5.0 (not a language feature)

f2py ran on NumPy 2.5.4 with GNU Fortran 16.2; the C and Python parts of Flang were not re-run for f2py

std::mdspan, std::layout_left

C++23

built with g++ -std=c++23 (GCC 16.2)

Not checked here, because no compiler for them is installed: Intel ifx, NAG, HPE Cray and NVIDIA HPC compilers. The runtime library names and the mangling rules above are those of GNU Fortran 16.2 and LLVM Flang 23; other compilers differ, which is the reason bind(c) exists. See History, Standards and Compilers for the per-compiler support table.

Pitfalls

  • No -lgfortran (or the Flang runtime) when the link driver is C. The link fails on _gfortran_* or _FortranAio* symbols, often only after someone adds -fcheck=all, a print or an allocate. Link with the Fortran driver, or name the runtime.

  • Libraries before objects on the link line. undefined reference to heatflow_step with the archive on the command line is almost always an order problem.

  • No extern "C" in the header. C++ asks for a mangled name. Guard every public header with #ifdef __cplusplus.

  • Arguments without value are pointers. A C caller that passes 3 instead of &n to a legacy routine crashes, not fails.

  • Layouts are silently wrong. Row-major data read as column-major produces a transposed or scrambled matrix, and sums and norms hide it. Test with a non-square, asymmetric array.

  • real(dp) in an f2py kernel without an f2cmap. Corrupts memory with no build-time warning. Use real(c_double) or a .f2py_f2cmap.

  • intent(inout) with a C-ordered array raises in f2py; intent(in,out) copies and returns a new array, so the argument is not updated.

  • Copies in a loop. An f2py call with a C-ordered or integer array copies the whole field every time. Allocate once, in Fortran order, in the right dtype.

  • stop and error stop in a library. They terminate the calling process, including a Python interpreter or a C++ application. Return a status code.

  • I/O units are the runtime’s, not the caller’s. Fortran print buffers independently of C stdout and Python’s sys.stdout, so the order of lines from the two sides may not be the order of the calls. Flush on the Fortran side before returning, or avoid printing from library code.

  • The shared library needs the runtime at run time. libgfortran.so.5 was a dependency of the library built here; check with ldd on the deployment machine.

  • Mixing compilers. gfortran objects and flang objects in one program need both runtimes and separate module directories.

  • Threads and the GIL. f2py holds the Python global interpreter lock during the call unless the extension is built with --freethreading-compatible, which the f2py documentation says to use only after inspecting the Fortran code for thread safety.

See also Best Practices and Style and Common Gotchas and Pitfalls.

References