Final Exam

Introduction to Fortran Programming: The Language of Supercomputers

Time: 3 hours. Total: 100 points.

Instructions. This is a comprehensive final. It draws on the whole book but is weighted toward Parts V–IX (Chapters 20–38): floating point, numerical methods, data and performance, and parallelism. Answer all four parts. A basic (non-programmable) calculator is permitted but not required; no computer, compiler, or interpreter may be used — in the spirit of the book, every numeric answer is to be hand-computed, and partial credit is given for correct reasoning even when the arithmetic slips. Show your work. All code you write should be modern Fortran (free-form, implicit none, real(dp), intent on arguments) that would compile with gfortran -std=f2018 -Wall.

Throughout, the running project is the 2D heat-equation solver you built across the book: $\dfrac{\partial u}{\partial t} = \alpha\nabla^2 u$ on a square plate, discretized with the five-point stencil and marched with the explicit FTCS scheme, with diffusion number $r \equiv \alpha\,\Delta t / h^2$.

Part A — Concepts (25 points)

Answer in a few sentences each. Precision of claim is graded, not length.

A1. Floating point (4 pts). (a) Explain why, in real(dp) arithmetic, 0.1_dp + 0.2_dp == 0.3_dp evaluates to .false.. Your answer must mention how 0.1 is stored. (b) State the value of epsilon(1.0_dp) (machine epsilon for double precision) as a power of two and as an approximate decimal, and say in one sentence what it measures.

A2. Finite differences and the stencil (3 pts). (a) Write the central-difference approximation to $f'(x)$ and state its order of accuracy in $h$; write the three-point approximation to $f''(x)$ and state its order. (b) The 2D five-point stencil for $\nabla^2 u$ on a square grid of spacing $h$ divides the neighbor-sum-minus- four-times-center by $h^2$. In one sentence, why is dividing by $h^2$ (rather than dropping it, or using $h$) not optional?

A3. CFL stability (3 pts). The explicit FTCS scheme for the 2D heat equation is stable only when $r = \alpha\Delta t/h^2 \le \tfrac14$. (a) Feeding the worst-case "checkerboard" mode through the update gives an amplification factor $G = 1 - 8r$. Show how the stability requirement $|G| \le 1$ yields $r \le \tfrac14$. (b) If you halve the grid spacing $h$ to sharpen the picture but leave $\Delta t$ and $\alpha$ unchanged, what happens to $r$, and why does the simulation blow up?

A4. Amdahl's Law (3 pts). (a) State Amdahl's Law for the speedup $S(N)$ on $N$ processors when a fraction $p$ of the work is parallelizable, and give the ceiling as $N \to \infty$. (b) A program is 80% parallelizable. What is the largest speedup it can ever achieve, on any number of cores? (c) In one sentence, what different question does Gustafson's Law answer that lets real supercomputers earn their keep?

A5. The no-aliasing advantage (3 pts). State "the no-aliasing advantage" from the compiler's point of view — what does the Fortran standard promise about a procedure's arguments, and what optimization does that promise unlock? Name the C99 keyword a C programmer must add, per pointer, to recover the same freedom.

A6. Column-major performance (3 pts). For a large real(dp) :: a(n,n) swept in a nested do loop, which index belongs on the inner loop for speed, and why? Name the hardware unit whose size (64 bytes = 8 doubles) drives the effect, and state whether the two loop orders compute the same result.

A7. Coarrays, OpenMP, MPI (3 pts). Place each of the three parallel tools on the memory-model map: for OpenMP, MPI, and coarrays, say whether the model is shared-memory or distributed-memory (coarrays: say what is unusual about it), and how the parallel workers coordinate (shared variables vs. explicit messages).

A8. Conditioning vs. stability (3 pts). Chapter 20 insists these are different. Define conditioning and numerical stability, say which is a property of the problem and which of the algorithm, and state which one you can fix by choosing a better method.

Part B — Numerical Methods: Read and Analyze the Code (20 points)

Every numeric answer here is hand-computable. Show the arithmetic.

B1. An Euler step and an RK4 step (6 pts). Consider the initial-value problem $\dfrac{dy}{dt} = -y$, $y(0) = 1$, whose exact solution is $y(t) = e^{-t}$ (so $y(1) = e^{-1} \approx 0.3679$).

(a) Compute one explicit Euler step of size $h = 1$ to estimate $y(1)$: $y_{1} = y_0 + h\,f(t_0, y_0)$.

