The starting point for a KF track finding algorithm is a track “seed”, an initial guess at the track state. For the full mkFit algorithm, the input track seeds are built in an upstream algorithm using three or four hits from the innermost layers of the CMS detector. For simplicity, the p2z and p2r benchmarks use an artificial standalone input consisting of a single track seed. The parameters of the initial seed are smeared according to a Gaussian distribution, in order to prevent the algorithm from performing identical numerical operations for each track.

The tracks are built by first propagating the initial seed parameters to the next layer. In the full combinatorial mkFit application, the next step is to search for compatible hits on that layer. In the p2z and p2r benchmarks, only one hit per layer is considered, similar to what is required for track fitting. The hit parameters—the hit locations and uncertainties on each layer—are smeared per track according to the same procedure as the input track parameters.

The total number of tracks generated with this procedure is split into “events” with a fixed number of tracks per event.

For the GPU implementations of the benchmark, the track propagation and Kalman update steps are run as a single GPU kernel. The data for all tracks and hits are prepared on the CPU, transferred to the GPU for the computations in two bulk transfers, and the output data are transferred back to the CPU.

The “propagation to z” or p2z benchmark uses the expected helical trajectory of the charged particle to calculate the track parameters and the covariance matrix on endcap layer N + 1. In contrast, the “propagation to r” or p2r benchmark uses an iterative approach to propagation, advancing the track state from the initial barrel radius to the final radius in discrete steps. In practice, both approaches involve a series of operations involving sin() and cos() functions, as well as matrix multiplication of up to 6x6 matrices. The two benchmarks are expected to be similar in terms of arithmetic intensity.

The second task in both benchmarks is the Kalman update step, which “updates” the track state using the parameters of a hit on that layer. The hit parameters include three coordinates for the hit position and a 3x3 covariance matrix. Similar to the propagation step, the update step involves small matrix multiplication and matrix inverse operations as well as trigonometric functions.

Both benchmarks employ an Array-Of-Structures-Of-Arrays (AOSOA) data structure, as illustrated in Figure 1. Each benchmark runs over a fixed number of events (nevts) and processes a fixed number of tracks in each event (ntrks). Within one event, tracks are grouped into batches of bsize tracks, and each batch of tracks is put into a Structure-Of-Arrays construct referred to as an MPTRK. Here, as in the full mkFit algorithm, the goal of organizing tracks into different batches is to enable SIMD operations across batch elements. The value of bsize can be optimized for different platforms; for example, on a GPU it might be the NVIDIA warp size of 32, while on a CPU it might be a multiple of the AVX-512 vector width of 16; for consistency we use 32 everywhere.

To maximize the opportunities for SIMD operations, the storage order in the AOSOA follows the same general scheme as in mkFit. Within an MPTRK, the first elements of bsize different vectors (or matrices) get stored in adjacent locations, followed by the second elements of bsize different vectors (or matrices), and so on, until the full structure representing a batch of bsize tracks is completely populated. Then the next MPTRK structure is written into memory, and the next, until all ntrks tracks (8192 for p2r and 9600 for p2z) for the first event are present in memory. This first event corresponds to the first row of SOAs in Figure 1. It is followed by a second row of SOAs for the second event, etc. A similar memory layout applies to the hit data.

Our reference CPU implementation is based on the oneAPI Threading Building Blocks (or oneTBB, or just TBB) library. TBB is a template library originated by Intel for parallel programming on multi-core processors that simplifies the definition and optimization of threaded applications. For a given executable, TBB collects all tasks that are specified to run in parallel and the library manages and schedules threads to execute them. TBB is used as the CPU thread manager in the software framework of the CMS experiment.

In our code, nested TBB parallel_for loops are over events and batches of tracks. Each batch of tracks is then vectorized, so that tracks in the batch are processed in a SIMD fashion. Vectorization is implemented following the approach used in mkFit, where groups of matrices are processed concurrently and loops over the matrix index in the group are decorated with omp simd pragmas. These pragmas were activated by the -fopenmp compiler option.

