Exercises: Procedures

These exercises drill the one habit and the one tool that make Fortran code robust: an intent on every argument, and arrays passed as assumed-shape. Type and compile everything — the compiler is a teacher here, and several problems are designed to be rejected by it so you learn to read its messages.

Difficulty: ⭐ warm-up · ⭐⭐ standard · ⭐⭐⭐ deeper. Solutions: worked solutions to the daggered (†) and odd-numbered problems are in appendices/answers-to-selected.md; the computational ones are also provided as compilable code in code/exercise-solutions.f90. Predict every output before you compile.

Throughout, assume use, intrinsic :: iso_fortran_env, only: dp => real64 and implicit none unless a snippet shows otherwise.

Part A — Warm-ups ⭐

6.1 † For each task, say whether you would write a subroutine or a function, and why: (a) compute the mean of a vector; (b) swap two variables; (c) return the area of a circle; (d) write a field to a file; (e) solve $A\mathbf{x}=\mathbf{b}$, returning the solution vector x.

6.2 Give the correct intent for each argument: (a) subroutine normalize(v, s) divides v by the scalar s in place and leaves s unchanged; (b) subroutine minmax(x, lo, hi) reads array x and returns its extremes; (c) subroutine accumulate(total, x) adds x to a running total.

6.3 Predict the exact printed line, then compile to check:

print '(f6.2)', twice(triple(2.0_dp))
! ...
pure function triple(x) result(y)
  real(dp), intent(in) :: x
  real(dp) :: y
  y = 3.0_dp * x
end function triple
pure function twice(x) result(y)
  real(dp), intent(in) :: x
  real(dp) :: y
  y = 2.0_dp * x
end function twice

6.4 The function skeleton below is missing its result clause and an intent. Add both so it compiles.

function hypotenuse(a, b)
  real(dp) :: a, b
  hypotenuse = sqrt(a**2 + b**2)
end function hypotenuse

Part B — Functions and Subroutines ⭐⭐

6.5 † Write a pure function norm2_of(x) result(nrm) that returns the Euclidean norm $\sqrt{\sum_i x_i^2}$ of an assumed-shape vector x(:). Test it on [3.0_dp, 4.0_dp] and predict the result. (Code solution provided.)

6.6 A colleague writes a function that "returns the mean and the standard deviation." Explain why a Fortran function cannot do this directly, and rewrite the interface as a subroutine.

6.7 † Predict the output, then compile:

real(dp) :: lo, hi
call bounds([4.0_dp, -1.0_dp, 7.0_dp, 2.0_dp], lo, hi)
print '(2f7.2)', lo, hi
! ...
subroutine bounds(x, mn, mx)
  real(dp), intent(in)  :: x(:)
  real(dp), intent(out) :: mn, mx
  mn = minval(x)
  mx = maxval(x)
end subroutine bounds

Part C — intent, the Contract ⭐⭐

6.8 † Find the bug. This will not compile. What does gfortran say, and what is the one-word fix?

subroutine normalize(v, s)
  real(dp), intent(in) :: v(:)
  real(dp), intent(in) :: s
  v = v / s
end subroutine normalize

6.9 Find the bug. This compiles but gives a wrong (and unpredictable) answer. Diagnose it and fix it without changing any intent.

subroutine add_all(x, total)
  real(dp), intent(in)  :: x(:)
  real(dp), intent(out) :: total
  integer :: i
  do i = 1, size(x)
     total = total + x(i)
  end do
end subroutine add_all

6.10 † Explain, in terms of "definable entities," why call swap(1.0_dp, 2.0_dp) fails to compile when swap's arguments are intent(inout), while call swap(a, b) (with variables a, b) succeeds.

Part D — Optional and Keyword Arguments ⭐⭐

6.11 † Write a pure function clamp(x, lo, hi) result(y) that returns x limited to the range [lo, hi], where hi is optional — if omitted, there is no upper limit. Predict the results of clamp(5.0_dp, 0.0_dp, 3.0_dp), clamp(-1.0_dp, 0.0_dp, 3.0_dp), and clamp(2.0_dp, 0.0_dp). (Code solution provided.)

6.12 A subroutine has the interface subroutine plot(x, y, color, width) where color and width are both optional. Write a call that supplies x, y, and width but not color. Why is a keyword mandatory here?

6.13 † Find the bug. Why is this dangerous, and what is the fix?

subroutine scale(x, factor)
  real(dp), intent(inout)        :: x(:)
  real(dp), intent(in), optional :: factor
  x = x * factor
