Case Study 1: Modernizing a COMMON-Block Rod-Conduction Kernel

"First make the change easy (warning: this may be hard), then make the easy change." — Kent Beck

Executive Summary

The PLATE migration in the chapter was a set piece; this is the everyday reality. You inherit RODEQ, a small but thoroughly legacy FORTRAN 77 program that solves one-dimensional steady-state heat conduction along a rod by the same Jacobi relaxation — fixed-form, IMPLICIT DOUBLE PRECISION, a COMMON /ROD/ block, a DATA statement, a statement function, and a GOTO loop. Your job is not to admire it but to modernize it, and to prove the modern version computes the identical answer. We work the eight-step recipe on real code, capture a regression reference first, and verify numerical equivalence against a profile we can check by hand. The rod is simpler than the plate on purpose: it lets you see the whole method end to end without the 2-D bookkeeping getting in the way.

Skills applied: reading legacy code and identifying its constructs (§18.1, and Chapter 17); the eight-step recipe (§18.1); incremental modernization and characterization testing (§18.2); regression testing and numerical equivalence (§18.3); hand-verification of a relaxation result (§18.4).

Background

Here is the code exactly as you found it. Read it first as archaeology — name every legacy trait before you change anything.

      PROGRAM RODEQ
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MX = 101)
      COMMON /ROD/ U(MX), M
      M = 5
      CALL INIT
      CALL SOLVE
      CALL SHOW
      STOP
      END
C
      SUBROUTINE INIT
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MX = 101)
      COMMON /ROD/ U(MX), M
      DATA UL, UR / 0.0D0, 100.0D0 /
      DO 10 I = 1, M
        U(I) = 0.0D0
   10 CONTINUE
      U(1) = UL
      U(M) = UR
      RETURN
      END
C
      SUBROUTINE SOLVE
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MX = 101)
      COMMON /ROD/ U(MX), M
      DOUBLE PRECISION UN(MX)
      DATA TOL, MAXIT / 1.0D-6, 10000 /
      HALF(A,B) = 0.5D0 * (A + B)
      IT = 0
   20 IT = IT + 1
      DM = 0.0D0
      DO 30 I = 2, M-1
        UN(I) = HALF(U(I-1), U(I+1))
        D = ABS(UN(I) - U(I))
        IF (D .GT. DM) DM = D
   30 CONTINUE
      DO 40 I = 2, M-1
        U(I) = UN(I)
   40 CONTINUE
      IF (DM .GT. TOL .AND. IT .LT. MAXIT) GO TO 20
      WRITE (6,'(A,I5)') ' sweeps = ', IT
      RETURN
      END
C
      SUBROUTINE SHOW
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MX = 101)
      COMMON /ROD/ U(MX), M
      WRITE (6,'(5F8.3)') (U(I), I = 1, M)
      RETURN
      END

The traits: fixed columns; implicit double-precision typing (so U, UN, TOL, DM, D, HALF are double, while M, MX, I, IT, MAXIT are integers); a COMMON /ROD/ block repeated in all four units; a DATA statement for the end temperatures and another for the tolerance; a statement function HALF; and a GOTO-based convergence loop. Physically it is a rod with its left end held at 0° and its right end at 100°, relaxed to steady state.

Phase 1 — Capture the Reference (Before Touching Anything)

The first rule of a safe migration: build and run the trusted code once, and freeze its output.

$ gfortran -std=legacy rodeq.f -o rodeq && ./rodeq > reference.txt
$ cat reference.txt
 sweeps =    40
    0.000  25.000  50.000  75.000 100.000

Now verify that reference by hand, so you trust it independently of the program. At steady state each interior point is the average of its two neighbors, $u_i = \tfrac12(u_{i-1} + u_{i+1})$ — the discrete statement that the second derivative is zero, whose solution is a straight line. With ends at 0° and 100° across five equally spaced points, the line is $u_i = 25(i-1)$: 0, 25, 50, 75, 100. That is exactly the printed profile, and the values are exact (25, 50, 75 are representable exactly), so we will be able to demand bit-for-bit agreement from the modern version.

The reasoning that matters: the exact sweep count (40 here, at tol $10^{-6}$) depends on the relaxation's convergence rate and is not something we need to hand-derive — what we hand-derive is the answer the sweeps march toward. The regression test will confirm the modern code reproduces both the profile and the count; the hand-check confirms the profile is physically right in the first place.

Phase 2 — The Mechanical Steps (1, 2, 6)

Apply the safe, local, form-only changes first — free-form, implicit none, and GOTO → do/exit — one routine at a time, rebuilding and re-checking after each. These touch no arithmetic, so the reference must not move. The heart of SOLVE becomes:

