Chapter 15: Fortran-Python Interoperability — f2py, ctypes, and the Best of Both Worlds

"Make the common case convenient and the hot case fast: script in Python, compute in Fortran." — a working principle of scientific computing, not any one person's words

Overview

Here is the moment this book has been promising you since Chapter 1. You have a Python program. It is pleasant to write, it reads like the mathematics, and it does everything you need — right up until the profiler tells you that ninety-five percent of its runtime is spent in one numerical loop, an honest for loop over a million grid points that Python executes one slow, interpreted step at a time. You could rewrite the whole program in a faster language and lose everything you like about Python. Or you could do what the professionals do: leave the program in Python, lift out that one hot loop, rewrite it in Fortran, and call the Fortran from Python as if it were an ordinary function. The rest of your program never knows the difference. The loop that took a minute now takes a second.

That is the entire promise of this chapter, and it is not a trick — it is the standard architecture of modern scientific software. NumPy itself is built this way: a thin, friendly Python skin over compiled C and Fortran, so that when you write numpy.linalg.solve you are calling down into exactly the kind of Fortran this book teaches you to write. In this chapter you stop being a consumer of that arrangement and become a producer of it. You will take a Fortran routine, run one command, and get back something Python can import and call. Then you will do it to the heat solver you have been building — wrap its update kernel, drive the whole simulation from a Python script, and plot the result with matplotlib — and you will have, in your own hands, the two-language workflow that runs a large fraction of computational science.

In this chapter, you will learn to:

• Split a numerical program cleanly into a Fortran compute kernel and a Python orchestration layer, and recognize which parts of a real code belong on each side.
• Wrap a Fortran procedure with f2py — the tool that ships inside NumPy — into an extension module you can import like any .py file, and control the Python-facing signature with intent and !f2py directives.
• Pass NumPy arrays to Fortran correctly: match the dtype, understand C-contiguous versus F-contiguous memory, and use order='F' to avoid the silent copy that the column-major/row-major mismatch otherwise forces on every call.
• Reach for ctypes or cffi when f2py is not the right fit, and see why both depend on the bind(c) interface you built in Chapter 14.
• Measure the speedup yourself — port a slow pure-Python loop to Fortran and benchmark honestly, seeing the 10-to-100× that makes the whole exercise worthwhile.

Learning Paths

How to read this chapter by track. - 🔬 Scientist ("my Python is too slow") — this is your chapter. Read every section; §15.5 and the Project Checkpoint are the payoff you came for. Do the benchmark on your own machine. - 📖 Standard — read straight through. f2py is not part of the Fortran standard (it is a NumPy tool), so this chapter is about a workflow, not language rules; §15.3 is where the standard's column-major guarantee (Chapter 5) meets NumPy's world. - 🔧 Legacy — §15.2 and §15.4 matter most: f2py and ctypes are how you put a Python-friendly face on an old Fortran library without touching its internals. - ⚡ HPC — skim §15.2, read §15.3 closely (the copy on every call is a real bottleneck at scale), and note the honest benchmarking discipline in §15.5 that Chapter 28 makes rigorous.

15.1 The Two-Language Workflow

Before any tools, the idea. A great deal of scientific computing is done in what people loosely call "Python," but almost none of the arithmetic is. What actually happens is a division of labor so common it deserves a name.

Definition (the two-language workflow). A program structured so that a high-level, interactive language handles everything except the numerical hot spot — reading inputs, orchestrating the run, making decisions, plotting, saving results — while a compiled language handles the small, arithmetic-heavy kernel where nearly all the time is spent. In scientific computing the two languages are, over and over, Python for orchestration and Fortran (or C) for the kernel.

The reason this works is a fact about where time goes, which you will make rigorous in Chapter 28 but can accept now: in a typical numerical program, the overwhelming majority of the runtime is spent in a tiny fraction of the code. One stencil sweep, one matrix factorization, one inner loop. That hot fraction is where a compiled language earns its 10-to-100× advantage. The other ninety-plus percent of the lines — argument parsing, file handling, setting up the problem, drawing the figure — runs a handful of times and could be written in anything; you may as well write it in the language that makes it pleasant.

🚪 Threshold Concept. The speed of your whole program is set by the language of its hot loop, not the language of its bulk. Once you internalize that, the "Python versus Fortran" argument dissolves. It was never a choice. You write the ten percent that runs a billion times in Fortran, and the ninety percent that runs once in Python, and you get Python's convenience at Fortran's speed. This is the sixth recurring theme of the book made concrete: Fortran and Python are better together.

Think about who does what in a real run of the heat solver you are building:

  PYTHON (orchestration)                     FORTRAN (the kernel)
  ---------------------------------          ------------------------------
  read the config / CLI arguments
  set up the initial temperature field  ---> step(u, u_new, alpha, dt)   <-- the hot loop,
  for n in range(nsteps):                     the five-point stencil          called nsteps
      u = step(u, ...)   ------------------>  swept over every cell            times
  every 100 steps: plot / save
  make the final figure with matplotlib
  write a summary report

