Evaluating OpenMP Offloading for Intra-node Multi-GPU Programming across NVIDIA, AMD, and Intel Architectures: A 3D Heat Transfer Case Study

Evaluating OpenMP Offloading for Intra-node Multi-GPU Programming across NVIDIA, AMD, and Intel Architectures: A 3D Heat Transfer Case Study Ezhilmathi Krishnasamy 1

∗

1,2,3,∗

Rudolfovo Science and Technology, Novo Mesto, Slovenia 2 University of Ljubljana, Ljubljana, Slovenia 3 University of Luxembourg, Belval, Luxembourg

Corresponding author: ezhilmathi.krishnasamy@rudolfovo.eu

arXiv:2608.11882v1 [cs.DC] 12 Aug 2026

Abstract Currently, most supercomputers are equipped with GPUs from manufacturers such as NVIDIA, AMD, or Intel, which provide substantial parallelism and high throughput. It is common for a single compute node (intra-node) to host multiple GPUs, typically four or more. Therefore, effectively leveraging all these GPUs within a single compute node is essential for applications in scientific and engineering domains. However, several factors must be considered before utilizing these GPUs for scientific computing, including the implementation of data communication, the programming models available for use across these GPUs, and the level of performance that can be achieved with a single codebase across different GPU architectures and configurations within a single compute node. OpenMP Offloading is a prominent directive-based programming model that can be executed on all three GPU types: NVIDIA, AMD, and Intel. In this research, we present an analysis of the benefits and performance challenges of using OpenMP Offloading to address the 3D heat equations, which involve both primary computation, as well as halo computation and communication. We investigate how performance varies in relation to native GPU programming models—CUDA for NVIDIA, HIP for AMD, and SYCL for Intel. Furthermore, we demonstrate that OpenMP Offloading can achieve performance improvements of approximately 2x for 2 GPUs and around 4x for 4 GPUs when compared to single-GPU OpenMP Offloading implementations across all three GPU types. This analysis is conducted systematically through various OpenMP Offloading implementations that utilize different low-level APIs for memory allocation, memory transfer options (synchronous, asynchronous, and peer-to-peer), and other native GPU programming models such as CUDA (NVIDIA), HIP (AMD), and SYCL (Intel).

Keywords: Multi-GPU; OpenMP Offloading; CUDA; HIP; SYCL; 3D Heat Equation; Stencil Computation; AMD; NVIDIA; Intel GPU

1

Introduction

In recent decades, it has become increasingly common for supercomputers and small clusters equipped with Graphics Processing Units (GPUs) to accelerate computational applications. Typically, GPUs are employed to enhance the speed of these applications, a process often referred to as Offloading computation to GPUs. However, the diversity of hardware and the range of programming options available can present developers with numerous challenges. While this variety offers benefits, it can also necessitate significant adjustments or even complete rewrites of the original code to accommodate different GPU architectures using languages such as CUDA [1], SYCL [2], and OpenCL [3]. For instance, code designed for execution on an NVIDIA GPU using CUDA cannot be directly run on an AMD GPU; however, it can be relatively easily converted into a HIP [4] version for execution. Nevertheless, questions continue to arise about whether similar performance levels can be achieved on AMD GPUs as on NVIDIA GPUs, primarily due to the architectural differences between the two. 1

For example, in the CUDA framework, threads are organized into groups called warps, consisting of 32 threads, whereas in AMD’s HIP, 64 threads are grouped into wavefronts. Consequently, when converting code from CUDA to HIP, one must manually adjust the thread blocks to align with their respective architectures to achieve optimal performance. Moreover, CUDA programming requires a deep understanding of thread organization and other CUDA-specific syntax. Should we consider Intel GPUs, utilizing their architecture necessitates either SYCL or OpenMP Offloading. Thus, it is essential to explore alternative programming models that can be easily executed across NVIDIA, AMD, and Intel GPUs. Fortunately, such programming models do indeed exist. Low-level programming models, like SYCL and OpenCL, as well as directive-based programming models such as OpenMP Offloading, can be employed on these GPUs. However, both SYCL and OpenCL require extensive syntax learning and a deeper understanding of thread management for GPU kernels. In contrast, OpenMP Offloading integrates seamlessly into existing code with minimal changes, facilitating execution across all these GPUs. Therefore, in this research, we investigate the OpenMP Offloading [5, 6] programming model on intra-node (single compute node) GPUs, using the 3D steady-state heat equation as a test case. The 3D heat equation demands intensive kernel computation and halo exchanges, which are typical in many scientific computing routines.

1.1

Research Objectives

Directive-based programming models have existed even before GPUs became a significant component of scientific computing. When GPUs began to play a major role in this field, there was an initial surge in demand for learning low-level programming models such as CUDA and OpenCL. However, over time, other programming models have emerged to support GPU programming, including OpenACC [5, 7], OpenMP Offloading, and SYCL. While SYCL, CUDA, and OpenCL require a comprehensive understanding of thread organization to execute kernels or primary computations on the GPU, OpenACC and OpenMP Offloading simplify this process through their directive-based programming models. These models extend OpenMP syntax, providing a more straightforward programming approach and reducing complexity. Generally, GPU programming involves mapping parallel threads across the available GPU cores, based on the Single Instruction Multiple Threads (SIMT) [8] architecture. CUDA, HIP, and SYCL effectively utilize this concept by employing threads and thread blocks, which include warps and 1D, 2D, and 3D thread block configurations. In contrast, directive-based programming models streamline this process with simplified syntax, significantly diminishing the amount of code required compared to programming models such as CUDA, OpenCL, and SYCL. Although OpenACC is a directive-based programming model, it unfortunately suffers from limited compiler support for AMD and Intel GPUs. It can accommodate FORTRAN programming models on AMD but has limited support for C/C++ [9, 10]. Therefore, this research aims to examine the feasibility of OpenMP Offloading across various well-established GPUs available in the market, including NVIDIA, AMD, and Intel GPUs within a single compute node. We will evaluate whether there are any performance variations among different GPUs when using the same source code with OpenMP Offloading. To this end, we focus on a key application: 3D heat transfer, which is commonly encountered in compute-intensive applications. In multi-GPU scenarios, this would necessitate complete computations, along with main and halo computations, as well as halo data exchange.

1.2

Contributions

This research makes several significant contributions to the field, which are outlined below: • A comprehensive analysis of single-node multi-GPU OpenMP Offloading implementations, specifically focusing on three-dimensional heat equation solvers. • A detailed examination of data transfer methodologies between GPUs within a single compute node, utilizing benchmark cases proposed for evaluation.

2

• The development of a performance model designed to analyze single-node multi-GPU systems across various hardware platforms, including NVIDIA, AMD, and Intel architectures. • The provision of source code and research artifacts to enhance the reproducibility of the findings; detailed instructions for accessing and reproducing these artifacts are provided in Artifact Description Appendix A.

1.3

Organization of the Paper

This paper is organized as follows. Section 2 reviews several significant research contributions related to OpenMP Offloading. In Section 4, we present the mathematical details of the 3D heat equation, which serves as the use case for this research. Section 5 offers a comprehensive overview of the hardware utilized, including details about the compiler and its relevant flags, ensuring reproducibility of our results. Section 6 provides an in-depth discussion of key proposals that enhance the feasibility of OpenMP Offloading across various GPUs in multi-GPU contexts. This section includes illustrative code snippets demonstrating different OpenMP Offloading techniques, along with the derivation of performance models. Finally, Section 7 delivers a thorough analysis of multiple implementations of OpenMP Offloading, assessing their feasibility across various GPU architectures, including those from NVIDIA, AMD, and Intel. This section also features an extensive discussion of the performance models. The paper concludes with a summary of the principal findings derived from this research.

2

Related Work

Numerous studies have investigated OpenMP Offloading, focusing on single GPU configurations, multiple GPUs within a single compute node (intra-node), and multiple nodes (inter-node). This section reviews a selection of the most pertinent research in this area. E. Krishnasamy et al. [11] examined the single GPU OpenMP Offloading option for parallelizing the Tsunami scientific code, originally written in FORTRAN. Their findings indicated that implementing OpenMP Offloading results in a significant speedup compared to traditional CPU-based multi-core OpenMP implementations when utilizing an NVIDIA GPU. Fridman et al. [12] explored the application of OpenMP Offloading in conjunction with OpenMP Validation and Verification (OMPVV) as a test case on both NVIDIA and Intel GPUs, specifically focusing on single GPU scenarios. Furthermore, E. Krishnasamy et al. [10] assessed OpenMP Offloading across various GPUs in the context of matrix-vector multiplication. This research utilized low-level Application Programming Interfaces (APIs) for memory allocation and synchronized memory transfers, considering both GPU-aware and non-GPU-aware methodologies. Additional research has validated the advantages of utilizing low-level APIs for memory transfer and computational kernel optimizations, with significant investigations conducted by E. Krishnasamy et al. [13, 14] on NVIDIA and AMD GPUs. Moreover, Ferat et al. [15] evaluated hidden helper thread configurations within heterogeneous MPI+OpenMP applications, while Tian et al. [16] investigated the concurrent execution of deferred OpenMP target tasks. However, no research has systematically explored the utilization of multiple GPUs within a single compute node using OpenMP Offloading—especially utilizing low-level APIs—across three distinct platforms: NVIDIA, AMD, and Intel. This research aims to address this gap by investigating various approaches related to low-level API memory management, diverse parallelization strategies, performance models, data transfer between the CPU and GPU, as well as inter-GPU communication, all within the context of typical scientific applications such as 3D heat transfer.

3

GPU Programming Models

GPU programming models facilitate significant computational speedups by leveraging the capabilities of GPUs. These models can be broadly categorized into two types: low-level programming models and high-level programming models. Notable examples of low-level programming models include CUDA, HIP, OpenCL, and SYCL. In contrast, OpenACC and OpenMP Offloading are classified as high-level 3

or directive-based programming models. Both OpenACC and OpenMP Offloading support hybrid programming approaches, which enable effective targeting of both CPUs and GPUs. Like other programming models, OpenACC and OpenMP Offloading continuously improve their APIs with each release, introducing new features related to compute kernels, memory allocation, and data transfer [5–7]. These models provide three types of memory transfer between CPUs and GPUs: Structured High-Level API, Unstructured High-Level API, and Low-Level API for Data Movement [13, 14, 17]. • Structured High-Level API: This API is user-friendly, though it may be insufficient for larger computations. • Unstructured High-Level API: While this API is also easy to use, its performance has some limitations. • Low-Level API for Data Movement: This API delivers optimal performance, as it is closely aligned with low-level GPU programming models. Data movement remains one of the primary bottlenecks in scientific computing. Computational workloads are often classified as either compute-bound or memory-bound. Compute-bound problems are primarily limited by the available computational throughput of the hardware and the arithmetic requirements of the task. In contrast, memory-bound problems are restricted mainly by the speed at which data can be transferred through the memory hierarchy. Their performance depends on various factors, including the nature of the problem, memory access patterns, hardware bandwidth and latency, cache utilization, and the efficiency of the software interfaces and implementations used for data transfer. Previous research by E. Krishnasamy et al. [13, 14, 17] demonstrated that utilizing low-level APIs for memory allocation and testing various combinations of optimized compute kernel options revealed that the default parallel constructs and clauses offer the best computational efficiency for Basic Linear Algebra Subroutines (BLAS). Consequently, this work will primarily focus on low-level APIs for memory allocation alongside default parallel constructs and clauses for compute kernels.

4

Mathematical Details

To test our intra-node multi-GPU parallelization using OpenMP Offloading, we aimed to select an application that demands an intensive computational kernel and involves communication between the CPU and GPU, as well as among GPUs, such as in halo exchanges. A prime example of this need is the requirement for halo exchanges during each iteration across a 3D domain, a common scenario in scientific computing simulations. Given these considerations, we chose the 3D heat equation. Below, we present the mathematical details of the 3D model that we considered in this work: ∂u = ∇ · (κ(x, y, z)∇u) + f (x, y, z), (1) ∂t where u(x, y, z, t) is the scalar field and κ(x, y, z) is the spatially varying diffusion coefficient. In this simulation, the spatially varying coefficient reads κ(x, y, z) = 0.5 + 0.45 sin(πx) sin(πy) sin(πz),

(2)

and the source term is defined by f (x, y, z) = 1 + cos(πx) cos(πy) cos(πz).

(3)

The equation is solved on a uniform grid with spacing h and time step ∆t. Let um i,j,k denote the numerical solution at time level tm . Using explicit time stepping, the update scheme for interior grid points is

4

∆t (Dx + Dy + Dz ) , (4) 2h2 Boundary condition values are fixed to zero (Dirichlet boundary conditions), and the initial condition m um+1 i,j,k = ui,j,k + ∆t fi,j,k +

is u(x, y, z, 0) = sin(πx) sin(πy) sin(πz).

(5)

The solution is advanced in time using an explicit time-stepping scheme until the final simulation time T . For numerical stability, the time step must satisfy the CFL-type condition ∆t < 0.5. h2 Pseudocode 1: OpenMP parallel explicit solver for the 3D variable-coefficient heat equation. Input: Grid size n and final time T Output: Numerical solution u 1 h ← 1/(n − 1), κmax ← 0.95, ∆t ← h2 /(12κmax ), β ← ∆t/(2h2 ) // Initialize the solution and coefficient arrays 2 #pragma omp parallel for collapse(3) 3 foreach grid node (i, j, k) do 4 fi,j,k ← ∆t [1 + cos(πxi ) cos(πyj ) cos(πzk )] 5 κi,j,k ← 0.5 + 0.45 sin(πxi ) sin(πyj ) sin(πzk ) 6 if (i, j, k) is a boundary node then 7 ui,j,k ← 0 8 else 9 ui,j,k ← sin(πxi ) sin(πyj ) sin(πzk ) 10 end 11 end 12 t ← 0 13 #pragma

omp parallel

14 { 15

16 17 18 19

20 21 22 23 24 25

while t < T do // Distribute the stencil update among the existing threads #pragma omp for collapse(3) foreach interior node (i, j, k) do unew i,j,k ← ui,j,k + fi,j,k + βL(u)i,j,k end // One thread advances the shared time level #pragma omp single { swap(u, unew ) t ← t + ∆t } end

26 } 27 return u

5

(6)

Table 1: Comparison of technical specifications of H100, MI250X, and Intel Data Center GPU Max 1550 accelerators. Property

