Final Exam — Solutions

Introduction to Fortran Programming: The Language of Supercomputers

Full worked answers with point breakdowns. Every numeric result is hand-computed (no code was run to produce it) and stated exactly. Illustrative performance figures are labeled as such. Instructors: award partial credit generously for correct reasoning; the book's ethic is that a right argument with a slipped digit beats a right number with no reasoning.

Part A — Concepts (25 points)

A1. Floating point (4 pts)

(a) Why 0.1_dp + 0.2_dp /= 0.3_dp (2 pts). In base 2, 0.1 is a repeating fraction ($0.000110011001\ldots_2$), so it cannot be stored exactly in a finite 53-bit significand; it is rounded to the nearest representable grid point, a value very slightly larger than one-tenth. The same happens to 0.2 and to 0.3. When the two rounded values 0.1 and 0.2 are added, their small excesses add too, and the result lands one grid point (one ULP) away from the stored value of 0.3, which itself rounded the other way. Asked whether the two distinct grid points are equal, Fortran correctly answers .false.

Grading: 1 pt for "0.1 is a repeating binary fraction / not exactly representable," 1 pt for "the rounded sum lands one ULP off the rounded 0.3." Do not accept "floating point is imprecise" alone.

(b) Machine epsilon (2 pts). epsilon(1.0_dp) $= 2^{-52} \approx 2.22\times10^{-16}$. It is the gap between $1.0$ and the next larger representable double — the relative resolution of the format (equivalently, why double carries about 15–16 significant decimal digits). (1 pt value, 1 pt meaning. Accept the note that only the argument's kind matters, not its value.)

A2. Finite differences and the stencil (3 pts)

(a) (2 pts). Central difference: $f'(x) \approx \dfrac{f(x+h) - f(x-h)}{2h}$, which is second order, $O(h^2)$. Three-point second difference: $f''(x) \approx \dfrac{f(x+h) - 2f(x) + f(x-h)}{h^2}$, also $O(h^2)$. (½ pt each formula, ½ pt each order.)

(b) (1 pt). The bare neighbor-sum-minus-4-center is just a comparison of neighbors — a pure number with no units. Dividing by $h^2$ is what converts it into an actual second derivative (temperature per length squared). Drop the scaling, or use $h$ instead of $h^2$, and every term is off by a mesh-dependent factor, so the answer would change when you refine the grid — the signature of a broken discretization.

A3. CFL stability (3 pts)

(a) (2 pts). Stability of a single Fourier mode requires the per-step amplification factor to not grow: $|G| \le 1$. With $G = 1 - 8r$ and $r > 0$, the upper bound $1 - 8r \le 1$ is automatic; the binding constraint is the lower bound

$$ 1 - 8r \ge -1 \;\Longrightarrow\; 8r \le 2 \;\Longrightarrow\; r \le \tfrac14. $$

(b) (1 pt). $r = \alpha\Delta t/h^2$ depends on $h^2$, so halving $h$ quadruples $r$. If the old run sat between $\tfrac{1}{16}$ and $\tfrac14$, the finer grid pushes $r$ past $\tfrac14$; then $G = 1 - 8r < -1$, so the checkerboard mode (seeded by round-off) flips sign and grows every step, doubling without bound until the numbers overflow to Inf/NaN.

A4. Amdahl's Law (3 pts)

(a) (1 pt). $S(N) = \dfrac{1}{(1-p) + p/N}$; ceiling as $N \to \infty$ is $S_{\max} = \dfrac{1}{1-p}$.

(b) (1 pt). $p = 0.80 \Rightarrow S_{\max} = \dfrac{1}{1 - 0.80} = \dfrac{1}{0.20} = \mathbf{5\times}$ — no number of cores beats it.

(c) (1 pt). Gustafson asks the scaled question: "in fixed time on more cores, how much bigger a problem can I solve?" As the problem grows, the parallel work grows while the fixed serial part shrinks as a fraction, so speedup climbs nearly linearly with $N$ instead of hitting Amdahl's ceiling.

A5. The no-aliasing advantage (3 pts)

The Fortran standard forbids a procedure from having a written (defined) dummy argument aliased to any other argument or entity it also reads. Because keeping arguments distinct is the programmer's responsibility, the compiler may assume distinctness for free and reorder/vectorize loads and stores without any run-time overlap check — it knows writing one array can never change another. C pointers may alias by default, so the C programmer must add the C99 keyword restrict, per pointer, to promise the same thing. (1 pt the guarantee, 1 pt the optimization it unlocks, 1 pt restrict.)

A6. Column-major performance (3 pts)

The first index belongs on the inner loop. Fortran stores arrays column-major — the first index varies fastest — so a(i,j) and a(i+1,j) are adjacent in memory; sweeping the inner loop over i walks contiguous addresses, so each cache line (64 bytes = 8 doubles) fetched is used fully before the next is needed. The other order strides a whole column per step, touching a new cache line almost every access and wasting most of each one (numerical code is memory-bound). Both orders compute the identical result — only the speed differs. (1 pt inner index = first, 1 pt cache line + why, 1 pt "same result.")

A7. Coarrays, OpenMP, MPI (3 pts)

• OpenMP — shared memory: threads share one address space and coordinate through common variables (no messages). Confined to one node's cores.
• MPI — distributed memory: each process has private memory and they coordinate by explicit messages (send/receive). Scales across a whole cluster.
• Coarrays — a native, standardized Fortran parallel model (since Fortran 2008) that presents one notation across both shared and distributed memory: you reach a neighbor's data through a coindexed reference u(:)[q], and the runtime turns it into a local read or a network transfer. That single model spanning both worlds is what is unusual about it.

(1 pt per tool correctly placed.)

A8. Conditioning vs. stability (3 pts)

Conditioning measures how sensitive a problem's answer is to small changes in its inputs — a property of the problem (the mathematics), independent of any algorithm; an ill-conditioned problem loses digits no method can recover. Numerical stability measures how much extra error an algorithm adds on top of what conditioning already forces — a property of the algorithm. You can fix stability by choosing a better-arranged method (e.g., the conjugate trick, or LAPACK's pivoting); you cannot beat bad conditioning by changing the algorithm. (1 pt each definition + correct problem/algorithm attribution, 1 pt "stability is the one you can fix.")

