In the three parts of this exercise, we use the example of matrix multiplication to demonstrate three different strategies for migrating code to the GPU. The code that we obtain in each case is portable, because it is based on a platform-independent standard, and it does not depend on any single vendor's proprietary interface. Here are the strategies we'll explore:

• OpenMP offload, which is a cross-language standard that we illustrate in C, though it also applies to C++ and Fortran.
• Parallel execution policies ("stdpar") in C++, which make use of specific features of the C++17 language standard.
• DO CONCURRENT loops in Fortran, which likewise depend on specific features of the Fortran 2018 language standard.

Comparison of portability strategies in the 3 codes

All of these methods have the advantage that the same source code can be compiled so that it will run on the host's CPUs with no GPU offloading, when the compiler is given different options.

Why choose matrix multiplication, with N = 16384?

Matrix-matrix multiplication is a standard operation in dense linear algebra that offers opportunities for parallelization. The codes create and multiply two N by N matrices, where the constant N is set to 16384. At this size there are ~268 million (i, j) pairs and approximately 8.8 TFLOP of work. A straightforward implementation generates 268 million concurrent threads on NVIDIA GPUs, ensuring that every SM stays saturated on recent NVIDIA models: the Blackwell B200 has 26,624 FP32 CUDA cores across 208 SMs, while the Hopper H100 has 16,896 FP32 CUDA cores across 132 SMs. Reducing N significantly would leave cores idle at the start and end of the kernel, and it might cause the codes to run too quickly to time accurately.

C / OpenMP

OpenMP's target directives may be used together with the usual OpenMP directives to parallelize ordinary for loops, thereby preserving the typical 2D loop structure used in C-style matrix multiplication. For programmers coming from a C or Fortran background this style is more natural and readable. Optimized compilations for both GPU and CPU can proceed with no difficulty.

Portability

OpenMP target offload is supported by NVIDIA NVHPC compilers (nvc, nvc++, nvfortran) on NVIDIA platforms, as well as GCC and Clang/LLVM on the majority of platforms. With GCC, you have the ability to choose different GPU targets with -foffload=. Thus, the same source code can target AMD GPUs by changing the compiler target from nvptx-none to amdgcn-amdhsa and recompiling. Clang has analgous options, but they can vary depending on vendor customizations to Clang. The flag in basic Clang (from LLVM Project) is -fopemmp-targets=<triple>. NVHPC compilers of course assume the target is NVIDIA when given -mp=gpu, and specific NVIDIA architectures can be selected with the -gpu= flag.

Example Code

Open and download the file: matmul_omp.c. Compile it with any of the following commands, then run the executable on the host. Certain parts of the code will be offloaded to the NVIDIA GPU target:

Syntax Notes

OpenMP loop directives, target teams distribute parallel for collapse(2)

The full directive that controls the nested for loops is the OpenMP equivalent of a CUDA kernel launch:

• teams distribute partitions the iteration space across the grid of thread blocks (teams);
• parallel for distributes iterations across threads within each block;
• collapse(2) merges the i and j loops into a single 67-million-iteration range before distributing, maximizing the parallelism visible to the runtime.

OpenMP data directives, target data map(alloc:)

Wrapping both kernels inside a single #pragma omp target data region keeps the arrays resident on the GPU between the two kernel launches. Without this, each target directive would independently map and unmap the arrays, transferring data unnecessarily at the boundary between the two kernels.

In OpenMP, the programmer may explicitly map memory through the map clauses (as is described in more detail in the Data Management topic of OpenMP Offload to GPUs):

Using map(alloc:) for all three arrays avoids unnecessary data transfers between the host and the device. Matrices A and B can exist entirely on the device, where they are written by the initialization kernel and consumed by the multiply kernel. Likewise, values of the product matrix C do not need to be transferred back to the host, as they do not matter for the purposes of this exercise.

C++ / stdpar

C++ Standard Parallelism (often called stdpar) enables developers to express parallelism directly through the C++ Standard Library, primarily using parallel algorithms introduced in C++17. By specifying an execution policy such as std::execution::par or std::execution::par_unseq, programmers indicate that operations like sorting, transforming, reducing, or searching may be executed concurrently across multiple CPU cores or on accelerator hardware. The key advantage of stdpar is that it allows existing sequential code to be parallelized with minimal changes, without requiring explicit thread management, synchronization primitives, or device-specific programming models.

Portability

Modern compiler implementations, such as NVIDIA's nvc++, can offload stdpar algorithms to GPUs, providing performance acceleration while preserving standard, portable C++ source code. As a result, stdpar provides a high-level, standards-based approach to parallel programming that improves code maintainability and portability while enabling scalable performance on contemporary computing systems.

Why stdpar does not suit all use cases

The standard parallel algorithms express only a limited set of patterns (map, reduce, scan, sort, etc.). Work that requires explicit synchronization between threads, shared memory tiling, or fine-grained control over a thread hierarchy cannot be expressed with execution policies and requires either OpenMP target directives or direct CUDA/HIP.

Example Code

Open and download the file: matmul_stdpar.cpp. Compile it so that the stdpar sections are offloaded to the NVIDIA GPU, then run the executable:

Syntax Notes

Why par_unseq and not par

The par_unseq policy means parallel and unsequenced. It implies that iterations may run in any order, concurrently, and with SIMD. This provides the necessary conditions for a given set of iterations to be expressed as a GPU kernel by the compiler. In particular, if a code is compiled with nvc++ -stdpar=gpu, then a code block enclosed by a loop like the following becomes a CUDA kernel:

In contrast, std::execution::par enables parallel execution but not vectorization within a thread. For this reason, par_unseq is the required policy for triggering CUDA code generation with nvc++.

Why a flat index vector rather than a 2-D range

C++ parallel algorithms operate over a single iterator range. There is no standard way to express a 2D iterator range, so the (i, j) index space is flattened to a 1D integer range, and the row/column are recovered inside the lambda with integer division and modulo:

Why data pointers are captured rather than vectors

A std::vector object is not trivially copyable—it holds a pointer, a size, and a capacity. When nvc++ turns a par_unseq call into a CUDA kernel, the lambda's capture list becomes the kernel's parameter block, which must consist entirely of trivially copyable values. Capturing a std::vector by reference across a kernel boundary is undefined behavior. Extracting raw pointers with .data() and capturing those by value is the correct pattern.

Blocking behavior and timing accuracy

Unlike CUDA kernel calls, stdpar algorithm calls are blocking from the host's perspective—they do not return until the GPU has finished. This means a std::chrono::high_resolution_clock measurement placed immediately before and after the std::for_each call accurately captures GPU wall-clock time, without requiring the equivalent of cudaDeviceSynchronize.

C++ functor as an alternative to lambda

A lambda in C++ is just a convenient, inline way of defining a functor, i.e., a named struct or class that overloads operator(). In other words, a lambda is simply "syntactic sugar"; it is completely equivalent to write the lambda's callable as a functor and pass it to std::for_each that way. The functor approach can be preferable when the kernel is large enough to warrant a named type, or when the callable needs to be reused in multiple call sites.

If you want to try the alternative of calling a functor, download matmul_cpp_functor.h and follow the instructions in matmul_stdpar.cpp. Then compile and run the code as before. In principle, the generated code should be identical to the lambda version; in practice, nvc++ appears to treat the two cases differently, and the functor version runs faster.

Fortran / DO CONCURRENT

DO CONCURRENT is a loop construct introduced in Fortran 2008/2018 to express iterations that are logically independent and can be executed in any order. Unlike a traditional DO loop, DO CONCURRENT tells the compiler that there are no data dependencies between iterations, allowing it to apply optimizations such as vectorization, threading, or execution on accelerators. Often it is paired with the LOCAL clause (from the Fortran 2018 standard, Section 11.1.7.5) to identify the variables that are local to each thread. The combined construct enables programmers to write clear, high-level, parallelizable code while preserving the mathematical intent of the algorithm.

Portability

Because parallelism is expressed as part of the Fortran language standard, rather than through platform-specific directives or APIs, the same source code can be compiled and run on a wide range of systems, including multicore CPUs, clusters, and GPUs. Different compilers can choose the best execution strategy for the target hardware, allowing applications to benefit from new architectures without requiring significant code changes.

Example Code

Open and download the file: matmul_concurrent.f90. Compile it for offloading to an NVIDIA GPU, and run it as shown.

The key construct is the main DO CONCURRENT loop that performs the matrix multiplication:

Syntax Notes

Swapping index order in DO CONCURRENT

You may have noticed that the loop indices are ordered (j=1:N, i=1:N) rather than (i=1:N, j=1:N) as they were in C/C++. The reason for this is performance. The nvfortran compiler maps the last declared index in DO CONCURRENT to the fast CUDA thread dimension (threadIdx.x). With i last, consecutive threads get consecutive i values. Since Fortran stores arrays in column-major order (i varies fastest in memory), A(i,k) and C(i,j) become coalesced accesses. Without the swap, the stride for threads in a warp would be N=16384 floats when accessing B and C. That means each thread would hit a different cache line, which is the worst possible access pattern on a GPU.

The locality clause: LOCAL(k, acc)

The LOCAL(k, acc) clause tells the compiler that k and acc are private to each iteration. This is what makes the loop parallelizable; otherwise, there would be a race condition in which all the threads attempt to write to the same variables. The k and acc variables must be declared in the host scope prior to the DO CONCURRENT construct, because the LOCAL clause by itself can only reference existing variable names; it does not declare new ones.

Inner k-loop inside DO CONCURRENT

The regular DO loop nested inside DO CONCURRENT is permitted by the standard. The outer (i, j) pairs are parallel; the inner k loop is sequential within each thread. The compiler maps the outer index space to the GPU thread grid and compiles the inner loop as ordinary sequential code inside each thread. If the same loop were compiled for a CPU, the k loop would be vectorized.

Initialization as a second DO CONCURRENT

The initialization loop is also a DO CONCURRENT loop. With -stdpar=gpu, the compiler allocates all arrays in CUDA unified memory and runs the initialization kernel on the device. This means the data is already resident on the GPU before the multiply kernel launches—no host-to-device transfer occurs between the two kernels.

What nvfort does under the hood

1. Allocates A, B, C in CUDA unified memory (because -stdpar=gpu uses cudaMallocManaged by default).
2. Converts the initialization DO CONCURRENT to a kernel that writes to device memory.
3. Expresses the multiply DO CONCURRENT as a second CUDA kernel; no explicit host-to-device transfers are needed because the data is already on the GPU from step 2.
4. Inserts a cudaDeviceSynchronize before system_clock so the timing is accurate.

This exercise was constructed with the assitance of Claude Code and Microsoft Copilot.