H100 [18, 19]

Compute & memory (per accelerator package) Compute units (vendor 132 SMs terminology) Compute dies/tiles 1 Peak FP64 vector (GFLOP/s) 34,000 Memory bandwidth (GB/s) 3,350

MI250X [20, 21]

Intel Max 1550 [22, 23]

220 CUs

128 Xe-cores

2 GCDs 47,900 3,276.8

2 tiles ∼52,000 3,276.8

GPU–GPU interconnect Interconnect technology Links/ports per accelerator pair Total P2P links/ports per accelerator

NVLink 4 NV6 18

Infinity Fabric 2 IF links 5 P2P links

Xe Link 4 Xe Link ports 12 active ports

Host interconnect (GPU–CPU) Host link Host-link bandwidth, 1 dir. (GB/s)

PCIe Gen5 ×16 ∼64

Infinity Fabric ∼36

PCIe Gen5 ×16 ∼64

halo exchange

interior plane

GPU3

(±z between GPUs)

halo plane

ghost plane increasing memory address GPU2

ghost

halo

recv. halo send to GPU−1 from GPU−1 ghost

GPU1

z(k)

halo

ghost

recv. halo send to GPU+1 from GPU+1 ghost

Inside one halo plane (fixed k): the (i, j) cells are contiguous flatten

j y(j)

interior points

k = 0 … N−1

GPU0

Nᵢ·Nⱼ contiguous cells (i fastest)

i x(i)

Figure 1: Schematic overview of the data partitioning in a 3D domain, illustrating halo data transfer during halo and interior-point computations.

6

Infinity Fabric - 2 links

4 x IF = 400 GB/s/ bi-dir.

GPU0

GPU3

NVLink 4: 25 GB/s/dir per link (per direction) with NV6

GPU1

tile 0/1

ic

br

a yF

(IF

):

1

/s GB

ir.

i-d

/b

00

it fin

In

tile 3/0

tile 2/0 GPU2

GPU3

tile 3/1

tile 1/0

MDF

tile 0/1

tile 1/1

tile 3/0 GPU2

~180 GB/s unidirec. meas.

2 × XL8 routes per physicalGPU pair = 4 ports = 16 lanes

MDF tile 1/1

MDF tile 3/1

tile 2/1

53.125 Gbit/s per lane 106.25 GB/s/dir (212.5 GB/s bidir.)

tile 2/0 12 active Xe Link ports

MDF

GPU2

G1 6×25 = 300 GB/s/bi-dir

1 * XL8 route = 2 ports × four-lane Xe Link ports = 8 lanes

tile 0/0

tile 1/0

GPU1

GPU1 G0

GCD (Graphics Compute Die)

Σ 18 links = 450 GB/s/dir

tile 0/0

GPU3

G1 6×25 = 150 GB/s/dir

2 X IF: 200 GB/s/bi-dir.

G0 GPU0

PCIe Gen5 X 16 (Per GPU) x16 lanes = 64 GB/s each way

NUMA X :36 + 36 GB/s/bi-dir. to each GPUX

GPU0

PCIe Gen5 X 16 (Per GPU) x16 lanes = 64 GB/s each way

tile 2/1

MDF: Connected with Multi-Die Fabric Interface

(a) H100 (BSC).

(b) AMD MI250X (LUMI-G).

(c) Intel Max 1550 (SuperMUCNG).

Figure 2: Comprehensive comparison of the three GPUs, including theoretical memory bandwidth, CPU–GPU connectivity, and inter-GPU communication.

5

Testing Platforms and Experimental Setups

In this research, we examined three categories of GPUs: the NVIDIA H100 (BSC, Spain), the AMD MI250X (LUMI-G, Finland), and the Intel Data Center 1550 (SuperMUC-NG, Germany). Each GPU features a distinct architecture that provides varying levels of computational performance and memory bandwidth, as well as end-to-end network connectivity between the CPU and GPU, and among the GPUs on the NUMA node. Figure 2 illustrates the data transfer network connectivity and theoretical bandwidths of these GPUs. Key specifications for each of the GPUs are presented in Table 1. For details regarding the compiler suite utilized for various tasks, please refer to Table 2. To evaluate the efficacy of our proposed methodologies, we implemented the 3D heat equation along with its multi-GPU parallelization strategies, which are described in Sections 4 and 6. This implementation was used to test both kernel and halo computations, as well as the necessary halo communication. We maintained a consistent time step of 500 while exploring various 3D domain sizes, ranging from 5123 to 12803 , with increments of 128. This parameterization enabled us to assess the runtime across the different GPU programming models and their respective implementations examined in this research.

6

Methodology

The primary objective of this research is to investigate the feasibility of utilizing OpenMP Offloading as a unified source code across the three GPUs examined in this research. Specifically, we aim to assess the performance enhancements achieved when employing one, two, and four GPUs to solve scientific computing problems, particularly in the context of three-dimensional (3D) heat transfer. Additionally, we seek to develop a performance model that can predict the runtime requirements for the current problem as well as for similar computational challenges. To address these objectives, we propose the following key approaches. First, we will implement various versions of multi-GPU OpenMP Offloading to solve the heat equations on a single compute node (intra-node), which will allow us to identify the most effective optimization strategies. Second, we will explore additional optimization opportunities provided by OpenMP, which can enhance asynchronous operations related to multiple kernels and data transfer tasks when utilizing multi-GPU configurations within a single compute node. Third, we will present detailed benchmark cases for data transfer that will enable us to evaluate performance under realistic conditions. Finally, we aim to develop a performance model that considers various data transfer scenarios and diverse OpenMP multi-GPU implementations. This model will assist in predicting runtimes for 3D heat transfer and similar computational problems.

7

Table 2: Hardware and software environment per platform and programming models. Platform

NVIDIA H100

AMD MI250X

GPU/arch

Compiler

Toolkit Stack & Key Flags

NVHPC 25.3 (nvcc)

CUDA: NVHPC/25.3-CUDA-12.8.0; -O2 -arch=sm_90 --fmad=false --prec-div=true --prec-sqrt=true --ftz=false

Clang 18.1.8

OpenMP: clang/18.1.8-cuda12.8, CUDA 12.8, GCC 13.2.0; -O3 -fopenmp --offload-arch=sm_90 -ffp-contract=off

AMDClang (ROCm 6.4.4)

HIP & OpenMP: LUMI/25.09, PrgEnv-amd, craype-accel-amd-gfx90a, ROCm 6.4.4; OpenMP -O3 -fopenmp -ffp-contract=off (via cc)

Intel oneAPI icpx (SYCL)

SYCL: intel-toolkit (+intel-dpct); -fsycl -fiopenmp -ffp-contract=off; SPIR-V JIT; ZE_FLAT_DEVICE_HIERARCHY=COMPOSITE

Intel oneAPI icx (OpenMP)

OpenMP: intel-toolkit; -fiopenmp -fopenmp-targets=spir64 -ffp-contract=off; COMPOSITE

Hopper, sm_90

CDNA2, gfx90a (2 GCDs/card)

Intel Max 1550 Xe-HPC, 2 stacks

6.1

Multi-GPU OpenMP Offloading Versions

We initialize the 3D computational domain in the order of k, j, i and subsequently partition the domain along the k direction. Algorithm 1 presents an OpenMP implementation of the 3D heat equations, where we delineate the primary computational segment designated for offloading to the GPUs. Figure 1 illustrates the discretization of the 3D domain and clarifies how such a partitioning can enhance computational efficiency in a multi-GPU environment. It is crucial to note that omitting the split at the final index prevents the acquisition of flattened halo data. The core principle of this approach is to facilitate efficient halo computations while simply transferring contiguous data values for subsequent iterations. To achieve efficient OpenMP Offloading, we draw from previous implementations of OpenMP Offloading [13, 14], which include low-level APIs specifically designed for memory allocation and data transfer between the Host and Device (HtoD) and Device and Host (DtoH). Furthermore, OpenMP Offloading provides options for synchronous, asynchronous, and Peer-to-Peer (P2P) data transfer, as detailed in Table 4. To evaluate the performance portability of OpenMP Offloading, we developed nine distinct variants: one version for single GPU configurations, four versions tailored for dual GPU setups, and four versions designed for quad GPU configurations: • OpenMP Offloading 1-GPU Baseline: A baseline implementation utilizing a single GPU. • OpenMP Offloading {2,4}-GPU Version 1: A basic multi-GPU implementation. – A single host thread oversees all program calls, which encompass the execution of interior point computational kernels, halo computations, memory calls, and data synchronization associated with the copy APIs. • OpenMP Offloading {2,4}-GPU Version 2: An asynchronous operation that employs streams, directed by a single host thread. – This version uses a single host thread along with a streaming concept for asynchronous operations related to both interior and halo computations, as well as for data transfer. • OpenMP Offloading {2,4}-GPU Version 3: Utilization of multiple host threads to manage the GPUs in conjunction with streams, where each thread is responsible for an individual GPU. 8

– Multiple host threads each govern a GPU, managing its respective kernel execution and data transfer while employing asynchronous data copying APIs. • OpenMP Offloading {2,4}-GPU Version 4: Facilitating P2P data transfer between GPUs using streams, managed by multiple host threads. – A scheme in which multiple host threads supervise the GPUs, with each thread responsible for kernel execution and data transfer of a single GPU, utilizing asynchronous data copying with P2P functionality APIs. To facilitate a comparative analysis, we also developed equivalent implementations using CUDA, HIP, and SYCL, which will be benchmarked against the OpenMP Offloading versions. While one could theorize that the P2P implementation may yield superior performance, it is essential to evaluate it against the other three versions (i.e., Version 1: synchronous, Version 2: asynchronous, and Version 3: multiple host threads with asynchronous execution). To ensure a thorough comparison, we have also implemented three additional variants alongside Version 4 (P2P). Table 3 provides a straightforward interpretation of these versions. Table 3: Summary of the comparison of the 2-GPU and 4-GPU parallel strategy implementations of OpenMP Offloading {2,4}-GPU Versions. Version OpenMP Offloading {2,4}-GPU Version 1 OpenMP Offloading {2,4}-GPU Version 2 OpenMP Offloading {2,4}-GPU Version 3 OpenMP Offloading {2,4}-GPU Version 4

6.1.1

Host Threads

Halo-Transfer Path

Overlap with Computation

1 1 2, 4 2, 4

GPU → host → GPU GPU → host → GPU GPU → host → GPU GPU → GPU

None Limited Yes, via host Yes, via P2P

OpenMP Offloading 1-GPU Baseline

This implementation does not incorporate data partitioning, which eliminates the need for halo computations and halo data communication. OpenMP Offloading allows for comparable control over thread block configuration in targeted regions through the ‘num_teams‘ and ‘thread_limit‘ clauses, similar to CUDA’s grid and block dimensions (refer to Table 4 in [14]). Furthermore, OpenMP Offloading offers limited combinations of compute kernels with parallel constructs and clauses that, in theory, should enhance performance. This interplay has been examined in previous work by E. Krishnasamy [13, 14], as shown in Table 7 of [13]. Notably, our investigations into various configurations for compute kernels revealed that using the default compute constructs and clauses resulted in superior performance compared to manually adjusting the number of threads and block configurations. This finding aligns with earlier studies conducted by E. Krishnasamy [13, 14]. Consequently, this research refrains from engaging in manual thread and grid block configuration. For compute kernel directives, we employed the ‘collapse(3)‘ directive to facilitate interior-point computation on the 3D dataset. For data transfer between the CPU and GPU, we utilized the low-level API, which has been shown to be more efficient than high-level APIs for data transfers [13, 14]. Detailed specifications for both main computation and kernel computation are clearly outlined in Pseudocode 2. 6.1.2

OpenMP Offloading {2,4}-GPU Version 1

In this approach, we partitioned the data blocks when utilizing multiple GPUs. For the case of two GPUs, we implement two distinct kernels: one for interior computations and another for halo computations. In the four-GPU scenario, we use a single kernel for interior computations and two separate kernels for halo computations, designated for the top and bottom layers, respectively, on the 2nd and 3rd GPUs. It is important to note that synchronization of kernel execution and data transfer operations across the entire process can be costly in terms of performance. This synchronization can 9

Pseudocode 2: OpenMP Offloading 1-GPU Baseline: allocation, host–device transfers, single-GPU computation, and deallocation. /* –– Setup ––

*/

1 read n and T from command line

κmax ← 0.95, ∆t ← h2 /(12κmax ), 3 allocate host arrays uold , unew , rhs, κ 4 initialise uold , unew , rhs, κ 5 bytes ← n3 · sizeof(double) 6 host_id ← omp_get_initial_device() 7 device_id ← omp_get_default_device()

2 h ← 1/(n − 1),

β ← ∆t/(2h2 )

// CPU // GPU

/* –– GPU memory allocation –– 8 d_uold ← omp_target_alloc(bytes, device_id) 9 d_unew ← omp_target_alloc(bytes, device_id) 10 d_rhs ← omp_target_alloc(bytes, device_id) 11 d_κ ← omp_target_alloc(bytes, device_id)

*/

/* –– Host → device transfer –– 12 omp_target_memcpy( d_uold , uold , bytes, 0, 0, device_id, host_id) 13 omp_target_memcpy( d_unew , unew , bytes, 0, 0, device_id, host_id) 14 omp_target_memcpy( d_rhs, rhs, bytes, 0, 0, device_id, host_id) 15 omp_target_memcpy( d_κ, κ, bytes, 0, 0, device_id, host_id)

*/

/* –– Time-stepping loop ––

*/

16 start ← omp_get_wtime() 17 while t < T do 18 19 20

COMPUTE_SINGLE_GPU( d_unew , d_uold , d_rhs, d_κ, device_id) swap(d_uold , d_unew ) // device-pointer swap; no data movement t ← t + dt

21 elapsed ← omp_get_wtime() − start

/* –– Device → host transfer –– 22 omp_target_memcpy( uold , d_uold , bytes, 0, 0, host_id, device_id)

*/

/* –– GPU memory deallocation –– 23 omp_target_free(d_uold , device_id) omp_target_free(d_unew , device_id) 24 omp_target_free(d_rhs, device_id) omp_target_free(d_κ, device_id) 25 free host arrays uold , unew , rhs, κ

*/

