Chapter 35: GPU Computing — OpenACC and CUDA Fortran for Graphics Processors

"More science, less programming." — motto of the OpenACC initiative

Overview

In Chapter 31 we drew a map of parallelism with three territories: shared memory (the cores of one node), distributed memory (the nodes of a cluster), and a third that we sketched and then set aside — the graphics processor, a different kind of chip entirely, with thousands of small cores built to do the same arithmetic to enormous quantities of data at once. Chapters 32 through 34 conquered the first two territories with coarrays, OpenMP, and MPI. This chapter is the third, and it closes Part VIII.

A modern GPU is the most parallel device most programmers will ever touch — a single card can hold many thousands of arithmetic units and deliver a peak throughput that dwarfs the CPU beside it. For the right computation — a stencil swept over a million-cell grid, a matrix multiplied, a field of particles advanced — it is transformative, and it is exactly the shape of computation Fortran was born for. But the GPU comes with a catch that will dominate this chapter and every real offload you ever write: it has its own memory, separate from the CPU's, and getting your data across the narrow bridge between the two is slow. Master that one fact and the GPU is a gift; ignore it and your "accelerated" code will run slower than the plain CPU version it replaced. We will make that failure happen on purpose so you never make it by accident.

We take two roads onto the device. The first is OpenACC — directives you sprinkle onto ordinary Fortran, exactly like the OpenMP of Chapter 33, that ask the compiler to move the work to the GPU for you. It is the gentle path: portable, incremental, and often enough. The second is CUDA Fortran — explicit GPU kernels you write and launch yourself, NVIDIA-only, more work and more control. Most scientists should start, and often finish, with OpenACC; but knowing what a kernel is demystifies the whole machine, so we build one by hand too.

Everything here obeys an honesty this book insists on and this chapter must underline heavily. We have no GPU and no nvfortran to compile against, so none of this code was run, and every speedup you read is illustrative — a plausible order of magnitude, never a promise. The OpenACC and CUDA Fortran syntax is written to be exactly correct, but it is not machine-checked here; where a detail is compiler-specific, we say so. Treat the numbers as arithmetic to reason with, not benchmarks to quote.

In this chapter, you will learn to:

• Describe what a GPU is — a throughput machine of thousands of simple cores — and the host/device model that governs every program that uses one.
• Offload a data-parallel loop to the GPU with OpenACC (!$acc parallel loop), and keep data on the device across many steps with a !$acc data region so you do not pay to move it every step.
• Write an explicit CUDA Fortran kernel — attributes(global), device arrays, a <<<grid,block>>> launch — and get the one-based thread indexing right.
• See why host–device transfer is almost always the bottleneck, and structure your code to minimize it.
• Judge, before you invest a week, whether your computation will win on a GPU at all.

Learning Paths

How to read this chapter by track. - ⚡ HPC ("I need the accelerator") — this is your chapter; read all of it. §35.2 (OpenACC) and §35.4 (data movement) are what you will use on Monday; the Project Checkpoint is the pattern you will copy. - 🔬 Scientist ("my solver is too slow") — read §35.1 for the model and §35.2 for the easy path, then §35.5 to decide whether your problem suits a GPU at all before you spend effort. You can skim §35.3. - 📖 Standard — read straight through; this chapter completes the parallel-programming toolkit and the running project's journey across every model. - 🔧 Legacy ("I inherited GPU code") — §35.1 and §35.3 tell you what attributes(global), <<<>>>, and device arrays mean when you meet them in an old CUDA Fortran file; §35.4 explains the data-region scaffolding around them.

35.1 GPUs as Massively Parallel Processors; the Host/Device Model

Start with the hardware, because the GPU is not "a faster CPU" and treating it like one is the root of most GPU disappointment. A CPU and a GPU are optimized for opposite goals.

A CPU is a latency machine. It has a handful of large, sophisticated cores — a few, or a few dozen — each engineered to finish one thread of instructions as fast as possible: deep pipelines, out-of-order execution, branch prediction, and large caches, all spent making a single stream of work race to the finish. It is a small team of brilliant generalists.

A GPU is a throughput machine. It spends its transistors not on making one core clever but on packing in thousands of small, simple ones, each individually slower and dumber than a CPU core, but collectively able to perform a staggering number of arithmetic operations per second — provided they are all doing essentially the same thing at the same time. It is an army of thousands of laborers who thrive on identical, repetitive work and stall the moment you ask them each to do something different. The GPU executes threads in lock-step groups (NVIDIA calls a group of 32 a warp): all 32 threads run the same instruction at the same time on different data. Give them a million grid cells to update with the same stencil and they are unstoppable; give them a snarl of branches where each thread takes a different path and most of them sit idle while the others work.

Definition (GPU; accelerator). A graphics processing unit (GPU) is a processor built around a very large number of simple cores that execute the same operation on many data elements simultaneously — a throughput-optimized, data-parallel machine, originally designed to shade millions of pixels and now used for general numerical computation. An accelerator is the general term for any such specialized device attached to a CPU to speed up a particular kind of work; a GPU is by far the most common accelerator in scientific computing, and in this chapter the two words are effectively interchangeable.

The GPU does not work alone. It is a device that plugs into a host computer, and the two form a partnership with a strict division of labour. The CPU — the host — runs your program, controls the flow, does the serial work, and issues commands to the GPU. The GPU — the device — sits idle until the host hands it a chunk of parallel work, does that work at enormous speed, and hands the results back. Crucially, they have separate memories. The host's arrays live in the computer's main RAM; the device's arrays live in the GPU's own on-board memory. Neither can directly read the other's.