Our reference GPU implementation is based on the CUDA programming model, which is a multi-threaded SIMD model for general purpose GPU programming, introduced by NVIDIA. The CUDA implementation is ported from the TBB version, which shares the same AOSOA data structure for input data. The main difference in CUDA is that each MPTRK is processed by a block of GPU threads, and each thread processes the computation for one track in the MPTRK. Since the computation of each track is independent of the others, we find keeping the intermediate results in the local registers to have the most efficient memory access. We explored using shared memory to store the intermediate results within an MPTRK for all the threads in the block, but it was shown to have significantly lower memory throughput in a detailed profiling study.

In relation to other portability technologies, CUDA provides a level of abstraction similar to the general accelerator execution model and memory model that is also employed by general GPU programming models such as OpenCL and SYCL. As a proprietary NVIDIA GPU programming model, however, it exposes several NVIDIA-GPU-specific features, which allows available architecture-specific features to be fully exploited, but it also means that the code is not portable across heterogeneous accelerators.

HIP is the vendor-supported, native programming model for AMD GPUs and is designed to be portable across NVIDIA and AMD GPUs. It is also designed to be syntactically similar to CUDA so that most API calls can be simple translations of names. In the case of p2z and p2r, the kernels only rely on the core functionalities of CUDA, such as memory allocation and kernel dispatches, which are supported in HIP. This leads to a straightforward port to the HIP version starting from the CUDA version.

Directive-based, high-level programming models such as OpenMP and OpenACC use a set of directives (a special type of comments that a compiler can understand) that allow a programmer to provide the compilers with important characteristics of an application, such as available parallelism and data sharing/mapping rules, so that much of the low-level programming and optimization burdens are automatically handled by the directive compilers.

The biggest advantage of the directive programming model is that it allows incremental porting of existing sequential CPU applications written in standard programming languages such as C, C++, and Fortran to parallel versions that can offload work to heterogeneous accelerators, without requiring major changes in the existing code structures. The initial OpenMP version was created by converting the reference TBB CPU implementation into an OpenMP CPU implementation, which is relatively straightforward due to the similarity between the TBB parallel syntax (using lambdas) and OpenMP parallel constructs. Then, the GPU offloading version was created by extending the OpenMP CPU implementation with OpenMP target-offloading directives and data mapping directives. The CPU version and the GPU version have different parallelism and data mapping strategies. For instance, on GPUs, the OpenMP target and data mapping directives are essential, but when targeting CPUs, they are unnecessary since the original host data and corresponding device data will share storage; in the latter case, how the OpenMP compiler implements those unnecessary directives on the CPUs is implementation-defined. Another issue is that on CPUs, team-level parallelism may or may not be ignored depending on the compiler implementations, possibly resulting in additional overheads due to loop serialization.

The initial conversion from the OpenMP GPU implementation into the OpenACC implementation was straightforward since both OpenMP and OpenACC are directive-based accelerator programming models and provide very similar execution and memory models. However, the main issue in converting between OpenMP and OpenACC was that different OpenMP/OpenACC compilers may choose different parallelism mapping strategies and vary in terms of maturity and supported features (see Section 5.1.1 for results and a more detailed discussion).

Alpaka (Matthes et al., 2017) is a single-source, header-only C++ parallelization library. The API level is similar to CUDA, with an abstraction layer added between the application and the vendor-specific programming models to achieve portability. For example, the kernel functions in Alpaka are templated with an accelerator type, which is resolved at compile time for different execution backends of the kernel. One difference between Alpaka and CUDA is that Alpaka has an additional abstraction level called elements in the parallel hierarchy model, where multiple elements can be processed in a thread. Having an explicit level allows compilers to take advantage of the SIMD vector registers when compiling for the CPU backends. Each SOA (MPTRK) is processed by a block, and the thread/element level is mapped differently between CPU and GPU backends to take full advantage of the parallel hierarchy. For GPU backends, blocks of bsize threads are assigned to process each MPTRK, whereas one thread with bsize elements is assigned to process each MPTRK for CPU backends. In each case, enough threads are defined so that all ntrks tracks are processed for an event. Profiling results confirm the use of vectorized instructions in the regions where the original mkFit implementation was also able to vectorize.

