Case Study 1: Wrapping a C Library

"Most interoperability bugs are not in the code. They are in the contract between two pieces of code that each believe they are right."

Executive Summary

A colleague hands you dsp.c, a small, fast, well-tested C library of three signal-processing routines, and asks you to call it from your Fortran analysis program. Nobody will rewrite the C — it works — so your job is purely to wrap it: to write the Fortran interface blocks that let Fortran call each routine correctly. This is the most common interoperability task in practice, and it is deceptively easy to get wrong, because a wrong interface usually compiles and then misbehaves silently. We will read the C header, classify every argument by how it is passed, write a first set of interfaces, watch two of them fail, diagnose the failures, and ship the corrected wrapper with a hand-checked driver.

Skills applied: reading a C signature and mapping each type to a C-interoperable kind (§14.1); bind(c) linkage and the name= clause (§14.2); the by-value-versus-by-reference decision — the crux of §14.3–§14.4; diagnosing the two silent bugs that a missing or misplaced value attribute cause.

Background

Here is the library. Three routines: a scalar interpolator, a whole-array reduction, and an in-place mutator.

/* dsp.c -- a small C signal-processing library to be called from Fortran. */
#include <math.h>

/* Linear interpolation between a and b at fraction t. All by value. */
double lerp(double a, double b, double t) {
    return a + t * (b - a);
}

/* Root-mean-square of x[0..n-1]. Reads the array; does not modify it. */
double rms(int n, const double *x) {
    double s = 0.0;
    for (int i = 0; i < n; i++) s += x[i] * x[i];
    return sqrt(s / n);
}

/* Normalize x in place: divide every element by the maximum. Modifies x. */
void normalize(int n, double *x) {
    double m = x[0];
    for (int i = 1; i < n; i++) if (x[i] > m) m = x[i];
    for (int i = 0; i < n; i++) x[i] /= m;
}

The corresponding header a C programmer would #include is just the three prototypes:

/* dsp.h */
double lerp(double a, double b, double t);
double rms(int n, const double *x);
void   normalize(int n, double *x);

Your Fortran program has an array of samples and needs all three. Nothing about the C changes; everything happens in how you describe it to Fortran.

Phase 1 — Read the Header, Classify Every Argument

The entire art is in one question asked of each argument: is it passed by value or by reference, and (if a pointer) is it modified? C's rules make this readable straight off the prototype.

Two tells are worth naming. First, a plain scalar (double, int) is always by value in C, so it always needs value in the Fortran interface. Second, a pointer argument is by reference; the const qualifier tells you whether the routine writes through it — const double *x in rms is read-only (intent(in)), while the bare double *x in normalize is written (intent(inout)). Reading const correctly is how you choose the right intent.

Phase 2 — Write the Interfaces (the First Attempt)

A first pass, written a little too quickly:

! dsp_wrap.f90 -- FIRST ATTEMPT (contains two bugs).
module dsp
  use, intrinsic :: iso_c_binding, only: c_int, c_double
  implicit none
  interface

    function lerp(a, b, t) bind(c, name="lerp") result(r)
      import :: c_double
      real(c_double) :: a, b, t          ! BUG 1: no  value
      real(c_double) :: r
    end function lerp

    function rms(n, x) bind(c, name="rms") result(r)
      import :: c_int, c_double
      integer(c_int), value :: n
      real(c_double), intent(in) :: x(n)
      real(c_double) :: r
    end function rms

    subroutine normalize(n, x) bind(c, name="normalize")
      import :: c_int, c_double
      integer(c_int), value :: n
      real(c_double), value :: x(n)       ! BUG 2: value on an array
    end subroutine normalize

  end interface
end module dsp

This module compiles without error. That is exactly the danger: an interoperability bug lives in the agreement between languages, and neither compiler sees both sides.

Phase 3 — Two Silent Failures, Diagnosed

Compile the module with a driver and the two bugs surface as wrong answers, not error messages.

Bug 1 — lerp without value. The C lerp receives its arguments by value, but the Fortran interface, lacking value, passes addresses. C reads those three pointers as double bit patterns and interpolates between two astronomically large or tiny numbers. The symptom is a wildly wrong result or a crash. The rule from §14.3: every scalar C passes by value needs value in Fortran. Fix: add it.