(b) Compute one classical RK4 step of size $h = 1$ to estimate $y(1)$. Show all four stages $k_1, k_2, k_3, k_4$ and the final combination $y_1 = y_0 + \tfrac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)$.

(c) Compare both estimates to the exact $e^{-1}$. In one sentence, why is RK4 so much closer despite Euler's step being "simpler"?

B2. Trapezoid vs. Simpson (5 pts). Estimate $\displaystyle\int_0^2 x^2\,dx$ (whose exact value is $\tfrac{8}{3} \approx 2.6667$) using two panels ($n = 2$, so $h = 1$; nodes at $x = 0, 1, 2$).

(a) Apply the composite trapezoidal rule $T_2 = h\big[\tfrac12 f_0 + f_1 + \tfrac12 f_2\big]$.

(b) Apply Simpson's rule $S_2 = \tfrac{h}{3}\big[f_0 + 4f_1 + f_2\big]$.

(c) One of your two answers is exact. Which one, and why is that not a coincidence?

B3. A LAPACK dgesv call (5 pts). Read the following program fragment.

integer,  parameter :: n = 2
real(dp) :: a(n,n), b(n)
integer  :: ipiv(n), info

a = reshape([ 2.0_dp, 1.0_dp,   &
              1.0_dp, 3.0_dp ], [n, n], order=[2, 1])
b = [ 5.0_dp, 10.0_dp ]

call dgesv(n, 1, a, n, ipiv, b, n, info)

(a) Decode the routine name dgesv letter by letter, and state in one sentence what the call computes. (b) After the call returns with info == 0, what does the array b hold? What has happened to a? Compute the numeric contents of b on exit. (The system is $2x_1 + x_2 = 5$ and $x_1 + 3x_2 = 10$.) (c) A colleague never checks info. Give one value of info that signals their bug (a bad argument) and one that signals a numerical failure, and state what info == 0 does — and does not — guarantee.

B4. A stencil update (4 pts). The explicit FTCS interior update on a square grid ($r = \alpha\Delta t/h^2$) is:

u_new(i,j) = u(i,j) + r*( u(i-1,j) + u(i+1,j) + u(i,j-1) + u(i,j+1) - 4.0_dp*u(i,j) )

Take a $4\times4$ plate with the top row (i = 1) held at 100 and every other cell 0 (Dirichlet edges), $r = 0.2$. Only the four interior cells $(2,2), (2,3), (3,2), (3,3)$ update; edges are held fixed.

(a) Compute the full field after one step. (Hint: the interior of row 2 sees the hot edge; the interior of row 3 does not, yet.) (b) Compute the value of cell $(2,2)$ after a second step. You must read the neighbor values from the step-1 field, not update in place — explain in one sentence why that distinction matters.

Part C — Performance and Parallelism (25 points)

C1. An Amdahl calculation (8 pts). Profiling shows the heat solver spends 96% of its run time in the parallelizable stencil sweep and 4% in irreducibly serial work (setup, the time-loop bookkeeping, periodic output). So $p = 0.96$.

(a) Using $S(N) = \dfrac{1}{(1-p) + p/N}$, compute the exact speedup on $N = 4$, $N = 8$, and $N = 16$ processors. Show each denominator. (b) What is the hard ceiling as $N \to \infty$? (c) Compute the parallel efficiency $E = S/N$ at $N = 16$. (d) You run it on 16 real cores and measure only $7\times$, not the ideal you found in (a). Give two reasons (named in Chapters 31/33/38) why a measured speedup falls below the Amdahl ceiling.

C2. A slow loop, and the fix (6 pts). This subroutine sums a large matrix and is measured to be several times slower than it should be.

pure function total(a) result(s)
  real(dp), intent(in) :: a(:,:)
  real(dp)             :: s
  integer              :: i, j
  s = 0.0_dp
  do i = 1, size(a,1)          ! outer over rows
    do j = 1, size(a,2)        ! inner over columns
      s = s + a(i,j)
    end do
  end do
end function total

(a) Explain precisely why this loop nest is slow on a large array — name the memory layout and what happens to each cache line under this access pattern. (b) Rewrite the loop nest so it runs with the grain of memory. (You may change only the loop structure.) (c) Does your fix change the value returned? Justify your answer.

C3. Add correct OpenMP scoping (6 pts). Parallelize this serial reduction with OpenMP so it runs correctly on any number of threads. Write the full directive using default(none), and then justify the data-sharing attribute you gave to each variable (a, n, i, j, s).

s = 0.0_dp
do j = 1, n
  do i = 1, n
    s = s + a(i,j)**2
  end do