it = 0
do
  it = it + 1
  dm = 0.0_dp
  do i = 2, m-1
    un(i) = 0.5_dp * (u(i-1) + u(i+1))     ! HALF inlined (used once)
    dm = max(dm, abs(un(i) - u(i)))
  end do
  u(2:m-1) = un(2:m-1)
  if (dm <= tol .or. it >= maxit) exit      ! negation of the legacy GO TO test
end do

Two things to notice. The statement function HALF was used exactly once, so we inlined it rather than making an internal procedure — the simplest correct choice. And the exit condition is the exact logical negation of DM .GT. TOL .AND. IT .LT. MAXIT; getting that negation right is precisely what the regression test guards.

Phase 3 — The Structural Steps (3, 4, 5)

Now the larger moves. The COMMON /ROD/ block becomes a module — and then, better still, the rod is passed as an assumed-shape intent(inout) argument so the solver has no global state at all:

module rod
  implicit none
  integer, parameter :: dp = selected_real_kind(15, 307)
contains
  subroutine solve(u, tol, maxit, sweeps)
    real(dp), intent(inout) :: u(:)          ! assumed-shape: the rod knows its own length
    real(dp), intent(in)    :: tol
    integer,  intent(in)    :: maxit
    integer,  intent(out)   :: sweeps
    real(dp) :: un(size(u)), dm
    integer  :: i, m
    m = size(u)
    un = u
    sweeps = 0
    do
      sweeps = sweeps + 1
      dm = 0.0_dp
      do i = 2, m-1
        un(i) = 0.5_dp * (u(i-1) + u(i+1))
        dm = max(dm, abs(un(i) - u(i)))
      end do
      u = un
      if (dm <= tol .or. sweeps >= maxit) exit
    end do
  end subroutine solve
end module rod

The DATA-initialized end temperatures move into the driver as ordinary assignments, and the rod is set up with plain array syntax: u = 0.0_dp; u(1) = 0.0_dp; u(m) = 100.0_dp. There is no COMMON, no statement function, no GOTO, and no implicit typing left.

Phase 4 — Error Handling (Step 8) and the Final Check

Guard the precondition the arithmetic depends on — a rod needs at least three points to have an interior — and rebuild against the reference one last time:

if (m < 3) error stop 'solve: rod needs at least 3 points'

$ gfortran -std=f2018 -Wall -O2 rod.f90 -o rod && ./rod
sweeps = 40
   0.000  25.000  50.000  75.000 100.000

Compare to reference.txt. The profile is identical to the last digit and, as we argued, to the last bit: same operations, same order, same double precision. The sweep count is identical too, because the exit condition preserved the stopping rule exactly. The migration is numerically equivalent.

Phase 5 — What We Gained, Concretely

The modern rod.f90 is not merely tidier; the difference is checkable. implicit none guarantees the mistyped-variable class of bug cannot exist. The module and argument passing mean no two routines can disagree about the rod's storage. intent(inout) on u and intent(in) on tol mean the solver physically cannot corrupt the tolerance. Assumed-shape means the rod length travels with the array, so a mismatch is caught by -fcheck=all rather than silently walking off the end. And solve is now reusable: because it takes the rod as an argument rather than reading a global, you can relax a rod of any length — including the 4-, 3-, or 101-point cases — without editing the routine.

Discussion Questions

1. The statement function HALF was inlined because it was used once. If it had been used in five places, would you still inline it? What does an internal procedure give you that five copies of an expression do not?
2. We claimed bit-for-bit equivalence. Which single edit, made carelessly during Phase 2 or 3, would most plausibly have broken it while still producing a "reasonable-looking" profile? How would the regression test have caught it?
3. The legacy code over-dimensions U(MX) with MX = 101 but uses only M = 5. How does the modern assumed-shape version make the MX parameter unnecessary, and why is that safer?

Your Turn: Extensions

• Option A. Actually carry out the migration: type in rodeq.f, capture its reference, and modernize it one step at a time, rebuilding and diffing (numerically!) after each. Keep a commit per step.
• Option B. Add a second boundary condition option — an insulated right end (a Neumann condition, $u_M = u_{M-1}$) — to the modern solve, and hand-derive the new steady profile to check it.
• Option C. Instrument solve to print the max change each sweep, and confirm empirically that it decays geometrically. Estimate the ratio; relate it to why a 101-point rod needs far more sweeps than a 5-point one.

Key Takeaways

• The eight-step recipe is not just for the textbook's set piece; it is a repeatable procedure you run on real inherited code, mechanical steps first, structural steps second, guarded by a reference captured before you start.
• Hand-verify the answer the code converges to (here, the linear profile 0/25/50/75/100), not the incidental iteration count — that is what makes the regression reference trustworthy.
• Because the migration preserved the arithmetic, we could demand bit-for-bit agreement and get it. The proof that you modernized only the engineering is that the numbers did not move.