Definition (host; device). In accelerator programming, the host is the CPU and its main memory, which runs the main program and orchestrates the computation; the device is the accelerator (the GPU) and its separate on-board memory, which executes the parallel work the host sends it. Host and device memories are distinct address spaces: a value the host computes is not visible to the device until it is explicitly copied across, and vice versa. Keeping straight which memory a given array lives in — and when it must be moved — is the central discipline of every GPU program.

That separation forces a three-step rhythm on every GPU computation, and the whole chapter is really about performing this dance well:

        HOST (CPU + RAM)                          DEVICE (GPU + its own memory)
   +--------------------------+                 +-------------------------------+
   |  program flow, serial    |   1. copy in    |  thousands of cores,          |
   |  work, I/O               | ==============> |  huge memory bandwidth        |
   |                          |   (host->dev)   |                               |
   |  arrays in main RAM      |                 |  arrays in device memory      |
   |                          |   3. copy back  |                               |
   |                          | <============== |  2. COMPUTE (the parallel     |
   |                          |   (dev->host)   |     kernel runs here)         |
   +--------------------------+                 +-------------------------------+
        ^                                                   ^
    latency-optimized:                              throughput-optimized:
    few clever cores                                thousands of simple cores

Offloading is the name for handing a piece of computation to the device.

Definition (offload). To offload a computation is to move a portion of a program — typically a data-parallel loop or kernel and the data it needs — from the host to the device for execution, then bring the results back. A program that uses a GPU is a host program that offloads its hot, parallel regions to the device while keeping its serial control flow, setup, and I/O on the host. The art is choosing what to offload (enough parallel work to be worth it) and how to manage the data it touches (as few crossings of the host–device bridge as possible).

💡 Intuition: Think of the host as a head chef and the GPU as a warehouse full of a thousand line cooks across the street. The line cooks are fast, but there is one narrow door between the buildings. If the chef sends a huge tray of identical prep work across once, lets the thousand cooks tear through it, and receives one tray back, the arrangement is a triumph. If instead the chef runs across the street for every single onion — carry it over, have it chopped, carry it back, repeat — the thousand cooks spend their lives waiting at the door and the whole scheme is slower than chopping onions in the main kitchen. The cooking is fast; the doorway is the problem. Hold that image: §35.4 is entirely about the doorway.

🔗 Connection. This is the "GPU / accelerator" row of the taxonomy table from Chapter 31, now made concrete. There we noted a GPU is "dramatically better at wide, regular, number-crunching work and dramatically worse at branchy, irregular, one-thing-at-a-time logic," and that "using one means shipping data across to the device, computing there, and shipping results back — and that transfer is very often the bottleneck." Every sentence of that preview is a section of this chapter. The GPU is the one model in that taxonomy with its own memory, and that is what makes it both powerful and treacherous.

🔄 Check Your Understanding. 1. In one sentence each, contrast what a CPU core and a GPU core are each optimized for. 2. What does it mean that the host and device have "separate memory spaces," and what consequence does that force on every GPU program? 3. What is "offloading," and which parts of a program do you typically not offload?

Answers

(1) A CPU core is latency-optimized — few large cores that finish a single instruction stream as fast as possible; a GPU core is one of thousands of throughput-optimized simple cores that are fast only in lock-step on the same operation across much data. (2) Host RAM and device memory are distinct; neither can read the other's arrays directly, so data must be explicitly copied across the bus before the device can compute on it and after, to retrieve results. (3) Offloading is handing a parallel region (and its data) to the device to execute. You typically keep serial control flow, setup, and I/O on the host and offload only the hot, data-parallel kernels.

35.2 OpenACC: Directive-Based Offload — the Easiest Path

The gentlest way onto the GPU is to not write GPU code at all — to keep writing ordinary Fortran and add directives that ask the compiler to do the offloading for you. This is OpenACC, and if you learned OpenMP in Chapter 33, you already know the shape of it.

Definition (OpenACC). OpenACC (Open Accelerators) is an open standard for directive-based parallel programming of accelerators. You annotate ordinary Fortran (or C/C++) with structured comments beginning !$acc, and an OpenACC-aware compiler generates the code that moves data to the device, runs the annotated loops there in parallel, and moves results back. Like OpenMP directives, they are comments to a compiler that does not understand them, so the same source still compiles and runs correctly (just on the CPU) without OpenACC support — the offload is additive and portable.

The parallel with OpenMP is exact and worth making explicit, because it is your fastest way in. In Chapter 33 you wrote !$omp parallel do above a loop to spread its iterations across CPU threads. In OpenACC you write !$acc parallel loop above a loop to spread its iterations across GPU cores. Same idea, different target.

🔗 Connection — OpenMP directives → OpenACC directives. OpenACC was deliberately designed to feel like OpenMP, and the two share a heritage; OpenMP itself later grew its own accelerator-offload directives (!$omp target). The mental model transfers wholesale: a directive is a comment that asks the compiler to parallelize the following loop, data-scoping clauses say what is shared and what is private, and the underlying loop must have independent iterations for it to be correct. If you internalized the fork–join discipline of Chapter 33, OpenACC is that discipline pointed at a different piece of hardware. The one genuinely new worry the GPU adds — because it has separate memory — is data movement, which OpenMP on a shared-memory node never had to think about. That new worry is §35.4.

Here is the simplest possible offload: a SAXPY, the "scaled vector add" $y \leftarrow a\,x + y$ that is the "hello, world" of numerical kernels. Each element is independent, so the loop is perfectly data-parallel.