The overall conversion from CUDA to Alpaka is relatively smooth due to the similarity between the two programming models, except that the heavy use of templating often leads to a more verbose coding style and convoluted error messages during debugging. Nevertheless, Alpaka versions are often able to produce close-to-native performance after some effort of optimization. A particular relevant note for HEP experiment is that CMS has chosen to use Alpaka as its supported portability solution for the GPU usage in LHC Run 3 (Kortelainen et al., 2021; Bocci et al., 2023).

Similar to Alpaka, Kokkos (Edwards et al., 2014; Trott et al., 2022) serves as a single-source C++ template metaprogramming (TMP) library, intended to achieve architecture agnosticism and alleviate programmers from the complexities of vendor- or target-specific programming paradigms and heterogeneous hardware architectures. By embracing TMP methodologies, Kokkos facilitates device-specific code generation and optimizations via template specialization. To cater to diverse computing environments, Kokkos offers multiple device-specific backends, implemented as template libraries atop various HPC programming models like CUDA, HIP, OpenMP, HPX, SYCL, and OpenACC. These backends are tailored to adhere to advancements in the C++ standard, ensuring compatibility and efficacy.

A notable departure from Alpaka is Kokkos' emphasis on descriptive rather than prescriptive parallelism. Kokkos prompts developers to articulate algorithms in general parallel programming concepts, which are subsequently mapped to hardware by the Kokkos framework. The Kokkos programming model revolves around two fundamental abstractions: the first being the user data abstraction (Kokkos::View), a template library facilitating the representation of multidimensional arrays while managing efficient data layout for both CPU and GPU. The second abstraction revolves around parallel execution patterns (parallel_for, parallel_reduce, and parallel_scan), which can be executed under three distinct execution policies: RangePolicy for mapping single parallel loops, MDRangePolicy for mapping directly nested parallel loops, and TeamPolicy for hierarchical mapping of multiple nested parallel loops.

In this study, Kokkos versions were developed by translating CUDA code to Kokkos using the Kokkos View and parallel dispatch abstractions. The initial translation process was relatively straightforward, owing to the striking similarities between the Kokkos and CUDA execution models and memory models, except for the added complexity due to Kokkos-specific restrictions on C++ template programming. Efficient execution on both CPU and GPU took further optimization and was achieved by configuring the teamSize and vectorSize in the TeamPolicy used in the parallel_for execution pattern, as the defaults were found not to be optimal.

The C++ programming language is often the preferred choice for implementing high performance scientific applications. The recent revisions of the ISO C++ standard introduced a suite of algorithms capable of being executed on accelerators. Although this approach may not yield the best performance, it can present a viable balance between code productivity and computational efficiency. Numerous production-grade compilers are available, such as clang and its variants from various providers, or the recently released NVHPC from NVIDIA.

With the introduction of the C++17 standard, the C++ Standard Template Library (STL) underwent a substantial overhaul of its suite of algorithms, now updated with execution policies to adapt across various computing architectures, including multi-core x86 systems and GPUs. These parallel algorithms extended most of the existing STL algorithms with an additional argument, which is an execution policy. The policy enables programmers to specify the intended parallelism of an algorithm, which can result in performance improvements for the computational task. In particular, the execution policies in C++17 include:

The first option forces the algorithm to run sequentially, while the remaining three options allow the algorithm to be vectorized or run in parallel (with additional vectorization). Both the p2z and p2r stdpar implementations use the std::execution::par_unseq execution policy. Currently, only the nvc++ compiler offers support for stdpar algorithms to be offloaded on NVIDIA GPUs. It leverages CUDA Unified Memory to handle automatic data movement between CPU and GPU. On the systems that do not support Heterogeneous Memory Management, only data that is dynamically allocated in CPU code compiled by nvc++ is automatically managed, whereas memory allocated in GPU code is exclusively for GPU use and remains unmanaged. Thus, on such systems, CPU and GPU stack memory, along with global objects, are outside nvc++'s automatic management scope. Even data allocated on the CPU heap outside units compiled by nvc++ is not managed. When dealing with parallel algorithm invocations, pointers and objects must refer to data on the managed CPU heap to avoid errors. Any dereferencing of pointers to CPU stack or global objects in GPU code can lead to memory violations. These aspects encapsulate the nvc++ compiler's precise approach to memory management across CPU and GPU, emphasizing careful allocation and reference handling to ensure efficient operations. In our experience, developing code for this application is largely similar to standard C++ programming, with the primary distinction being the need to consider the previously mentioned limitations.