Pseudocode 3: COMPUTE_SINGLE_GPU: OpenMP target computation over the interior grid points. 1 procedure COMPUTE_SINLGE_GPU(unew , uold , rhs, κ, device_id) {

8

#pragma omp target is_device_ptr(unew , uold , rhs, κ) device(device_id) { #pragma omp teams distribute parallel for collapse(3) for k ← 1 to nk − 2 do for j ← 1 to ny − 2 do for i ← 1 to nx − 2 do unew [k, j, i] ← STENCIL(uold , rhs, κ, k, j, i)

9

}

2 3 4 5 6 7

10 }

10

pose a significant issue in both the 2-GPU and 4-GPU configurations, with the 4-GPU case exhibiting greater performance degradation compared to the 2-GPU setup. In contrast, CUDA inherently supports asynchronous kernel execution, allowing host threads to proceed with subsequent instructions after initiating kernel computations. Unfortunately, the OpenMP Offloading model differs in this respect. To address these limitations, we implemented an asynchronous computation approach for both the interior and halo point kernels using the nowait directive. The pseudocode illustrating this algorithm for a 2-GPU configuration is presented in 5, accompanied by a schematic workflow depicted in Appendix Figure 9. Pseudocode 4: Common preparation and two-GPU domain decomposition used by Pseudocodes 5–11. 1 read n and T ; h ← 1/(n − 1), κmax ← 0.95, ∆t ← h2 /(12κmax ), β ← ∆t/(2h2 ) 2 decompose the domain along the k-axis into data blocks D0 and D1 , including one ghost plane per interface 3 initialize host arrays ug , rhs g , and κg for each GPU g ∈ {0, 1} 4 foreach GPU g ∈ {0, 1} do 5 6 7 8

// Allocate device-resident arrays dg uold ← omp_target_alloc(size g , g) dg unew ← omp_target_alloc(size g , g) dg rhs ← omp_target_alloc(size g , g) dg κ ← omp_target_alloc(size g , g)

12

// Transfer initialized data from the host to GPU g omp_target_memcpy (dg uold , ug , size g , . . . , GPUg , host) omp_target_memcpy (dg unew , ug , size g , . . . , GPUg , host) omp_target_memcpy (dg rhs, rhs g , size g , . . . , GPUg , host) omp_target_memcpy (dg κ, κg , size g , . . . , GPUg , host)

6.1.3

OpenMP Offloading {2,4}-GPU Version 2

9 10 11

In this scenario, we utilize asynchronous operations for both computation and data transfer, specifically between the Host and Device, as well as between Device and Host. The interior and halo computations are initialized on each GPU block, while the halo computations are linked to a depobj ID to facilitate synchronization for the subsequent Device-to-Host (DtoH) copy operation. Once the DtoH operation is complete, the updated values of u are swapped. To guarantee that the DtoH operation fully completes before the swapping occurs, we introduce the directive #pragma omp taskwait. The pseudocode for this algorithm, demonstrated in a 2-GPU scenario, is presented in 8, and the corresponding schematic workflow is illustrated in Appendix Figure 10. 6.1.4

OpenMP Offloading {2,4}-GPU Version 3

This version closely resembles the previously described OpenMP Offloading {2,4}-GPU Version 2, but it introduces multiple threads for managing GPU control. Each GPU operation, including both kernel execution and memory copying, is handled by its respective host thread. As in the prior version, each kernel execution and data transfer is identified by a unique depobj ID, ensuring that computation and communication can overlap for each device. To mitigate the risk of data races, #pragma omp barrier constructs are strategically placed to guarantee numerical correctness before any data swapping occurs. The pseudocode for this algorithm, demonstrated in a 2-GPU scenario, is shown in 10, and the corresponding workflow is illustrated in Appendix Figure 11.

11

Pseudocode 5: OpenMP Offloading 2-GPU Version 1: synchronous two-GPU version with host-staged halo exchange. 1 setup + host initialization as in Pseudocode 4 2 allocate device arrays and transfer D0 to GPU 0 and D1 to GPU 1 3 allocate host halo buffers H0 and H1 4 for step ← 0 to Nsteps − 1 do 5 6

// Launch computation on GPU 0 Interior(D0 , GPU0 ) Boundary(B0 , GPU0 )

8

// Launch computation on GPU 1 Interior(D1 , GPU1 ) Boundary(B1 , GPU1 )

9

// Complete computation on both GPUs #pragma omp taskwait

7

10 11

12 13

14 15

// Halo block 0: GPU 0 → host → GPU 1 omp_target_memcpy (H0 , B0 , halo_bytes, . . . , host, GPU0 ) omp_target_memcpy (G1 , H0 , halo_bytes, . . . , GPU1 , host) // Halo block 1: GPU 1 → host → GPU 0 omp_target_memcpy (H1 , B1 , halo_bytes, . . . , host, GPU1 ) omp_target_memcpy (G0 , H1 , halo_bytes, . . . , GPU0 , host) // Advance to the next time level new SWAP(uold 0 , u0 ) old SWAP(u1 , unew 1 )

16 copy the local solutions back to the host and reconstruct the global domain

Pseudocode 6: OpenMP target kernel for the interior data block. 1 procedure INTERIOR(unew , uold , rhs, κ, kstart , kend , device_id) {

8

#pragma omp target is_device_ptr(unew , uold , rhs, κ) device(device_id) nowait { #pragma omp teams distribute parallel for collapse(3) for k ← kstart to kend do for j ← 1 to ny − 2 do for i ← 1 to nx − 2 do unew [k, j, i] ← STENCIL(uold , rhs, κ, k, j, i)

9

}

2 3 4 5 6 7

10 }

12

Pseudocode 7: OpenMP target kernel for an interface boundary block. 1 procedure BOUNDARY(unew , uold , rhs, κ, kboundary , device_id) {

7

#pragma omp target is_device_ptr(unew , uold , rhs, κ) device(device_id) nowait { #pragma omp teams distribute parallel for collapse(2) for j ← 1 to ny − 2 do for i ← 1 to nx − 2 do unew [kboundary , j, i] ← STENCIL(uold , rhs, κ, kboundary , j, i)

8

}

2 3 4 5 6

9 }

Pseudocode 8: OpenMP Offloading 2-GPU Version 2: asynchronous two-GPU stream version with host-staged halo exchange. 1 setup + host initialization as in Pseudocode 4 2 allocate device arrays and transfer D0 to GPU 0 and D1 to GPU 1 3 allocate host halo buffers H0 and H1 4 initialize dependency tokens tk0→1 and tk1→0 5 #pragma

omp depobj(dep0→1 ) depend(inout: 6 #pragma omp depobj(dep1→0 ) depend(inout:

tk0→1 ) tk1→0 )

7 for step ← 0 to Nsteps − 1 do 8 9

10 11

12 13

// Launch computation on GPU 0 Interior(D0 , GPU0 ) Boundary_async(B0 , GPU0 , tk0→1 ) // Launch computation on GPU 1 Interior(D1 , GPU1 ) Boundary_async(B1 , GPU1 , tk1→0 ) // Asynchronous halo chain: GPU 0 → host → GPU 1 omp_target_memcpy_async (H0 , B0 , halo_bytes, . . . , host, GPU0 , 1, &dep0→1 ) omp_target_memcpy_async (G1 , H0 , halo_bytes, . . . , GPU1 , host, 1, &dep0→1 )

15

// Asynchronous halo chain: GPU 1 → host → GPU 0 omp_target_memcpy_async (H1 , B1 , halo_bytes, . . . , host, GPU1 , 1, &dep1→0 ) omp_target_memcpy_async (G0 , H1 , halo_bytes, . . . , GPU0 , host, 1, &dep1→0 )

16

// Complete both interior kernels and halo chains #pragma omp taskwait

14

17 18

// Advance to the next time level new SWAP(uold 0 , u0 ) new SWAP(uold 1 , u1 )

19 #pragma

omp depobj(dep0→1 ) destroy omp depobj(dep1→0 ) destroy 21 copy the local solutions back to the host and reconstruct the global domain

20 #pragma

13

Pseudocode 9: OpenMP asynchronous target kernel for an interface boundary block. 1 procedure BOUNDARY_ASYNC(unew , uold , rhs, κ, kboundary , device_id, dep_token) {

7

#pragma omp target is_device_ptr(unew , uold , rhs, κ) device(device_id) nowait depend(out: dep_token[0]) { #pragma omp teams distribute parallel for collapse(2) for j ← 1 to ny − 2 do for i ← 1 to nx − 2 do unew [kboundary , j, i] ← STENCIL(uold , rhs, κ, kboundary , j, i)

8

}

2

3 4 5 6

9 }

6.1.5

OpenMP Offloading {2,4}-GPU Version 4

This iteration builds upon the previous version, OpenMP Offloading {2,4}-GPU Version 3. It eliminates the host-stage data transfer of halo data by leveraging P2P functionality. This approach negates the need for data transfers between the device and host, as detailed in the pseudocode for the 2-GPU scenario presented in 11. The corresponding schematic workflow is illustrated in Appendix Figure 12. 6.1.6

CUDA, HIP, and SYCL

In this research, we utilize baseline versions of CUDA, HIP, and SYCL alongside OpenMP Offloading implementations. Our primary objective is to avoid the redundant task of developing three separate versions for CUDA, HIP, and SYCL from the ground up. Instead, we will create the CUDA version first, which can then be converted to HIP using the command hipify-perl mycode.cu > mycode.hip and to SYCL with dpct mycode.cu –out-root=sycl_out. However, it is important to note that some fine-tuning will be necessary after conversion to achieve optimal performance; otherwise, the performance may fall significantly short compared to the OpenMP Offloading implementations. Leveraging our extensive expertise in CUDA programming, particularly within single-GPU and multi-GPU environments [24, 25], we have devised an efficient strategy for both the CUDA single-GPU and multi-GPU implementations that parallels our approach with OpenMP Offloading. This strategy includes effective thread block management, memory coalescing, and enhancements to communication and memory transfer efficiency. 6.1.7

Numerical Correctness

The numerical double-precision correctness of various implementations of the GPU programming model was compared element-wise with a serial CPU reference using an absolute tolerance of 10−6 and a relative tolerance of 10−5 . A result was accepted when every value satisfied uGPU − uCPU ≤ 10−6 + 10−5 uCPU , i i i with matching array sizes and no NaN or infinite values. All evaluated GPU implementations are satisfied this criterion.

6.2

Evaluation of Hidden Helper Thread Configurations

OpenMP provides a mechanism to utilize hidden helper threads, which can be controlled via the environment variable LIBOMP_NUM_HIDDEN_HELPER_THREADS=N. These threads enhance the asynchronous behavior of compute kernels that utilize the nowait configuration during GPU Offloading and various data copy operations. By default, eight hidden helper threads are employed. Therefore, it is crucial to evaluate the impact of these threads within the context of Multi-GPU OpenMP Offloading and their

14

Pseudocode 10: OpenMP Offloading 2-GPU Version 3: asynchronous two-GPU version with two host threads and host-staged halo exchange. 1 setup + host initialization as in Pseudocode 4 2 allocate device arrays and transfer D0 to GPU 0 and D1 to GPU 1 3 allocate shared host halo buffers H0 and H1 4 initialize dependency tokens tk0→1 and tk1→0 5 #pragma

omp depobj(dep0→1 ) depend(inout: 6 #pragma omp depobj(dep1→0 ) depend(inout:

tk0→1 ) tk1→0 )

omp parallel num_threads(2) default(shared) { tid ← omp_get_thread_num() device_id ← tid

7 #pragma 8 9 10 11

12 13

14 15 16

17 18

for step ← 0 to Nsteps − 1 do if tid = 0 then // Thread 0 controls GPU 0 Interior(D0 , GPU0 ) Boundary_async(B0 , GPU0 , tk0→1 ) // Push halo block 0: GPU 0 → host → GPU 1 omp_target_memcpy_async (H0 , B0 , halo_bytes, . . . , host, GPU0 , 1, &dep0→1 ) omp_target_memcpy_async (G1 , H0 , halo_bytes, . . . , GPU1 , host, 1, &dep0→1 ) else // Thread 1 controls GPU 1 Interior(D1 , GPU1 ) Boundary_async(B1 , GPU1 , tk1→0 ) // Push halo block 1: GPU 1 → host → GPU 0 omp_target_memcpy_async (H1 , B1 , halo_bytes, . . . , host, GPU1 , 1, &dep1→0 ) omp_target_memcpy_async (G0 , H1 , halo_bytes, . . . , GPU0 , host, 1, &dep1→0 )

19 20

21

22 23 24 25 26

// Complete both compute and communication chains #pragma omp barrier // One host thread advances both GPU time levels #pragma omp single { new SWAP(uold 0 , u0 ) old SWAP(u1 , unew 1 ) }

27 } 28 #pragma

omp depobj(dep0→1 ) destroy omp depobj(dep1→0 ) destroy 30 copy the local solutions back to the host and reconstruct the global domain

29 #pragma

15

Pseudocode 11: OpenMP Offloading 2-GPU Version 4: asynchronous two-GPU version with two host threads and direct P2P halo exchange. 1 setup + host initialization as in Pseudocode 4 2 allocate device arrays and transfer D0 to GPU 0 and D1 to GPU 1

// No host halo buffers are required 3 initialize dependency tokens tk0→1 and tk1→0 4 #pragma 5 #pragma

omp depobj(dep0→1 ) depend(inout: omp depobj(dep1→0 ) depend(inout:

tk0→1 ) tk1→0 )

omp parallel num_threads(2) default(shared) { tid ← omp_get_thread_num() device_id ← tid

6 #pragma 7 8 9 10

11 12

for step ← 0 to Nsteps − 1 do if tid = 0 then // Thread 0 controls GPU 0 Interior(D0 , GPU0 ) Boundary_async(B0 , GPU0 , tk0→1 )

16

// Direct halo push: GPU 0 → GPU 1 omp_target_memcpy_async (G1 , B0 , halo_bytes, . . . , GPU1 , GPU0 , 1, &dep0→1 ) else // Thread 1 controls GPU 1 Interior(D1 , GPU1 ) Boundary_async(B1 , GPU1 , tk1→0 )

17

// Direct halo push: GPU 1 → GPU 0 omp_target_memcpy_async (G0 , B1 , halo_bytes, . . . , GPU0 , GPU1 , 1, &dep1→0 )

13 14

15

18

19 20 21 22 23

// Complete both compute and peer-copy chains #pragma omp barrier // One host thread advances both GPU time levels #pragma omp single { new SWAP(uold 0 , u0 ) new SWAP(uold 1 , u1 ) }

24 } 25 #pragma