program acc_saxpy
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  integer,  parameter :: n = 8
  real(dp) :: x(n), y(n), a
  integer  :: i

  x = [(real(i, dp), i = 1, n)]      ! x = 1, 2, 3, ..., 8
  y = 1.0_dp                          ! y = 1 everywhere
  a = 2.0_dp

  !$acc parallel loop copyin(x) copy(y)   ! offload: copy x and y over, run on the GPU
  do i = 1, n
    y(i) = a*x(i) + y(i)              ! SAXPY, one independent element per GPU thread
  end do

  print '(a, 8f6.1)', 'y = a*x + y = ', y
end program acc_saxpy

The directive !$acc parallel loop does three things at once: it copies the arrays named in its data clauses to the device, runs the loop's iterations in parallel across GPU cores, and (via copy) brings y back afterward. Because every y(i) is computed from its own x(i) and y(i) with no cross-talk, the compiler can hand each iteration to a different thread. The result is arithmetic you can check by hand — y(i) = 2i + 1:

$ nvfortran -acc -Minfo=accel example-01-acc-saxpy.f90 -o saxpy && ./saxpy
y = a*x + y =    3.0   5.0   7.0   9.0  11.0  13.0  15.0  17.0

Two things about that compile command matter. The -acc flag turns on OpenACC and, with nvfortran (the NVIDIA HPC compiler you met as an option in Chapter 30), targets the GPU by default. The -Minfo=accel flag is one you will use constantly: it makes the compiler tell you what it offloaded — which loops became GPU kernels, what data it moved, how it mapped the parallelism — so you are never guessing. And a portability note the honesty of this book requires: gfortran also implements OpenACC, via the -fopenacc flag, though its GPU-offload maturity trails nvfortran's; OpenACC can even target the CPU's own cores with nvfortran -acc=multicore, so a single directive-annotated source can run three ways — serial, multicore CPU, or GPU — from the same file. That portability is much of OpenACC's appeal.

Data clauses and the data region

The clauses copyin, copyout, and copy are how you tell OpenACC which arrays to move and in which direction. They are the most important words in the language, because — as §35.4 will hammer — the movement, not the compute, is usually what costs you.

Putting data clauses directly on !$acc parallel loop, as the SAXPY did, is fine for a single offloaded loop. But real programs offload the same data to many kernels in a row — most of all, a time-stepping loop that runs the same update thousands of times. If each step re-copies the field to and from the device, you pay the transfer cost thousands of times over. The cure is to hoist the data movement out of the loop into a data region that surrounds it.

Definition (data region). A data region is a block of code, delimited in Fortran by !$acc data and !$acc end data, over which named arrays are kept resident in device memory. The data clauses on the !$acc data directive move the arrays once — into the device at region entry, back out at region exit — and every !$acc parallel loop inside the region reuses that resident copy without moving it again (referring to it with present). A data region is how you separate the decision to move data from the decision to compute on it, so that data crosses the host–device bridge as few times as possible — ideally once in and once out for an entire multi-step computation.

We will build the canonical example — a data region wrapped around the heat solver's whole time loop — in the Project Checkpoint, because it is the single most important pattern in the chapter. For now, hold the principle: move the data once, compute on it many times.

OpenACC also handles the other everyday parallel pattern, a reduction — combining many values into one, a sum or a max — with a reduction clause identical in spirit to OpenMP's from Chapter 33: !$acc parallel loop reduction(+:total). The GPU computes the partial sums across thousands of threads and combines them for you; doing that correctly by hand is surprisingly involved, which is one more reason the directive path earns its keep (Case Study 2 pulls that thread).

🐍 Python Comparison. Python's GPU story is a family of libraries that mirror the two roads of this chapter. CuPy gives you a NumPy-shaped array that lives on the GPU, so a + b on CuPy arrays runs on the device — pleasant, and, like OpenACC, it hides the kernels from you. Numba can compile a decorated Python function into a CUDA kernel, closer to the explicit CUDA Fortran of §35.3. Both, underneath, call the very same NVIDIA GPU runtime that nvfortran targets, and both face the identical host–device transfer tax: a CuPy program that shuttles arrays between NumPy (host) and CuPy (device) every iteration is slow for exactly the reason a badly-structured OpenACC program is. The lesson is language-independent. Where Fortran pulls ahead is the same place it always does — when the kernel is a tight, compiled, array-shaped numerical loop, nvfortran generates device code as good as anything, with none of the interpreter overhead, and the !$acc directives sit directly on the Fortran you already wrote.

🔄 Check Your Understanding. 1. What does the directive !$acc parallel loop ask the compiler to do, and what must be true of the loop for it to be correct? 2. Which data clause would you use for an array the device only reads as input? Which for scratch that never needs to leave the device? 3. Why does wrapping a time loop in a !$acc data region beat putting copy on each step's kernel?

Answers

(1) It asks the compiler to run the following loop's iterations in parallel on the accelerator (moving any data named in clauses); the loop's iterations must be independent — no iteration may depend on another's result — or the parallel result is wrong. (2) copyin for a read-only input; create for device-only scratch (allocated on the device, never copied either way). (3) A !$acc data region moves the data once (in at entry, out at exit) and every kernel inside reuses the resident copy; putting copy on each step's kernel re-transfers the whole field every step, paying the expensive host–device crossing thousands of times.

35.3 CUDA Fortran: Explicit Kernels — More Control, NVIDIA-Only

