Chapter 35 — Key Takeaways (GPU Computing)

A one-page reference to offloading Fortran to a GPU with OpenACC and CUDA Fortran — and the data-movement discipline that decides whether it is worth it.

The vocabulary

Term Meaning
GPU / accelerator a throughput processor of thousands of simple cores; fast only on wide, regular, data-parallel work
host / device host = CPU + its RAM (control, serial, I/O); device = GPU + its own separate memory (the kernels)
offload hand a parallel region + its data to the device; keep serial flow and I/O on the host
OpenACC directive-based (!$acc), portable, incremental accelerator programming — the easy path
CUDA Fortran explicit GPU kernels in Fortran; NVIDIA-only, more control
kernel a procedure run on the device by many threads at once; CUDA Fortran: attributes(global)
data region !$acc data … end data; keeps arrays resident on the device so they transfer once, not per step
host–device transfer copying across the CPU–GPU bus; ~10–30× slower than device memory — usually the bottleneck

OpenACC — the essential directives

!$acc parallel loop        offload the following loop; run its iterations on the device
!$acc parallel loop collapse(2)          fold nested loops into one big iteration space
!$acc parallel loop reduction(+:total)   safe parallel sum/max/etc. (private partials, tree combine)
!$acc data  copy(u) create(tmp)   ... !$acc end data     resident data region (move once)
!$acc update self(u)       refresh the HOST copy from the device (a.k.a. update host)
!$acc update device(u)     refresh the DEVICE copy from the host
Data clause Entry Exit Use for
copyin(a) host → device — read-only inputs
copyout(a) allocate device → host write-only outputs
copy(a) host → device device → host read and written
create(a) allocate — device-only scratch
present(a) assert already there — data an enclosing region moved

CUDA Fortran — the explicit kernel

attributes(global) subroutine k(a, x, y, n)   ! kernel: called from host, runs on device
  real(dp), value :: a                          ! scalars pass by VALUE
  integer,  value :: n
  real(dp) :: x(n), y(n)                         ! dummy arrays live in DEVICE memory
  integer  :: i
  i = (blockIdx%x - 1)*blockDim%x + threadIdx%x  ! ONE-BASED global index (Fortran!)
  if (i <= n) y(i) = a*x(i) + y(i)               ! guard the tail
end subroutine k
! host side:
real(dp), device :: x_d(n), y_d(n)               ! device arrays
x_d = x                                          ! assignment IS the host->device copy
call k<<<nblocks, tpb>>>(a, x_d, y_d, n)         ! launch: <<<blocks, threads-per-block>>>
y = y_d                                          ! device->host copy (waits for the kernel)
  • Kernels must be module procedures (not internal contains procedures).
  • nblocks = (n + tpb - 1)/tpb covers n with tpb-thread blocks (256 is a sane default).
  • ONE-based indices: threadIdx%x, blockIdx%x count from 1 (unlike CUDA C) — subtract 1 from blockIdx%x.

Compile commands introduced

$ nvfortran -acc -Minfo=accel prog.f90 -o prog    # OpenACC on the GPU; -Minfo=accel reports offloads
$ nvfortran -acc=multicore prog.f90 -o prog       # OpenACC on the CPU's cores (portability)
$ gfortran -fopenacc prog.f90 -o prog             # gfortran also does OpenACC (GPU maturity varies)
$ nvfortran -cuda prog.f90 -o prog                # CUDA Fortran (or name the file prog.cuf)

The one commandment: minimize host–device transfer

WRONG (per-step transfer -> can be SLOWER than the CPU):
  do step = 1, nsteps
    !$acc parallel loop copy(u)     ! copies the whole field in AND out, every step
    ...
  end do

RIGHT (resident data -> transfer once):
  !$acc data copy(u) create(u_new)
  do step = 1, nsteps
    !$acc parallel loop present(u, u_new)   ! data already on the device
    ...
  end do
  !$acc end data                            ! copied back once, here

Move the data once; compute on it many times. For a memory-bound kernel (low flops per byte), per-step transfers cost more than all the computation — the difference between a GPU win and a GPU slowdown.

When a GPU helps — and when it doesn't

GPU wins GPU loses
massively data-parallel (a stencil over millions of cells) small (too little to cover overhead)
large, and reused on the device across many steps branchy / irregular (warp divergence, scattered access)
regular, branch-free arithmetic transfer-bound (low intensity, no reuse)
high arithmetic intensity (dense matmul) serial / dependency-heavy

Pitfalls

  • Per-step transfers — a copy on the loop-body kernel instead of a !$acc data region around the loop. The classic way an offload runs slower than the CPU.
  • CUDA C's zero-based index in CUDA Fortran — blockIdx%x*blockDim%x + threadIdx%x is off by one every block; CUDA Fortran needs (blockIdx%x - 1)*blockDim%x.
  • A reduction as a plain parallel loop — thousands of threads writing one accumulator is a data race; use reduction(+:var).
  • Expecting bit-for-bit CPU/GPU agreement — a parallel float reduction reorders the sum (not associative); compare by tolerance, not bits.
  • Forgetting Amdahl — offloading speeds only the offloaded part; the rest plus the transfers is the new serial fraction and caps the whole-program speedup.

Numbers and rules worth carrying

  • Host–device bus ≈ 10–30× slower than device memory — transfers dominate; minimize them.
  • Transfer time ≈ bytes ÷ bandwidth; a 32 MB field over 16 GB/s is ~2 ms one way, 4 ms round trip.
  • Arithmetic intensity = flops ÷ bytes moved; a stencil is low (~0.1–0.4) → wins only via on-device reuse.
  • <<<blocks, threads>>>; 256 threads/block is a reasonable default; always guard if (i <= n).

Project piece added this chapter (optional, advanced)

The heat solver's stencil update, offloaded to a GPU with OpenACC: a single !$acc data copy(u) create(u_new) region around the whole time loop keeps the field resident (transfer once in, once out), and each step runs !$acc parallel loop collapse(2) present(u, u_new) over the interior. The result is identical to Chapter 24's CPU solver — a correct offload reproduces the serial physics exactly. It pays only on a large grid swept many times; on the 5×5 toy grid the GPU is slower (too little work to fill it). To scale beyond one GPU, combine with MPI (Chapter 34): one rank per GPU, halos exchanged between devices — the multi-GPU configuration the Chapter 38 capstone can present. The solver has now run on every model in the taxonomy: images, threads, processes, and a graphics processor.