SYCL represents a cross-platform abstraction layer that enables code for heterogeneous processors to be written in a “single-source” style using completely standard C++. This approach aims to enhance the efficiency and accessibility of programming for a variety of compute architectures. While the SYCL programming language was promoted by the Khronos Group, it is predominantly advocated for by Intel, so the primary focus is on the optimization for Intel GPUs.

One of the key advantages of SYCL is the ability to handle regular C++ code for the host CPU and a subset of C++ for the device code within the same source file. This ability paves the way for an integrated and simplified development process. It also enables parallelism and the usage of memory hierarchies through a class template library. This effectively allows the expression of parallel STL, thus further integrating SYCL with standard C++ features. Since it was designed to be fully compatible with standard C++, it allows developers to utilize any C++ library within a SYCL application. This compatibility with standard C++ makes SYCL a versatile tool for developers. Moreover, SYCL's design prioritizes performance portability, aiming to provide high performance across a wide range of hardware architectures. Its abstractions are constructed to allow optimization but do not require a particular architecture or kernel language.

For the benchmarks described in this paper, the programming approach for SYCL is nearly a direct replication of the CUDA approach. In our implementation, we utilized SYCL's Unified Shared Memory feature for data management. To compile SYCL for NVIDIA and AMD GPUs, we compiled the dpc++ tool chain following the instructions in Intel's open-source llvm-sycl repository.¹

The throughput of the p2z and p2r benchmarks were measured using two different systems with two different NVIDIA GPUs. To test the p2r implementations, we used the Joint Laboratory for System Evaluation (JLSE), a collection of HPC testbeds hosted at Argonne National Lab. The NVIDIA GPU that was used for testing is an A100 GPU with an AMD 7532 CPU as the host machine. For the p2z implementations, the measurements were performed on a test node for the Summit HPC system at Oak Ridge National Laboratory. Each Summit node includes six NVIDIA V100 GPUs and two IBM Power9 CPUs, although only one GPU was utilized in our tests.

The throughput measurements on NVIDIA GPUs are shown in Figure 2 for both the p2z and p2r benchmarks. In both cases, the native CUDA implementation is used as the reference implementation, and the measurement time includes the kernel execution only. Several different compilers are used for the different implementations: the CUDA, Alpaka, and Kokkos implementations are compiled with the nvcc compiler; OpenMP and OpenACC are compiled with the OpenARC (Lee and Vetter, 2014) compiler; the stdpar versions are compiled with nvc++; and for SYCL, the intel/llvm-sycl repository is used to compile dpc++ tool chains to be compatible with NVIDIA/AMD GPUs. For evaluation, whenever possible, we use the same launch parameters, including number of blocks and number of threads per block. In the OpenACC and OpenMP versions, different compilers varied in how they enforced the user-specified configurations (see Section 5.1.1). The launch parameters need to be manually specified for Alpaka and Kokkos, otherwise the libraries choose suboptimal values and the performance is about 30% worse, taking those versions even further from the native CUDA version. Setting the number of registers per thread is another approximately 10% effect. These parameters cannot be manually specified in the stdpar implementation. In the SYCL implementations, we specified the execution volume manually but relied on the defaults for the other launch parameters.

In general, most of the different portability solutions managed to produce close-to-native performance. The stdpar implementation did not perform well for either benchmark, mainly because the stdpar implemented in the nvc++ compiler relies on CUDA Unified Memory for all data movements between CPU and GPU memory, which fetches data needed by the GPU kernel on demand, exposing memory transfer overheads to the kernel execution. CUDA Unified Memory provides APIs to prefetch or migrate data between CPU and GPU memory to hide or reduce the memory transfer overheads, but the current stdpar does not include such functionalities. In order to mitigate the effects of data transfers while measuring the kernel execution time, we introduced implicit prefetching using the parallel std::copy algorithm. For the p2r benchmark, the worst performing version is the SYCL implementation. Detailed profiling using NVIDIA NSight Compute shows significant branching when using SYCL, but preliminary investigations have not revealed an obvious explanation for the branching.