omp depobj(dep0→1 ) destroy 26 #pragma omp depobj(dep1→0 ) destroy 27 copy the local solutions back to the host and reconstruct the global domain

16

60

11523

5123

140 120

1.98×

80 1.98×

7683

100 80 60

11523

4.12×

0

12803

5123

50

6403

1.89×

7683 8963 Grid Size (N 3).

3.97×

12803

30 20 10

10243

11523

Intel 1550 - 4 GPUs (composite) (~500 time steps). SYCL 1-GPU Baseline OpenMP Off. 4-GPU Version 2 OpenMP Off. 1-GPU Baseline SYCL 4-GPU Version 3 SYCL 4-GPU Version 1 OpenMP Off. 4-GPU Version 3 OpenMP Off. 4-GPU Version 1 SYCL 4-GPU Version 4 SYCL 4-GPU Version 2 OpenMP Off. 4-GPU Version 4

2.12×

1.84×

1.79×

5123

1.91×

20

0

8963 10243 Grid Size (N 3).

3.62×

1.93×

Solver Time in Seconds.

40

30

10

7683

(d) AMD MI250X, 4 GPUs.

Intel 1550 - 2 GPUs (composite) (~500 time steps). SYCL 1-GPU Baseline OpenMP Off. 2-GPU Version 2 OpenMP Off. 1-GPU Baseline SYCL 2-GPU Version 3 SYCL 2-GPU Version 1 OpenMP Off. 2-GPU Version 3 OpenMP Off. 2-GPU Version 1 SYCL 2-GPU Version 4 SYCL 2-GPU Version 2 OpenMP Off. 2-GPU Version 4

1.85×

Solver Time in Seconds.

40

6403

3.99×

2.06×

40

(c) AMD MI250X, 2 GPUs. 50

4.39×

AMD MI250X - 4 GPUs (1 GCD per module) (~500 time steps). HIP 1-GPU Baseline OpenMP Off. 4-GPU Version 2 OpenMP Off. 1-GPU Baseline HIP 4-GPU Version 3 HIP 4-GPU Version 1 OpenMP Off. 4-GPU Version 3 OpenMP Off. 4-GPU Version 1 HIP 4-GPU Version 4 HIP 4-GPU Version 2 OpenMP Off. 4-GPU Version 4

20

8963 10243 Grid Size (N 3).

12803

0

12803

(e) Intel 1550, 2 GPUs.

5123

6403

3.42×

6403

11523

3.28×

5123

2.01×

2.00×

2.00×

20 0

Solver Time in Seconds.

100

40

8963 10243 Grid Size (N 3).

3.94×

AMD MI250X - 2 GPUs (1 GCD per module) (~500 time steps). HIP 1-GPU Baseline OpenMP Off. 2-GPU Version 2 OpenMP Off. 1-GPU Baseline HIP 2-GPU Version 3 HIP 2-GPU Version 1 OpenMP Off. 2-GPU Version 3 OpenMP Off. 2-GPU Version 1 HIP 2-GPU Version 4 HIP 2-GPU Version 2 OpenMP Off. 2-GPU Version 4

60

7683

(b) NVIDIA H100, 4 GPUs.

1.97×

Solver Time in Seconds.

120

4.29×

6403

(a) NVIDIA H100, 2 GPUs. 140

5.07×

0

12803

4.12×

8963 10243 Grid Size (N 3).

4.13×

7683

3.07×

6403

10

4.16×

5123

20

2.75×

2.10×

2.35×

10 0

2.11×

20

30

4.16×

2.07×

2.03×

30

40

4.25×

2.07×

40

4.61×

Solver Time in Seconds.

50

2.03×

Solver Time in Seconds.

50

NVIDIA H100 - 4 GPUs (~500 time steps). CUDA 1-GPU Baseline OpenMP Off. 4-GPU Version 2 OpenMP Off. 1-GPU Baseline CUDA 4-GPU Version 3 CUDA 4-GPU Version 1 OpenMP Off. 4-GPU Version 3 OpenMP Off. 4-GPU Version 1 CUDA 4-GPU Version 4 CUDA 4-GPU Version 2 OpenMP Off. 4-GPU Version 4

3.82×

60

NVIDIA H100 - 2 GPUs (~500 time steps). CUDA 1-GPU Baseline OpenMP Off. 2-GPU Version 2 OpenMP Off. 1-GPU Baseline CUDA 2-GPU Version 3 CUDA 2-GPU Version 1 OpenMP Off. 2-GPU Version 3 OpenMP Off. 2-GPU Version 1 CUDA 2-GPU Version 4 CUDA 2-GPU Version 2 OpenMP Off. 2-GPU Version 4

7683 8963 Grid Size (N 3).

10243

12803

(f) Intel 1550, 4 GPUs.

Figure 3: Comparison of runtime performance across different GPU programming models along with various GPU architectures for different grid sizes N 3 , utilizing configurations with 2 and 4 GPUs. The speedup values indicate the performance of OpenMP Offloading {2,4}-GPU version 4 relative to the OpenMP Offloading 1-GPU Baseline.

17

20

33.3

39.7

7.52

Transferred size (GB, type: double)

13.9

20.4

(a)

100

HIP (2 GPUs) OpenMP Off.(2 GPUs)

60

40 20

33.3

39.7

0

50

30

7.52

13.9

20.4

HIP (2 GPUs) OpenMP Off.(2 GPUs)

60

40

33.3

39.7

0

50

30

7.52

13.9

20.4

33.3

0

39.7

60 40

7.52

13.9

20.4

Theoretical BW, 1 GPU (64 GB/s)

60 Mean absolute DtoH gap CUDA vs. OpenMP Off. (1 GPU): 4.9 GB/s CUDA vs. OpenMP Off. (2 GPUs): 5.5 GB/s

40

Bandwidth (GB/s)

80

20

33.3

39.7

0

60

40 30 20

Mean absolute DtoH gap HIP vs. OpenMP Off. (1 GPU): 0.4 GB/s HIP vs. OpenMP Off. (2 GPUs): 1.2 GB/s

20.4

33.3

0

39.7

80 60 40

7.52

13.9

20.4

NVIDIA H100 - P2P

200 Theoretical BW, 1 GPU (150 GB/s)

150 100

Mean absolute P2P gap CUDA vs. OpenMP Off. (1 GPU): 2.1 GB/s CUDA vs. OpenMP Off. (2 GPUs): 5.3 GB/s

50

4.16

Transferred Size (GB, type: double)

9.92

0

HIP (1 GPU) OpenMP Off (1 GPU)

Bandwidth (GB/s)

250

OpenMP Off. (2 GPUs)

200

150 125 Theoretical BW, 1 GPU (100 GB/s)

100 75

Mean absolute P2P gap HIP vs. OpenMP Off. (1 GPU): 0.1 GB/s HIP vs. OpenMP Off. (2 GPUs): 10.1 GB/s

1.74

2.55

OpenMP Off. (1 GPU) SYCL (1 GPU)

175

HIP (2 GPUs) OpenMP Off (2 GPUs)

0.940

33.3

39.7

Bandwidth (GB/s)

0

Intel 1550 (composite) - P2P

Theoretical BW, 2 GPUs (200 GB/s)

300

26.8

(l)

AMD MI250X (1 GCD) - P2P

Theoretical BW, 2 GPUs (300 GB/s)

20

Transferred size (GB, type: double)

(k)

CUDA (1 GPU) OpenMP Off. (1 GPU) CUDA (2 GPUs)

100

Theoretical BW, 1 GPU (64 GB/s)

Transferred size (GB, type: double)

(j)

0

120

Mean absolute DtoH gap SYCL vs. OpenMP Off. (1 GPU): 1.6 GB/s SYCL vs. OpenMP Off. (2 GPUs): 2.3 GB/s

10 26.8

39.7

Theoretical BW, 2 GPUs (128 GB/s)

Theoretical BW, 1 GPU (36 GB/s)

13.9

33.3

SYCL (1 GPU) OpenMP Off. (1 GPU) SYCL (2 GPUs) OpenMP Off. (2 GPUs)

70

50

7.52

26.8

(i)

OpenMP Off.(2 GPUs)

Transferred size (GB, type: double)

20

Intel 1550 (composite) - DtoH: asynchronous

Theoretical BW, 2 GPUs (72 GB/s)

100

100

Transferred size (GB, type: double)

HIP (1 GPU) OpenMP Off.(1 GPU) HIP (2 GPUs)

120

120

80

(h)

Theoretical BW, 2 GPUs (128 GB/s)

0

Theoretical BW, 1 GPU (64 GB/s)

AMD MI250X (1 GCD) - DtoH: asynchronous

CUDA (1 GPU) OpenMP Off. (1 GPU) CUDA (2 GPUs) OpenMP Off. (2 GPUs)

39.7

Mean absolute HtoD gap SYCL vs. OpenMP Off. (1 GPU): 3.4 GB/s SYCL vs. OpenMP Off. (2 GPUs): 1.2 GB/s

10 26.8

33.3

SYCL (1 GPU) OpenMP Off. (1 GPU) SYCL (2 GPUs) OpenMP Off. (2 GPUs)

Transferred size (GB, type: double)

(g)

(m)

20

Mean absolute HtoD gap HIP vs. OpenMP Off. (1 GPU): 1.5 GB/s HIP vs. OpenMP Off. (2 GPUs): 2.2 GB/s

NVIDIA H100 - DtoH: asynchronous

3.35

40

Theoretical BW, 1 GPU (36 GB/s)

Transferred size (GB, type: double)

26.8

Theoretical BW, 2 GPUs (128 GB/s)

100

20

2.55

20.4

Intel 1550 (composite) - HtoD: asynchronous 70

Bandwidth (GB/s)

Mean absolute HtoD gap CUDA vs. OpenMP Off. (1 GPU): 3.3 GB/s CUDA vs. OpenMP Off. (2 GPUs): 6.0 GB/s

1.74

13.9

Theoretical BW, 2 GPUs (72 GB/s)

60

0.94

7.52

20

(f)

HIP (1 GPU) OpenMP Off.(1 GPU)

Theoretical BW, 1 GPU (64 GB/s)

26.8

40

Transferred size (GB, type: double)

120

80

20.4

60

AMD MI250X (1 GCD) - HtoD: asynchronous

Theoretical BW, 2 GPUs (128 GB/s)

13.9

0

39.7

80 Theoretical BW, 1 GPU (64 GB/s)

(e)

NVIDIA H100 - HtoD: asynchronous

7.52

33.3

100

Mean absolute DtoH gap SYCL vs. OpenMP Off. (1 GPU): 1.0 GB/s SYCL vs. OpenMP Off. (2 GPUs): 1.6 GB/s

10 26.8

0

120

OpenMP Off. (2 GPUs)

Transferred size (GB, type: double)

CUDA (1 GPU) OpenMP Off. (1 GPU) CUDA (2 GPUs) OpenMP Off. (2 GPUs)

26.8

20

Mean absolute DtoH gap HIP vs. OpenMP Off. (1 GPU): 0.2 GB/s HIP vs. OpenMP Off. (2 GPUs): 0.6 GB/s

(d)

20.4

40

Theoretical BW, 1 GPU (36 GB/s)

Transferred size (GB, type: double)

39.7

SYCL (1 GPU) OpenMP Off. (1 GPU) SYCL (2 GPUs)

Bandwidth (GB/s)

Mean absolute DtoH gap CUDA vs. OpenMP Off. (1 GPU): 0.9 GB/s CUDA vs. OpenMP Off. (2 GPUs): 8.6 GB/s

33.3

Theoretical BW, 2 GPUs (128 GB/s)

70

60

26.8

(c)

Theoretical BW, 2 GPUs (72 GB/s)

Theoretical BW, 1 GPU (64 GB/s)

13.9

20.4

Intel 1550 (composite) - DtoH: synchronous

HIP (1 GPU) OpenMP Off.(1 GPU)

80

7.52

13.9

20

Transferred size (GB, type: double)

120

Bandwidth (GB/s)

OpenMP Off. (2 GPUs)

26.8

7.52

(b)

Theoretical BW, 2 GPUs (128 GB/s)

20.4

39.7

AMD MI250X (1 GCD) - DtoH: synchronous

CUDA (1 GPU) OpenMP Off. (1 GPU) CUDA (2 GPUs)

13.9

0

Transferred size (GB, type: double)

NVIDIA H100 - DtoH: synchronous

7.52

33.3

40 Mean absolute HtoD gap SYCL vs. OpenMP Off. (1 GPU): 1.3 GB/s SYCL vs. OpenMP Off. (2 GPUs): 0.9 GB/s

10 26.8

60

Bandwidth (GB/s)

26.8

20

Mean absolute HtoD gap HIP vs. OpenMP Off. (1 GPU): 1.3 GB/s HIP vs. OpenMP Off. (2 GPUs): 2.2 GB/s

80 Theoretical BW, 1 GPU (64 GB/s)

Bandwidth (GB/s)

20.4

30

50

Bandwidth (GB/s)

13.9

0

40

Theoretical BW, 1 GPU (36 GB/s)

100

Bandwidth (GB/s)

40

50

120

OpenMP Off. (2 GPUs)

Bandwidth (GB/s)

60

60

Bandwidth (GB/s)

Theoretical BW, 1 GPU (64 GB/s)

SYCL (1 GPU) OpenMP Off. (1 GPU) SYCL (2 GPUs)

70

Bandwidth (GB/s)

80

Bandwidth (GB/s)

100

Mean absolute HtoD gap CUDA vs. OpenMP Off. (1 GPU): 0.0 GB/s CUDA vs. OpenMP Off. (2 GPUs): 6.5 GB/s

Theoretical BW, 2 GPUs (128 GB/s)

Theoretical BW, 2 GPUs (72 GB/s) HIP (1 GPU) OpenMP Off.(1 GPU) HIP (2 GPUs) OpenMP Off.(2 GPUs)

120

OpenMP Off. (2 GPUs)

7.52

Intel 1550 (composite) - HtoD: synchronous

AMD MI250X (1 GCD) - HtoD: synchronous

Theoretical BW, 2 GPUs (128 GB/s)

Theoretical BW, 2 GPUs (212 GB/s)

200

OpenMP Off. (2 GPUs) SYCL (2 GPUs)