Part B — Numerical Methods (20 points)

B1. An Euler step and an RK4 step (6 pts)

Here $f(t, y) = -y$, $y_0 = 1$, $h = 1$.

(a) Euler (1 pt). $$ y_1 = y_0 + h\,f(t_0, y_0) = 1 + 1\cdot(-1) = \mathbf{0}. $$ Euler overshoots the decay all the way to zero.

(b) RK4 (3 pts). The four stages: $$ \begin{aligned} k_1 &= f(0,\,1) = -1 \\ k_2 &= f\!\left(0.5,\; 1 + \tfrac12(1)(-1)\right) = f(0.5,\,0.5) = -0.5 \\ k_3 &= f\!\left(0.5,\; 1 + \tfrac12(1)(-0.5)\right) = f(0.5,\,0.75) = -0.75 \\ k_4 &= f\!\left(1,\; 1 + (1)(-0.75)\right) = f(1,\,0.25) = -0.25 \end{aligned} $$ Combine with weights $\tfrac16(1, 2, 2, 1)$: $$ y_1 = 1 + \tfrac{1}{6}\big(k_1 + 2k_2 + 2k_3 + k_4\big) = 1 + \tfrac{1}{6}\big(-1 - 1 - 1.5 - 0.25\big) = 1 + \tfrac{1}{6}(-3.75) = 1 - 0.625 = \mathbf{0.375}. $$