The step kernel — a double loop over the grid, run once per timestep, perhaps hundreds of thousands of times over a simulation — is exactly the code that must be fast, and exactly the code Fortran is best at. Everything around it is glue, and glue is where Python shines: a few lines of matplotlib give you a figure that would be a chapter's worth of work in Fortran, and Python's ecosystem (SciPy, pandas, scikit-learn, Jupyter) is right there for the analysis.

🐍 Python Comparison — "then why not just use NumPy for the kernel too?" Often you should! If your hot loop can be written as a handful of whole-array NumPy operations — u[1:-1,1:-1] = ... on entire array slices — then NumPy is calling compiled C underneath and will be fast, and you need no Fortran at all. The trouble is that many real kernels cannot be vectorized that way. The instant your loop carries a dependency — each timestep needs the array the previous timestep produced, each iteration reads the value the last one wrote — you are forced back into an explicit Python loop, and an explicit Python loop over a million elements is where the interpreter's per-element overhead crushes you, routinely 50-to-100× slower than compiled code. The heat time-loop is precisely such a case: step n+1 needs the field from step n. NumPy can vectorize the space sweep inside one step, but not the time loop across steps — and that is the moment you reach for Fortran. We measure exactly this gap in §15.5.

So the plan for the rest of the chapter is: (§15.2) learn the mechanical act of turning a Fortran routine into something Python can call; (§15.3) get the data across the boundary correctly, which is entirely a question of memory layout; (§15.4) learn the fallbacks for when the easy path does not fit; and (§15.5) prove the payoff with a measurement. Then the Project Checkpoint makes it real on the solver.

15.2 f2py: Wrapping a Fortran Routine for Python

The tool that makes this easy is called f2py, and you very likely already have it — it ships as part of NumPy.

Definition (f2py). f2py — "Fortran to Python interface generator" — is a program, distributed with NumPy, that reads Fortran source and automatically generates and compiles a Python extension module exposing your Fortran procedures as callable Python functions. It handles the tedious middle: parsing the Fortran, generating the C glue code that bridges the Python and Fortran calling conventions, mapping types and array shapes, and invoking the compilers to build the result. You write Fortran and Python; f2py writes and compiles the layer between them.

Definition (extension module). An extension module is a compiled shared library — .so on Linux/macOS, .pyd on Windows — that Python can import exactly like a pure-Python .py file, but whose code is native machine code rather than interpreted Python. This is not exotic: NumPy, SciPy, and most of the scientific Python stack are extension modules. f2py simply lets you build your own from Fortran.

The smallest possible example

Start with a kernel small enough to check by hand: sum the squares of an array. Put it in a module (our house style — Chapter 8):

! example-01-sum-squares.f90
module fastmath
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
contains

  subroutine sum_squares(x, s, n)
    integer,  intent(in)  :: n
    real(dp), intent(in)  :: x(n)
    real(dp), intent(out) :: s
!f2py intent(hide), depend(x) :: n = shape(x, 0)
    integer :: i
    s = 0.0_dp
    do i = 1, n
      s = s + x(i) * x(i)
    end do
  end subroutine sum_squares

end module fastmath

Two things in that source are for f2py, not for Fortran. The argument s is declared intent(out): Fortran already knows this means "I produce this value," and f2py reads the same declaration to decide that s should be a Python return value, not something the caller passes in. And the comment line beginning !f2py is a directive that f2py obeys and the Fortran compiler ignores (it is just a comment). Here it says "hide the argument n from Python and compute it from the shape of x" — because Python already knows how long the array is; making the user pass the length by hand would be un-Pythonic and error-prone. Everything else is ordinary modern Fortran that compiles on its own:

$ gfortran -std=f2018 -Wall -c example-01-sum-squares.f90

That command just checks it compiles (it produces an object file and a .mod); f2py will do the real build.

The build command

The one command you must remember:

$ f2py -c -m flib example-01-sum-squares.f90

Read it as three pieces. -c means "compile and build the extension module" (as opposed to only generating a signature file). -m flib names the module flib — this is the name you will import in Python, and it is yours to choose; it need not match the filename or the Fortran module name. The final argument is the source. f2py parses the Fortran, generates the C wrapper, calls gfortran under the hood, and leaves a file like flib.cpython-311-x86_64-linux-gnu.so (the exact suffix encodes your Python version and platform) in the current directory. That file is the extension module.

⚠️ Common Pitfall — the build backend on new Python. For decades f2py built through distutils, which was removed from the standard library in Python 3.12. On recent NumPy and Python, f2py builds with Meson instead, which you must have installed (pip install meson ninja). The command is the same; if you hit a build error mentioning distutils or a missing backend, add --backend meson: f2py -c --backend meson -m flib example-01-sum-squares.f90. This is a moving target across versions — if a build fails, check which NumPy you have (python -c "import numpy; print(numpy.__version__)") and read that version's f2py notes rather than guessing.