150 Mean absolute P2P gap SYCL vs. OpenMP Off. (1 GPU): 44.4 GB/s Theoretical BW, 1 GPU (106 GB/s) SYCL vs. OpenMP Off. (2 GPUs): 19.5 GB/s

100

50

Bandwidth (GB/s)

NVIDIA H100 - HtoD: synchronous CUDA (1 GPU) OpenMP Off. (1 GPU) CUDA (2 GPUs)

25 3.35

4.16

Transferred Size (GB, type: double)

(n)

9.92

0

0.940

1.74

2.55

3.35

4.16

9.92

0

Transferred Size (GB, type: double)

(o)

Figure 4: Theoretical and measured intra-node transfer bandwidth on three GPU architectures (columns: NVIDIA H100, AMD MI250X, Intel Max 1550) for five transfer modes (rows: synchronous HtoD and DtoH, asynchronous HtoD and DtoH, and P2P). In the host-transfer rows, a single GPU is driven by one host thread and two GPUs by two host threads using OpenMP; the P2P row uses two GPUs and two host threads throughout.

18

Table 4: Comparison of selected low-level API features across CUDA, HIP, SYCL, and OpenMP Offloading. Category

CUDA

HIP

SYCL

OpenMP Offloading

Memory management

cudaMalloc() cudaFree()

hipMalloc() hipFree()

malloc_device() free()

omp_target_alloc() omp_target_free() Supported through device tasks/nowait

Streams/QueuescudaStreamCreate() hipStreamCreate() queue Synchronous cudaMemcpy() memory copy

hipMemcpy()

queue.memcpy().wait() omp_target_memcpy()

Asynchronous cudaMemcpyAsync() memory copy

hipMemcpyAsync()

queue.memcpy()

omp_target_memcpy_async()

omp_target_memcpy() Device-tocudaMemcpyPeer() hipMemcpyPeer() queue.memcpy().wait() with device IDs device omp_target_memcpy_async() hipMemcpyPeerAsync() queue.memcpy() memory copy cudaMemcpyPeerAsync() with device IDs

overall effect on performance. To conduct this investigation, we established three categories of hidden helper thread configurations, as detailed in Table 5. Table 5: Host-thread and hidden-helper-thread configurations used for the evaluting OpenMP Offloading {2,4}-GPU Versions. Implementation

Host threads

Hidden Helper Threads (HHT) HHT-0

HHT-GPU

HHT-8

OpenMP Offloading 2-GPU Version 1 OpenMP Offloading 2-GPU Version 2 OpenMP Offloading 2-GPU Version 3 OpenMP Offloading 2-GPU Version 4

1 1 2 2

0 0 0 0

2 2 2 2

8 8 8 8

OpenMP Offloading 4-GPU Version 1 OpenMP Offloading 4-GPU Version 2 OpenMP Offloading 4-GPU Version 3 OpenMP Offloading 4-GPU Version 4

1 1 4 4

0 0 0 0

4 4 4 4

8 8 8 8

HHT-0 disables hidden helper threads; HHT-GPU uses one hidden helper thread per GPU; and HHT-8 uses a fixed total of eight hidden helper threads per GPU.

6.3

Bandwidth Measurement on NVIDIA H100, AMD MI250, and Intel 1550

To fully understand the capabilities of the GPU bandwidth network between GPUs and the CPU, it is essential to analyze the performance of memory transfer APIs across various GPU programming models and architectures. We have established a set of benchmarks designed to measure the performance of GPU programming model APIs for data transfer between the CPU and GPU, as well as between GPUs (P2P). For instance, we utilized all the functionalities of the memory transfer APIs outlined in Table 4. This methodology expands upon the framework we established in our earlier work on OpenACC [26]. In the context of OpenMP Offloading, we developed several benchmark cases to evaluate the actual bandwidth offered by the OpenMP Offloading API for communication between the CPU and GPU, as well as between multiple GPUs. • HtoD (Host to Device): – 1 GPU: Using omp_target_memcpy() and omp_target_memcpy_async() with a single host thread, data is transferred from the host to the device, as depicted in Figure 5a. 19

– 2 GPUs: Using omp_target_memcpy() and omp_target_memcpy_async() with two host threads, two sets of data can be sent concurrently from the host to both devices, as illustrated in Figure 5b. • DtoH (Device to Host): – 1 GPU: Using omp_target_memcpy() and omp_target_memcpy_async() with a single host thread, data is transferred from the device to the host, as shown in Figure 5a. – 2 GPUs: Using omp_target_memcpy() and omp_target_memcpy_async() with two host threads, two sets of data can be sent concurrently from the device to the host, as illustrated in Figure 5b. • P2P (Peer-to-Peer): – 1 GPU: Using omp_target_memcpy_async() with a single host thread, data is transferred from one GPU to another, as depicted in Figure 5a. – 2 GPUs: Using omp_target_memcpy_async() with two host threads, data can be transferred concurrently between the two GPUs, as illustrated in Figure 5b. In this research on solving the 3D heat equation, the estimates for halo data movement are as follows: For double-precision floating-point values (b = 8 bytes) and a one-layer halo (h = 1), each GPU is responsible for sending and receiving one N 2 plane. Consequently, the communication volume per GPU can be calculated as: VGPU = 2hN 2 b = 2 × 1 × N 2 × 8 = 16N 2 bytes. For global domain sizes N 3 in the set {5123 , 6403 , 7683 , 8963 , 10243 , 11523 , 12803 }, the corresponding send-and-receive communication volumes per GPU are 0.004194, 0.006554, 0.009437, 0.012845, 0.016777, 0.021234, and 0.026214 GB, respectively. This volume of data is essential for developing the performance model for multi-GPU systems, which is further discussed in Subsections 6.4 and 7.3. Host / NUMA node

2 threads (omp sections)

DtoH · thread 0

d_A=h_A buffer 0=d_B

System memory

PCIe/IF run concurrently

CPU Host / NUMA node CPU multicore

System memory host DRAM

HtoD PCIe/IF DtoH

GPU 0 Device mem VRAM

d_A

d_A=h_A

NV6 IF Xe Link P2P

GPU 1 Device mem VRAM

d_B

d_B=d_A

h_A

(a)

h_A

h_B

host DRAM

host DRAM

HtoD · thread 1 DtoH · thread 1

P2P

thd.0

GPU 0 Device VRAM d_A

thd.1

HtoD · thread 0

NV6 / IF / Xe Link

GPU 1 Device VRAM d_B d_B=h_B buffer 1=d_A

(b)

Figure 5: Comparison of bandwidth measurement scenarios: (a) 1 GPU: Schematic overview of measurements for HtoD, DtoH, and P2P between the host and GPU, as well as between GPUs using 1 host thread; (b) 2 GPUs: Schematic overview of measurements for HtoD, DtoH, and P2P between the host and GPU, and between GPUs using 2 host threads concurrently.

6.4

Performance Model

Here, we derive the performance model to gain a better understanding of the parallel implementations of OpenMP Offloading, along with hardware bandwidth measurements. This will help predetermine the expected runtime cost and, overall, draw the performance model. For this, we assume that the computational cost is determined by the number of points that need to be calculated in the 3D domain. Specifically, (N − 2)3 points need to be computed, where N represents the total number of points in the grid, excluding the boundary layers. Therefore, the computation time per step can be modeled as: 20

tstep = k(N − 2)3 ,

(7)

where k denotes the overhead cost, which may arise from latency or other machine-related factors, and we assume it remains constant. Let nsteps denote the total number of time steps. In our measurements, the reported runtime Ttotal refers to the execution of nsteps tstep . We can express this relationship as follows: Ttotal = nsteps tstep . Substituting the expression for tstep from Eq. (7), we obtain: Ttotal = nsteps k(N − 2)3 . Rearranging this equation to solve for k, we arrive at: k=

Ttotal (N − 2)3 . nsteps

Since the grid size changes, we need to take the average of k. To do this, we model k̄ as follows: k̄ =

M 1 X k(Nj ). M j=1

(8)

Here, Nj represents the j-th grid size, and M is the total number of grid sizes considered. For instance, k(Nj ) is evaluated at Nj = 512, 640, . . . , 1280. Now we can model the estimated total computation time for the multi-GPU version. For the multi-GPU setup, we must consider the data movement involving HtoD, DtoH, and P2P transfers. Each operation is modeled as follows: HtoD =

N 2 bytes , Meas.BWHtoD

DtoH =

N 2 bytes , Meas.BWDtoH

P2P =

N 2 bytes . Meas.BWP2P

Since we have two approaches for data transfer, the total time for data movement via CPU and between GPUs through P2P can be represented as: Tstage = nsteps (HtoD + DtoH) TP2P = nsteps (P2P) The total computation time for the multi-GPU setup can be calculated as follows: For overlapping (concurrent) transfers: " stage Ttotal_G = nsteps

where G represents the number of GPUs. For P2P transfers: P2P Ttotal_G = nsteps

#

k̄(N − 2)3 + Tstage , G

"

(9)

#

k̄(N − 2)3 + TP2P . G

(10)

If communication is fully overlapped with computation, the equations referenced in Eq. 9 and Eq. 10 become: " stage,overlap Ttotal,G = nsteps max

21

#

k̄(N − 2)3 , Tstage , G

(11)

and

" P2P,overlap Ttotal,G = nsteps max

#

k̄(N − 2)3 , TP2P . G

(12)

These models assume ideal computational load balancing and neglect additional costs such as kernel launch overhead, synchronization, communication latency, and bandwidth contention.

7

Results and Discussion

In this section, we present a discussion of the results obtained from the experiment based on the methodology described in Section 6. This discussion includes an examination of various versions of multi-GPU OpenMP Offloading, along with comparisons to CUDA, HIP, and SYCL. We will also analyze the impact of HHT and discuss the proposed performance model using the measured bandwidth. The remainder of this section will elaborate on these findings in detail.

7.1

Multi-GPU intra-node OpenMP Offloading Performance Analysis

This subsection focuses on the performance analysis of solving 3D heat transfer problems using various versions of OpenMP Offloading in an intra-node setting. Figure 3 illustrates the overall runtime for different grid sizes of N , comparing the performance of CUDA, HIP, and SYCL. We will further analyze these approaches based on the results we obtained. 7.1.1

NVIDIA H100

In our tests on the NVIDIA H100, we examined the implementations of CUDA and OpenMP Offloading using single, 2, and 4 GPU configurations. As shown in Figures 3a and 3b, the performance of OpenMP Offloading {2,4}-GPU Version 4 in both the 2 (2x) and 4 (4x) GPU scenarios is competitive with the single GPU OpenMP Offloading implementation. Notably, the 4 GPU implementations (Versions 3 and 4) outperformed the single GPU implementation of the CUDA version. Below, we summarize the key findings from our analysis of the NVIDIA H100 regarding the 4 GPU implementation: • The 2 GPU implementation is less efficient than the 4 GPU implementation primarily because, while the workload can be distributed across two GPUs, it remains insufficient for optimal performance. Additionally, the overhead from communication contributes to the reduced efficiency of the 2 GPU implementations compared to their 4 GPU counterparts. • In contrast, all 4 GPU implementations allow for a more effectively distributed workload, making computations less expensive compared to those in single and 2 GPU implementations. Although the increased number of GPUs results in more communication, the overall computation cost is still lower. This combination leads to a significant speedup, particularly when compared to all versions of the 2 GPU implementations and the single GPU implementation. • The OpenMP Offloading {2,4}-GPU Version 4 not only surpasses the single GPU CUDA implementations but also outperforms other variations of OpenMP Offloading. Its enhanced speed is primarily achieved by concurrently launching GPU kernels with multiple host threads. Furthermore, the H100’s P2P bandwidth, as illustrated in Figure 2a, provides higher bandwidth, which significantly contributes to the speedup by avoiding data transfer routing through the CPU. 7.1.2

AMD MI250X

Within the AMD MI250X, we evaluate the HIP implementations as a reference to compare them against OpenMP Offloading implementations. As illustrated in Figure 2b, the bandwidth between the Graphics Compute Dies (GCDs) of the GPUs is modest, resulting in only a slight difference in bandwidth between the CPU and GPU. However, the AMD MI250X features a unique design in which each module 22

contains two GCDs. The memory bandwidth between these two GCDs is significantly higher than the data transfer rates between the modules of the MI250X. Notably, all multi-GPU implementations of HIP and OpenMP Offloading utilize a single GCD from each module of the MI250X, a design choice that impacts the communication dynamics between multiple GPUs compared to single-GPU implementations. Despite this characteristic, the overarching conclusion remains valid: the implementation utilizing four GPUs in OpenMP Offloading performs competitively against the two-GPU implementations and the single-GPU implementation of the HIP version. Below, we summarize key insights regarding the multi-GPU implementations of HIP and OpenMP Offloading. • Ideally, the OpenMP Offloading {2,4}-GPU Version 4 is expected to outperform all other variants. However, due to the limited computational power of the AMD MI250X (1 GCD), we do not observe a significant difference between the OpenMP Offloading 2-GPU Version 3 and OpenMP Offloading 2-GPU Version 4. In scenarios where computational time is extensive, the communication costs become negligible, leading to minimal differences in performance between Versions 3 and 4. Interestingly, 2-GPU Version 2 performs slightly better than Versions 3 and 4, indicating that the latter two achieve a better overlap of computation compared to 2-GPU Versions 1 and 2. • Additionally, in the context of the 4-GPU implementations, both OpenMP Offloading {2,4}-GPU Version 3 and OpenMP Offloading {2,4}-GPU Version 4 demonstrate significant speedup compared to the other versions within the 4-GPU category. As the workload is distributed across the four GPUs, we observe a remarkable speedup in 2- and 4-GPU Version 4, achieving nearly 2x and 4x the performance compared to a single GPU implementation. 7.1.3

Intel 1550