Bug 2 — value on normalize's array. Here the interface over-applies value. Marking x(n) as value tells Fortran to pass a copy of the array. C dutifully normalizes the copy — and the copy is discarded on return, so the caller's array is unchanged. The routine appears to do nothing. Worse, copying a large array by value on every call is a performance disaster (§14.3's performance note). The rule: an array that C takes as double * is passed by reference — never value. Its intent (in vs inout) follows the const qualifier. Fix: replace value with intent(inout).

A useful mnemonic falls out of these two bugs: scalars from C want value; pointers from C do not. Getting either backwards gives a program that builds cleanly and lies to you.

Phase 4 — The Corrected Wrapper and a Driver

! dsp_wrap.f90 -- CORRECTED.
module dsp
  use, intrinsic :: iso_c_binding, only: c_int, c_double
  implicit none
  interface

    function lerp(a, b, t) bind(c, name="lerp") result(r)
      import :: c_double
      real(c_double), value :: a, b, t        ! fixed
      real(c_double) :: r
    end function lerp

    function rms(n, x) bind(c, name="rms") result(r)
      import :: c_int, c_double
      integer(c_int), value :: n
      real(c_double), intent(in) :: x(n)
      real(c_double) :: r
    end function rms

    subroutine normalize(n, x) bind(c, name="normalize")
      import :: c_int, c_double
      integer(c_int), value :: n
      real(c_double), intent(inout) :: x(n)   ! fixed
    end subroutine normalize

  end interface
end module dsp

! dsp_demo.f90 -- exercise the wrapped library.
program dsp_demo
  use, intrinsic :: iso_c_binding, only: c_int, c_double
  use dsp, only: lerp, rms, normalize
  implicit none
  real(c_double) :: a(2) = [1.0_c_double, 7.0_c_double]
  real(c_double) :: b(3) = [2.0_c_double, 4.0_c_double, 8.0_c_double]

  print '(a, f6.2)', 'lerp(10,20,0.5) = ', lerp(10.0_c_double, 20.0_c_double, 0.5_c_double)
  print '(a, f6.2)', 'rms([1,7])      = ', rms(2_c_int, a)
  call normalize(3_c_int, b)
  print '(a, 3f6.2)', 'normalize(b)    = ', b
end program dsp_demo

The three results are hand-computable:

• lerp(10, 20, 0.5) = 10 + 0.5*(20-10) = 15.00.
• rms([1,7]) = sqrt((1 + 49)/2) = sqrt(25) = 5.00.
• normalize([2,4,8]): the max is 8, so the array becomes [0.25, 0.50, 1.00].

Phase 5 — Build and Verify

Compile the C library with the C compiler, the Fortran with gfortran, and link with gfortran (the C uses sqrt, so add -lm):

$ gcc -c dsp.c
$ gfortran -std=f2018 -Wall -c dsp_wrap.f90
$ gfortran -std=f2018 -Wall dsp_demo.f90 dsp_wrap.o dsp.o -o dsp_demo -lm
$ ./dsp_demo
lerp(10,20,0.5) =  15.00
rms([1,7])      =   5.00
normalize(b)    =   0.25  0.50  1.00

A final sanity check confirms the fixes did what we claimed. Deliberately revert Bug 2 (put value back on normalize's array) and the last line becomes 2.00 4.00 8.00 — the array untouched, because C normalized a throwaway copy. Revert Bug 1 and lerp prints nonsense. Seeing the failures and the fixes side by side is how the by-value/by-reference rule stops being a rule you memorize and becomes one you feel.

Discussion Questions

1. rms takes const double *x and normalize takes double *x. Explain how the presence or absence of const in the C prototype maps to the choice between intent(in) and intent(inout) on the Fortran side. What would go wrong if you declared normalize's array intent(in)?
2. Both bugs compiled cleanly. Articulate why the compiler cannot catch either, and what that implies about how you should test a freshly written interface.
3. The library uses int for the array length. On a machine with more than two billion samples, why might c_int be the wrong choice, and what interoperable kind would you use instead?

Your Turn: Extensions

• Option A. Add a fourth routine to dsp.c: double dot(int n, const double *a, const double *b);. Write its interface, add it to the driver on a=[1,2,3], b=[4,5,6], and predict the printed value before compiling ($1\cdot4 + 2\cdot5 + 3\cdot6 = 32$).
• Option B. The real dsp.h your colleague ships may declare normalize to return the maximum it divided by (double normalize(int n, double *x);). Rewrite the interface as a function and adjust the driver. What is the returned value for [2,4,8]?
• Option C. Take any small real C library you actually use (a hashing routine, a random-number generator, a units converter) and wrap one function of it from Fortran. Document, for each argument, the by-value/by-reference decision and the const-to-intent mapping before you write a line.

Key Takeaways

• Wrapping a C library is a reading exercise first: classify each argument as by-value scalar or by-reference pointer straight from the prototype, and let const choose the intent.
• The two failure modes are mirror images: a missing value on a scalar passes an address as a value (garbage); a spurious value on an array passes a copy the caller never sees (no effect, and slow).
• Interoperability bugs compile. The contract lives between the two languages, where no single compiler can see it — so a new interface is not "done" until you have run it against a hand-computed result.
• The corrected wrapper is a pure interface module: the C stays exactly as tested, and Fortran gains three fast routines at the cost of getting one attribute right per argument.