Calling it from Python

Now the Python side. Because we put sum_squares inside a Fortran module named fastmath, f2py exposes it as a sub-namespace of the extension module: flib.fastmath.sum_squares. (Had we written a bare subroutine outside any module, it would be flib.sum_squares directly. The extra level is the price of the tidy module organization, and it is worth naming so it does not surprise you.)

# example-01-driver.py
import numpy as np
import flib                       # the extension module f2py just built

x = np.array([1.0, 2.0, 3.0, 4.0, 5.0], dtype=np.float64)
s = flib.fastmath.sum_squares(x)  # n is hidden; s comes back as the return value
print("sum of squares =", s)

# Expected output:
# sum of squares = 55.0

$ python example-01-driver.py
sum of squares = 55.0

Hand-check: $1^2 + 2^2 + 3^2 + 4^2 + 5^2 = 1 + 4 + 9 + 16 + 25 = 55$. Notice what f2py did for you. You did not pass n — it was inferred from x. You did not pass s — it came back as the function's return value, the natural Python way. The Fortran signature sum_squares(x, s, n) became the Python signature sum_squares(x) -> s. That reshaping, driven by the intent and !f2py directives, is most of what f2py does, and learning to control it is most of learning f2py.

When a wrapper does not behave the way you expected — an argument you thought was hidden is still required, or a return value is missing — do not guess. f2py writes the exact Python signature it generated into the function's docstring, so ask it:

print(flib.fastmath.sum_squares.__doc__)

sum_squares(x)

Wrapper for ``sum_squares``.

Parameters
----------
x : input rank-1 array('d') with bounds (n)

Returns
-------
s : float

That array('d') is f2py telling you it expects a float64 ('d' = C double) rank-1 array; Returns: s confirms the intent(out) became the return value; and n is absent, confirming it was hidden. Reading this docstring is the fastest way to see how your intent and !f2py directives actually landed — make it your first move whenever a call raises a TypeError about arguments.

📜 From History. f2py was written by Pearu Peterson around the turn of the millennium and folded into what became NumPy; it has been the standard bridge from Python to Fortran for over two decades. Its longevity is itself a small piece of evidence for a theme of this book: Fortran is not dead — a tool whose entire job is to call Fortran from the world's most popular scientific-scripting language would make no sense for a language nobody used.

💡 Intuition — intent is the contract, in both directions. You met intent in Chapter 6 as a safety feature: it tells the compiler who may write what. f2py reads the very same declarations to build the Python interface: intent(in) becomes a Python argument, intent(out) becomes a return value, intent(inout) becomes an argument that is modified in place. The habit the book has drilled since Chapter 6 — an intent on every dummy argument — turns out to be exactly the information f2py needs. Careful Fortran wraps itself.

🔄 Check Your Understanding. 1. In f2py -c -m mykernel foo.f90, which part names the thing you will import in Python? 2. Why did the Python caller not have to pass n (the array length)? 3. A Fortran subroutine argument is declared intent(out). How does f2py present it to Python?

Answers

(1) -m mykernel — you would write import mykernel. The .f90 filename and the Fortran module name are independent of it. (2) The !f2py intent(hide), depend(x) :: n = shape(x,0) directive told f2py to hide n and compute it from x's shape, which Python already knows. (3) As a return value — f2py removes it from the input arguments and hands it back as (part of) the function's result.

15.3 NumPy Arrays Meet Fortran Arrays: dtype and the order='F' Issue

The one-dimensional example hid the single most important detail in Fortran↔Python interop, because a 1-D array has nothing to hide. As soon as you pass a 2-D array — a grid, a matrix, an image, the temperature field of your solver — you collide with the fact that has haunted this entire book: Fortran stores arrays column-major, and NumPy stores them row-major. This section is that collision, and how to survive it.

Two ways an array must match: dtype and layout

For f2py to pass a NumPy array to Fortran without trouble, two properties of the array must line up with the Fortran declaration.

The first is the dtype — the element type. This one is straightforward: the NumPy dtype must match the Fortran kind, byte for byte.

Definition (dtype). A NumPy array's dtype is the type of its elements — float64, float32, int32, complex128, and so on. It must match the Fortran declaration's kind: real(real64) ↔ float64, real(real32) ↔ float32, integer(int32) ↔ int32, complex(real64) ↔ complex128. If they disagree, f2py must convert — which means copying the whole array on every call, and possibly losing precision.

Because the book computes in double precision (dp = real64, from Chapter 3), your NumPy arrays should be dtype=np.float64 — which is NumPy's default for floating-point, so you often get it for free. The trap is integers and single precision: a Python list of whole numbers becomes an int64 array, and a literal like np.array([1, 2, 3]) is int64, not float64. Be explicit.

The second property is memory layout, and this is the one that bites.

