Exercises: GPU Computing

These exercises exercise the two roads onto the device — OpenACC directives and CUDA Fortran kernels — and, above all, the discipline that decides whether an offload is worth doing: managing host–device data movement. Because this book has no GPU to compile against, the "type, compile, run" problems ask you to predict the result (a correct offload reproduces the CPU answer exactly) and, where you have access to a GPU and nvfortran, to confirm it yourself. The back-of-the-envelope problems are the most valuable here: transfer arithmetic tells you, before you write a line of kernel code, whether the GPU will help at all.

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 worked as runnable, hand-checked code in code/exercise-solutions.f90 (that file is plain Fortran and compiles with gfortran -std=f2018 -Wall — it is the arithmetic about GPU programs, not GPU code). Try every problem before you look. OpenACC snippets assume nvfortran -acc; CUDA Fortran assumes nvfortran -cuda.

Part A — The Host/Device Model ⭐

35.1 † In one sentence each, define host and device, and state the single most important consequence of their having separate memory spaces.

35.2 A CPU has a few large cores; a GPU has thousands of small ones. For which of these workloads does the GPU's design win, and for which does it lose? (a) applying the same five-point stencil to ten million grid cells; (b) a recursive tree traversal with data-dependent branching at every node; (c) multiplying two large dense matrices; (d) a loop of 20 iterations, each depending on the previous.

35.3 † Explain "offloading" in one sentence, and name two parts of a typical program you would not offload to the device.

35.4 True or false, with a one-sentence justification: "A GPU is simply a faster CPU, so any program will run faster if you move it to the GPU."

Part B — OpenACC Directives ⭐⭐

35.5 † Predict the exact printed output of this program (a correct offload gives the CPU answer):

integer, parameter :: n = 6
real(dp) :: v(n)
integer  :: i
v = [(real(i, dp), i = 1, n)]
!$acc parallel loop copy(v)
do i = 1, n
  v(i) = v(i)*v(i) + 1.0_dp
end do
print '(6f7.1)', v

35.6 Name the OpenACC data clause you would use for each array: (a) an input the kernel only reads; (b) an output the kernel only writes; (c) an array the kernel reads and writes; (d) scratch that never leaves the device; (e) an array an enclosing !$acc data region already placed on the device.

35.7 † Rewrite this loop with an OpenACC directive so it offloads to the GPU, moving x and y in and the result z out — and no more than necessary:

do i = 1, n
  z(i) = x(i) + 2.0_dp*y(i)
end do

35.8 (Port it.) Here is an OpenMP loop from Chapter 33. Translate it to the OpenACC equivalent that offloads to a GPU, and name the one new concern OpenACC has that the OpenMP version did not.

!$omp parallel do
do i = 1, n
  c(i) = a(i)*b(i)
end do
!$omp end parallel do

35.9 † What does the -Minfo=accel flag do when you compile with nvfortran -acc, and why is it one of the first flags you should reach for when an offload is not behaving?

Part C — CUDA Fortran Kernels ⭐⭐

35.10 † This kernel is meant to compute c = a + b elementwise, but its thread-index formula was copied from a CUDA C tutorial and is wrong in CUDA Fortran. Find the bug and fix it. (Hint: with a single 256-thread block and n <= 256 it writes nothing at all; in general it is off by a full block.)

attributes(global) subroutine add_kernel(a, b, c, n)
  real(dp) :: a(n), b(n), c(n)
  integer, value :: n
  integer :: i
  i = blockIdx%x*blockDim%x + threadIdx%x
  if (i <= n) c(i) = a(i) + b(i)
end subroutine add_kernel

35.11 You launch a kernel over n = 1000 elements with 256 threads per block. (a) How many blocks do you need? (b) How many threads are launched in total? (c) Why does the kernel need an if (i <= n) guard?

35.12 † In CUDA Fortran, why must a kernel (attributes(global)) be a module procedure, and why are scalar arguments like a and n given the value attribute?

35.13 Name three differences between the OpenACC SAXPY of §35.2 and the CUDA Fortran SAXPY of §35.3, and state, for a scientist writing new code, which you would reach for first and why.

Part D — Data Movement and Residency ⭐⭐

35.14 † A time loop runs 10,000 steps, each offloaded with !$acc parallel loop copy(u) on the per-step kernel. (a) How many host–device transfers of u does the whole run perform? (b) Rewrite the structure so the field crosses the bridge only twice for the entire run, and give the transfer count.

35.15 In the residency pattern, you still need to write a VTK frame to disk every 500 steps. How do you get the current field from the device to the host for output without ending the data region? Give the exact directive.