Concerning Intel 1550, we do the comparison of the reference GPU programming model SYCL against OpenMP Offloading on Intel GPUs. Figures 3e and 3f illustrate the runtime implementations of SYCL and OpenMP Offloading with 2 GPUs and 4 GPUs, respectively. The performance characteristics of the Intel 1550 GPU closely resemble those of the NVIDIA H100. Below, we summarize the key insights derived from the OpenMP Offloading implementation: • The data transfer between GPUs is slightly higher than the data transfer via the CPU, making OpenMP Offloading {2,4}-GPU Version 4 an improvement over OpenMP Offloading {2,4}-GPU Version 3. • As an end user, programmers can opt to use either a single tile or the entire composite GPU of the Intel Max 1550. This flexibility minimizes overall computational time compared to the other GPUs evaluated in this research. For specifications of the Intel GPU, please refer to Table 1. • The P2P data transfer performance between SYCL and OpenMP Offloading does not appear particularly impressive. The measured bandwidth gap is slightly greater than what we observed with the NVIDIA H100 and the AMD MI250X. This discrepancy is primarily due to OpenMP Offloading not delivering performance on par with SYCL when the Intel 1550 is considered as a composite. • The ideal speedup of 4x is not achieved when using 4 GPUs with the OpenMP Offloading 4-GPU Version 4. This shortfall is attributed to an increase in communication overhead, kernel launches, synchronization, and other factors, which ultimately raise the overall computational time. As a result, employing 4 GPUs only yields approximately a 3.5x speedup with the OpenMP Offloading 4-GPU version. In contrast, utilizing 2 GPUs achieves nearly a 2x speedup with the OpenMP Offloading 2-GPU Version 4. This is due to the substantial computational workload, which is sufficient to mask the communication overhead.

23

7.2

Performance Impact of Hidden Helper Threads

The impact of Hidden Helper Threads (HHT) has shown a notable performance gain across the three GPUs tested in this research. Figure 6 provides an overview of the speedup achieved through the three methodologies considered, as summarized in Table 5. In every case, HHT-8 demonstrated superior performance compared to both HTT-0 and HTT-GPU. It is evident that utilizing a greater number of HHTs significantly enhances the concurrency of asynchronous operations within the application. To further illustrate this scenario, we plotted another model using the geometric mean (gmean). As shown in Figure 7, we divided the execution time by the number of time steps (approximately 500) for each GPU architecture, implementation version, and grid size. The geometric mean of the normalized times for HHT-0, HHT-GPU, and HHT-8 served as a neutral reference. The gain was then calculated by dividing this reference time by the normalized time of each configuration. A gain above 1 indicates better performance, a gain of 1 indicates equal performance, and a gain below 1 indicates lower performance. For the overall comparison, gains across the selected implementation versions were also combined using the geometric mean. We did not experiment with more than 8 threads because, in the current configuration involving four GPU cases, there are approximately 10 asynchronous operations (comprising 4 interior kernels and 6 boundary kernels) that must be managed by the host. Additionally, each GPU is assigned a corresponding OpenMP thread to facilitate the streaming of operations. We hypothesize that increasing the number of threads beyond 8 would not yield further speedups within this problem setting. In terms of performance efficiency among the three GPUs—NVIDIA H100, AMD MI250X, and Intel 1550—the Intel 1550 exhibited slightly lower speedup compared to the NVIDIA H100 and AMD MI250X. This is attributed to the Intel 1550’s superior computational power, which allows it to complete computational tasks more quickly but incurs a penalty in communication time. Conversely, when computation takes longer, as observed with the other two GPUs (NVIDIA H100 and AMD MI250X), there is sufficient time for halo computation and communication to be completed. Consequently, we observe nearly perfect speedups with the NVIDIA H100 and AMD MI250X. This observation is further supported by the lower computational power of the NVIDIA H100 and AMD MI250X compared to the Intel 1550. The production analysis plots in Figure 3 are based on the HHT-8 configuration (with HHT set to 8).

7.3

Performance Model Analysis

To evaluate the performance efficiency of our parallelized OpenMP implementations, we compared the estimated times against the measured times. Specifically, we utilized the equation from Eq. 10 in Subsection 6.4. While we proposed several estimated time models, we are primarily focusing on the P2P concept (OpenMP Offloading {2,4}-GPU Version 4), which demonstrates the most effective parallelization among the other OpenMP Offloading models. To analyze OpenMP Offloading {2,4}-GPU Version 4, we calculated the average value, denoted as k̄, using Eq. 8, along with P2P values gleaned from Figure 4. This data is summarized in Table 6. By substituting these values into Eq. 10, we generated the model plot illustrated in Figure 8. Upon examining Figure 8, it is evident that our estimated values closely align with the measured values across all scenarios for OpenMP Offloading {2,4}-GPU Version 4. However, we noted some fluctuations in the results for the NVIDIA H100 and Intel 1550 GPUs. To enhance clarity, we also plotted the mean absolute error (MAE), including bias, alongside the relative error. Notably, for smaller grid sizes on the NVIDIA H100, we observed significant deviations. In contrast, the performance measurements for the AMD MI250X closely corresponded with our estimated values, although some discrepancies also appeared with the Intel 1550. 6.4042e-11 On the AMD MI250X, the deviation between estimated and measured values was minimal. This is primarily because we utilized only 1-GCD on the AMD MI250X, which offers lower computational power compared to the other GPUs in this research. Consequently, we experienced higher overall computational times, effectively masking the communication overhead in the OpenMP Offloading {2,4}-GPU Versions, resulting in a mean absolute percentage error (MAPE) of less than 3%. Similarly, while the NVIDIA H100 showed slightly higher deviations than the AMD MI250X, these deviations

24

N = 512³

N = 640³

1 GPU: H100 3.7 s PVC 3.5 s MI250X 8.6 s

1 GPU: H100 6.2 s PVC 6.6 s MI250X 17.1 s

5 4×

4

4×

3 2

2×

2×

1 0 2 Ver. 1

2 Ver. 2

2 Ver. 3

2 Ver. 4

4 Ver. 1

4 Ver. 2

4 Ver. 3

4 Ver. 4

2 Ver. 1

2 Ver. 2

2 Ver. 3

2 Ver. 4

4 Ver. 1

4 Ver. 2

4 Ver. 3

N = 768³

N = 896³

1 GPU: H100 13.7 s PVC 10.9 s MI250X 30.0 s

1 GPU: H100 21.6 s PVC 16.9 s MI250X 48.5 s

4 Ver. 4

5 4×

4

4×

speedup over 1 GPU (t_1GPU / t).

3 2

2×

2×

1 0 2 Ver. 1

2 Ver. 2

2 Ver. 3

2 Ver. 4

4 Ver. 1

4 Ver. 2

4 Ver. 3

4 Ver. 4

2 Ver. 1

2 Ver. 2

2 Ver. 3

2 Ver. 4

4 Ver. 1

4 Ver. 2

4 Ver. 3

N = 1024³

N = 1152³

1 GPU: H100 34.7 s PVC 24.6 s MI250X 73.0 s

1 GPU: H100 49.6 s PVC 35.8 s MI250X 103.8 s

4 Ver. 4

5 4×

4

4×

3 2

2×

2×

1 0 2 Ver. 1

2 Ver. 2

2 Ver. 3

2 Ver. 4

4 Ver. 1

4 Ver. 2

4 Ver. 3

4 Ver. 4

2 Ver. 1

2 Ver. 2

2 Ver. 3

2 Ver. 4

4 Ver. 1

4 Ver. 2

4 Ver. 3

4 Ver. 4

N = 1280³ 1 GPU: H100 65.7 s PVC 48.9 s MI250X 142.5 s

5 4×

4

H100-HHT-0

Intel-1550-HHT-0

AMD-MI250X-HHT-0

H100-HHT-GPU

Intel-1550-HHT-GPU

AMD-MI250X-HHT-GPU

H100-HHT-8

Intel-1550-HHT-8

AMD-MI250X-HHT-8

1 GPU (no gain) perfect scaling (2×, 4× w.r.t. each GPU group) 2 Ver. 1 = OpenMP Off. 2-GPU Version 1

3

2 Ver. 2 = OpenMP Off. 2-GPU Version 2

2

2 Ver. 3 = OpenMP Off. 2-GPU Version 3

2×

2 Ver. 4 = OpenMP Off. 2-GPU Version 4 4 Ver. 1 = OpenMP Off. 4-GPU Version 1

1

4 Ver. 2 = OpenMP Off. 4-GPU Version 2 4 Ver. 3 = OpenMP Off. 4-GPU Version 3 4 Ver. 4 = OpenMP Off. 4-GPU Version 4

0 2 Ver. 1

2 Ver. 2

2 Ver. 3

2 Ver. 4

4 Ver. 1

4 Ver. 2

4 Ver. 3

4 Ver. 4

OpenMP Offloading {2,4}-GPU Versions on H100, MI250X, and Intel-1550.

Figure 6: Speedup of the OpenMP Offloading {2,4}-GPU Versions implementations over each machine’s own OpenMP Offloading 1-GPU Baseline single-GPU run, for problem sizes N = 512 to N = 1280. Variant names are abbreviated: 2 Ver. 1 = OpenMP Offloading 2-GPU Version 1; and so on.

25

OpenMP Offloading {2,4}-GPU Versions gain vs gmean (geo-mean over 9 variant(s))

H100

Intel-1550

MI250X

1.2 1.1 1.0

HHT-0 HHT-GPU HHT-8

0.9 0.8 0.7 5123

6403

7683

8963

10243

Problem size, N 3

11523

12803

5123

6403

7683

8963

10243

Problem size, N 3

11523

12803

5123

6403

7683

8963

10243

Problem size, N 3

11523

12803

Figure 7: Geometric-mean performance gain of the HHT-0, HHT-GPU, and HHT-8 approaches for the OpenMP Offloading {2,4}-GPU Versions on the NVIDIA H100, Intel-1550, and AMD MI250X. The dashed horizontal line represents a gain of one relative to the geometric-mean reference. HHT-8 remains above unity for all evaluated problem sizes and GPU architectures. Table 6: Measured values of βP2P and the average computational constant ratio k̄ (see Eq. 8) for the evaluated GPU architectures, which are used to model performance (see Figure 8). Architecture

GPUs

k̄ (×10−11 )

βP2P (GB/s)

NVIDIA H100 AMD MI250X Intel 1550

2, 4 2, 4 2, 4

6.4042 13.584 4.6676

253.11 61.28 62.85

were most pronounced with smaller problem sizes. As the problem size increased, the estimated and measured values demonstrated a closer agreement. Overall, for both GPU cases, the NVIDIA H100 and AMD MI250X, we must consider the overhead costs associated with kernel launches, synchronization, and communication contention. On the Intel 1550 GPU, our estimated model tended to underestimate the measured wall-clock time. For example, in the 2-GPU case, the MAPE was around 9.0%, indicating an almost optimal scaling model. However, in the 4-GPU wall-clock scenario involving the Intel 1550, we found a mean absolute error of 23.6%, coupled with a slight negative bias of -23.6%. This negative bias reflects inherent execution overheads, including kernel launch delays, synchronization delays, communication latencies, and bandwidth contention, which our idealized estimated overlap model did not account for. It is important to note that the Intel 1550 (composite stack of GPUs) possesses greater computational power than the other two GPUs examined in this research. This increased power results in reduced computation time, consequently creating additional communication overhead. Therefore, we observe larger deviations at smaller grid sizes, such as N=512. However, this gap tends to decrease as the computational domain expands, indicating that as computation time increases, the significance of communication overhead diminishes. In summary, we conclude that our proposed model and methodology provide an optimized framework for multi-GPU computations on a single compute node. This model can be confidently applied to similar methodologies and to various problem-domain sizes to test parallel efficiency.

7.4

Summary of Key Findings

We summarize the key findings of this research as follows: • By default, OpenMP Offloading for kernel computations is generally less efficient than native GPU programming models, which are specifically designed for their respective GPU architectures. While OpenMP Offloading provides specialized clauses and constructs similar to those found in 26

Wall-clock time (s)

P2P, overlap Estimated Ttotal , 2, GPUs P2P, overlap Measured Ttotal , 2, GPUs

101 6 × 100 4 × 100 3 × 100 2 × 100

MAPE = 16.3% bias = +16.3%

6403

7683

6 × 100 4 × 100 3 × 100 2 × 100

8963 10243 11523 12803

50 0 −50

5123

Grid Size (N 3)

(a) NVIDIA H100, 2 GPUs.

Wall-clock time (s)

101 6 × 100 4 × 100

MAPE = 2.8% bias = +2.0%

20 0

−20

5123

6403

7683

4 × 101 3 × 101 2 × 101

8963 10243 11523 12803

MAPE = 9.0% bias = −9.0%

20 0

6403

7683

MAPE = 4.5% bias = +4.1%

0

−20

5123

6403

7683

8963 10243 11523 12803

Grid Size (N 3)

8963 10243 11523 12803

Grid Size (N 3)

(e) Intel 1550 (composite), 2 GPUs.

P2P, overlap Estimated Ttotal , 4, GPUs P2P, overlap Measured Ttotal, 4, GPUs

101 6 × 100 4 × 100 3 × 100 2 × 100

MAPE = 23.6% bias = −23.6%

100

Rel. err. (%)

Rel. err. (%)

Wall-clock time (s)

Wall-clock time (s)

P2P, overlap Estimated Ttotal , 2, GPUs P2P, overlap Measured Ttotal, 2, GPUs

5123

P2P, overlap Estimated Ttotal , 4, GPUs P2P, overlap Measured Ttotal , 4, GPUs

(d) AMD MI250X (1 GCD), 4 GPUs.

101

−20

Grid Size (N 3)

20

Grid Size (N 3)

6 × 100 4 × 100 3 × 100 2 × 100

8963 10243 11523 12803

6 × 100 4 × 100 3 × 100 2 × 100

(c) AMD MI250X (1 GCD), 2 GPUs.

2 × 101

7683

101

Rel. err. (%)

Rel. err. (%)

Wall-clock time (s)

4 × 101 3 × 101 2 × 101

6403

(b) NVIDIA H100, 4 GPUs.

P2P, overlap Estimated Ttotal , 2, GPUs P2P, overlap Measured Ttotal , 2, GPUs

6 × 101

MAPE = 21.2% bias = +21.2%

100

50 0 −50

5123

P2P, overlap Estimated Ttotal , 4, GPUs P2P, overlap Measured Ttotal , 4, GPUs

101

Rel. err. (%)

Rel. err. (%)

Wall-clock time (s)

3 × 101 2 × 101

50 0

−50

5123

6403

7683

8963 10243 11523 12803

Grid Size (N 3)

(f) Intel 1550 (composite), 4 GPUs.

Figure 8: Comparison of measured vs. estimated wall-clock time for OpenMP Offloading {2,4}-GPU Version 4 across three GPU architectures: NVIDIA H100, AMD MI250X, and Intel 1550: model estimation (dashed) against measurement (solid). Lower panels: relative error (Testi − Tmeas )/Tmeas per grid size, with ±10% and ±20% guides; note the wider vertical scale in panel (f). Mean absolute percentage error and signed bias are annotated in each upper panel. 27