(Grading: 2 pts for the four stages correct, 1 pt for the weighted combination = 0.375. A single mistyped stage should collapse the answer — that is the book's point about RK4 coefficients being unforgiving.)

(c) Comparison (2 pts). Exact $e^{-1} \approx 0.3679$. Euler gives $0$ (error $\approx 0.37$); RK4 gives $0.375$ (error $\approx 0.007$, about fifty times smaller). RK4 samples the slope at four points within the step and takes a weighted average, matching the true solution's Taylor expansion through the $h^4$ term, where Euler follows a single left-endpoint slope and matches only through $h^1$ — so four extra function evaluations buy three extra orders of accuracy.

B2. Trapezoid vs. Simpson (5 pts)

Nodes $x = 0, 1, 2$; $f(x) = x^2$, so $f_0 = 0,\ f_1 = 1,\ f_2 = 4$; $h = 1$. Exact value $\int_0^2 x^2\,dx = \left[\tfrac{x^3}{3}\right]_0^2 = \tfrac83 \approx 2.6667$.

(a) Trapezoid (2 pts). $$ T_2 = h\big[\tfrac12 f_0 + f_1 + \tfrac12 f_2\big] = 1\big[\tfrac12(0) + 1 + \tfrac12(4)\big] = 0 + 1 + 2 = \mathbf{3}. $$

(b) Simpson (2 pts). $$ S_2 = \tfrac{h}{3}\big[f_0 + 4f_1 + f_2\big] = \tfrac{1}{3}\big[0 + 4(1) + 4\big] = \tfrac{1}{3}(8) = \tfrac{8}{3} \approx \mathbf{2.6667}. $$

(c) Which is exact (1 pt). Simpson returns the exact $\tfrac83$. This is not luck: Simpson fits a parabola through the three points, and the integrand $x^2$ is a parabola, so the fit is exact — Simpson integrates any quadratic (indeed any cubic) exactly. The trapezoid fits straight-line chords, which lie above this convex curve, so it overestimates ($3 > 2.6667$).

B3. A LAPACK dgesv call (5 pts)

(a) Decode and purpose (2 pts). dgesv = d (double precision) + ge (general matrix) + sv (solve a linear system). The call solves the dense general system $A\mathbf{x} = \mathbf{b}$ in double precision — here a $2\times2$ system.

(b) Contents on exit (2 pts). On return, b holds the solution $\mathbf{x}$ (it is overwritten in place), and a has been overwritten by the $L$ and $U$ factors of its LU decomposition (the original matrix is gone unless copied). Solving $$ \begin{bmatrix} 2 & 1 \\ 1 & 3 \end{bmatrix}\begin{bmatrix} x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} 5 \\ 10 \end{bmatrix}, \qquad \det = 2\cdot3 - 1\cdot1 = 5, $$ gives $x_1 = \dfrac{5\cdot3 - 1\cdot10}{5} = \dfrac{5}{5} = 1$ and $x_2 = \dfrac{2\cdot10 - 1\cdot5}{5} = \dfrac{15}{5} = 3$. Check: $2(1)+1(3) = 5$ ✓ and $1(1)+3(3) = 10$ ✓. So on exit b = [1.0, 3.0]. (The ipiv array holds the pivot row-interchanges from partial pivoting.)

(c) The info flag (1 pt). info < 0 signals a bad argument in your call — e.g., info = -4 means the 4th argument (lda) was illegal (your bug). info > 0 signals a numerical failure — for dgesv, info = i means $U_{ii}$ came out exactly zero, so $A$ is singular and no solution was computed. info == 0 guarantees the routine finished — not that the answer is accurate: an ill-conditioned $A$ can return info = 0 and a badly wrong $\mathbf{x}$ (use dgecon to estimate the condition number when it matters).

B4. A stencil update (4 pts)

Grid $4\times4$, top row $i=1$ held at $100$, all else $0$, $r = 0.2$. Interior cells are $(2,2),(2,3),(3,2), (3,3)$; edges are fixed.

(a) After one step (2 pts). For each interior cell, $u_{\text{new}} = u + r(\text{sum of 4 neighbors} - 4u)$:

• $(2,2)$: neighbors $u(1,2){=}100,\ u(3,2){=}0,\ u(2,1){=}0,\ u(2,3){=}0$; center $0$. $0 + 0.2(100 + 0 + 0 + 0 - 0) = 0.2(100) = 20$.
• $(2,3)$: neighbors $u(1,3){=}100,\ u(3,3){=}0,\ u(2,2){=}0,\ u(2,4){=}0$; center $0$. $0 + 0.2(100) = 20$.
• $(3,2)$: all four neighbors $0$; $0 + 0.2(0) = 0$.
• $(3,3)$: all four neighbors $0$; $0$.

after step 1:
  100  100  100  100
    0   20   20    0
    0    0    0    0
    0    0    0    0

(b) Cell $(2,2)$ after step 2 (2 pts). Read neighbors from the step-1 field: center $u(2,2) = 20$, neighbors $u(1,2){=}100$ (north, still hot), $u(3,2){=}0$ (south), $u(2,1){=}0$ (west), $u(2,3){=}20$ (east). $$ u_{\text{new}}(2,2) = 20 + 0.2\big(100 + 0 + 0 + 20 - 4\cdot20\big) = 20 + 0.2\big(120 - 80\big) = 20 + 0.2(40) = \mathbf{28}. $$ The snapshot matters because FTCS is defined in terms of the old neighbor values $u^n$. If you overwrote u in place, cell $(2,2)$ would read the already-updated $(2,1)$ or $(3,2)$, silently computing a different (Gauss–Seidel) scheme and the wrong physics.

(For reference, the full step-2 interior is row 2 = 28 28, row 3 = 4 4 — but the question only asks for $(2,2)$.)

Part C — Performance and Parallelism (25 points)

C1. An Amdahl calculation (8 pts)

$p = 0.96$, so $1 - p = 0.04$, and $S(N) = \dfrac{1}{0.04 + 0.96/N}$.

(a) (4 pts). $$ \begin{aligned} S(4) &= \frac{1}{0.04 + 0.96/4} = \frac{1}{0.04 + 0.24} = \frac{1}{0.28} \approx \mathbf{3.57\times} \\ S(8) &= \frac{1}{0.04 + 0.96/8} = \frac{1}{0.04 + 0.12} = \frac{1}{0.16} = \mathbf{6.25\times} \\ S(16) &= \frac{1}{0.04 + 0.96/16} = \frac{1}{0.04 + 0.06} = \frac{1}{0.10} = \mathbf{10.0\times} \end{aligned} $$ (Grading: each denominator shown correctly is worth the point even if the final division is slightly off.)

(b) Ceiling (2 pts). $S_{\max} = \dfrac{1}{1 - 0.96} = \dfrac{1}{0.04} = \mathbf{25\times}$. No number of cores beats it.

(c) Efficiency at 16 (1 pt). $E = S/N = \dfrac{10.0}{16} = \mathbf{0.625}$, i.e. 62.5% — over a third of the machine is already wasted on this fixed problem.