35.16 † Explain, in terms of arithmetic intensity (flops per byte moved), why a single five-point stencil sweep is a poor GPU offload on its own but a heat solver that sweeps a resident field 10,000 times is a good one.

35.17 True or false, with justification: "If my kernel runs 20× faster than the CPU loop, my program will run about 20× faster."

Part E — Find the Bug ⭐⭐

Each snippet is flawed. Diagnose it and give the corrected version.

35.18 † A colleague's offloaded heat solver is 3× slower than the CPU version. The time loop looks like this. What is the cause, and what is the structural fix?

do step = 1, nsteps
  !$acc parallel loop collapse(2) copy(u) create(u_new)
  do j = 2, ny-1
    do i = 2, nx-1
      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))
    end do
  end do
  u(2:nx-1,2:ny-1) = u_new(2:nx-1,2:ny-1)
end do

35.19 This OpenACC loop gives the wrong answer on the GPU. Why? (Hint: is each iteration independent?)

!$acc parallel loop copy(a)
do i = 2, n
  a(i) = a(i) + a(i-1)          ! running partial sum
end do

35.20 † A CUDA Fortran programmer ports an index formula straight from a CUDA C tutorial: i = blockIdx%x*blockDim%x + threadIdx%x, then guards with if (i >= 1 .and. i <= n). Their results are wrong — even with a single block, elements go untouched. Explain why the formula is wrong in CUDA Fortran, and give the correct index.

35.21 A developer writes real(dp) :: a_d(n) (no device attribute) and then a_d = a before launching a kernel that expects a device array. Explain what is wrong and give the one-word fix.

Part F — Design It and Back of the Envelope ⭐⭐⭐

Confirm the daggered computational answers against code/exercise-solutions.f90.

35.22 † (Back of the envelope — transfer time.) A $1000\times1000$ field of real(dp) is offloaded over a bus with 16 GB/s bandwidth (take 1 GB $= 10^9$ bytes). (a) How many bytes is the field? (b) How long does one host→device copy take? (c) A round trip? (d) If the GPU sweeps the stencil in 0.05 ms but the host→device and device→host copies happen every step, is the GPU faster or slower per step than a CPU that sweeps in 1.0 ms? Show the arithmetic.

35.23 (Back of the envelope — the crossover.) Same field. With the residency pattern, the field is transferred once in and once out (round trip 1.0 ms total), and each GPU sweep costs 0.05 ms; the CPU sweeps in 1.0 ms. After how many steps $K$ does the resident GPU version become faster than the CPU? At $K = 10{,}000$, what is the speedup?

35.24 † (Back of the envelope — arithmetic intensity.) A five-point stencil update reads 5 values and writes 1 (real(dp), 8 bytes each) and does about 6 floating-point operations per cell. Estimate its arithmetic intensity in flops per byte, counting (a) all 6 memory touches, then (b) only the 2 that must reach main memory in a well-cached sweep (one read miss + one write). Why is a low intensity the warning sign of a transfer-bound kernel?

35.25 † (Design it — solver.) Sketch, in directives only (no full program), how you would offload the heat solver's time loop to the GPU so that: the field is resident across all steps; each step's stencil update runs on the device; and a frame is written to disk every 100 steps. Mark where each transfer happens.

35.26 (Back of the envelope — Amdahl with transfer.) A solver is 98% stencil (offloadable) and 2% serial setup/I-O. You make the stencil 25× faster on the GPU, but the one-time transfer adds an amount equal to 1% of the original run time to the serial part. Using Amdahl's Law with the transfer folded into the serial fraction, estimate the whole-program speedup. Compare to the no-transfer ideal.

Part G — Interleaved (Chapters 5 and 33) ⭐⭐

35.27 † (Ch. 5.) The offloaded stencil sweeps u(i,j) with collapse(2). Recalling column-major storage from Chapter 5, which index varies fastest in memory, and why does mapping that dimension across neighbouring GPU threads matter (name the GPU term)?

35.28 (Ch. 33.) You have the same 98%-parallel stencil running two ways: OpenMP on 16 CPU cores, and OpenACC on a GPU. Both obey Amdahl's Law. Name the one extra cost that enters the GPU version's serial fraction but not the OpenMP version's, and say why the shared-memory model is spared it.

Solutions to the daggered and odd-numbered problems are in appendices/answers-to-selected.md; the computational ones (35.5, 35.22, 35.23, 35.24, 35.26) are worked as runnable, hand-checked code in code/exercise-solutions.f90. The design problems (35.25) give a model answer; a plan that keeps the field resident, offloads the per-step stencil, and uses !$acc update self only for output is on the right track — the transfer count is the score.