Definition (C-contiguous vs F-contiguous). A 2-D (or higher) array is C-contiguous (row-major) when its last index varies fastest in memory — rows stored one after another. It is F-contiguous (column-major, "F" for Fortran) when its first index varies fastest — columns stored one after another. NumPy arrays are C-contiguous by default; Fortran expects F-contiguous. You can check with a.flags['F_CONTIGUOUS'] (or np.isfortran(a)), and create an F-contiguous array with order='F' or np.asfortranarray(a). A 1-D array is both — which is why §15.2 sailed through.

This is the same column-major/row-major fact you first met in Chapter 5, §5.6 — Fortran's "first index varies fastest" — viewed from the Python side. Nothing has changed about Fortran; NumPy simply made the opposite default, and the boundary between them is where the two conventions meet.

What f2py actually does about it — and the cost

Here is the reassuring part and the expensive part, together. f2py preserves your logical indexing. It does not silently transpose your array. If you pass a NumPy array a and read a[i, j] in Python, the Fortran code sees the corresponding element as u(i+1, j+1) (plus one for Chapter 5's 1-based indexing). Element $(i,j)$ maps to element $(i,j)$. You do not have to mentally flip your grid.

But it can only keep that promise one of two ways. If the array you pass is already F-contiguous with the right dtype, f2py hands Fortran a pointer straight to your data — zero copy, instant. If the array is C-contiguous (the default!), f2py must make an F-contiguous copy of the entire array before the call, pass the copy to Fortran, and — for output arrays — copy the result back afterward. The values are correct either way. But the copy is $O(\text{array size})$ on every single call, and if you are calling the kernel hundreds of thousands of times in a time loop, that copy can cost more than the computation you came to accelerate.

⚡ Performance Note. The silent C-to-F copy is the classic way a Fortran-accelerated Python program ends up no faster than before. Picture a $1000 \times 1000$ field — eight megabytes — passed to step once per timestep for 100,000 steps. If every call copies the array in and the result out, you have added roughly $2 \times 8\,\text{MB} \times 10^{5} \approx 1.6$ terabytes of pure memory-shuffling to a program whose actual arithmetic you were trying to speed up. Allocate the field order='F' once, and the copy vanishes. This is not a micro-optimization; it is the difference between the workflow working and not.

The fix is a single keyword. When you create the array on the Python side, ask for Fortran order:

u = np.zeros((n, n), dtype=np.float64, order='F')   # F-contiguous from birth
# ... or convert an existing one:
u = np.asfortranarray(u)                            # copies once, then reuse

Do it once, up front, and reuse the same F-contiguous array for the whole run. Then every call is zero-copy.

Seeing it on the heat kernel

Let us wrap a genuine 2-D kernel: one explicit step of the heat update. This is the same physics-placeholder step you built in Chapter 6 — each interior cell moves toward the sum of its four neighbors — but written in the explicit-shape form f2py prefers, and producing a fresh array rather than updating in place:

! example-02-heat-step.f90
module heat_kernel
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
contains

  subroutine step(u, u_new, alpha, dt, n, m)
    integer,  intent(in)  :: n, m
    real(dp), intent(in)  :: u(n, m)
    real(dp), intent(out) :: u_new(n, m)
    real(dp), intent(in)  :: alpha, dt
!f2py intent(hide), depend(u) :: n = shape(u, 0), m = shape(u, 1)
    integer  :: i, j
    real(dp) :: lap

    u_new = u                          ! copy boundaries; interior overwritten below
    do j = 2, m - 1
      do i = 2, n - 1                  ! inner loop over first index — column-major (Ch. 5)
        lap = u(i-1, j) + u(i+1, j) + u(i, j-1) + u(i, j+1) - 4.0_dp * u(i, j)
        u_new(i, j) = u(i, j) + alpha * dt * lap
      end do
    end do
  end subroutine step

end module heat_kernel

Build it, naming the extension module heatlib:

$ f2py -c -m heatlib example-02-heat-step.f90

Now the Python driver, doing the layout correctly:

# example-02-driver.py
import numpy as np
import heatlib

n = 4
u = np.zeros((n, n), dtype=np.float64, order='F')   # F-contiguous, float64
u[0, :] = 100.0                                      # hold the top edge hot

print("F-contiguous?", np.isfortran(u))              # -> True: no copy on the call

u_new = heatlib.heat_kernel.step(u, alpha=1.0, dt=0.1)
print(u_new)

# Expected output:
# F-contiguous? True
# [[100. 100. 100. 100.]
#  [  0.  10.  10.   0.]
#  [  0.   0.   0.   0.]
#  [  0.   0.   0.   0.]]

Hand-check the one interesting row. With $\alpha\,\Delta t = 1.0 \times 0.1 = 0.1$ and the top row held at 100, interior cell u[1,1] (Python 0-based; Fortran u(2,2)) has neighbors up $=100$, down $=0$, left $=0$, right $=0$, center $=0$, so $\text{lap} = 100 + 0 + 0 + 0 - 4\cdot 0 = 100$ and the new value is $0 + 0.1 \times 100 = 10$. Cell u[1,2] is identical by symmetry. Every other interior cell has all-zero neighbors and stays $0$; the boundary cells were copied straight through. These are the same numbers the native Fortran driver produced in Chapter 6 — because it is the same kernel. You have run Fortran from Python and gotten Fortran's answer.

⚠️ Common Pitfall — intent(inout) demands F-contiguity. The kernel above takes u as intent(in) and returns a new array (intent(out)), which is forgiving: f2py may copy the input if it must. But if you write a kernel that updates its array in place — intent(inout) :: u — f2py cannot fake it with a copy, because your in-place changes would be made to the copy and thrown away. So f2py refuses: pass a C-contiguous (or wrong-dtype) array to an intent(inout) argument and it raises a ValueError about the array not being Fortran-contiguous, rather than silently doing the wrong thing. The cure is the same order='F'. In-place update is the faster pattern (no per-call allocation) and the one you will want in a tight time loop — but it makes the F-contiguity requirement a hard error instead of a silent cost, which is arguably a feature.

🐛 Find the Bug. A colleague reports their f2py-wrapped solver is "somehow slower than the pure-Python version." Their setup line is u = np.zeros((2000, 2000)) and they call step 50,000 times in a loop. What is wrong, and what is the one-word-per-array fix?

Answer

np.zeros((2000, 2000)) is C-contiguous (NumPy's default) — so every one of the 50,000 calls forces f2py to copy the whole 32 MB array into Fortran order and copy the result back. The copying dwarfs the arithmetic. Fix: create it F-contiguous, np.zeros((2000, 2000), order='F'), once, and reuse it. (Also confirm it is float64; np.zeros defaults to float64, so that part is fine here.)

🔗 Connection. This is the Python-facing shadow of the most important performance idea in the book. Chapter 5, §5.6 taught you to loop down columns because Fortran is column-major; here the same column-major layout is what your NumPy array must adopt (order='F') to cross the boundary for free. And it is the mirror image of the row/column-major gotcha you handled going the other way, from Fortran into C, in Chapter 14, §14.5. One layout fact, three appearances — respect the memory order and it repays you everywhere.

15.4 When f2py Isn't Enough: ctypes and cffi

f2py is the right first choice for numerical Fortran, because it understands arrays and dtypes and builds the module for you. But it is not the only door, and sometimes not the right one. Two alternatives from the Python standard library and ecosystem are worth knowing.

ctypes is a foreign-function interface built into Python's standard library. Instead of generating and compiling a wrapper from source, you compile your Fortran into a plain shared library yourself and load it at runtime:

$ gfortran -std=f2018 -Wall -shared -fPIC example-04-cfuncs.f90 -o libheat.so

import ctypes
import numpy as np

lib = ctypes.CDLL("./libheat.so")           # load the shared library
# You must declare the C-level argument and return types yourself:
lib.c_sum_squares.restype = ctypes.c_double
lib.c_sum_squares.argtypes = [
    np.ctypeslib.ndpointer(dtype=np.float64, flags="F_CONTIGUOUS"),
    ctypes.POINTER(ctypes.c_int),
]

Notice what ctypes does not do for you: it does not read your Fortran, so it cannot infer types, hide dimension arguments, or manage array layout. You declare every type by hand, you pass array lengths yourself, and — the crucial part — you must give the routine a predictable C name and calling convention, which is exactly what bind(c) from Chapter 14 provides:

! example-04-cfuncs.f90  — a bind(c) routine ctypes can find and call
function c_sum_squares(x, n) result(s) bind(c, name="c_sum_squares")
  use, intrinsic :: iso_c_binding, only: c_double, c_int
  implicit none
  integer(c_int), intent(in) :: n
  real(c_double), intent(in) :: x(n)
  real(c_double) :: s
  integer :: i
  s = 0.0_c_double
  do i = 1, n
    s = s + x(i) * x(i)
  end do
end function c_sum_squares

Without bind(c), gfortran would mangle the exported symbol — typically lowercasing it and appending an underscore, so c_sum_squares becomes c_sum_squares_ — and, worse, the details of how the value is returned and how the hidden array length is passed would follow Fortran's private conventions rather than the C ABI ctypes speaks. bind(c, name="c_sum_squares") pins the symbol name and forces the C calling convention, so lib.c_sum_squares resolves and behaves. This is the payoff of Chapter 14: ctypes and cffi can only call Fortran that has been given a C face.

💡 Intuition — the C ABI is the universal adapter. Fortran has no standard binary interface of its own, but essentially every language can call C. So the portable way for language X to call Fortran is: Fortran wears a C mask (bind(c)), and X calls it as if it were C. f2py builds that mask and the Python side automatically; ctypes and cffi make you supply the mask (via bind(c)) and describe the C signature by hand. Same strategy, different amount of hand-work.

cffi — the C Foreign Function Interface — is a third-party package (widely used, originally from the PyPy project) with the same goal as ctypes but an interface many prefer: you paste the C declarations (double c_sum_squares(double *x, int *n);) and cffi compiles and manages the binding. Like ctypes, it targets the C ABI, so it too calls Fortran through bind(c).

When should you use which? A practical rule:

⚠️ Common Pitfall — passing by reference. Fortran passes arguments by reference (it hands over the address, not the value). C, and therefore ctypes, passes scalars by value unless you say otherwise. That is why the n in the ctypes example above is typed ctypes.POINTER(ctypes.c_int) and you pass ctypes.byref(...), not a bare Python int. Forget the pointer and you pass Fortran a garbage address for the array length — usually a crash, sometimes worse. f2py hides all of this; ctypes and cffi make it your responsibility. It is one more reason to prefer f2py for numerical work and keep ctypes for the cases that need it.

🔄 Check Your Understanding. 1. Why can neither ctypes nor cffi call a Fortran subroutine that lacks bind(c) (reliably)? 2. Give one situation where ctypes is a better choice than f2py.

Answers

(1) Without bind(c), the compiler mangles the symbol name (e.g., appends an underscore) and uses Fortran's private calling conventions; ctypes/cffi speak the C ABI and look up a C symbol name, so they cannot reliably find or correctly call the routine. bind(c) gives it a stable C name and the C calling convention. (2) When you have a prebuilt shared library you cannot rebuild from source, or you want to avoid a NumPy-based build step entirely — ctypes just loads the .so/.dll at runtime.

15.5 The Payoff: Port, Wrap, and Measure the Speedup

Everything so far has been machinery. This section is the reason the machinery exists — and it is the anchor this book has foreshadowed since Chapter 1: take a numerical function that is slow in pure Python, rewrite it in Fortran, call it from Python, and measure how much faster it got. You will do this yourself, on your own machine, because a speedup you measured is worth a hundred you were told about.

The honest-benchmarking caveat, up front

This book never runs code while it is being written, and speedups are the most abused numbers in computing, so read this carefully: the numbers below are illustrative orders of magnitude, not measurements we made. The real number depends on your CPU, your compiler and its flags, the array size, your NumPy version, and a dozen other things. What is robust — reported for decades, and easy to reproduce — is the order of magnitude: a pure-Python element-by-element loop is typically 10 to 100× slower than the equivalent compiled Fortran, and can be far worse for large sizes. Your job in this section is not to trust "50×"; it is to run the benchmark and find your number. Chapter 28 makes benchmarking rigorous — warm-up runs, repetition, variance, wall-versus-CPU time; here we do the honest-but-simple version.

The function to port

A clean, self-contained example: a one-dimensional three-point smoother — replace each interior element with the average of itself and its two neighbors. Simple enough to verify by hand, and a stand-in for the kind of neighbor-averaging that appears everywhere in numerical code.

! example-03-smooth.f90
module smoother
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
contains

  subroutine smooth(x, y, n)
    integer,  intent(in)  :: n
    real(dp), intent(in)  :: x(n)
    real(dp), intent(out) :: y(n)
!f2py intent(hide), depend(x) :: n = shape(x, 0)
    integer :: i
    y(1) = x(1)                          ! endpoints copied through
    y(n) = x(n)
    do i = 2, n - 1
      y(i) = (x(i-1) + x(i) + x(i+1)) / 3.0_dp
    end do
  end subroutine smooth

end module smoother

$ f2py -c -m smoothlib example-03-smooth.f90

Verify it on a tiny input before trusting it on a big one:

import numpy as np
import smoothlib

x = np.array([0.0, 0.0, 9.0, 0.0, 0.0])
print(smoothlib.smoother.smooth(x))

# Expected output:
# [0. 3. 3. 3. 0.]

Hand-check: endpoints copy, so y[0]=0, y[4]=0. Interior: y[1]=(0+0+9)/3=3, y[2]=(0+9+0)/3=3, y[3]=(9+0+0)/3=3. The spike spread into a small plateau — the smoother works.

The pure-Python version and the benchmark

Now the same computation written the way you would before you knew any better — an explicit Python loop — and a harness to time both:

# benchmark.py
import time
import numpy as np
import smoothlib

def smooth_py(x):
    """Pure-Python three-point smoother — the slow baseline."""
    n = len(x)
    y = np.empty(n)
    y[0], y[-1] = x[0], x[-1]
    for i in range(1, n - 1):            # one interpreted step per element
        y[i] = (x[i-1] + x[i] + x[i+1]) / 3.0
    return y

n = 5_000_000
x = np.random.rand(n)                    # 1-D: C- and F-contiguous both, no order issue

t0 = time.perf_counter()
y_py = smooth_py(x)
t1 = time.perf_counter()
y_f  = smoothlib.smoother.smooth(x)
t2 = time.perf_counter()

# Same answer? (guards against a wrapping mistake)
assert np.allclose(y_py, y_f)

py_time = t1 - t0
f_time  = t2 - t1
print(f"pure Python : {py_time:8.4f} s")
print(f"Fortran     : {f_time:8.4f} s")
print(f"speedup     : {py_time / f_time:6.1f}x")

# Expected output (ILLUSTRATIVE — your numbers will differ, possibly a lot):
# pure Python :   ~1.5    s
# Fortran     :   ~0.01   s
# speedup     :   ~100x

The structure is the honest part; treat the numbers as placeholders you will overwrite with your own. Two disciplines are already built in and worth keeping every time you benchmark: the assert np.allclose(...) that proves the fast version computes the same thing as the slow one (a speedup that gets the wrong answer is not a speedup), and timing each version separately with time.perf_counter(), the right clock for wall-time measurement.

⚡ Performance Note — where the gap comes from. The Fortran loop and the Python loop do the same arithmetic, the same $3n$ additions and $n$ divisions. The difference is per-iteration overhead. Each pass of the Python for loop dispatches bytecode, checks types, boxes and unboxes floats, and indexes Python objects — tens of machine instructions of interpreter bookkeeping wrapped around one useful division. The Fortran loop, compiled, is a handful of machine instructions per element with no bookkeeping, and the compiler may vectorize several elements at once. Multiply that per-element overhead by five million elements and the "10-to-100×" stops being folklore and becomes arithmetic.

🐍 Python Comparison — "but NumPy could vectorize this smoother." True, and important to say honestly: this particular one-pass smoother can be written as one NumPy slice expression, y[1:-1] = (x[:-2] + x[1:-1] + x[2:]) / 3, which is fast because it runs in compiled C — and if that is all you need, use it and skip Fortran. The Fortran wins decisively in the cases NumPy cannot vectorize: when the loop carries a dependency the array form cannot express. The heat solver is exactly that case — step n+1 needs the field step n produced, so you cannot collapse the time loop into one array expression; you must iterate, and iterating in Python is slow. That is why the Project Checkpoint, not this smoother, is the real demonstration — and why "port the hot loop to Fortran" remains the answer even in a NumPy-fluent shop.

🧩 Try It Yourself. Before you read the checkpoint, actually build smoothlib and run benchmark.py. Write down the three numbers you get. Then change n from five million to fifty thousand and run again. Does the speedup shrink? (It should — at small n, the fixed cost of the Python-to-Fortran call and the array handling is a larger share of the total, so the advantage of the fast inner loop is diluted. The lesson: the two-language workflow pays off precisely when the kernel does enough work per call to dwarf the crossing cost.)

🔄 Check Your Understanding. 1. Why does this chapter refuse to print a specific measured speedup and ask you to measure your own? 2. What does the assert np.allclose(y_py, y_f) line protect you from?

Answers

(1) Because no code in the book is executed while it is written, and because a real speedup depends on hardware, compiler, flags, and size — a single number would be dishonest and probably wrong on your machine. The robust, reproducible claim is the order of magnitude (10-to-100×). (2) From "accelerating" your way to the wrong answer: it verifies the fast Fortran path computes the same result as the trusted slow Python one, catching a wrapping bug (transposed indices, a dtype mishap, an off-by-one at the boundary) before you rely on the speed.

Project Checkpoint

This is the one the whole book has been pointing at. You are going to wrap your heat solver's step kernel with f2py, drive the entire simulation loop from Python, and plot the evolving field with matplotlib — the moment Fortran justifies itself to a Python user, in your own code, on your own machine.

The kernel. Save the explicit-shape heat_kernel module from §15.3 as heat-solver/heat_kernel.f90. It is the f2py-friendly face of the step you have carried since Chapter 6: the same five-point update, but with the grid dimensions passed explicitly (and hidden from Python) instead of assumed-shape, and producing a new field rather than mutating in place. Build it:

$ f2py -c -m heatlib heat_kernel.f90

The driver. The simulation now lives in Python — setup, the time loop, and the figure — with the hot kernel called across the boundary each step:

# project-checkpoint-driver.py — drive the Fortran heat kernel from Python
import numpy as np
import matplotlib.pyplot as plt
import heatlib

n, nsteps, alpha, dt = 50, 500, 1.0, 0.1
u = np.zeros((n, n), dtype=np.float64, order='F')   # F-contiguous: zero-copy calls
u[0, :] = 100.0                                      # hot top edge; other edges cold

for step_no in range(nsteps):
    u = heatlib.heat_kernel.step(u, alpha, dt)       # the hot loop, in Fortran
    # u is returned F-contiguous, so the next call is zero-copy too.

plt.imshow(u.T, origin='lower', cmap='inferno')      # .T for the usual x-right/y-up view
plt.colorbar(label='temperature')
plt.title(f'heat field after {nsteps} steps')
plt.savefig('heat_field.png', dpi=150)
print('wrote heat_field.png; final max temperature =', u.max())

The column-major NOTE — read this, it is the point of the checkpoint. The single line that makes this work correctly and quickly is order='F'. The field u is created F-contiguous and every call to step returns an F-contiguous array, so no copy ever happens at the boundary — the 8-byte doubles NumPy holds are handed to Fortran exactly as they lie in memory, in the column-major order Fortran has assumed since 1957 (Chapter 5). Drop order='F' and the program still gives the right picture — f2py preserves your indexing — but it silently copies the whole field in and out on all 500 steps, and at production sizes that copy is the bottleneck. The habit to build: the array that crosses into Fortran is born order='F' and dtype=np.float64, once, and is reused.

A note on the dp kind and f2py. Our kernel declares real(dp) with dp => real64. Recent f2py resolves real64 cleanly, but a custom kind name like dp can confuse older versions, which do not know what dp maps to. The bulletproof fix is a one-line file named .f2py_f2cmap in the build directory:

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

It tells f2py "Fortran real(dp) is a C double" — i.e., NumPy float64. (Equivalently, declare the boundary arguments real(real64) directly, or compile the project's real kinds.f90 alongside the kernel: f2py -c -m heatlib kinds.f90 heat_kernel.f90.) This is the kind of small, real friction that separates a tutorial from production; now you know the fix.

How it feeds the capstone. You now have the architecture of a real scientific code: a fast Fortran kernel with a stable interface, driven and visualized from Python. The physics inside step is still the Chapter 6 placeholder — the true finite-difference stencil, the CFL-stable timestep, and the boundary conditions are the heart of Chapter 24, and because the interface of step will not change when its body becomes real, this driver will keep working. The Chapter 38 capstone presents the finished solver; this checkpoint is where it learned to speak Python.

Summary

This chapter turned Fortran into something Python can call, and made the two-language workflow concrete on your own solver.

The two things to remember: first, the boundary is a layout question — match the dtype and create your arrays order='F', and data crosses for free; forget it and you either copy on every call or (for intent(inout)) get an error. Second, you write the hot ten percent in Fortran and the other ninety in Python, and the payoff — the 10-to-100× on the loop that dominates your runtime — is the whole reason the workflow exists.

Spaced Review

Retrieval practice on the two chapters this one leans on hardest: arrays and memory layout (Chapter 5) and C interoperability (Chapter 14).

1. (Ch. 5) In Fortran's column-major order, which index of a 2-D array varies fastest through memory, and what is the corresponding NumPy layout called?
   

   Answer

   The **first** index varies fastest (all of column 1 contiguously, then column 2, …). That layout is **F-contiguous** in NumPy; NumPy's *default*, the opposite, is **C-contiguous** (row-major, last index fastest).

2. (Ch. 5) You have a NumPy array a of shape (1000, 1000) created with plain np.zeros((1000,1000)), and you pass it to an f2py-wrapped Fortran routine 10,000 times in a loop. Why might this be slow, and what is the fix?
   

   Answer

   `np.zeros` produces a **C-contiguous** array, so f2py copies all 1,000,000 elements into Fortran (column-major) order on every one of the 10,000 calls. Create it `np.zeros((1000,1000), order='F')` once and reuse it — then the calls are zero-copy.

3. (Ch. 14) What does bind(c) do to a Fortran procedure, and why do ctypes and cffi require it?
   

   Answer

   `bind(c)` gives the procedure the **C calling convention** and a **stable, unmangled symbol name** (optionally set with `name="..."`). Without it the compiler mangles the name (e.g., appends an underscore) and uses Fortran's private conventions. ctypes and cffi look up a C symbol and call through the C ABI, so they can only reliably reach a procedure that has been given a C face with `bind(c)`.

4. (Ch. 14) Which iso_c_binding kind matches a NumPy float64 element, and which Fortran kind does the book normally use for the same 8-byte double?
   

   Answer

   `real(c_double)` from `iso_c_binding` matches NumPy `float64`. The book normally writes the same 8-byte IEEE double as `real(dp)` with `dp = selected_real_kind(15, 307)`, which equals `real64` from `iso_fortran_env` — all three name the same type.

5. (Ch. 5 + 14, synthesis) True or false: because Fortran is column-major and C is row-major, f2py transposes your array's indices, so a[i, j] in Python becomes u(j, i) in Fortran.
   

   Answer

   **False.** f2py **preserves logical indexing**: `a[i, j]` maps to `u(i+1, j+1)` (the +1 is 1-based indexing, not a transpose). It reconciles the layouts by *copying* when the array is not already F-contiguous — never by silently transposing your indices.

What's Next

You have now connected Fortran outward in both directions this part covers — to C in Chapter 14, to Python here — and in doing so you kept reaching for tools (NumPy, matplotlib, and f2py itself) that live in an ecosystem larger than the language. Chapter 16 maps that ecosystem: LAPACK and BLAS (the Fortran libraries under NumPy you have been calling all along, met properly at last), FFTW, NetCDF and HDF5, and the modern toolchain — the Fortran Package Manager fpm, the standard library, documentation and testing tools — that make Fortran genuinely pleasant to build and ship in 2026. You have made Fortran talk to Python; next you will see the neighborhood Fortran lives in, and put your solver into a real project structure with fpm.