For the p2z benchmark, most of the different implementations achieved similar performance except for the stdpar version. While similar to other p2z implementations, the relative performance of p2z Kokkos version is lower than that of the p2r Kokkos version, when compared to their corresponding CUDA versions. Detailed profiling shows that the Kokkos version of p2z uses more registers than the CUDA version, while the Kokkos version of p2r uses a similar number of registers compared to the CUDA version. Further investigations of the branching of the p2r SYCL version and of the difference in register usage in the p2z and p2r Kokkos version are beyond the scope of the present paper.

Figure 3 shows the throughput of the p2z implementations on the NVIDIA V100 GPUs, this time including both the kernel execution time and the memory transfer times. The transfer times are generally 2 to 5 times larger then the kernel execution times, which means that much of the variability between implementations is concealed. With the exception of stdpar, all implementations have close to identical performance because they provide memory management features to facilitate explicit data transfers, and in some cases memory prepinning as discussed in Section 5.1.2.

Support for other GPU architectures is in general less mature than the support for NVIDIA GPUs, although this is an area of rapid expansion. We explored the preliminary performance of the p2r benchmark on both AMD and Intel GPUs, and the kernel-only throughput measurements are shown in Figure 6. This is an out-of-the-box comparison of each tool's portability; no dedicated effort was made to optimize for AMD or Intel GPU architectures. For example, the HIP implementation for the AMD GPU is a carbon copy of the native CUDA version. For both Alpaka and Kokkos, switching backends is relatively seamless and does not require any code changes.

The results on the AMD GPU are shown in the left plot of Figure 6. The AMD GPU tests were performed on the JLSE testbed, which includes two AMD EPYC 7543 32c (Milan) CPUs and four AMD MI100 32GB GPUs. Only one GPU was used to perform the measurements. Both Kokkos and Alpaka include HIP backends which achieve reasonable performance: Alpaka actually outperforms the HIP version that was ported from CUDA, and Kokkos is within a factor of about 2. The same launch parameters as the CUDA implementation are used for the native HIP, Alpaka:HIP, and Kokkos:HIP measurements.

The throughput measurements on an Intel A770 GPU are shown in the right plot of Figure 6. Since the A770 is not an HPC-class GPU, all calculations were converted to single-precision operations. Relying on double-precision emulation results in performance that is 3 to 30 times slower, depending on the implementation. The SYCL backend for Alpaka has only experimental support (introduced in v0.9.0) and was not tested here. The Kokkos SYCL backend is still under active development: we observed a factor of 2 improvement in the throughput of the Kokkos:SYCL implementation when we updated from Kokkos 3.6.1 to Kokkos 4.1.0.

The original mkFit application is parallelized using the Threading Building Blocks (TBB) library from Intel, so we also used TBB as the reference native implementation on the CPU. This implementation is multithreaded and vectorized. It is worth noting that the original version was initially developed based on the “classic” Intel C++ Compiler (icc version 19), which led to improved vectorization performance compared to more recent compilers, resulting in p2z execution times that were approximately 2.7x faster. However, since this version is not supported anymore, we choose not to include it in our main results. Figure 7 shows the throughput performance of the p2z and p2r benchmarks on a two-socket system equipped with Intel Xeon Gold 6248 CPUs. All implementations were compiled with gcc except stdpar, which was compiled with nvc++, and SYCL, which was compiled with dpcpp. Nonstandard options such as -ffast-math are not included in the compilation since, while they may help vectorizing loops including trigonometric functions, they do not guarantee numerical reproducibility.

Using the portability layers, we are able to achieve throughput equal to or better than 70% of the original, native performance on the CPU for most implementations. The Alpaka implementation of the p2z benchmark actually outperformed the TBB reference implementation. The reason for the better performance of p2z Alpaka is not currently known, and the results might be platform dependent. Note that the p2z Alpaka implementation uses the OpenMP execution backend, while the p2r Alpaka implementation relies on the TBB execution backend. The p2r SYCL implementation only achieves 27% of the reference implementation, but unlike many of the other implementations, SYCL is a language extension and depends heavily on compiler optimizations rather than a performance-tuned library. When optimizing these implementations, considerations included deciding on an optimal data layout that works for both CPUs and GPUs and ensuring that the loops are properly vectorized when run on the CPU.