Exercises: Fortran-Python Interoperability

These exercises make you do the two-language workflow, not just read about it. Several ask you to build a real extension module with f2py and call it from Python; a few ask you to measure a speedup on your own machine. You will need NumPy installed (f2py ships with it) and a working gfortran — the same compiler you have used all book.

Difficulty: ⭐ warm-up · ⭐⭐ standard · ⭐⭐⭐ deeper. Solutions: worked solutions to the daggered (†) and odd-numbered problems are in appendices/answers-to-selected.md; the compilable ones are in code/exercise-solutions.f90. Try every problem before you look. Because benchmark numbers depend on your hardware, the timing exercises give a model result and ask for your number.

Part A — Warm-ups ⭐

15.1 † Explain, flag by flag, what f2py -c -m mymod kernel.f90 does. Which part names the thing you will import in Python?

15.2 What is an extension module? Name two extension modules you almost certainly already use without having thought of them that way.

15.3 † For a Fortran dummy argument, state how f2py presents each of these to Python: intent(in), intent(out), intent(inout).

15.4 Give the matching NumPy dtype for each: real(real64), real(real32), integer(int32), complex(real64). Why is np.array([1, 2, 3]) a trap when the Fortran expects real(dp)?

Part B — Wrap It ⭐⭐

Type, build, and run. Predict the Python signature before you build.

15.5 † Here is a Fortran function that returns the Euclidean norm of a vector:

module veclib
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
contains
  function norm2_kernel(x, n) result(r)
    integer,  intent(in) :: n
    real(dp), intent(in) :: x(n)
    real(dp) :: r
    r = sqrt(sum(x**2))
  end function norm2_kernel
end module veclib

Add the one !f2py directive that hides n, write the f2py build command naming the module veclib, and write the two-line Python driver that calls it on [3.0, 4.0]. What does it print, and why?

15.6 A subroutine has three arguments: intent(in) :: a(n), intent(out) :: total, and intent(out) :: mean, with n hidden. After wrapping, what does the Python call look like, and what does it return? (Hint: two intent(out) arguments.)

15.7 † Wrap the smooth kernel from code/example-03-smooth.f90 and write a Python driver that smooths [0, 0, 0, 12, 0, 0, 0]. Predict the output array by hand before running.

Part C — The order='F' Issue ⭐⭐

15.8 † A colleague's f2py-accelerated code is slower than their pure-NumPy version. The relevant lines are field = np.zeros((4000, 4000)) and a loop that calls step(field, ...) 20,000 times. Diagnose the problem in one sentence and give the one-keyword fix.

15.9 True or false, and explain: "Because Fortran is column-major and NumPy is row-major, f2py silently transposes your array, so a[i, j] in Python is u(j, i) in Fortran."

15.10 † You pass a default np.zeros((100, 100)) array to an f2py routine whose argument is intent(inout). What happens, and why is it different from passing the same array to an intent(in) argument?

15.11 Given a 2-D array a, write the one line that tells you whether it is F-contiguous, and the one line that returns an F-contiguous version (copying only if necessary).

Part D — Port It ⭐⭐

Translate the Python to Fortran, wrap it, and compare speed. Show your reasoning; exact times will differ.

15.12 † Port this pure-Python running-difference loop to a wrapped Fortran kernel, then benchmark both on an array of five million elements:

def diffs_py(x):
    n = len(x)
    d = np.empty(n)
    d[0] = 0.0
    for i in range(1, n):
        d[i] = x[i] - x[i-1]
    return d

Write the Fortran, the build command, and the benchmark harness (with an assert np.allclose check). Report the speedup you measure and explain where the gap comes from.

15.13 Port a 2-D four-neighbor averaging smoother (out(i,j) = average of the four edge-neighbors) to Fortran and wrap it. What must the Python driver do to the input array to make every call zero-copy?

15.14 † Back of the envelope. A kernel is called 100,000 times on a 2000 × 2000 float64 field. If the field is C-contiguous, f2py copies it in and (for the output) back on every call. Estimate the total bytes moved by copying alone. How does order='F' change that number?

Part E — Alternatives: ctypes and cffi ⭐⭐⭐

15.15 † Explain why ctypes cannot reliably call a plain Fortran subroutine that lacks bind(c). What symbol name would gfortran typically export for a module-less subroutine foo, and why does that defeat ctypes?

15.16 Given the bind(c) function in code/example-04-cfuncs.f90, write the ctypes restype and argtypes declarations and the call. Why must the length argument be passed with ctypes.byref?

15.17 † Give two concrete situations in which you would choose ctypes (or cffi) over f2py, and one in which f2py is clearly the better choice. Justify each.

Part F — Design It (extend the solver) ⭐⭐⭐

15.18 † Modify the heat kernel so a single Python call advances the field by k steps inside Fortran (loop k times internally, ping-ponging between two buffers), returning only the final field. Why does this amortize the boundary-crossing cost, and when does it matter most?

15.19 Write an intent(inout) in-place variant of step (update u directly, no u_new). What two properties must the Python caller guarantee about the array it passes, and what does f2py do if they are not met?

15.20 † Design (in prose + skeleton code) a Python driver that runs the solver and saves a PNG frame every 50 steps, so the frames can be assembled into an animation. Which part is Python's job and which is Fortran's? Tie your answer to the two-language workflow of §15.1.

Part G — Back of the Envelope & Interleaved ⭐⭐⭐

15.21 † (looks ahead to Chapter 31) Your kernel is 95% of the program's runtime. You make just the kernel 50× faster by moving it to Fortran; the other 5% (Python orchestration) is unchanged. What is the speedup of the whole program? (This is Amdahl's Law in disguise — set it up and compute it.)

15.22 Estimate the memory footprint of one C-to-F copy of a 5000 × 5000 float64 array. If you make that copy on every one of 10,000 steps, how much total copying have you done, and what is the fix?

15.23 † (interleaved Chapter 5) In a Fortran double loop over a 2-D array, why should the inner loop run over the first index? Connect your answer to why a NumPy array bound for Fortran should be created order='F'.

15.24 (interleaved Chapter 14) Which iso_c_binding kind matches a NumPy float64? Write the C prototype for a bind(c) function double dot(double *x, double *y, int *n) and the matching Fortran interface.

15.25 † (interleaved Chapter 6) You wrote intent on every dummy argument since Chapter 6 as a safety habit. Explain how that same habit now, in Chapter 15, does double duty by shaping the Python interface f2py generates.

Solutions to the daggered and odd-numbered problems are in appendices/answers-to-selected.md; the compilable Fortran solutions (15.5, 15.7, 15.12, 15.13, 15.19, 15.24) are in code/exercise-solutions.f90. The timing exercises (15.12, 15.14, 15.22) give a model result — your measured numbers depend on your machine, and finding them is the point.