OpenACC asks the compiler to generate the GPU code. CUDA Fortran lets you write it yourself. It is more work and less portable, but it hands you direct control over the device — exactly how the parallel threads are organized, what lives in which kind of device memory, how the work is launched — and there are kernels you can express and tune in CUDA Fortran that a directive cannot reach. For most scientific work OpenACC is the right first (and often last) tool; but seeing a kernel built by hand removes the mystery from the whole machine, and you will meet CUDA Fortran in existing NVIDIA-tuned codes.

Definition (CUDA Fortran). CUDA Fortran is a set of extensions to Fortran, defined by NVIDIA (from the PGI compiler, now nvfortran), that let you write explicit GPU kernels in Fortran and launch them from host code. It exposes NVIDIA's CUDA programming model directly: you mark a procedure as running on the device, declare arrays that live in device memory, and launch the kernel across a grid of thousands of threads with a special syntax. It is NVIDIA-specific — it runs only on NVIDIA GPUs, compiled only by nvfortran — which is the price of its control, and the contrast with portable OpenACC.

The central new object is the kernel.

Definition (kernel). A kernel is a procedure that runs on the device, executed simultaneously by many GPU threads, one thread per data element. In CUDA Fortran a kernel is a subroutine marked attributes(global) — "global" meaning it is called from the host but runs on the device. Every one of the thousands of launched threads executes the same kernel body; what makes them do different work is that each thread computes its own unique index from built-in variables and operates on its own slice of the data.

Here is the same SAXPY as §35.2, written the explicit way. Read it against the OpenACC version and the trade — control for verbosity — is on the page:

module saxpy_mod
  use, intrinsic :: iso_fortran_env, only: dp => real64
  use cudafor
  implicit none
contains
  attributes(global) subroutine saxpy_kernel(a, x, y, n)
    real(dp), value :: a               ! scalar arguments pass by VALUE from the host
    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   ! this thread's global index (1-based!)
    if (i <= n) y(i) = a*x(i) + y(i)   ! guard: threads past n do nothing
  end subroutine saxpy_kernel
end module saxpy_mod

program cuda_saxpy
  use, intrinsic :: iso_fortran_env, only: dp => real64
  use cudafor
  use saxpy_mod
  implicit none
  integer,  parameter :: n = 8, tpb = 256    ! tpb = threads per block
  real(dp) :: x(n), y(n), a
  real(dp), device :: x_d(n), y_d(n)         ! _d arrays live on the device
  integer  :: i, nblocks

  x = [(real(i, dp), i = 1, n)]
  y = 1.0_dp
  a = 2.0_dp

  x_d = x                                     ! host -> device: array assignment copies
  y_d = y
  nblocks = (n + tpb - 1)/tpb                 ! enough blocks to cover n (=1 here)

  call saxpy_kernel<<<nblocks, tpb>>>(a, x_d, y_d, n)   ! LAUNCH across the grid

  y = y_d                                     ! device -> host (waits for the kernel)
  print '(a, 8f6.1)', 'y = a*x + y = ', y
end program cuda_saxpy

$ nvfortran -cuda example-02-cuda-saxpy.f90 -o saxpy && ./saxpy
y = a*x + y =    3.0   5.0   7.0   9.0  11.0  13.0  15.0  17.0

Identical output to the OpenACC SAXPY — the same physics, computed the hard way. Now walk the machinery, because every piece is a genuine CUDA Fortran feature:

• attributes(global) subroutine marks the kernel: called from the host, runs on the device. Kernels must be module procedures (they cannot be internal contains procedures of the program), which is why the kernel lives in saxpy_mod.
• real(dp), device :: x_d(n) declares device arrays. The device attribute puts them in GPU memory. The plain assignment x_d = x triggers the host-to-device copy (it compiles to a cudaMemcpy); y = y_d copies back. Assignment is the transfer — clean, but easy to do too often.
• <<<nblocks, tpb>>> is the launch configuration: it says "run this kernel on nblocks blocks of tpb threads each." The GPU organizes its thousands of threads into a grid of equal-sized blocks; you choose the block size (256 is a common, sane default) and compute how many blocks cover your data.
• i = (blockIdx%x - 1)*blockDim%x + threadIdx%x is how each thread finds its own element. threadIdx%x is the thread's position within its block, blockIdx%x its block's position in the grid, and blockDim%x the block size — combine them and you get a unique global index per thread.
• if (i <= n) is the guard. We launched 256 threads to cover 8 elements, so threads 9–256 have i > n and must do nothing; without the guard they would write out of bounds. This "launch a round number of threads, guard the tail" pattern is universal in GPU code.

The single most important — and most Fortran-specific — subtlety is hiding in that index arithmetic.

⚠️ Common Pitfall — CUDA Fortran thread indices are ONE-based. In CUDA C, threadIdx.x and blockIdx.x count from 0, and the global index is blockIdx.x*blockDim.x + threadIdx.x. In CUDA Fortran, matching the language's one-based arrays, threadIdx%x runs from 1 to blockDim%x and blockIdx%x from 1 to gridDim%x, so the correct global index is (blockIdx%x - 1)*blockDim%x + threadIdx%x. Copy a formula out of a CUDA C tutorial verbatim and you will be off by one at every block boundary — thread 1 of block 2 will compute the wrong element, silently corrupting your results past the first block. This is the classic trap when porting CUDA C knowledge to CUDA Fortran; write the -1 on blockIdx%x by reflex. (blockDim and gridDim are ordinary counts, not indices, so they are not shifted — only the two Idx variables are one-based.)

The compile flag is -cuda (older PGI/nvfortran versions used -Mcuda); naming the source file with a .cuf extension enables CUDA Fortran automatically, without the flag. And use cudafor brings in the CUDA runtime and the built-in kernel variables. There is no gfortran equivalent — CUDA Fortran is nvfortran only, the sharp edge of the portability trade you accept for the control.