CUDA—such as thread blocks and grids—the performance does not match that of CUDA on NVIDIA, HIP on AMD, or SYCL on Intel. This performance discrepancy can be attributed to several factors: – CUDA, HIP, and SYCL effectively utilize the concept of SIMT, which optimizes GPU resource usage. – In contrast, we hypothesize that OpenMP Offloading primarily adheres to Single Instruction Multiple Data (SIMD) [27], a method that may not achieve optimal performance on GPUs. Although there is an option to manually set the thread grid and block in OpenMP Offloading, the current compilers associated with the OpenMP Offloading API do not execute this task as efficiently as those in CUDA, HIP, or SYCL for GPUs. • This performance gap is clearly illustrated in Figure 3, which shows that CUDA, HIP, and SYCL demonstrate significantly superior performance compared to OpenMP Offloading. Additionally, as depicted in Figure 8, CUDA, HIP, and SYCL also exhibit slightly higher measured bandwidth between the CPU and GPU, as well as among GPUs. This enhanced bandwidth contributes to their overall better performance relative to OpenMP Offloading. • The use of HHT significantly enhances the overall performance of multi-GPU OpenMP Offloading implementations. The LLVM backend enables asynchronous operations of OpenMP Offloading to be executed concurrently with the assistance of HHT threads. In this research, setting HHT to 8 resulted in improved performance; however, the HHT thread size can be adjusted according to the specific requirements of the problem at hand. • Our performance model scales well with the measured time; nevertheless, it does not account for penalties associated with kernel launches, synchronization, and communication latencies, resulting in some discrepancies. Despite this, our model provides reasonable scaling predictions across all three GPUs and can be applied to other similar problems. • For larger or more complex codebases, utilizing CUDA may prove advantageous, although it typically comes at the cost of increased development time. Furthermore, even if a conversion to another programming model is attempted, it may not yield significant performance benefits. OpenMP Offloading, on the other hand, offers a more user-friendly approach, making it easier for individuals with backgrounds in physics, mathematics, or engineering to adopt. As such, the same codebase can be executed without further modifications. • However, we believe that further compiler optimizations for kernel computations in OpenMP Offloading could enable performance levels comparable to those of CUDA on NVIDIA, HIP on AMD, and SYCL on Intel. Achieving such optimizations could significantly contribute to effective exascale computing that rivals other low-level GPU programming models.

8

Conclusions

In this research, we investigated four distinct OpenMP Offloading programming models targeting 2 and 4 GPUs, specifically those developed for NVIDIA, AMD, and Intel architectures. We compared these four OpenMP Offloading variants, operating on 2 and 4 GPUs, against native GPU programming models, including CUDA, HIP, and SYCL. Among the multi-GPU OpenMP Offloading configurations, the OpenMP Offloading {2,4}-GPU Version 4 exhibited the best performance compared to the other three variants. This performance enhancement can largely be attributed to several factors: the concurrent execution of GPU kernels using multiple host threads, the utilization of P2P data copies for halo data exchanges, and efficient workload distribution across the available GPUs. Notably, the OpenMP Offloading {2,4}-GPU Version 4 achieved approximately 2x speedup with 2 GPUs and around 4x speedup with 4 GPUs compared to the single GPU implementation of the OpenMP Offloading 1-GPU Baseline. The Intel configuration with 4 GPUs, however, achieved slightly less than the 4x speedup, at about 3.5x. 28

The proposed parallelization strategy achieved better performance under configurations with 2 and 4 GPUs, demonstrating effective communication overlap. This overlap improves the efficiency of our parallel implementation. Utilizing 8 HHT appears to be optimal; however, we anticipate that this number can be increased based on the number of asynchronous operations executed on a single GPU. Additionally, our bandwidth analysis informs the development of a performance model to evaluate the efficiency of our parallel implementations. While OpenMP Offloading did not achieve performance levels competitive with established native programming models such as NVIDIA’s CUDA, AMD’s HIP, and Intel’s SYCL, the significant advantage lies in having a single source code that is compatible across various GPU vendors. Adapting extensive scientific codes to different architectures can be a complex task. Even after conversion, the resulting code often requires substantial modifications to achieve satisfactory performance. Therefore, this work does not advocate for OpenMP Offloading as a replacement for CUDA, HIP, or SYCL. Instead, it outlines potential performance expectations for OpenMP Offloading in scientific domains, exemplified through a 3D heat transfer problem. Our parallelization methodology, enhanced by the use of multiple host threads, a performance model, and bandwidth measurements, provides a robust framework for achieving optimal performance in OpenMP Offloading for multi-GPU computations on a single compute node. This research offers clear insights and a deeper understanding of the application of OpenMP Offloading for multi-GPU programming on single compute nodes, paving the way for broader adoption of OpenMP Offloading in various scientific computing applications, independent of specific GPU architectures.

Author Contributions Conceptualization, methodology, software, validation, formal analysis, investigation, resources, data curation, writing—original draft preparation, writing—review and editing, visualization, supervision, project administration, and funding acquisition, Ezhilmathi Krishnasamy. The author has read and agreed to the published version of the manuscript.

Funding This research was carried out within the QCLS project—Quantum Computing for Scientific Computing (ID OP27.01207)—co-funded by the European Union through the European Regional Development Fund (ERDF) and by the Republic of Slovenia through the Ministry of Higher Education, Science and Youth, under the JR RZK call implemented by the Slovenian Research and Innovation Agency (ARIS).

Institutional Review Board Statement Not applicable.

Informed Consent Statement Not applicable.

Data Availability Statement The source code and supporting materials underlying the results presented in this research are publicly available on Zenodo and GitHub. The archival version corresponding to release v1.0.0 is available on Zenodo at https://doi.org/10.5281/zenodo.21878593, and the actively maintained version is available on GitHub at https://github.com/ezhilmathik/OpenMP-Offloading-Intranode-Artifact/.

29

Acknowledgments The author gratefully acknowledges the Leibniz Supercomputing Centre (LRZ) for providing computing resources on the SuperMUC-NG supercomputer. The author also gratefully acknowledges the EuroHPC Joint Undertaking (EuroHPC JU) for providing access to the MareNostrum 5 supercomputer at the Barcelona Supercomputing Center (BSC) and to the LUMI supercomputer through a Development Access call.

Conflicts of Interest The author declares no conflicts of interest. The funder had no role in the design of the research; in the collection, analysis, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Abbreviations The following abbreviations are used in this manuscript: API BSC CFD CPU GPU CU CUDA FLOP/s FP64 GB/s GFLOP/s GPU HIP HPC IF LUMI LUMI-G OpenACC OpenCL OpenMP P2P PCIe SM SYCL HtoD DtoH Meas.BW GCDs SIMT SIMD HHT k̄ βP2P

Application Programming Interface Barcelona Supercomputing Center Computational Fluid Dynamics Central Processing Unit Graphics Processing Unit Compute Unit Compute Unified Device Architecture Floating-point operations per second Double-precision, 64-bit floating-point format Gigabytes per second Billions of floating-point operations per second Graphics Processing Unit Heterogeneous-compute Interface for Portability High-Performance Computing Infinity Fabric Large Unified Modern Infrastructure GPU partition of the LUMI supercomputer Open Accelerators Open Computing Language Open Multi-Processing Peer-to-Peer Peripheral Component Interconnect Express Streaming Multiprocessor SYCL heterogeneous programming model Host-to-Device Device-to-Host Measured Bandwidth Graphics Compute Dies Single Instruction Multiple Threads Single Instruction Multiple Data Hidden Helper Threads Average Computational Constant Ratio P2P Measured Bandwidth

30

A

Artifact Description

A.1

Availability

All source code, run scripts, raw measurement data, and figure-generation scripts are publicly available at https://github.com/ezhilmathik/OpenMP-Offloading-Intranode-Artifact and permanently archived on Zenodo (DOI: 10.5281/zenodo.21878593, release tag: v1.0.0, commit: e41a873563d1). Source code is licensed under the European Union Public Licence v1.2 (EUPL-1.2); measurement data and generated figures are licensed under the Creative Commons Attribution 4.0 International licence (CC BY 4.0). Every directory ships both the sources and the results they produced. All figures and tables in this paper can be regenerated from the recorded data without access to the GPUs; re-running the measurements requires the hardware listed in Tables 2 and 7.

A.2

Repository structure

OpenMP_Single_Source/

nine OpenMP sources + serial CPU reference, compiled unmodified on all three vendors 3D_Heat_Benchmark/ main results; one directory per machine <MACHINE>/Source-Results/ OpenMP/ build/run/verify/aggregate pipeline + raw results HIP/ CUDA/ SYCL/ native-model counterpart, same pipeline model/ analytical performance model plots/ final figures HHT/ hidden-helper-thread sensitivity study <MACHINE>-{01,02,03}/ one full pipeline copy per HHT configuration gain.py cross-machine comparison and figures Bandwidth_Measurements/ intra-node bandwidth and latency microbenchmarks

Each directory carries its own README.md with the machine-specific modules, compiler flags, device masking and submission commands.

A.3

Experimental environment

Table 7: Hardware and software environment. One OpenMP device is one GCD on MI250X, one GPU on H100, and one whole card (two stacks, COMPOSITE hierarchy) on Max 1550. Compiler / runtime

Device selection

NVIDIA H100 (sm_90)

Cray cc (PrgEnv-amd), ROCm 6.4.4, LUMI/25.09 clang 18.1.8, CUDA 12.8, GCC 13.2.0

Intel Max 1550 (PVC)

icx (oneAPI)

ROCR_VISIBLE_DEVICES = 0 / 0,2 / 0,2,4,6 CUDA_VISIBLE_DEVICES = 0 / 0,1 / 0,1,2,3 ZE_AFFINITY_MASK = 0 / 0,1 / 0,1,2,3

AMD MI250X (gfx90a)

A.4

Problem parameters

Explicit finite-difference solve of ∂u/∂t = ∇ · (κ(x, y, z) ∇u) on an N 3 grid, decomposed in the slowestvarying dimension, with one halo-face exchange per internal boundary per timestep. The timestep is derived identically in submit.sh and verify.sh: h=

1 , N −1

∆t =

h2 , 12 κmax

T = Ns · ∆t,

(13)

with κmax = 0.95. The production sweep uses N ∈ {512, 640, 768, 896, 1024, 1152, 1280}, Ns = 500 timesteps, and five independent repetitions per grid size, each submitted as a separate Slurm job. Reported times are medians over the repetitions, taken from the solver’s own timer rather than the harness wall clock.

31

A.5

Reproducing the measurements

The pipeline is identical on every machine and in every model directory: cd 3 D _ H e a t _ B e nc h m a r k / < MACHINE >/ Source - Results / OpenMP ./ submit . sh 5 # one build job + 5 repetitions x 7 grid sizes squeue -u $USER # wait for completion ./ aggregate . sh # merge results_ */ into summary . csv ( medians ) ./ plot . sh -a # terminal summary of the aggregated data

The sweep can be narrowed with NS, FAMILIES, BASELINE and STEPS; SKIP_BUILD=1 reuses existing binaries. The first argument of submit.sh is the repetition count, not a grid size. Correctness is checked separately, because timing runs never write solution files: sbatch -- export = ALL , N =128 , STEPS =50 verify . sh

This builds the serial CPU reference and compares every variant against it element-wise. Timings produced by verify.sh are not meaningful, as file I/O falls inside the measured region.

A.6

Mapping results to the artifact Table 8: Source files associated with each reported result.

Result

Source

Figure 3 (OpenMP Off. 2,4-GPU Versions) Figure 4 (Bandwidth benchmark) Figure 6 (Machines, HHTs, and N 3 ) Figure 7 (HHT recommendation) Figure 8 (Model versus measurements)