(d) Why measured < ideal (1 pt). Any two of: parallel overhead — thread fork/join cost, or MPI message latency, that did not exist in the serial program and grows with core count; memory-bandwidth saturation — the stencil is memory-bound, so once enough cores saturate the shared memory bus, more cores queue for starved memory (a second wall below Amdahl's); Amdahl's formula is the optimistic ceiling that assumes overhead-free parallelism. (7/16 ≈ 44% efficiency — the gap from the ideal $10\times$ is exactly this overhead.)

C2. A slow loop, and the fix (6 pts)

(a) Why it is slow (2 pts). Fortran is column-major: a(i,j) and a(i+1,j) are adjacent in memory, but a(i,j) and a(i,j+1) are a whole column (size(a,1) elements) apart. The inner loop here runs over j (the second index), so consecutive reads jump a full column and land on a different cache line almost every time — the processor fetches a 64-byte line, uses 8 bytes, and discards the rest, then does it again. It moves several times the necessary memory and stalls on it (the loop is memory-bound).

(b) The fix (3 pts). Swap the loops so the inner one runs over the first index i:

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

(Accept also s = sum(a), the whole-array intrinsic, which sweeps in memory order — an equally correct and idiomatic fix. Full marks for either.)

(c) Same value? (1 pt). The returned value is the same sum of the same elements (to within floating-point rounding, which can differ in the last bit because summation order changes and floating-point addition is not associative). For a plain sum this is negligible; the two loop orders visit the same numbers, so the result is the same to any meaningful precision.

C3. Add correct OpenMP scoping (6 pts)

s = 0.0_dp
!$omp parallel do default(none) shared(a, n) private(i, j) reduction(+:s)
do j = 1, n
  do i = 1, n
    s = s + a(i,j)**2
  end do
end do
!$omp end parallel do

Justification of each attribute (1 pt each, up to 5; +1 pt for a correct, compiling directive with default(none) and reduction):

• a — shared: read-only; every thread safely reads the same one copy.
• n — shared: read-only scalar bound.
• j — private: it is the !$omp do loop index, made private automatically, but under default(none) you must name it. Each thread needs its own.
• i — private: the inner loop index is not scoped automatically; left shared it would be a data race, with all threads clobbering one i. This is the single most common nested-loop scoping mistake.
• s — reduction(+:s): it accumulates a sum across all iterations. A plain shared s would be raced (lost updates, wrong and nondeterministic); a plain private s would discard each thread's partial. The reduction gives each thread a private accumulator and combines them once at the end — correct and race-free.

C4. Spot the race (5 pts)

(a) The raced variable (2 pts). The inner loop index j is not scoped, so it defaults to shared. All threads execute their inner loop over the same shared j, incrementing and reading it concurrently — a data race on the loop counter. Threads clobber each other's j, so iterations are skipped, repeated, or run with corrupt bounds, and the accumulated y differs from run to run.

(b) Why y(i) is not raced (1 pt). The outer index i is private and the !$omp do partitions the i-range across threads, so each thread owns a disjoint set of rows i and writes only its own y(i) elements. No two threads write the same y(i), so y itself is not raced — the bug is entirely the shared j.

(c) The fix (2 pts).

!$omp parallel do default(none) shared(a, x, y, n) private(i, j)
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

Add j to private. The one clause that would have caught this at compile time is default(none) — it refuses to guess an attribute for any variable, so the unscoped j becomes a compile error instead of a silent runtime race.

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

D1. Choose a CFL-safe timestep (8 pts)

$\alpha = 0.5$, square grid $h = 0.1$.

(a) Limit and largest stable $\Delta t$ (2 pts). In 2D, $r = \alpha\Delta t/h^2 \le \tfrac14$, i.e. $$ \Delta t \le \frac{h^2}{4\alpha} = \frac{(0.1)^2}{4(0.5)} = \frac{0.01}{2} = \mathbf{0.005}. $$

(b) 90%-safe $\Delta t$ and its $r$ (2 pts). $$ \Delta t = 0.9 \times 0.005 = \mathbf{0.0045}, \qquad r = \frac{\alpha\,\Delta t}{h^2} = \frac{0.5\times0.0045}{0.01} = \frac{0.00225}{0.01} = 0.225 = 0.9\times\tfrac14, $$ which is safely under $\tfrac14$. ✓

(c) Refine to $h = 0.05$ (2 pts). $$ \Delta t_{\text{limit}} = \frac{(0.05)^2}{4(0.5)} = \frac{0.0025}{2} = 0.00125, \qquad \Delta t = 0.9\times0.00125 = \mathbf{0.001125}. $$ That is exactly one quarter of the previous $0.0045$ ($0.0045/4 = 0.001125$). Halving $h$ quarters $\Delta t$ because $\Delta t \propto h^2$ — so you also take $4\times$ as many steps to reach the same final time. This is the "$\Delta t \sim h^2$ tax": refining an explicit diffusion solver costs $4\times$ the steps (and, with $4\times$ the cells in 2D, $16\times$ the work) per halving.

(d) stable_dt (2 pts).

pure function stable_dt(alpha, h) result(dt)
  real(dp), intent(in) :: alpha, h
  real(dp)             :: dt
  real(dp), parameter  :: safety = 0.9_dp
  dt = safety * h**2 / (4.0_dp * alpha)    ! 2D CFL limit, at 90% margin
end function stable_dt

Hard-coding $\Delta t$ is the most common way these solvers blow up because the moment anyone changes the resolution, a fixed $\Delta t$ that was safe at the old $h$ can silently violate $r \le \tfrac14$ at the new $h$ (recall halving $h$ quadruples $r$) — deriving $\Delta t$ from the grid keeps the run stable by construction.

(Accept the general form dt = safety / (2*alpha*(1/dx**2 + 1/dy**2)) for a non-square grid; it reduces to the above when dx = dy = h.)

D2. Decompose the domain for MPI (8 pts)

Global plate $400 \times 400$, $P = 8$ ranks, 1D row strips.

(a) Owned rows and ghosts (2 pts). Each rank owns $400 / 8 = \mathbf{50}$ rows. In u(nx, 0:nloc+1) the local rows 1..nloc (here nloc = 50) are owned; local rows 0 and nloc+1 are ghost cells (halo cells) — stored copies of the neighboring rank's edge rows (or the fixed physical boundary), not owned by this rank.

(b) Halo exchange, unchanged stencil (2 pts). Refreshing the ghost rows each step is the halo exchange (halo swap). Once the halo holds the neighbors' current edge values, a rank's interior update reads u(i, k-1) and u(i, k+1) uniformly — some of those are ghosts, but the stencil neither knows nor cares — so the exact serial Chapter-24 stencil runs on the owned rows with no modification. All the distributed-memory complexity is quarantined in the one exchange routine.

(c) Exchange volume and MPI_PROC_NULL (2 pts). An interior rank has a neighbor above and below, so it sends/receives two full rows per step: $2 \times 400 = \mathbf{800}$ real(dp) values. Rank 0 (top strip) has no neighbor above — that side is the physical top edge — so it exchanges only its one bottom row: $1 \times 400 = \mathbf{400}$ values. MPI_PROC_NULL is a null rank: a send/recv to it is a no-op that returns immediately, so setting the missing neighbor of the top and bottom strips to MPI_PROC_NULL lets the identical exchange code run on every rank with no boundary if, leaving the physical-edge ghost rows (the fixed Dirichlet values) untouched.

(d) Deadlock (2 pts). The failure is a deadlock: if every rank calls mpi_send first, each blocks inside its send waiting for a matching receive that the partner cannot post because it too is stuck in mpi_send. It "works" in small tests because mpi_send in standard mode may buffer small messages and return immediately (eager protocol); a large enough row (real problems have them) exceeds the buffer threshold, mpi_send switches to blocking rendezvous mode, and both ranks hang forever. The one-call fix is mpi_sendrecv, which bundles the send and receive into a single call that MPI guarantees will not deadlock.

D3. Validate against an analytical solution (8 pts)

(a) Verification vs. validation (2 pts). Verification asks "are we solving the equations right?" — does the code correctly solve the intended mathematical model (compare to an exact solution of the same PDE)? Validation asks "are we solving the right equations?" — does the model match physical reality (compare to experiment)? Comparing the solver to an exact solution of the same heat equation is verification. (Full validation would need a real heated plate in a lab, which this study does not have — and saying so is part of an honest write-up.)

(b) Substitution (3 pts). Let $u = \sin(\pi x)\sin(\pi y)\,e^{-2\alpha\pi^2 t}$.

Space — each second derivative brings down $-\pi^2$: $$ \frac{\partial^2 u}{\partial x^2} = -\pi^2 u, \qquad \frac{\partial^2 u}{\partial y^2} = -\pi^2 u \;\Longrightarrow\; \nabla^2 u = -2\pi^2 u, \quad\text{so}\quad \alpha\nabla^2 u = -2\alpha\pi^2 u. $$

Time — differentiating the exponential brings down its rate: $$ \frac{\partial u}{\partial t} = \big(-2\alpha\pi^2\big)\,u. $$

The two sides are equal, $\dfrac{\partial u}{\partial t} = \alpha\nabla^2 u$ ✓, so $u$ is an exact solution. (It also vanishes on all four edges — each sine is zero at $0$ and $1$ — matching the zero-Dirichlet boundary.)

(c) Observed order (2 pts). With errors $E_{\text{coarse}} = 1.0\times10^{-2}$ and $E_{\text{fine}} = 2.5\times10^{-3}$ at $h$ and $h/2$: $$ \text{ratio} = \frac{E_{\text{coarse}}}{E_{\text{fine}}} = \frac{1.0\times10^{-2}}{2.5\times10^{-3}} = 4.0, \qquad \text{observed order} = \log_2(4.0) = \mathbf{2}. $$ The error falls by $4\times$ when $h$ is halved, so the method is second-order accurate in space — exactly what the five-point stencil's $O(h^2)$ truncation error predicts. That agreement is the verification.

(d) Order 1 where 2 was expected (1 pt). No — it is not "close enough"; it is a bug announcing its own address. A first-order result where the method is second-order almost always means a broken discretization: a dropped $1/h^2$ scaling, a lopsided (one-sided) difference in place of the centered one, or an in-place update that mixes time levels. The wrong order is a specific diagnostic, not a rounding issue.

D4. Reproducibility (6 pts)

(a) Four things to record (3 pts; any four, 0.75 pt each).

• The governing equation, discretization, and scheme stated precisely (heat equation, five-point stencil, FTCS, the CFL timestep rule) — enough for a stranger to rebuild the method.
• The exact parameters and inputs: $\alpha$, domain, boundary values, initial condition, grid sizes, $\Delta t$ (or the safety factor), final time, and any config file / namelist.
• The build configuration: compiler and version (gfortran), optimization flags (-O3 -march=native, -fopenmp), and libraries/versions — plus the build manifest (fpm.toml).
• The code itself, published, with the regression test that pins the analytical-solution error, so the convergence result can be regenerated and can't silently drift.

(Accept also: random seeds if any; the machine/environment. The point is enough to regenerate the number.)

(b) Why bit-identical determinism is required (2 pts). The physics is deterministic, so a correct parallel solver must reproduce the serial answer exactly — if it does not, you have a bug (a scoping race, an in-place update), not a faster answer; matching the serial result digit-for-digit is how you prove the parallelization is correct. A legitimate perturbation of the last digits, even with a correct algorithm, is -ffast-math / -Ofast (it lets the compiler reorder floating-point additions, which are not associative) — or, equivalently, a parallel reduction whose summation order differs from the serial one. (Either example earns the point.)

(c) What the speedup claim needs (1 pt). The baseline it is measured against (8× versus what — an untuned serial run, or the best tuned serial code?) and the core count plus problem size and build (8× on how many cores, on what grid, with which compiler/flags). Without both, "8×" is a boast, not a measurement — and a slow baseline inflates any speedup.

Point summary

Chapter coverage. A1/A8 → Ch. 20 (floating point, conditioning); A2/B2/B4/D1/D3 → Ch. 22, 24 (finite differences, stencil, CFL, verification); A3/D1 → Ch. 24 (CFL); A4/C1 → Ch. 31 (Amdahl); A5/A6/C2 → Ch. 27 (no-aliasing, column-major); A7/C3/C4/D2 → Ch. 31, 33, 34 (parallel models, OpenMP scoping/races, MPI halo); B1 → Ch. 23 (Euler/RK4); B3 → Ch. 21 (dgesv); D3/D4 → Ch. 38 (verification, reproducibility). The exam is comprehensive but centered, as intended, on Parts V–IX.