🔄 Check Your Understanding. 1. What does attributes(global) mean, and where must such a subroutine be defined? 2. In the launch call k<<<nblocks, tpb>>>(...), what do the two numbers control? 3. Why is the global index (blockIdx%x - 1)*blockDim%x + threadIdx%x in CUDA Fortran rather than the CUDA C form blockIdx*blockDim + threadIdx?

Answers

(1) It marks a kernel — a procedure called from the host that runs on the device, executed by many threads at once; it must be a module procedure (not an internal procedure). (2) The launch configuration: nblocks = number of thread blocks in the grid, tpb = threads per block; together they set how many threads run the kernel. (3) Because CUDA Fortran's threadIdx%x and blockIdx%x are one-based (to match Fortran's one-based arrays), unlike CUDA C's zero-based versions; you subtract 1 from blockIdx%x to get the correct offset, or you are off by one at every block boundary.

35.4 Host↔Device Data Movement — Usually the Bottleneck

Now the fact that governs everything. Return to the doorway between the two kitchens. The GPU's on-board memory is ferociously fast — its bandwidth, the rate at which its cores can stream data from its own memory, is far higher than a CPU's access to its RAM. That internal bandwidth is much of why the GPU is fast. But the bridge between host and device — historically the PCI Express bus, and even on the tightest modern interconnects — is dramatically slower than either processor's access to its own memory. Data crossing that bridge is the slowest step in most GPU programs by a wide margin.

Definition (host–device transfer). A host–device transfer is the copying of data across the bus that connects the CPU and the GPU — from host memory to device memory before a kernel can use it, and back after. Its bandwidth is much lower, and its per-transfer latency much higher, than either processor's access to its own memory, so in the great majority of GPU programs the transfers, not the on-device computation, dominate the run time. The first rule of GPU performance follows directly: minimize host–device transfers — move data to the device once, keep it resident while you compute on it repeatedly, and bring back only what you must, as rarely as you can.

Put rough numbers to it to feel the asymmetry (illustrative orders of magnitude, not a spec sheet — Tier 2). A GPU's internal memory bandwidth is commonly on the order of hundreds of gigabytes to a terabyte or more per second; the host–device bus delivers on the order of tens of gigabytes per second, often 10–30× slower than the GPU's own memory. So copying an array across the bridge can take longer than dozens of passes of computation over it on the device. If your kernel does only a little arithmetic per element — as a stencil does, a few adds and a multiply — then a single round-trip transfer of the field can cost more than all the computation you came to the GPU to accelerate. The computation was never the problem; the commute was.

This is why the data region of §35.2 is not a convenience but the whole game. The failure mode is so common, and so devastating, that it deserves to be seen happening.

⚠️ Common Pitfall — per-step transfers that quietly kill your performance. The most common way a GPU port ends up slower than the CPU code it replaced is transferring data every time step. It happens when a data clause sits on the per-step kernel instead of on a region around the whole time loop:

fortran ! WRONG: copies the ENTIRE field to the device and back on EVERY step. do step = 1, nsteps !$acc parallel loop collapse(2) copy(u) ! <-- copy in AND out, every step do j = 2, ny-1 do i = 2, nx-1 u_new(i,j) = ... ! a few flops per cell end do end do end do

For a memory-bound stencil (little arithmetic per value), the two transfers per step swamp the tiny computation, and the "accelerated" run can be several times slower than the plain CPU version — you have paid the commute ten thousand times to do ten thousand blinks of work. The fix is structural: hoist the data movement into a region around the entire loop, so the field crosses the bridge once in and once out, and every step computes on the resident copy:

fortran ! RIGHT: one transfer in, one transfer out, for the whole run. !$acc data copy(u) create(u_new) do step = 1, nsteps !$acc parallel loop collapse(2) present(u, u_new) ! data already on device do j = 2, ny-1 do i = 2, nx-1 u_new(i,j) = ... end do end do end do !$acc end data

Same kernel, same result — but the transfers went from 2*nsteps down to 2. That structural change, nothing else, is often the entire difference between a GPU win and a GPU embarrassment.

Here is a minimal, self-contained demonstration of the residency pattern: an array doubled three times inside a single data region. The point is not the arithmetic (trivially checkable) but the transfer count — one in, one out, for all three passes.

program acc_data_region
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  integer,  parameter :: n = 5, nsteps = 3
  real(dp) :: a(n)
  integer  :: i, step

  a = [(real(i, dp), i = 1, n)]        ! a = 1, 2, 3, 4, 5

  !$acc data copy(a)                    ! ONE copy in now; ONE copy out at end data
  do step = 1, nsteps
    !$acc parallel loop present(a)      ! a is already resident: NO transfer here
    do i = 1, n
      a(i) = 2.0_dp*a(i)               ! double every element, on the device
    end do
  end do
  !$acc end data                        ! results copied back to the host here, once

  print '(a, 5f7.1)', 'a after 3 doublings = ', a
end program acc_data_region

$ nvfortran -acc -Minfo=accel example-03-acc-data-region.f90 -o dregion && ./dregion
a after 3 doublings =     8.0   16.0   24.0   32.0   40.0

The field starts at $[1,2,3,4,5]$ and is doubled three times, so each element is multiplied by $2^3 = 8$: $[8, 16, 24, 32, 40]$. The !$acc data copy(a) moved a to the device once and retrieves it once at !$acc end data; the three parallel loop present(a) kernels ran on the resident copy with no transfer between them. Had we written !$acc parallel loop copy(a) on the loop instead, the same result would have cost six transfers instead of two.