end subroutine scale

Part E — pure and elemental ⭐⭐

6.14 † Write an elemental function sigmoid(x) result(s) computing $s = 1/(1 + e^{-x})$, and apply it to the array [-1.0_dp, 0.0_dp, 1.0_dp] in a single statement. Predict the three values to four decimals. (Code solution provided.)

6.15 For each, state whether it can be pure, and why: (a) a function returning sum(x) of its array argument; (b) a function that prints a debug line and returns 2*x; (c) a function that increments a module-level call counter and returns it; (d) a function returning x**2 + 1.0_dp. Which of the pure ones could also be elemental?

6.16 † Take the impure function from 6.15(b) and make it pure, or argue that it cannot be. If it can, show the corrected version and say where the debug output should go instead.

Part F — Recursion and Internal Procedures ⭐⭐⭐

6.17 † Write a recursive function gcd(a, b) result(g) implementing Euclid's algorithm (gcd(a, 0) = a, otherwise gcd(b, mod(a, b))). Predict gcd(48, 36) and gcd(1071, 462). (Code solution provided.)

6.18 Rewrite the recursive factorial from §6.5 as an iterative function using a do loop. Which version would you ship in a numerical library, and why?

6.19 † Find the hazard. The internal subroutine below shares the host's loop variable. Describe the bug this can cause and give the one-line fix.

program p
  implicit none
  integer :: i
  do i = 1, 3
     call greet()
  end do
contains
  subroutine greet()
    do i = 1, 5           ! uses the host's i
       ! ... some work ...
    end do
  end subroutine greet
end program p

Part G — Array Arguments, Design, and Estimation ⭐⭐/⭐⭐⭐

6.20 † Modernize it. Rewrite this FORTRAN 77-style routine in modern Fortran: assumed-shape, real(dp), intent, implicit none, and a whole-array operation. Note how many arguments disappear.

      SUBROUTINE DSCALE(A, N, S)
      INTEGER N, I
      DOUBLE PRECISION A(*), S
      DO 10 I = 1, N
         A(I) = A(I) * S
   10 CONTINUE
      END

6.21 Port it. Translate this NumPy-style Python function to an elemental Fortran function, then say why the Fortran version needs no loop even when applied to a large array.

def relu(x):
    return x if x > 0.0 else 0.0   # applied elementwise to a NumPy array

6.22 † Design it (heat solver). Extend step(field, alpha, dt) with an optional argument bc_value that, when present, resets the plate's four edges to that fixed temperature before the interior update. Write the new interface and the edge-setting code (use array sections). Why is optional the right choice rather than always requiring it?

6.23 Back of the envelope. In our step, allocate(old(nx,ny)); old = field copies the whole field each call. For a $1000 \times 1000$ grid of real(dp) run for 10,000 steps: (a) how many bytes is one copy? (b) how much data is copied over the whole run? (c) what cheap change removes the per-step allocate?

6.24 † Interleaved (Ch. 5). Write a pure function max_interior(f) result(mx) that returns the largest value among the interior cells of an assumed-shape 2-D array f(:,:) (exclude the first/last row and column) using a single maxval on an array section. Test it on the column-major matrix holding 1.0_dp through 16.0_dp in a $4\times4$ and predict the answer. (Code solution provided.)

6.25 Interleaved (Ch. 4 & 5). Write subroutine apply_dirichlet(field, edge) that sets all four boundaries of field(:,:) to the scalar edge using array-section assignments. Which intent does field need, and why not intent(out)?

6.26 † Find the bug. A program declares real(dp) :: g(5) and calls an external subroutine (defined in a separate file, not in a module or contains) whose dummy argument is x(:). The results are garbage. What is missing, and what are the two clean fixes?

6.27 Design it. Refactor the §6.1 describe subroutine into two pure functions — mean(x) and variance(x) — and rebuild describe on top of them. What did you gain, and what (if anything) did you lose?

6.28 † Back of the envelope. Someone proposes summing a $10^6$-element array with a recursive function that peels off one element per call. Estimate the recursion depth, explain why this risks a stack overflow where a do loop does not, and give the loop version's stack cost.

Solutions to the daggered and odd-numbered problems are in appendices/answers-to-selected.md; the computational ones (6.5, 6.11, 6.14, 6.17, 6.24) are compilable in code/exercise-solutions.f90. The "find the bug" problems are worth extra attention — being able to read a compiler's intent and interface errors fluently will save you more time than any other skill in this chapter.