GPU Programming
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.
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×".