When you do need the device's data on the host mid-computation — to write a checkpoint, print a diagnostic, or plot a frame — you do not tear down the data region; you use !$acc update self(a) (equivalently update host(a)) to refresh the host's copy of a from the device without ending residency, and !$acc update device(a) for the reverse. Periodic output from a resident field is exactly an update self every so often, and it is how the Project Checkpoint's solver would write a frame every hundred steps without surrendering the GPU.

🚪 Threshold Concept — on a GPU, optimize the data movement, not the arithmetic. Coming from the CPU world of Part VII, your instinct for "make it fast" is about arithmetic and cache: unroll the loop, vectorize it, respect column-major order. On the GPU those still matter, but they are the second question. The first question — the one that decides whether the GPU helps at all — is "how much data crosses the host–device bridge, and how often?" A kernel that computes twice as fast but transfers its data every step loses to a kernel half as clever that keeps its data resident. Once you internalize that the bridge is the bottleneck, you stop thinking "which loop do I offload?" and start thinking "what is the largest region I can keep resident on the device, crossing the bridge as few times as possible?" That reframing — from accelerating loops to managing residency — is what separates a GPU program that wins from one that merely compiles. It is the GPU's version of the theme that has run through this whole book: performance is not accidental. On an accelerator, it is a property of your data movement first and your arithmetic second.

⚡ Performance Note — arithmetic intensity decides. The quantity that predicts whether an offload is transfer-bound is its arithmetic intensity: flops performed per byte moved. A dense matrix multiply of $n \times n$ matrices does $O(n^3)$ arithmetic on $O(n^2)$ data — intensity grows with $n$, so for large matrices the compute dwarfs the transfer and the GPU shines. A single stencil sweep does a constant few flops per cell on every cell — low, fixed intensity — so a lone sweep is transfer-bound and a poor offload. What rescues the stencil is iteration: keep the field resident and sweep it ten thousand times, and the one-time transfer is amortized across ten thousand cheap sweeps, lifting the effective intensity of the whole run. That is precisely why the heat solver — many steps over resident data — is a good GPU candidate while one Laplacian is not, and why "move once, compute many times" is the rule that makes it so.

🔄 Check Your Understanding. 1. Why is host–device transfer, rather than on-device computation, usually the bottleneck in a GPU program? 2. A colleague's offloaded stencil is slower on the GPU than on the CPU. What is the single most likely cause, and the structural fix? 3. You need to write the field to disk every 100 steps without ending the device residency. What OpenACC directive do you use?

Answers

(1) Because the host–device bus bandwidth is far lower (often 10–30×) than either processor's access to its own memory; for kernels that do little arithmetic per element, moving the data costs more than computing on it. (2) Transferring the field every step (a data clause on the per-step kernel); the fix is to wrap the whole time loop in a !$acc data region so the field is moved once in and once out and each step uses the resident copy (present). (3) !$acc update self(field) (a.k.a. update host) — it refreshes the host copy from the device without ending the data region.

35.5 When GPUs Help, and When They Don't

A GPU is not a faster computer; it is a differently shaped one, and the whole skill of using it well is matching your computation to its shape. The honest answer to "should I use a GPU?" is "it depends on the computation," and this section is how to tell before you spend a week finding out.

A GPU rewards work that is:

• Massively data-parallel. Thousands of independent, identical operations — one per grid cell, per particle, per pixel. The GPU has thousands of cores; it needs thousands of independent work-items to fill them. This is the data parallelism of Chapter 31, and it is exactly the shape of a stencil sweep, a whole-array operation, a Monte-Carlo ensemble.
• Large. The one-time transfer and the kernel-launch overhead must be amortized over enough work to be worth it. A million-cell grid swept ten thousand times is ideal; a hundred-cell grid swept twice is not — the overheads swamp the work.
• Regular and branch-free. Because threads run in lock-step warps, they thrive when every thread takes the same path. Uniform arithmetic over a rectangular array is the sweet spot.
• High in arithmetic intensity, or reused on the device. Enough flops per byte moved (a dense matrix multiply), or enough reuse of resident data (a solver's many iterations) to hide the transfer.

A GPU punishes work that is:

• Small. Too little work to cover the transfer and launch overhead — you spend more setting up the GPU than computing.
• Branchy or irregular. Heavy data-dependent branching makes threads in a warp diverge — take different paths — and the hardware must run the paths serially, one after another, so a warp where every thread branches differently can run at a fraction of peak. Irregular memory access (chasing pointers, scattered indices, unstructured meshes) defeats the coalesced, streaming access the GPU is built for.
• Transfer-bound. Low arithmetic intensity with no reuse — one cheap pass over data you had to ship across the bridge. The commute costs more than the trip was worth.
• Serial, or riddled with dependencies. If each step depends tightly on the last with little parallel work inside, there is nothing for thousands of cores to do.

📜 From History — from pixels to physics. The GPU was not designed for science. It was designed to render graphics: transforming and shading millions of independent pixels and vertices per frame, which is the purest data parallelism imaginable — the same operations, on vast arrays, with no cross-talk. By the early 2000s researchers noticed that "shade a million pixels" and "update a million grid cells" are the same computation wearing different clothes, and began smuggling scientific work through the graphics pipeline. In 2007 NVIDIA released CUDA, a way to program the device for general computation directly, and general-purpose GPU computing (GPGPU) was born. CUDA Fortran and OpenACC are how that capability reaches the Fortran world. The lineage explains the shape: a GPU is spectacular at exactly the workloads that look like graphics — big, regular, embarrassingly parallel arithmetic over arrays — and mediocre at everything else, because that is the job it was built to do. Your stencil solver wins on a GPU because, to the hardware, a diffusing plate and a shaded surface are the same problem.

And Amdahl's Law has the last word, because the GPU does not repeal it — it sharpens it. Offloading accelerates only the part you offload; everything else, plus the transfers, becomes your new serial fraction. If the stencil is 98% of your run time and you make it 20× faster on the GPU while the serial 2% and the transfers are unchanged, Chapter 31's ceiling still bounds you: the un-accelerated remainder now dominates. The reflex from Chapter 31 is exactly right here — profile first, know your serial fraction, and estimate the ceiling before you port — with one GPU-specific addition: count the transfer time as part of the serial fraction, because it is work the accelerator adds, not removes.

🐛 Find the Bug. A team reports their GPU-accelerated particle code "gets no speedup, sometimes a slowdown." Profiling shows the kernel itself is 15× faster than the CPU loop, yet the program is not. The code offloads one short kernel per timestep inside a long time loop, with the particle array declared copy on each kernel, and the per-kernel arithmetic is a handful of flops per particle. What is happening, and what is the fix?

Answer

The kernel is fast but the program is transfer-bound: copy on each per-step kernel ships the entire particle array to the device and back every timestep, and for a low-arithmetic-intensity kernel those two transfers per step cost far more than the 15×-faster computation saves — sometimes more than the whole CPU version. The fix is structural, not a faster kernel: wrap the time loop in a !$acc data copy(particles) region so the array is resident for the whole run (transferred once in, once out), and mark the per-step kernels present. Use !$acc update self only when a frame must be written. This is the §35.4 pitfall in a different domain — the kernel was never the problem; the commute was.

🔄 Check Your Understanding. 1. List two properties of a computation that make it a good GPU candidate and two that make it a bad one. 2. What is warp divergence, and why does it hurt GPU performance? 3. How does Amdahl's Law apply to a GPU offload, and what GPU-specific cost must you count in the serial fraction?

Answers

(1) Good: massively data-parallel, large, regular/branch-free, high arithmetic intensity or heavy on-device reuse. Bad: small, branchy/irregular, transfer-bound (low intensity, no reuse), or serial/dependency-heavy. (2) Warp divergence is when threads in a lock-step group (a warp) take different branches; the hardware must execute the divergent paths serially, so a heavily branching warp runs at a fraction of peak. (3) Offloading speeds only the offloaded part; the rest plus the host–device transfers become the serial fraction, and Amdahl's ceiling $1/(1-p)$ still caps the whole-program speedup — so you count transfer time as serial and estimate the ceiling before porting.

Project Checkpoint

This checkpoint is optional and advanced — it needs a GPU and nvfortran, which most readers building the solver on a laptop will not have — but it is the natural climax of the data-region idea, and it completes the running project's tour across every parallel model. We offload the heat solver's stencil update to the GPU with OpenACC. The reader who has no GPU should read it for the pattern; the numbers below are hand-computed and identical to the CPU solver of Chapter 24, which is exactly the point — a correct offload reproduces the serial physics to the last digit.

The whole trick is one !$acc data region wrapped around the entire time loop, keeping the field resident on the device across all steps so the plate crosses the host–device bridge once in and once out — never per-step. Inside, each step is two !$acc parallel loop collapse(2) kernels over the interior: one that computes the FTCS update from the old field into u_new, one that commits it back into u. The collapse(2) folds the $i$ and $j$ loops into a single large iteration space so the GPU has the full interior — thousands of independent cells on a real grid — to spread across its cores.

r = alpha*dt/dx**2                     ! = 0.2 <= 1/4, CFL-stable (Chapter 24)
u = 0.0_dp;  u(1,:) = 100.0_dp         ! hot top edge (Dirichlet); other edges 0
u_new = u                              ! u_new carries the same fixed boundaries

!$acc data copy(u) create(u_new)       ! field RESIDENT on the device across ALL steps
do step = 1, nsteps
  !$acc parallel loop collapse(2) present(u, u_new)     ! stencil: read OLD u, write 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
  !$acc parallel loop collapse(2) present(u, u_new)     ! commit new interior into u
  do j = 2, ny-1
    do i = 2, nx-1
      u(i,j) = u_new(i,j)
    end do
  end do
end do
!$acc end data                         ! u copied back to the host here, ONCE

Compile it with nvfortran -acc -Minfo=accel project-checkpoint.f90 -o heat_gpu; the -Minfo=accel feedback will confirm the two kernels were offloaded and — the line to look for — that u and u_new are copied only at the data-region boundary, not inside the step loop. On the same $5\times5$ plate as Chapter 24 ($\alpha=1$, $\Delta x = 1$, $\Delta t = 0.2$ so $r = 0.2$, top edge held at $100$, two steps), the field after step 2 is identical to the CPU solver's — the GPU computes the same forward-Euler update, cell for cell:

after step 2:
  100.00  100.00  100.00  100.00  100.00
    0.00   28.00   32.00   28.00    0.00
    0.00    4.00    4.00    4.00    0.00
    0.00    0.00    0.00    0.00    0.00
    0.00    0.00    0.00    0.00    0.00

The hand-trace is Chapter 24's, unchanged, because the arithmetic is unchanged: interior cell $(2,3)$ goes $20 \to 20 + 0.2(100 - 80 + 20 + 20) = 32$, and $(3,j) \to 0 + 0.2(20) = 4$. That the answer matches to the last digit is the correctness check every offload must pass: the accelerator's job is to compute the same result faster, never a different one. On this toy grid the GPU would of course be far slower than the CPU — five interior cells cannot fill thousands of cores, and the launch overhead dwarfs the work — which is §35.5's lesson made personal: this pattern pays only when the grid is large (say $2000\times2000$) and the steps are many, so the resident field is swept thousands of times and the one-time transfer vanishes into the noise. The full program, with the two-buffer setup and a show routine, is in code/project-checkpoint.f90.

Two honest caveats close the checkpoint. First, the interior copy that commits u_new into u each step is the simple, correct choice; a production code swaps the two buffers by pointer instead, avoiding even that on-device copy — a refinement, not a correctness issue. Second, to scale beyond one GPU you combine this offload with the distributed-memory decomposition of Chapter 34: one MPI rank per GPU, each rank offloading its own tile and exchanging halos between GPUs — the multi-GPU pattern the Chapter 38 capstone can present as its most ambitious configuration. Your solver has now, in principle, run on every model in the taxonomy: coarray images, OpenMP threads, MPI processes, and a graphics processor.

Summary

This chapter put the running project — and your Fortran — onto the graphics processor, and drew the line between offloads that win and offloads that embarrass.

The two things to memorize. First, the host/device split with separate memories, and its one commandment: move data across the bridge as rarely as possible — once in, once out, compute many times in between. A data region around the time loop, not a copy on each step, is the difference between a GPU win and a GPU slowdown. Second, that the GPU rewards a specific shape of computation — large, regular, data-parallel — and punishes everything else, so you decide whether to offload (§35.5, Amdahl, arithmetic intensity) before you decide how.

Spaced Review

Reaching back to the array foundations (Chapter 5) and the OpenMP directives (Chapter 33) this chapter builds on. Answer before expanding.

1. (Ch. 5.) The stencil offload sweeps u(i,j) with !$acc parallel loop collapse(2). Recalling column-major storage, which index varies fastest in memory, and why does that still matter on a GPU?
   

   Answer

   The **first** index `i` varies fastest — Fortran is column-major, so `u(i,j)` and `u(i+1,j)` are adjacent in memory. It matters on the GPU because **coalesced** memory access — consecutive threads touching consecutive addresses — is what lets the device stream memory at full bandwidth, the exact analogue of the cache-friendly CPU access of [Chapter 27](../../part-07-performance/chapter-27-why-fortran-is-fast/index.md). Mapping the fastest-varying dimension across neighbouring threads keeps the GPU's memory system fed; the same column-major awareness that made the CPU sweep fast makes the GPU sweep fast.

2. (Ch. 5.) Why can the interior update be offloaded as a parallel loop at all — what property of the two-array (current/next) finite-difference structure makes the cells independent?
   

   Answer

   Each new value `u_new(i,j)` is computed purely from the **old** field `u`'s neighbours, and written to a *separate* array, so no cell's update reads another cell's *new* value — the interior updates of one step are mutually independent (embarrassingly data-parallel), which is exactly the independence `!$acc parallel loop` requires. It is the same two-array structure from Chapters 5 and 24 that made the update safe to vectorize and to parallelize with OpenMP.

3. (Ch. 33.) OpenACC's !$acc parallel loop` is deliberately shaped like OpenMP's `!$omp parallel do. Name one thing the GPU forces you to manage that a shared-memory OpenMP program never had to, and why.
   

   Answer

   **Host–device data movement.** In OpenMP on a shared-memory node, all threads see the same RAM, so there is nothing to copy — the data is simply *there*. A GPU has its own separate memory, so OpenACC must explicitly move arrays to and from the device (the `copy`/`copyin`/`data` clauses), and managing those transfers is the central new concern the accelerator adds.

4. (Ch. 33.) In OpenMP you used a reduction clause to sum across threads safely. OpenACC has the same reduction(+:total) clause. In one sentence, why is a reduction something the directive must handle specially rather than leaving it as an ordinary parallel loop?
   

   Answer

   Because every thread would otherwise be writing the *same* accumulator simultaneously — a data race that corrupts the result; the `reduction` clause gives each thread a private partial result and safely combines them at the end, on the GPU across thousands of threads (the subtlety Case Study 2 develops).

5. (Ch. 33 + this chapter.) A shared-memory OpenMP run and a GPU OpenACC run of the same 98%-parallel stencil both obey Amdahl's Law. What extra term enters the serial fraction for the GPU version?
   

   Answer

   The **host–device transfer time** — the cost of moving the field to the device and results back, which OpenMP (shared memory) never pays. It is work the GPU *adds*, so it counts toward the un-accelerated serial fraction and lowers the whole-program speedup ceiling; minimizing it (one data region, not per-step copies) is how you keep that extra serial term small.

What's Next

That closes Part VIII, and with it the parallel journey your solver began five chapters ago. It has now run — at least in principle — across coarray images, OpenMP threads, MPI processes, and a graphics processor, each model suited to a different machine and a different question. You have the full toolkit of high-performance Fortran: fast serial code (Part VII) made parallel four ways (Part VIII), all behind the same step interface you froze back in Chapter 6.

But a pile of fast, parallel kernels is not yet a scientific code. Real simulations are large, long-lived, multi-author software, and organizing them — so that a hundred thousand lines stay navigable, testable, and trustworthy — is its own discipline. Part IX turns from making code fast to making it real. Chapter 36 opens it by touring a genuine large scientific Fortran code — its directory layout, its module hierarchy, how to navigate a codebase far too big to hold in your head — and shows you how the driver/solver/physics/I-O structure you have been quietly building all along is exactly how the professionals organize a simulation. The performance is done; now we make it software.