3D_Heat_Benchmark/*/Source-Results/*/summary.csv via speedup.py Bandwidth_Measurements/ HHT/plots-approach/OpenMP/ allN_gain_OpenMP.pdf via gain.py HHT/plots-approach/OpenMP/ decision_OpenMP.pdf via gain.py 3D_Heat_Benchmark/*/ Source-Results/model/plot_p2p.py

The HHT figures are regenerated with: cd HHT python3 gain . py -- list python3 gain . py -- format pdf

A.7

# verify discovery of the nine directories

Resource requirements

A full sweep is 1 build job plus 35 run jobs per machine and model directory, each requesting one exclusive node for up to 90 minutes. Aggregation and plotting run on a laptop in seconds and require numpy, pandas and matplotlib; the terminal summaries produced by plot.sh.

B

Schematic Execution Workflows

To improve understanding of the pseudocode and parallelization strategies used in this research with various versions of OpenMP Offloading, the following figures provide a schematic overview of each method.

32

OpenMP Offloading 2-GPU Version 1 A single host thread with synchronized host-staged data copy. INITIALIZATION HOST SETUP Initialize global domain Split domain into D0 and D1

GPU 0 Receive D0

GPU 1 Receive D1 HOST Host allocates halo buffers H0 and H1

TIME-STEPPING LOOP For step = 0, . . . , Nsteps − 1

GPU 0

GPU 1

INTERIOR(D0 ) BOUNDARY(B0 )

INTERIOR(D1 ) BOUNDARY(B1 )

omp taskwait Wait for both GPUs

HALO EXCHANGE 0: GPU 0 → HOST → GPU 1

GPU 0 Boundary B0

device-to-host

HOST Buffer H0

host-to-device

GPU 1 Ghost G1

host staging

next step

HALO EXCHANGE 1: GPU 1 → HOST → GPU 0

GPU 1 Boundary B1

device-to-host

HOST Buffer H1

host-to-device

GPU 0 Ghost G0

host staging

ADVANCE TIME LEVEL GPU 0 new SWAP(uold 0 , u0 )

GPU 1 new SWAP(uold 1 , u1 ) Repeat next time step

FINALIZATION

after the final step

GPU 0 Local solution D0

GPU 1 Local solution D1

HOST Copy local solutions Reconstruct global domain

Figure 9: Synchronous OpenMP Offloading 2-GPU Version 1 with host-staged halo exchange.

33

OpenMP Offloading 2-GPU Version 2 A single host thread with asynchronous host-staged data copy. INITIALIZATION HOST SETUP Initialize global domain Split domain into D0 and D1

GPU 0 Receive D0

GPU 1 Receive D1 HOST Host allocates halo buffers H0 and H1

tk0→1 , dep0→1 depend(inout: token)

tk1→0 , dep1→0 depend(inout: token)

TIME-STEPPING LOOP For step = 0, . . . , Nsteps − 1

GPU 0

GPU 1

INTERIOR(D0 ) BOUNDARY_ASYNC( B0 , tk0→1 )

INTERIOR(D1 ) BOUNDARY_ASYNC( B1 , tk1→0 )

The host thread enqueues both GPU computations without a global wait.

ASYNCHRONOUS HALO CHAINS Dependency chain dep0→1

GPU 0 Boundary B0

device-to-host

HOST Buffer H0

host-to-device

GPU 1 Ghost G1

host staging

Dependency chain dep1→0 next step

GPU 1 Boundary B1

device-to-host

HOST Buffer H1

host-to-device

GPU 0 Ghost G0

host staging

omp taskwait Complete both interior kernels and both halo chains

ADVANCE TIME LEVEL GPU 0 new SWAP(uold 0 , u0 )

GPU 1 new SWAP(uold 1 , u1 ) Repeat next time step

FINALIZATION

after the final step

Destroy dep0→1

Destroy dep1→0

GPU 0 Local solution D0

GPU 1 Local solution D1

HOST Copy local solutions Reconstruct global domain

Figure 10: Asynchronous OpenMP Offloading 2-GPU Version 2 with one host thread, dependencycontrolled asynchronous data transfers, and host-staged halo exchange. Dashed arrows denote asynchronously enqueued operations or dependency-controlled execution.

34

OpenMP Offloading 2-GPU Version 3 Each host thread per GPU with asynchronous host-staged data copy. INITIALIZATION HOST SETUP Initialize global domain Split domain into D0 and D1

GPU 0 Receive D0

GPU 1 Receive D1 HOST Host allocates halo buffers H0 and H1

tk0→1 , dep0→1 depend(inout: token)

tk1→0 , dep1→0 depend(inout: token)

TWO-HOST-THREAD PARALLEL REGION omp parallel num_threads(2) For step = 0, . . . , Nsteps − 1

HOST THREAD 0 −→ GPU 0

HOST THREAD 1 −→ GPU 1

Host Thread 0 Enqueue GPU 0 compute and 0 → 1 halo chain

Host Thread 1 Enqueue GPU 1 compute and 1 → 0 halo chain

GPU 0 BOUNDARY_ASYNC( B0 , tk0→1 )

GPU 0 INTERIOR(D0 ) Interior computation continues while the boundary and halo chain execute

GPU 1 BOUNDARY_ASYNC( B1 , tk1→0 )

device-to-host

device-to-host

HOST Buffer H0

HOST Buffer H1

host-to-device

host-to-device

GPU 1 Ghost G1

GPU 0 Ghost G0

GPU 1 INTERIOR(D1 ) Interior computation continues while the boundary and halo chain execute next step

COMPUTE–COMMUNICATION OVERLAP WINDOW Independent host threads allow both halo chains to progress while the GPU interior kernels execute.

omp barrier Both host-thread lanes complete their compute and communication chains

ADVANCE TIME LEVEL omp single One host thread advances both GPU time levels GPU 0 new SWAP(uold 0 , u0 )

GPU 1 new SWAP(uold 1 , u1 ) Repeat next time step

FINALIZATION

after the final step

Destroy dep0→1

Destroy dep1→0

GPU 0 Local solution D0

GPU 1 Local solution D1

HOST Copy local solutions Reconstruct global domain

Figure 11: Asynchronous OpenMP Offloading 2-GPU Version 3 with two host threads and hoststaged halo data exchange. Each host thread independently controls one GPU and enqueues its local computation and outgoing halo data chain. This provides an overlap window in which boundary computation and halo data (communication) can be hidden behind the interior computation.

35

OpenMP Offloading 2-GPU Version 4 Each host thread per GPU with asynchronous P2P data copy between GPUs. INITIALIZATION HOST SETUP Initialize global domain Split domain into D0 and D1

GPU 0 Receive D0

GPU 1 Receive D1 DIRECT P2P EXCHANGE No host halo buffers are required

tk0→1 , dep0→1 depend(inout: token)

tk1→0 , dep1→0 depend(inout: token)

TWO-HOST-THREAD PARALLEL REGION omp parallel num_threads(2) For step = 0, . . . , Nsteps − 1

HOST THREAD 0 −→ GPU 0

HOST THREAD 1 −→ GPU 1

Host Thread 0 Host Thread 1 Enq. GPU 0 computation and 0 → 1 peer-copy chainEnq. GPU 1 computation and 1 → 0 peer-copy chain 0 → GPU 1 GPU 0 direct GPU dep0→1 GPU 1 BOUNDARY_ASYNC( Ghost region G1 B0 , tk0→1 )

GPU 0 INTERIOR(D0 )

GPU 1 INTERIOR(D1 )

PEER-TO-PEER DEVICE TRANSFERS Interior computation Interior computation Halo data moves directly between GPU address spaces; continues while the continues while the no device-to-host or host-to-device staging is performed. boundary and direct boundary and direct GPU 1 GPU 0 peer-copy chain BOUNDARY_ASYNC( peer-copy chain Ghost region Gdirect GPU 1 → GPU 0 0 execute execute B , tk ) dep1→0

1

1→0

next step COMPUTE–P2P COMMUNICATION OVERLAP WINDOW Independent host threads allow the boundary kernels and direct peer copies to progress while the interior kernels execute.

omp barrier Complete both interior computations, boundary computations, and peer-copy chains

ADVANCE TIME LEVEL omp single One host thread advances both GPU time levels GPU 0 new SWAP(uold 0 , u0 )

GPU 1 new SWAP(uold 1 , u1 ) Repeat next time step

FINALIZATION

after the final step

Destroy dep0→1

Destroy dep1→0

GPU 0 Local solution D0

GPU 1 Local solution D1

HOST Copy local solutions Reconstruct global domain

Figure 12: Asynchronous OpenMP Offloading 2-GPU Version 4 with two host threads and direct P2P halo exchange. Each host thread independently controls one GPU and enqueues its local computation and outgoing halo transfer. Boundary data is copied directly from the source GPU to the peer GPU’s ghost region, eliminating host staging. The P2P data copy and boundary computations may be hidden behind the corresponding interior computations.

36

References [1] NVIDIA Corporation. CUDA C++ Programming Guide. https://docs.nvidia.com/cuda/ cuda-c-programming-guide/, 2026. Release 13.3. Accessed: 2026-07-19. [2] Khronos SYCL Working Group. SYCL 2020 Specification (Revision 11). https://registry. khronos.org/SYCL/specs/sycl-2020/html/sycl-2020.html, 2025. Published 2025-11-07. Accessed: 2026-07-19. [3] Khronos OpenCL Working Group. The OpenCL Specification, Version 3.0. https://registry. khronos.org/OpenCL/specs/3.0-unified/html/OpenCL_API.html, 2025. Revision 3.0.19. Accessed: 2026-07-19. [4] Advanced Micro Devices, Inc. HIP Programming Guide (ROCm Documentation). https://rocm. docs.amd.com/en/latest/how-to/programming_guide.html, 2026. Heterogeneous-Compute Interface for Portability. Accessed: 2026-07-19. [5] OpenACC Organization. OpenACC Specification and Resources. https://www.openacc.org/ specification, 2025. Accessed: 2026-07-19. [6] OpenMP Architecture Review Board. OpenMP 6.0: Offloading and Device Features. https: //www.openmp.org/articles/openmp-6/, 2024. Accessed: 2026-07-19. [7] OpenACC Organization. The OpenACC Application Programming Interface, Version 3.4. https://www.openacc.org/sites/default/files/inline-images/Specification/ OpenACC-3.4.pdf, 2025. Accessed: 2026-07-19. [8] Erik Lindholm, John Nickolls, Stuart Oberman, and John Montrym. Nvidia tesla: A unified graphics and computing architecture. IEEE Micro, 28(2):39–55, 2008. doi: 10.1109/MM.2008.31. [9] Hewlett Packard Enterprise. Introduction to OpenACC. HPE Cray Programming Environment, version 24.03. Available online: https://cpe.ext.hpe.com/docs/24.03/cce/man7/intro_openacc. 7.html, 2024. (accessed on 14 July 2026). [10] Ezhilmathi Krishnasamy, James Trotter, Xing Cai, Dirk Pleiter, Leon Kos, Laura Saavedra, and Pascal Bouvry. Performance and programmability of MPI+X integration with CUDA, HIP, SYCL, OpenACC, and OpenMP offloading for supercomputing: A case study on dense matrix–vector multiplication. In Proceedings of the Supercomputing Asia and International Conference on High Performance Computing in Asia Pacific Region Workshops, pages 457–468, Osaka, Japan, 26–29 January 2026. doi: 10.1145/3784828.3786264. [11] Ezhilmathi Krishnasamy, Sven Harig, and Pascal Bouvry. Accelerating tsunami computation on a GPU using OpenMP offloading. In Computer and Information Science and Engineering, volume 1156, pages 163–178. Springer, Cham, Switzerland, 2024. doi: 10.1007/978-3-031-57037-7_12. [12] Yehonatan Fridman, Guy Tamir, and Gal Oren. Portability and scalability of OpenMP offloading on state-of-the-art accelerators. In High Performance Computing, volume 13999, pages 378–390, Hamburg, Germany, 21–25 May 2023. doi: 10.1007/978-3-031-40843-4_28. [13] Ezhilmathi Krishnasamy and Pascal Bouvry. OpenMP offloading on AMD and NVIDIA GPUs: Programmability and performance analysis. In Proceedings of the 2025 9th International Conference on High Performance Compilation, Computing and Communications, pages 44–56, Jinan, China, 27–29 August 2025. doi: 10.1145/3774949.3774956. [14] Ezhilmathi Krishnasamy and Pascal Bouvry. Comparative performance analysis of CUDA and OpenMP offloading for BLAS operations on GPU. In ICT for Intelligent Systems: Proceedings of ICTIS 2025, volume 126, pages 9–35, New York, NY, USA, 23–24 May 2025 2026. doi: 10.1007/978-981-95-1361-1_2. 37

[15] Manuel Ferat, Romain Pereira, Adrien Roussel, Patrick Carribault, Luiz-Angelo Steffenel, and Thierry Gautier. Enhancing MPI+OpenMP task based applications for heterogeneous architectures with GPU support. In Michael Klemm, Bronis R. de Supinski, Jannis Klinkenberg, and Brandon Neth, editors, OpenMP in a Modern World: From Multi-Device Support to Meta Programming (IWOMP 2022), volume 13527 of Lecture Notes in Computer Science, pages 3–16, Cham, Switzerland, 2022. Springer. doi: 10.1007/978-3-031-15922-0_1. [16] Shilei Tian, Johannes Doerfert, and Barbara Chapman. Concurrent execution of deferred OpenMP target tasks with hidden helper threads. In Barbara Chapman and José Moreira, editors, Languages and Compilers for Parallel Computing (LCPC 2020), volume 13149 of Lecture Notes in Computer Science, pages 41–56, Cham, Switzerland, 2022. Springer. doi: 10.1007/978-3-030-95953-1_4. [17] Ezhilmathi Krishnasamy and Pascal Bouvry. A study of GPU programming paradigms: OpenACC vs. CUDA in linear algebra computations. In Proceedings of the 2025 5th International Conference on Advances in Electrical, Electronics and Computing Technology (EECT), Guangzhou, China, 21–23 March 2025. doi: 10.1109/EECT64505.2025.10966957. [18] NVIDIA Corporation. NVIDIA H100 Tensor Core GPU. Available online: https://www.nvidia. com/en-us/data-center/h100/, 2026. (accessed on 7 July 2026). [19] Barcelona Supercomputing Center. MareNostrum 5. Available online: https://www.bsc.es/ marenostrum/marenostrum-5, 2026. (accessed on 7 July 2026). [20] Advanced Micro Devices, Inc. AMD Instinct MI250X Accelerator. Available online: https: //www.amd.com/en/products/accelerators/instinct/mi200/mi250x.html, 2026. (accessed on 7 July 2026). [21] LUMI Consortium. LUMI-G hardware. Available online: https://docs.lumi-supercomputer. eu/hardware/lumig/, 2026. (accessed on 7 July 2026). Intel Data Center GPU Max 1550—product specifications. Avail[22] Intel Corporation. able online: https://www.intel.com/content/www/us/en/products/sku/232873/ intel-data-center-gpu-max-1550/specifications.html, 2026. (accessed on 7 July 2026). [23] Leibniz Supercomputing Centre. Hardware of SuperMUC-NG phase 2. Available online: https: //doku.lrz.de/hardware-of-supermuc-ng-phase-2-222891050.html, 2026. (accessed on 7 July 2026). [24] Ezhilmathi Krishnasamy. Hybrid CPU–GPU parallel simulations of 3D front propagation. Master’s thesis, Linköping University, Linköping, Sweden, 2014. Available online: https://urn.kb.se/ resolve?urn=urn:nbn:se:liu:diva-114935 (accessed on 14 July 2026). [25] Ezhilmathi Krishnasamy, Mohammed Sourouri, and Xing Cai. Multi-GPU implementations of parallel 3D sweeping algorithms with application to geological folding. Procedia Comput. Sci., 51: 1494–1503, 2015. doi: 10.1016/j.procs.2015.05.339. [26] Ezhilmathi Krishnasamy, James D. Trotter, Xing Cai, and Pascal Bouvry. Performance evaluation of the openacc api targeting multi-gpu data transfers: A case study of porting parallel 3d sweeping algorithms for geological folding. Manuscript submitted for publication (under review), 2025. [27] Michael J. Flynn. Some computer organizations and their effectiveness. IEEE Transactions on Computers, C-21(9):948–960, September 1972. doi: 10.1109/TC.1972.5009071.

38

Related documents

• Accelerating Intra-Node GPU-to-GPU Communication Through Multi-Path Transfers with CUDA Graphs#134524
• ShuntServe: Cost-Efficient LLM Serving on Heterogeneous Spot GPU Clusters#287102
• PipeMax: Enhancing Offline LLM Inference on Commodity GPU Servers#155244
• Coral: Cost-Efficient Multi-LLM Serving over Heterogeneous Cloud GPUs#168288
• Performance Analysis of Low-Order, GPU-accelerated Finite Element Kernels using Kokkos#1049991