end do

C4. Spot the race (5 pts). A colleague parallelizes a matrix–vector product $\mathbf{y} = A\mathbf{x}$ and gets a different answer every run:

!$omp parallel do shared(a, x, y, n) private(i)
do i = 1, n
  do j = 1, n
    y(i) = y(i) + a(i,j) * x(j)
  end do
end do
!$omp end parallel do

(a) Identify the variable that is being raced and explain, in terms of what threads do to it, why the answer varies run to run. (b) Note that y(i) itself is not raced here even though it is written in the loop. Why not? (c) Give the corrected directive, and name the one clause that would have turned this silent runtime race into a compile-time error.

Part D — Design and Synthesis: The Heat Solver (30 points)

This part is one connected design problem on the running project. You are preparing the solver for a large parallel run and writing it up.

D1. Choose a CFL-safe timestep (8 pts). You will run the explicit 2D solver with thermal diffusivity $\alpha = 0.5$ on a square grid of spacing $h = 0.1$.

(a) State the 2D stability limit on $\Delta t$ in terms of $h$ and $\alpha$, and compute the largest stable $\Delta t$. (b) Your code chooses a timestep at 90% of the limit (a safety margin). Compute that $\Delta t$, and compute the resulting $r = \alpha\Delta t/h^2$ to confirm it is safely under $\tfrac14$. (c) You then refine the grid to $h = 0.05$ for more resolution. Compute the new safe $\Delta t$. By what factor did it change, and why is that the notorious "$\Delta t \sim h^2$ tax" of explicit diffusion? (d) Write a pure function stable_dt(alpha, h) that returns the 90%-safe timestep. It must be modern Fortran (kind-parameterized real, intent(in)), and it must derive the timestep from the grid, not hard-code it — say in one sentence why hard-coding $\Delta t$ is the single most common way these solvers blow up.

D2. Decompose the domain for MPI (8 pts). You will distribute a global plate of 400 rows × 400 columns across 8 MPI ranks using a 1D row-strip decomposition (rank 0 gets the top strip, rank 7 the bottom).

(a) How many interior rows does each rank own? Each rank stores its strip as u(nx, 0:nloc+1). What are the rows at local indices 0 and nloc+1 called, and what do they hold? (b) Before each update a rank must refresh those rows. Name this operation, and say why the serial Chapter-24 stencil then runs unchanged on each rank's owned rows. (c) How many real(dp) values does an interior rank exchange per step, and how many does the top rank (rank 0) exchange? Explain the role of MPI_PROC_NULL at the physical top and bottom edges. (d) The obvious exchange — "every rank calls mpi_send then mpi_recv" — passes small tests and then hangs in production. Name the failure, explain why it appears only at large message sizes, and give the one-call fix.

D3. Validate against an analytical solution (8 pts). To earn a reviewer's trust you verify the solver against the exact solution $u(x,y,t) = \sin(\pi x)\sin(\pi y)\,e^{-2\alpha\pi^2 t}$ on the unit square with all edges held at zero.

(a) Distinguish verification from validation, and state which one this comparison performs. (b) Show, by taking the derivatives, that this $u$ satisfies $\partial u/\partial t = \alpha\nabla^2 u$. (You need $\nabla^2 u$ and $\partial u/\partial t$; confirm they match.) (c) You run a convergence study at fixed $r$, halving $h$. The maximum error over the plate is $1.0\times10^{-2}$ on the coarse grid and $2.5\times10^{-3}$ on the grid with $h$ halved. Compute the observed order of accuracy and say what it confirms about the five-point stencil. (d) A different refactor of the code measures order 1 where you expected order 2. Is that "close enough"? State what it most likely means.

D4. Reproducibility (6 pts). Your write-up claims "the solver is second-order accurate and achieves an $8\times$ speedup." A reviewer will ask whether someone else could regenerate that result.

(a) List four distinct things you must record or publish so that the numerical result (the field, the convergence order) is reproducible by a stranger. (b) The parallel (OpenMP) solver must produce results bit-identical to the serial solver. Why is that determinism required rather than merely nice — and name one thing (a compiler flag or a coding choice) that can legitimately perturb a floating-point result in the last digits even when the algorithm is correct. (c) The speedup claim "$8\times$" is nearly meaningless as stated. Name the two pieces of information it must be reported with before a reviewer can judge it.

End of exam. Re-read Part D — the design questions reward a clear argument more than a long one. Correctness of every claim and number is what earns the marks, exactly as the book insists for the code itself.