GPU Programming

A GPU is a separate computer: thousands of simple cores with their own memory, connected to the host by a fast bus (PCIe/NVLink). Fortran reaches it three ways — OpenACC directives, OpenMP target offload, and CUDA Fortran. The industrial standard is the directive route: one source, portable across NVIDIA/AMD/Intel devices and CPUs, with the compiler managing the details.

The offload model

Three facts organize everything. (1) Data lives on both sides — the host cannot see device memory and vice versa, so arrays must be copied across the bus. (2) The bus is slow, the cores are fast — a kernel needs far more arithmetic than bytes copied, or the transfer time eats the win. (3) Loops, not programs, are offloaded — each parallel region is a loop whose iterations the device threads execute thousands at a time.

Host and device with data copies and offloaded loop

Fig. 1 — One directive names the loop; the runtime schedules threads and the data region manages the copies.

Your first OpenACC program

OpenACC (!$acc) is the directive language created for exactly this: annotate an existing loop, compile with gfortran -fopenacc (or NVHPC's nvc -acc), run on device or CPU. The two crucial clauses: parallel loop names the offloaded loop, and data copyin/copyout declares which arrays cross the bus and when:

program saxpy
  use iso_fortran_env, only: real64
  implicit none
  integer, parameter :: n = 100000000
  real(real64), allocatable :: x(:), y(:)
  integer :: i
  allocate (x(n), y(n))
  call random_number(x)
  y = 1.0_real64
  !$acc data copyin(x) copy(y)          ! bring arrays to the device once
  !$acc parallel loop                   ! thousands of threads, one per i
  do i = 1, n
     y(i) = 2.0_real64 * x(i) + y(i)    ! SAXPY: the BLAS-1 hello world
  end do
  !$acc end parallel loop
  !$acc end data                        ! copy(y) writes the result back
  print '(a,i0,a,f8.3)', 'elements ', n, ' last y = ', y(n)
end program saxpy

Reading the directive zoo

OpenACC organizes itself in three layers; OpenMP target offload maps 1:1 onto the same mental model. Regions: parallel (using gang/worker/vector to shape the thread hierarchy), kernels (let the compiler dissect several loops at once), and the combined parallel loop. Data: copyin (host → device), copyout (device → host), copy (both), present (already on the device, skip the copy), and the data region that keeps arrays resident across many kernels — the single most important performance knob, since bus traffic is the real cost. Reductions: the reduction(+:acc) clause on the parallel loop, exactly as in OpenMP.

OpenMP target: the converging standard

Since OMP 4.0 the same directives family governs offload: !$omp target, !$omp teams distribute parallel do, and map(to/from:) clauses replacing the acc data specifiers. OpenMP 5.x added deep copy of allocatable arrays and requires clauses, and both vendors (GCC with -fopenmp -foffload=nvptx-none, LLVM flang, NVHPC) now treat target as the long-term portable path. The choice between OpenACC and OpenMP target is mostly organizational — the compiler support table in your institution decides it faster than any committee argument.

! The same SAXPY in OpenMP target dialect.
!$omp target data map(to: x) map(tofrom: y)
!$omp target teams distribute parallel do
  do i = 1, n
     y(i) = 2.0_real64 * x(i) + y(i)
  end do
!$omp end target teams distribute parallel do
!$omp end target data

CUDA Fortran and the fine print

When directives cannot express the algorithm — custom thread blocks, shared memory tiles, atomics, streams — CUDA Fortran (PGI/NVHPC) exposes the CUDA execution model in Fortran syntax: attributes(global) device kernels, attributes(device) variables, cudaMemcpy data transfers wrapped as cudaMemcpy() runtime calls. It is NVIDIA-only by design, whereas the directive routes remain portable. For almost all scientific kernels the directive route delivers 95% of peak with 5% of the effort; CUDA Fortran earns its place for a small set of hand-tuned kernels.

Measuring before believing

GPUs promise speedups that do not survive contact with reality: a memory-bound kernel (SAXPY included) gains little because the bus saturates; the win lives in compute-bound, data-resident loops — matrix–matrix multiply, stencils with present data, FFTs. Measure with nvprof/nsys (NVIDIA) or rocm tools, and compare against an optimized CPU run (-O3 -march=native) before claiming victory: the honest question is "is my kernel compute-bound and my data resident?", not "does the brochure promise 100×".

If no GPU is installed: the compiler still accepts the directives and runs the loops on the CPU (serial or OpenMP if you combine flags) — the syntax is learnable anywhere; the FeiTeng/OpenMPI-free laptops and CI runners are exactly where the demos for this lesson are designed to run.