GPU Languages: Back to Basics
Jumanazarov Mardonbek
Plan:
This chapter focuses on lower-level languages for GPU processors. We call these native languages because they directly reflect the functionality of the target GPU hardware. We'll cover two such languages, CUDA and OpenCL, which have found widespread use. We'll also cover HIP, a new variant for AMD GPU processors. Unlike pragma-based implementations, these GPU languages are less dependent on the compiler. You should use these languages for more precise control over your program's performance. How do these languages differ from the languages presented in Chapter 11? We believe the difference is that those languages grew out of the characteristics of GPU and CPU hardware, whereas OpenACC and OpenMP began with high-level abstractions and relied on the compiler to map them to different hardware.
The set of native GPU languages, CUDA, OpenCL, and HIP, require the creation of separate source code for the GPU compute kernel. Separate source code often resembles CPU source code. The difficulty of maintaining two different source codes is a significant obstacle. If a native GPU language only supports one hardware type, then it may be necessary to support even more source code variants if you want to run on GPUs from multiple vendors. Some applications implement their algorithms in multiple GPU and CPU languages. Thus, the urgent need for more portable GPU programming languages becomes clear.
Fortunately, portability is increasingly being emphasized in some new GPU languages. OpenCL was the first open-standard language to run on a variety of GPU and even CPU hardware. After its initial popularity, OpenCL did not achieve the widespread adoption expected. Another language, HIP, was developed by AMD as a more portable version of CUDA that generates code for AMD GPUs. As part of AMD's portability initiative, support for GPUs from other manufacturers has been included.
The distinction between these native languages and higher-level languages is blurring as new languages emerge. SYCL, originally a C++ layer on top of OpenCL, is typical of these new, more portable languages. Along with Kokkos and RAJA, SYCL supports a single source code for both CPU and GPU hardware. We will discuss these languages at the end of the chapter. Figure 12.1 shows the current state of GPU language interoperability, which we will discuss in this chapter.
Figure 12.1: Mapping GPU language interoperability reveals an increasingly complex situation. Four GPU languages are shown at the top, and various hardware devices are shown at the bottom. Arrows indicate the code generation paths from the languages to the hardware. Dashed lines indicate hardware that is still under development.
The emphasis on language interoperability is gaining momentum as a wider variety of GPUs are deployed in major HPC installations. The Department of Energy's flagship HPC systems, Sierra and Summit, are powered by NVIDIA GPUs. In 2021, the Department of Energy's HPC roster will be expanded to include the Aurora system at Argonne, powered by Intel GPUs, and the Frontier system at Oak Ridge, powered by AMD GPUs. With the introduction of Aurora, SYCL has emerged from near obscurity to become a major player with numerous implementations. SYCL was originally developed to provide a more native C++ layer on top of OpenCL. Its sudden emergence stems from its adoption by Intel as part of the oneAPI programming model for Intel GPUs in Aurora. Because of its newfound importance, we discuss SYCL in Section 12.4. A similar growth in interest is being seen in other languages and libraries that support portability across the GPU landscape.
We conclude this chapter with a brief overview of two performance portability frameworks, Kokkos and RAJA, which were created to facilitate performance across a wide range of hardware, from CPUs to GPUs. They operate at a slightly higher level of abstraction but promise a single source code base that will run everywhere. Their development is the result of a major Department of Energy effort to support the portability of large scientific applications to newer hardware. The goal of RAJA and Kokkos is a one-time rewrite to create a single source code base that is portable and maintainable during periods of major hardware design changes.
1. Functionality of the native GPU programming language
A GPU programming language must have several basic features. It's useful to understand these features so that you can recognize them in each GPU language. Below, we briefly describe the necessary features of a GPU language.
1. Functionality of the native GPU programming language
When developing a language, one must also decide whether to place the host and design source code in the same file or in separate files. In either case, the compiler must distinguish between the host and design source code and must provide a way to generate the instruction set for different hardware. The compiler must even decide when to generate the instruction set. For example, OpenCL waits until the device is selected and then generates the instruction set using a just-in-time (JIT) compiler.
The mechanism for calling the device's computational cores from the host - is great, we now have the device source code—but we also need a way to call this code from the host. The syntax for this operation varies the most across languages. However, this mechanism is only slightly more complex than calling a standard routine.
1. Functionality of the native GPU programming language
This list isn't all that intimidating. Native GPU languages are, for the most part, not that different from current CPU source code. Furthermore, recognizing these commonalities in the functionality of a native GPU language helps you feel comfortable when switching from one language to another.
2. CUDA Languages and HIP GPU: Low-Level Performance Option
We'll begin by examining two low-level GPU languages, CUDA and HIP. These two GPU programming languages are the most widely used.
Compute Unified Device Architecture (CUDA) is NVIDIA's native language, running exclusively on their GPU processors. First released in 2008, it is currently the dominant native programming language for GPU processors. Over a decade of development, CUDA has acquired a rich set of features and performance-enhancing enhancements. CUDA closely mirrors NVIDIA's GPU architecture. It does not claim to be a universal accelerator language. However, the concepts of most accelerators are similar enough that the CUDA language design is applicable.
AMD (formerly ATI) GPU processors had a number of short-lived programming languages. They finally settled on a CUDA-like language that can be generated by "HIP-ifying" CUDA code using their HIP compiler. This functionality is part of the ROCm toolkit, which provides broad portability across GPU languages, including the OpenCL language for GPU (and CPU) processors described in Section 12.3.
2. CUDA Languages and HIP GPU: Low-Level Performance Option
Writing and Building Your First CUDA Application
We'll begin by answering the question of how to build and compile a simple CUDA application running on a GPU. We'll use the streaming triad example we've used throughout the book, which implements a loop for this calculation: C = A + scalar * B. The CUDA compiler splits the regular C++ code to feed the reference C++ compiler. It then compiles the remaining CUDA code. The source code from these two paths is linked together into a single executable.
To follow the example above, you may first need to install the CUDA software[1]. Each CUDA release works with a limited range of compiler versions. Starting with CUDA v10.2, GCC compilers up to version 8 are supported. If you work with multiple parallel languages and packages, this constant struggle with compiler versions is perhaps one of the most frustrating aspects of CUDA. On the plus side, you can use most of your regular toolchain and build systems with only version restrictions and a few special additions.
[1] For more information, see the CUDA Installation Guide (https://docs.nvidia.com/cuda/cuda-installation-guide-linux/).
2. CUDA Languages and HIP GPU: Low-Level Performance Option
We'll show three different approaches, starting with a simple makefile, and then two different ways to use CMake. We recommend you refer to the examples in this chapter, located at https://github.com/EssentialsofParallelComputing/Chapter12.
You can access this simple CUDA makefile by copying or linking it to Makefile, the default filename for make. The following listing shows the makefile itself.
1. To link to the file, type ln -s Makefile. simple Makefile.
2. Build the application using make.
3. Run the application using ./StreamTriad.
2. CUDA Languages and HIP GPU: Low-Level Performance Option
2. CUDA Languages and HIP GPU: Low-Level Performance Option
The key addition is the pattern rule in lines 9-10, which converts a file with the .cu suffix into an object file. We use the NVIDIA NVCC compiler for this operation. Next, we need to add the CUDA runtime library, CUDART, to the link string. You can use lines 4 and 5 to specify a specific NVIDIA GPU architecture and a special path to the CUDA libraries.
DEFINITION: A pattern rule is a description intended for the make utility, specifying a general rule for converting any file with one suffix pattern to a file with a different suffix pattern.
CUDA has extensive support in the CMake build system. Next, we'll look at both the old-style support and the new, modern CMake approach, which was introduced more recently. Listing 12.2 shows the old-style method. Its advantage is greater portability to systems with older versions of CMake and automatic detection of the NVIDIA GPU architecture. This latest hardware detection functionality is so convenient that it is now recommended to use the legacy CMake. To use this build system, link CMakeLists_old.txt to CMakeLists.txt:
ln -s CMakeLists_old.txt CMakeLists.txt mkdir build && cd build
cmake ..
make
2. CUDA Languages and HIP GPU: Low-Level Performance Option
2. CUDA Languages and HIP GPU: Low-Level Performance Option
2. CUDA Languages and HIP GPU: Low-Level Performance Option
Much of the CMake build system remains standard. The shared compilation attribute on line 11 is intended for a more robust build system focused on general development. It can then be disabled at a later stage to save a few registers in the CUDA compute cores, yielding minor optimizations in the generated code. The default CUDA values are intended for performance, not for a more general and robust build. The automatic detection of the NVIDIA GPU architecture on line 14 is a significant convenience, eliminating the need to manually modify the makefile.
With version 3.0, the CMake build system undergoes a significant overhaul, which, as they say, "modernizes" CMake. Key attributes of this style are a more integrated system and targeted use of attributes. This is nowhere more evident than in its CUDA support. Let's take a look at Listing 12.3 to see how to use it. To use this build system for CUDA support in the modern, new CMake style, you need to link CMakeLists_new.txt to CMakeLists.txt:
ln -s CMakeLists_new.txt CMakeLists.txt mkdir build && cd build
cmake ..
make
2. CUDA Languages and HIP GPU: Low-Level Performance Option
2. CUDA Languages and HIP GPU: Low-Level Performance Option
2. CUDA Languages and HIP GPU: Low-Level Performance Option
The first thing to note about using this modern CMake approach is how much simpler it is compared to the old style. The key is the inclusion of CUDA as the language in line 4. From here on, a little extra work is required.
As shown in lines 9-10, we can set flags to compile for a specific GPU architecture. However, in the modern CMake build system, we don't yet have an automatic way to detect the architecture. Without an architectural flag, the compiler generates code and optimizes for the sm_30 GPU device. The generated sm_30 source code runs on any device from the Kepler K40 or later, but it won't be optimized for newer architectures. Furthermore, it's possible to specify multiple architectures in a single compiler. Compilations will be slower, and the generated executable will be larger.
We can also set a separate compilation attribute for CUDA, but with a different syntax that applies it to a specific target. The optimization flag on line 10, -O3, is sent to the host compiler only for regular C++ source code. For CUDA source code, the default optimization level is -O3, which rarely needs to be modified.
Overall, building a CUDA program is simple and getting simpler. However, it's expected that build changes will continue. Clang is adding native support for compiling CUDA source code, giving you another option besides the NVIDIA compiler. Now let's move on to the source code. We'll start with the GPU compute kernel in the listing below.
2. CUDA Languages and HIP GPU: Low-Level Performance Option
2. CUDA Languages and HIP GPU: Low-Level Performance Option
As is typical with GPU compute kernels, we strip the for loop from the compute unit. This leaves the loop body on line 14. We need to add a conditional block on line 12 to prevent out-of-bounds access. Without this protection, compute kernels could crash silently. Then, on line 9, we obtain a global index from the block and thread variables defined by the CUDA runtime. Adding the _global_ attribute to the routine informs the compiler that this GPU compute kernel will be invoked from the host. Meanwhile, on the host side, we need to configure memory and execute the kernel invocation. This process is shown in the listing below.
2. CUDA Languages and HIP GPU: Low-Level Performance Option
2. CUDA Languages and HIP GPU: Low-Level Performance Option
2. CUDA Languages and HIP GPU: Low-Level Performance Option
First, we allocate memory on the host and initialize it in lines 31–39. We also need adequate memory space on the GPU to store the arrays while the GPU is accessing them. For this, we use the cudaMalloc routine in lines 43–45. Now we come to a few interesting lines (47–49) that are specific to the GPU. The block size is the workgroup size on the GPU. It is determined by the tile size, block size, or workgroup size, depending on the GPU programming language used (see Table 10.1). The next line, which calculates the grid size, is specific to the GPU code. We won't always have an array size that is an even integer multiple of the block size. This means we need an integer that is equal to or greater than the fractional number of blocks. Let's walk through the example step by step to understand what's going on.
Now all blocks except the last one have 512 values. The final block will be 512 in size but will only contain 488 data elements. The out-of-bounds check on line 12 of Listing 12.4 prevents problems with this partially filled block. The last few lines of Listing 12.5 free the device and host pointers. Remember to use cudaFree for device pointers and the C library function free for host pointers.
2. CUDA Languages and HIP GPU: Low-Level Performance Option
Example: Calculating the Block Size for a GPU
In line 3 of the listing below, we calculate the fractional number of blocks. In this example, with an array size of 1000, it equals 1.95 blocks. Instead of truncating it to 1, which would happen by default when using integer arithmetic, we need to round up to 2. If we simply calculated the array size divided by the block size, we would get integer truncation. Therefore, we must cast each of them to a floating-point value to obtain floating-point division. In reality, we only need to cast one of the values, and the C/C++ standard requires the compiler to provide the other elements. However, our programming conventions must explicitly mandate type conversion, otherwise it is a programming error. Compilers often don't signal these cases, but they can mask unintended situations.
The ceil function, used in lines 4 and 5 of the listing, rounds up to the next integer value equal to or greater than the floating-point number. We can achieve the same result using integer arithmetic by adding the block size minus one and then performing integer division with truncation, as done in line 6. We chose to use this version because the integer form doesn't require any floating-point operations and should be faster.
2. CUDA Languages and HIP GPU: Low-Level Performance Option
1. int stream_array_size = 1000
2. int blocksize = 512
3. float frac_blocks = (float)stream_array_size/(float)blocksize;
>>>frac_blocks = 1.95
4. int nblocks = ceil(frac_blocks);
>>> nblocks = 2
либо
5, int nblocks = ceil((float)stream_array_size/(float)blocksize);
либо
6. int nblocks = (stream_array_size + blocksize - 1)/blocksize;
2. CUDA Languages and HIP GPU: Low-Level Performance Option
All that's left to do is copy the memory to the GPU, call the GPU compute kernel, and copy the memory back. We do this in a timing loop (in Listing 12.6), which can be executed multiple times to obtain a more optimal measurement result. Sometimes the first GPU call will be slower due to initialization overhead. We can amortize this overhead by running multiple iterations. If this isn't enough, you can also discard the timing from the first iteration.
2. CUDA Languages and HIP GPU: Low-Level Performance Option
2. CUDA Languages and HIP GPU: Low-Level Performance Option
The timing loop template consists of the following steps:
1. Copy data to the GPU (lines 53-54).
2. Call the GPU compute kernel for array manipulation (line 59).
3. Copy the data back (line 64).
We add several synchronization and timer calls to obtain an accurate GPU compute kernel measurement. At the end of the loop, we then verify the correctness of the result. When deploying this source code to production, we can remove the timing, synchronization, and error checking. The GPU compute kernel call is easily identified by the three chevrons, or angle brackets. If you ignore the chevrons and the variables they contain, the line of code has the typical syntax of a C subroutine call:
StreamTriad(stream_array_size, scalar, a_d, b_d, c_d);
The values in parentheses are arguments to be passed to the GPU compute kernel. For example:
<<<grid size, block size>>>
2. CUDA Languages and HIP GPU: Low-Level Performance Option
So, what arguments are contained in the chevrons? These are arguments to the CUDA compiler about how to break the problem into blocks for the GPU. Furthermore, in lines 48–49 of Listing 12.2, we set the block size and calculated the number of blocks, or grid size, to fit all the data into the array. Here, the arguments are one-dimensional. We can also have two- or three-dimensional arrays by declaring and setting these arguments as values for an N×N matrix.
dim3 blocksize(16,16); dim3 blocksize(8,8,8);
dim3 gridsize( (N + blocksize.x - 1)/blocksize.x,
(N + blocksize.y - 1)/blocksize.y );
We can speed up memory transfers by eliminating data copying. This is possible thanks to a deeper understanding of how the operating system functions. Memory transferred over the network must be in a fixed location that cannot be moved during the operation. Regular memory allocations are placed in pageable memory, or memory that can be moved on demand. When transferring memory, the data must first be moved to pinned memory, or memory that cannot be moved. We first saw the use of pinned memory in Section 9.4.2 when comparing memory transfers over the PCI bus. We can eliminate memory copying by allocating our arrays in pinned memory rather than pageable memory. Figure 9.8 shows the performance difference we could achieve. Now the question arises: how can this be accomplished?
2. CUDA Languages and HIP GPU: Low-Level Performance Option
CUDA provides us with a function call, cudaHostMalloc, that does this for us. This function is a direct replacement for the regular system malloc routines, with a slight change in the arguments, where a pointer is returned as an argument, as shown below:
double *x_host = (double *)malloc(stream_array_size*sizeof(double)); cudaMallocHost((void**)&x_host, stream_array_size*sizeof(double));
Is there a downside to using pinned memory? The problem is that if you use a lot of pinned memory, there won't be room to bring in another application. Taking one application out of memory and bringing in another is a huge convenience for users. This process is called memory paging.
DEFINITION: Memory paging in multi-user, multi-application operating systems is the process of temporarily moving memory pages to disk to allow another process to run.
2. CUDA Languages and HIP GPU: Low-Level Performance Option
Memory paging is a major advancement in operating systems that allows you to appear to have more memory than you actually have. For example, it allows you to temporarily run Excel while working in Word without having to close the original application. This is accomplished by writing your data to disk and then reading it back when you return to Word. However, this operation is expensive, so in high-performance computing, we avoid memory paging due to the significant performance penalty it entails. Some heterogeneous computing systems, including both CPUs and GPUs, implement unified memory.
DEFINITION: Unified memory is memory that appears as a single address space to both the CPU and GPU.
By now, you've seen that managing separate memory spaces on the CPU and GPU significantly complicates writing GPU code. With unified memory, the GPU runtime handles this for you. You may still have two separate arrays, but the data is moved around automatically. On integrated GPUs, it's possible that memory may not need to be moved at all. However, it's recommended to write your programs with explicit memory copying so that they can be ported to systems without unified memory. Memory copying is skipped if the architecture doesn't require it.
2. CUDA Languages and HIP GPU: Low-Level Performance Option
CUDA's Reducer Compute Kernel: Things Get More Complex
When we require cooperation between GPU threads, things get more complicated in low-level, native GPU languages. We'll look at a simple summation example to understand how to handle this. The example requires two separate CUDA compute kernels and is shown in Listings 12.7–12.10. The following listing shows the first pass, in which we sum the values within a thread block and store the result back into the reducer scratchpad array, redscratch.
2. CUDA Languages and HIP GPU: Low-Level Performance Option
2. CUDA Languages and HIP GPU: Low-Level Performance Option
We begin the first pass by having all threads store their data in a scratchpad array in CUDA shared memory (lines 35–38). All threads in a block can access this shared memory. Shared memory can be accessed in one or two CPU cycles, compared to the hundreds required for main GPU memory. Shared memory can be treated as either a programmable cache or a scratchpad. To ensure that all threads complete their storage, we use a synchronization call in line 40.
Since the reduction sum within the block will be used in both reduction passes, we place the source code in a device routine and call it in line 42. A device routine is a routine that should be called not from the host, but from another device routine. After executing the routine, the resulting sum is stored back into the smaller data array, which we read during the second phase. We also save the result in line 47 in case the second pass might be skipped. Since we can't access values in other thread blocks, we must complete the operation in a second pass on another compute core. In this first pass, we reduced the data length by the size of our block.
Let's move on to the general device source code we mentioned in the first pass. We'll need summation reduction for the CUDA thread block in both passes, so we'll write it as a general device routine. The source code shown in the following listing can easily be modified for other reduction operators and requires only minor changes to accommodate HIP and OpenCL.
2. CUDA Languages and HIP GPU: Low-Level Performance Option
2. CUDA Languages and HIP GPU: Low-Level Performance Option
The shared device routine, which will be called from both passes, is defined on line 3. It performs summation reduction in the threaded block. The _device_ attribute before the routine specifies that it will be called from a GPU compute kernel. The basic concept of the routine is based on a pairwise reduction tree in O(log n) operations, as shown in Figure 12.2. The basic reduction tree from the figure is represented in the source code on lines 15–18. We make a few minor modifications when the working set exceeds the warp size on lines 8–13 and for the final pass level on line 19 to avoid unnecessary synchronization.
The same pairwise reduction concept is used for the full-threaded block, which can be as large as 1024 on most GPU devices, although 128 to 256 are more common. But what if your array size exceeds 1024? We add a second pass that uses only one flow block, as shown in the following listing.
2. CUDA Languages and HIP GPU: Low-Level Performance Option
Fig. 12.2 Pairwise reduction tree for a warp that sums values over log n steps
2. CUDA Languages and HIP GPU: Low-Level Performance Option
2. CUDA Languages and HIP GPU: Low-Level Performance Option
To avoid using more than two cores for larger arrays, we use a single thread block and loop in lines 67–69 to read and sum any additional data into a shared scratchpad array. We use a single-threaded block because synchronization can be performed within it, avoiding the need for another invocation of a compute core. If we use thread block sizes of 128 and have an array of a million elements, the loop will sum about 60 values at each location in shared memory (1,000,000 / 1282). The array size is reduced by 128 on the first pass, and then we sum into a scratchpad of size 128, yielding a division by 128 squared. If we used larger block sizes, such as 1024, we could reduce the loop from 60 iterations to a single read. Now we simply call the same shared thread block reduce we used before. The result will be the first value in the scratchpad array. The final part of the reduction is configuring and invoking these two compute cores from the host. We'll see how this is done in the listing below.
2. CUDA Languages and HIP GPU: Low-Level Performance Option
2. CUDA Languages and HIP GPU: Low-Level Performance Option
2. CUDA Languages and HIP GPU: Low-Level Performance Option
The host source code first calculates the sizes for the compute kernel calls in lines 100–103. Then we must allocate memory for the device arrays. For this operation, we need a scratchpad array in which we can store the sums for each block from the first compute kernel. We allocate it in line 108 with a size of gridsize because this number equals the number of blocks we have. We also need a shared-memory scratchpad array equal to the block size. We calculate this size in line 101 and pass it to the kernel in lines 112 and 115 as the third parameter to the chevron operator. The third parameter is optional; this is the first time we've seen its use. Look at Listing 12.9 (line 56) and Listing 12.7 (line 29) to see where the corresponding scratchpad source code is processed on the GPU device.
Trying to follow all the tangled loops can be difficult. Therefore, we created a version of the source code that executes the same cycles on the CPU and prints its values as it runs. This version is located in the CUDA/SumReductionRevealed directory at https://github.com/EssentialsofParallelComputing/Chapter12.
We don't have space to show the entire source code here, but you might find it useful to perform reconnaissance analysis of the values and print them as it runs. We'll show an edited version of the printout in the example below.
2. CUDA Languages and HIP GPU: Low-Level Performance Option
Example: CUDA/SumReductionRevealed
Calling first pass with gridsize 2 blocksize 128 blocksizebytes 1024
SYNCTHREADS after all values are in shared memory block Data count is 200
====== ITREE_LEVEL 1 offset 64 ntX is 128 MIN_REDUCE_SYNC_SIZE 32 ====
Data count is reduced to 128
Sync threads when larger than warp
====== ITREE_LEVEL 2 offset 32 ntX is 128 MIN_REDUCE_SYNC_SIZE 32 ====
Sync threads when smaller than warp Data count is reduced to 64
====== ITREE_LEVEL 3 offset 16 ntX is 128 MIN_REDUCE_SYNC_SIZE 32 ====
Sync threads when smaller than warp Data count is reduced to 32
====== ITREE_LEVEL 4 offset 8 ntX is 128 MIN_REDUCE_SYNC_SIZE 32 ====
Sync threads when smaller than warp Data count is reduced to 16
====== ITREE_LEVEL 5 offset 4 ntX is 128 MIN_REDUCE_SYNC_SIZE 32 ====
Sync threads when smaller than warp Data count is reduced to 8
====== ITREE_LEVEL 6 offset 2 ntX is 128 MIN_REDUCE_SYNC_SIZE 32 ====
Sync threads when smaller than warp Data count is reduced to 4
2. CUDA Languages and HIP GPU: Low-Level Performance Option
====== ITREE_LEVEL 7 offset 1 ntX is 128 MIN_REDUCE_SYNC_SIZE 32 ====
Data count is reduced to 2
Finished reduction sum within thread block End of first pass
Synchronization in second pass after loading data Data count is reduced to 2
====== ITREE_LEVEL 8 offset 1 ntX is 128 MIN_REDUCE_SYNC_SIZE 32 ====
Data count is reduced to 1
Finished reduction sum within thread block End of first pass
Synchronization in second pass after loading data Data count is reduced to 2
====== ITREE_LEVEL 8 offset 1 ntX is 128 MIN_REDUCE_SYNC_SIZE 32 ====
Data count is reduced to 1
Finished reduction sum within thread block
Synchronization in second pass after reduction sum Result -- total sum 19900
2. CUDA Languages and HIP GPU: Low-Level Performance Option
This example applies to an array of 200 integers, each element of which is initialized with an index value. We encourage you to consult the source code and Figure 12.1 to understand what's going on. The start and end of the first and second passes are printed. We see that the amount of data is reduced by half, until only two remain at the end of the first pass. The second pass quickly reduces them to a single value containing the sum.
We showed this reduction of thread blocks as a general introduction to computing kernels that require thread cooperation. It's easy to see how complex this is, especially compared to the single line required for an internal Fortran call. Along the way, we also achieved significant speedup on the CPU and stored the data for this operation on the GPU. This reduction algorithm can be further optimized, but you might also consider using some library services, such as CUDA UnBound (CUB), Thrust, or other GPU libraries.
2. CUDA Languages and HIP GPU: Low-Level Performance Option
HIPing the CUDA Source Code
The CUDA source code only runs on NVIDIA GPUs. However, AMD has implemented a similar GPU language called Heterogeneous Interface for Portability (HIP). It is part of AMD's ROCm (Radeon Open Compute Platform) toolchain. If you program in the HIP language, you can invoke the hipcc compiler, which uses NVCC on NVIDIA platforms and HCC on AMD GPUs.
To try these examples, you may need to install the ROCm software, information, and toolset. The installation process changes frequently, so please check the latest instructions. Some instructions are also included with the examples.
2. CUDA Languages and HIP GPU: Low-Level Performance Option
Example: A Simple makefile for HIPifying CUDA Source Code
There are two versions of the Makefile. One uses hipify-perl, and the other uses hipify-clang. hipify-perl is a simple Perl script. For a more accurate, syntax-aware translation, you can try hipify-clang. In any case, for more complex programs, you may need to make the final modifications manually. We'll be using the Perl version, so let's start by linking the Makefile.perl script shown in the following listing to the Makefile:
ln -s Makefile.perl Makefile
make
2. CUDA Languages and HIP GPU: Low-Level Performance Option
2. CUDA Languages and HIP GPU: Low-Level Performance Option
The only real addition to the standard makefile is replacing the compiler with hipcc and adding a template rule to convert the CUDA source code to HIP source code. We could perform the code conversion by simply calling the hipify-perl script manually and then applying the HIP version to both CUDA and AMD GPUs.
CMake also has good HIP support. This support has been available since CMake 2.8.3. A typical CMakeLists file for HIP is shown in the listing below.
2. CUDA Languages and HIP GPU: Low-Level Performance Option
2. CUDA Languages and HIP GPU: Low-Level Performance Option
In this listing, we first try specifying different paths for the possible HIP installation location, and then call find_package for HIP on line 15. Then, on line 22, we set the C++ compiler to hipcc. The HIP_ADD_EXECUTABLE command adds a build of our executable, and we finish the listing with settings for the HIP headers and libraries (lines 28–31). Now let's turn our attention to the HIP source code in Listing 12.12. We highlight the changes from the CUDA version of the source code, shown in Listings 12.5–12.6.
2. CUDA Languages and HIP GPU: Low-Level Performance Option
2. CUDA Languages and HIP GPU: Low-Level Performance Option
To convert from CUDA source code to HIP source code, we replace all instances of cuda in the source code with hip. The only significant change is in the execution call to the compute kernel, where HIP uses a more traditional syntax than the triple chevron used in CUDA. Surprisingly, the biggest change is the use of bilingual variable naming terminology.
3. OpenCL for Portable GPU Open Source Language
In response to the urgent need for portable GPU source code, a new GPU programming language called OpenCL emerged in 2008. OpenCL is an open standard GPU language that runs on NVIDIA and AMD/ATI graphics cards, as well as many other hardware devices. Apple led the development of the OpenCL standard, with contributions from many other organizations. One of the nice features of OpenCL is that virtually any C or even C++ compiler can be used for host source code. For GPU device source code, OpenCL was initially based on a subset of C99. Recently, OpenCL versions 2.1 and 2.2 added support for C++14, but implementations are still unavailable.
The OpenCL release started with great initial excitement. Finally, a way to write portable GPU code. For example, GIMP announced that it would support OpenCL as a GPU acceleration solution that would be available on many hardware platforms. The reality, however, has proven less convincing. Many believe that OpenCL is too low-level and verbose for widespread adoption. It's even possible that its ultimate role is to serve as a low-level portability layer for higher-level languages. But its value as a language for writing portable code across a wide range of hardware devices has been demonstrated by its availability in the embedded community for field-programmable gate arrays (FPGAs). One reason OpenCL is considered verbose is that device selection is more complex (and more powerful). You have to discover and select the device you'll be running on, which can lead to hundreds of lines of code right from the start.
3. OpenCL for Portable GPU Open Source Language
Almost everyone who uses OpenCL writes a library to handle low-level tasks. We're no exception. Our library is called EZCL. Almost every OpenCL call is wrapped with at least a lightweight layer to handle error conditions. Device detection, code compilation, and error handling consume many lines of code.
In our examples, we'll use a reduced version of our EZCL library, called EZCL_Lite, so we can see the actual OpenCL calls. The EZCL_Lite routines are used to select a device and configure it for the application, then compile the device code and handle errors. The source code for these operations is too long to show here, so please refer to the examples in the OpenCL directory at https://github.com/EssentialsofParallelComputing/Chapter12. This directory also contains the full EZCL library.
The EZCL routines provide detailed information about call errors and the line of source code where the error occurs. Before trying the OpenCL source code, make sure you have a properly configured environment and devices. You can use the clinfo command for this.
3. OpenCL for Portable GPU Open Source Language
Example: Obtaining OpenCL Installation Information
Run the OpenCL information command:
clinfo
If you receive the following output, OpenCL is not configured or you do not have the appropriate OpenCL device:
Number of platforms 0
If you do not have the clinfo command, try installing it using the appropriate command for your system. For Ubuntu, this is:
sudo apt install clinfo
The examples in this chapter provide some quick tips for installing OpenCL, but check for the latest information specific to your system. OpenCL has an extension that provides a detailed model of how each device should configure its driver in its Installable Client Driver (ICD) specification. This allows an application to have multiple OpenCL platforms and drivers.
3. OpenCL for Portable GPU Open Source Language
Writing and Building Your First OpenCL Application
Changing the standard makefile to embed OpenCL isn't too difficult. Typical changes are shown in Listing 12.13. To use the simple makefile for OpenCL, type:
ln -s Makefile.simple Makefile
Then build the application with make and run it with ./StreamTriad.
3. OpenCL for Portable GPU Open Source Language
3. OpenCL for Portable GPU Open Source Language
The makefile contains a method for setting the DEVICE_DETECT_DEBUG flag to print detailed information about available GPU devices. This flag enables more detailed information in the ezcl_lite.c source code. This is useful for troubleshooting device detection issues or for detecting the wrong device. Line 6 also adds a template rule that embeds the OpenCL source code into the program for use at runtime. On line 9, this Perl script converts the source code into a comment and as a dependency. It will be included in the StreamTriad.c file using the include statement.
We developed the embed_source.pl utility to bundle OpenCL source code directly with the executable. (The source code for this utility is provided in the examples in this chapter.) A common way to implement OpenCL code is to have a separate source file that must be localized at runtime, which is then compiled after the device is known. Using a separate file creates problems such as being unable to find it or getting the wrong version of the file. We strongly recommend embedding the source code in the executable to avoid these types of issues. As shown in the listing below, it is also possible to use CMake support for OpenCL in our build system.
3. OpenCL for Portable GPU Open Source Language
3. OpenCL for Portable GPU Open Source Language
3. OpenCL for Portable GPU Open Source Language
OpenCL support was added to CMake in version 3.1. We added this version requirement at the top of the CMakelists.txt file on line 1. A few other special points are worth noting. In this example, we used the -DDEVICE_DETECT_DEBUG=1 option in the CMake command to enable device detection verbosity. We also added a way to enable and disable OpenCL double-precision support. We used it in the EZCL_Lite source code to set the JIT (just-in-time) compilation flag for the OpenCL device source code. Finally, we added an application-specific command on lines 19–22 to embed the OpenCL device source code into the executable. The OpenCL compute kernel source code is located in a separate file called StreamTriad_kernel.cl, as shown in the following listing.
3. OpenCL for Portable GPU Open Source Language
3. OpenCL for Portable GPU Open Source Language
Compare this compute kernel source code with the CUDA compute kernel source code in Listing 12.4. The OpenCL source code is nearly identical, except that __kernel replaces __global__ in the routine declaration, the __global attribute is added to pointer arguments, and there is a different way to obtain the thread index. Furthermore, the CUDA compute kernel source code is in the same .cu file as the host source code, whereas the OpenCL source code is in a separate .cl file. We could separate the CUDA source code into its own .cu file and place the host source code in a standard C++ source file. This would be similar to the structure we use for our OpenCL application.
NOTE: Many of the differences between the CUDA and OpenCL compute kernel source codes are superficial.
So, how does the host-side OpenCL source code differ from the CUDA version? Let's look at the OpenCL version in Listing 12.16 and compare it to the source code in Listing 12.5. There are two versions of the OpenCL stream triad: StreamTriad_simple.c without error checking and StreamTriad.c with error checking. Error checking adds many lines of code that initially simply make it difficult to understand what's going on.
3. OpenCL for Portable GPU Open Source Language
3. OpenCL for Portable GPU Open Source Language
3. OpenCL for Portable GPU Open Source Language
At the beginning of the program, we encounter a few real differences in lines 34–37, where we must find our GPU device and compile the source code for our device. This is done for us behind the scenes by CUDA. Two lines of OpenCL code call our EZCL_Lite routines to detect the device and create a program object. We made these calls because the source code required for these functions is too large to show here. These routines represent hundreds of lines, with most of them devoted to error checking.
NOTE: The source code with the chapter's examples is available in the OpenCL/StreamTriad directory at https://github.com/EssentialsofParallelComputing/Chapter12. Some of the error checking source code was omitted from the short version of StreamTriad_simple.c, but it is contained in the long version of the code in StreamTriad.c. The rest of the setup and teardown code follows the same pattern we encountered in the CUDA code, with a little extra cleanup, again related to manipulating the device and program source code. Now, how does the section of source code that calls the OpenCL compute kernel in the timing loop from Listing 12.16 relate to the CUDA source code from Listing 12.6?
3. OpenCL for Portable GPU Open Source Language
3. OpenCL for Portable GPU Open Source Language
3. OpenCL for Portable GPU Open Source Language
What's going on in lines 57–61? OpenCL requires a separate call for each compute kernel argument. If we inspect the numeric code returned from each kernel, we'll get even more lines. This will be much more verbose than the single line 53 in Listing 12.6 in the CUDA version. But there's a direct correspondence between the two versions. OpenCL simply describes the argument-passing operations in more detail. With the exception of device detection and program compilation, both programs are similar in their operations. The biggest difference is the syntax used in the two languages.
In Listing 12.18, we show an approximate sequence of calls for device detection and creation. These routines are lengthy due to error checking and handling required for special cases. Good error handling is essential in these two functions. We need the compiler to report an error in the source code or that it received the wrong GPU device.
3. OpenCL for Portable GPU Open Source Language
3. OpenCL for Portable GPU Open Source Language
We conclude this presentation of OpenCL with a nod to the numerous language interfaces that have been created for it. There are versions in C++, Python, Perl, and Java. Each of these languages has a higher-level interface that hides some of the details found in the C version of OpenCL. We strongly recommend using our EZCL library or one of the many other OpenCL middleware libraries.
An unofficial C++ version has been available since OpenCL v1.2. Its implementation is merely a thin layer on top of the C version of OpenCL. Although it has not received approval from the standards committee, it is fully usable by developers. It is available at https://github.com/Khronos-Group/OpenCL-CLHPP. Official approval of C++ in OpenCL occurred only recently, but we are still awaiting implementations.
Reductions in OpenCL
Summation reduction in OpenCL is similar to reduction in CUDA. Instead of a step-by-step analysis of the source code, we'll simply look at the differences in the source code of the computing kernel. First, Figure 12.3 shows the differences in the sum_within_block procedure, which is common to both computing kernels, in a side-by-side comparison.
3. OpenCL for Portable GPU Open Source Language
Fig. 12.3 Comparison of reduction kernels in OpenCL and CUDA: sum_within_block
3. OpenCL for Portable GPU Open Source Language
The difference in this device compute kernel calling another compute kernel begins with the attributes in the declaration. CUDA requires the _device_ attribute in the declaration, while OpenCL does not. Arguments passed in a scratchpad array require the _local attribute, which CUDA does not need. Another difference is the syntax for defining the local thread index and block (tile) size (Figure 12.3, lines 5 and 6). The synchronization calls are also different. At the top of the routine, the warp size is defined by a macro, which helps ensure portability between NVIDIA and AMD GPUs. CUDA defines it as a warp size variable. In OpenCL, it is passed along with the compiler definition. We also change the terminology from block to tile in the actual source code to ensure consistency with the terminology of each language.
The following procedure represents the first of two passes of the compute kernel, called stage 1 of 2, in Figure 12.4. This kernel is invoked from the host. The __global__ attribute for CUDA becomes __kernel for OpenCL. We must also add the __global attribute to the pointer arguments for OpenCL.
3. OpenCL for Portable GPU Open Source Language
Fig. 12.4 Comparison for the first of two compute kernel passes for the OpenCL and CUDA reduction kernels
3. OpenCL for Portable GPU Open Source Language
The next difference is important and worth noting. In CUDA, we declare a shared-memory scratchpad as an extern shared variable in the body of the compute kernel. On the host side, the size of this shared memory space is specified as a number of bytes in the optional third argument enclosed in triple chevron brackets. In OpenCL, this is done differently. It is passed as the last argument in the argument list with the _local attribute. On the host side, the memory is specified in the call that sets the value of the fourth argument of the compute kernel:
clSetKernelArg(reduce_sum_1of2, 4,
local_work_size*sizeof(cl_double), NULL);
The size is the third argument in the call. The remaining changes concern the syntax for setting thread parameters and calling synchronization. The final part of the comparison is the second pass of the compute reduce kernel in Figure 12.5.
We have already encountered all the change patterns in the second kernel. We still have differences in the computation kernel declaration and arguments. The local scratchpad array also has the same differences as the first-pass computation kernel. Thread parameters and synchronization also have the expected differences.
Looking back at the three comparisons in Figures 12.3–12.5, we clearly see something we didn't have to note. The bodies of the computation kernels are essentially the same.
3. OpenCL for Portable GPU Open Source Language
Fig. 12.5 Comparison of the second pass of summation reduction
3. OpenCL for Portable GPU Open Source Language
The only difference is the synchronization call syntax. The host-side source code for summation reduction in OpenCL is shown in the listing below.
3. OpenCL for Portable GPU Open Source Language
3. OpenCL for Portable GPU Open Source Language
3. OpenCL for Portable GPU Open Source Language
The first pass of the compute kernel creates a local scratchpad array on line 46. Intermediate results are saved back to the redscratch array created on line 38. If more than one block exists, a second pass is required. The redscratch array is passed back to complete the reduction. Note that the compute kernel parameters in arguments 5 and 6 are set to local_work_size or one workgroup. This is done to ensure synchronization across all remaining data, and another pass is not required.
4. SYCL: Experimental C++ implementation becomes mainstream
SYCL was born in 2014 as an experimental C++ implementation on top of OpenCL. The goal of SYCL's developers is to provide a more natural extension to the C++ language than the C-based OpenCL implementation. It is being developed as a cross-platform abstraction layer that leverages the portability and efficiency of OpenCL. The focus of the experimental language suddenly changed when Intel selected it as one of its primary language paths for the Department of Energy's announced Aurora HPC system. The Aurora system will utilize Intel's new discrete GPUs, which are currently under development. Intel has proposed several extensions to the SYCL standard, which they prototyped in their Data Parallel C++ Compiler (DPCPP) in their open source oneAPI programming framework.
You can explore SYCL in several ways. Some even eliminate the need for software installation or the right hardware. You can try the following cloud systems first.
4. SYCL: Experimental C++ implementation becomes mainstream
You can also download and install versions of SYCL from the following websites:
We will be working with Intel's DPCPP version of SYCL. Instructions for setting up a oneAPI installation in VirtualBox with the examples included in this chapter are provided in the README.virtualbox file at https://github.com/EssentialsofParallelComputing/Chapter12. You should be able to run VirtualBox on almost any operating system. Let's start with a simple makefile for the DPCPP compiler, as shown in the following listing.
4. SYCL: Experimental C++ implementation becomes mainstream
Setting the C++ compiler to the Intel dpcpp compiler takes care of paths, libraries, and include files. The only other requirement is setting some flags for the C++ compiler. The following listing shows the SYCL source code for our example.
4. SYCL: Experimental C++ implementation becomes mainstream
4. SYCL: Experimental C++ implementation becomes mainstream
4. SYCL: Experimental C++ implementation becomes mainstream
4. SYCL: Experimental C++ implementation becomes mainstream
The first Sycl function selects a device and creates a queue to run on it. We request a CPU, although this code will also work on GPUs with unified memory.
Sycl::queue Queue(sycl::cpu_selector{});
We select the CPU for maximum portability, so that the code will run on most systems. For this code to run on GPUs without unified memory, we would need to explicitly copy data from one memory space to another. The default selector preferentially searches for the GPU but falls back to the CPU. If we want to select only the GPU or only the CPU, we can also specify other selectors, such as:
Sycl::queue Queue(sycl::default_selector{}); // takes the device used by default
//
Sycl::queue Queue(sycl::gpu_selector{}); // finds the GPU device Sycl::queue Queue(sycl::cpu_selector{}); // finds the CPU device Sycl::queue Queue(sycl::host_selector{}); // executes on the host (CPU)
4. SYCL: Experimental C++ implementation becomes mainstream
The last option means it will run on the host as if there was no SYCL or OpenCL source code. Setting up the device and queue is much simpler than what we did in OpenCL. Now we need to set up the device buffers using the SYCL buffer:
Sycl::buffer<double,1> dev_a { a.data(), Sycl::range<1>(a.size()) };
The first argument to the buffer is the data type, and the second is the data dimension. Then we give it a variable name, dev_a. The first argument to the variable is the host data array used to initialize the device array, and the second is the set of indices to use. In this case, we specify a one-dimensional range from 0 to the size of the variable a. On line 29, we encounter the first lambda for creating the command group handler for the queue:
Queue.submit([&](Sycl::handler& CommandGroup)
4. SYCL: Experimental C++ implementation becomes mainstream
We introduced lambdas in Section 10.2.1. The lambda capture expression, [&], specifies the capture of external variables used by reference in a routine. For this lambda, the capture gets nsize, scalar, dev_a, dev_b, and dev_c for use in the lambda. We could specify this with just the capture-by-reference [&] setting, or with the following form, where we specify each variable to be captured. Good programming practice would prefer the latter, but the lists can get long.
Queue.submit([&nsize, &scalar, &dev_a, &dev_b, &dev_c]
(Sycl::handler& CommandGroup)
In the lambda body, we access the device arrays and rename them for use within the device procedure. This is equivalent to the argument list for the command group handler. We then create the first operational task for the command group, parallel_for. parallel_for is also defined using a lambda.
CommandGroup.parallel_for<class StreamTriad>(Sycl::range<1>{nsize},[=]
(Sycl::id<1> it)
4. SYCL: Experimental C++ implementation becomes mainstream
The lambda is named StreamTriad. We then tell it that we will be operating on a one-dimensional range, ranging from 0 to nsize. The capture expression, [=], captures the variables a, b, and c by value. Determining whether the capture is by reference or by value is difficult. However, if the source code is hosted on the GPU, the initial reference may go out of scope and become invalid. Finally, we create a one-dimensional index variable, it, to loop through the range.
5. Higher-level languages for performance portability
By now, you've seen that the differences between CPU and GPU compute cores aren't all that significant. So why not generate each of them using C++ polymorphism and templates? This is precisely what several libraries developed by Department of Energy research labs have done. These projects were started to address the problem of porting large amounts of source code to new hardware architectures. Kokkos was created by Sandia National Laboratories and has become widely deployed. Lawrence Livermore National Laboratory has a similar project called RAJA. Both of these projects have already succeeded in achieving their goal of providing multiplatform capabilities with single-source support.
These two languages are in many ways similar to SYCL, which you saw in Section 12.4. In fact, they borrowed concepts from each other, aiming for performance portability. Each provides libraries that are fairly lightweight layers on top of lower-level parallel programming languages. We'll briefly cover each of them.
5. Higher-level languages for performance portability
Kokkos: A Performance Portability Ecosystem
Kokkos is a well-designed abstraction layer for languages like OpenMP and CUDA. It has been in development since 2011. Kokkos has the following named runtimes. They are enabled in the Kokkos build using the appropriate CMake flag (or the Spack build option). Some of these are more mature than others.
5. Higher-level languages for performance portability
Exercise: Threading Triad in Kokkos
For this exercise, we built Kokkos with an OpenMP backend, then built and executed the threading triad example. Start by typing:
git clone https://github.com/kokkos/kokkos
mkdir build && cd build
cmake ../kokkos -DKokkos_ENABLE_OPENMP=On
Then change to the Kokkos threading triad source directory and build out-of-tree with CMake:
mkdir build && cd build
export Kokkos_DIR=${HOME}/Kokkos/lib/cmake/Kokkos cmake ..
make
export OMP_PROC_BIND=true
export OMP_PLACES=threads
Kokkos, built with CMake, has been optimized to be simple, as shown in the following listing. The Kokkos_DIR variable should be set to the location of the CMake configuration file for Kokkos.
5. Higher-level languages for performance portability
Adding the CUDA option to the Kokkos build generates a version that runs on NVIDIA GPUs. Kokkos can use many other platforms and languages, and more platforms and languages are constantly being developed.
The Kokkos stream triad example in Listing 12.23 has some similarities with SYCL in that it uses C++ lambdas to inline functions for either CPUs or GPUs. Kokkos also supports functors for this mechanism, but lambdas are less verbose for practical use.
5. Higher-level languages for performance portability
5. Higher-level languages for performance portability
5. Higher-level languages for performance portability
A Kokkos program begins with Kokkos::initialize and Kokkos::finalize . These commands initiate things needed for the execution space, such as threads. Kokkos is unique in that it embeds flexible allocations of multidimensional arrays as data projections that can be switched depending on the target architecture. In other words, you can use a different data ordering on a CPU than on a GPU. We use Kokkos::View in lines 14–16, although this is only intended for single-dimensional arrays. The real value comes with multidimensional arrays. The general syntax for Kokkos::View is:
View < double *** , Layout , MemorySpace > name (...);
5. Higher-level languages for performance portability
Memory spaces are optional for the template, but by default they correspond to the execution space. Some memory spaces are:
A layout can also be specified, although it has a default value corresponding to the memory space:
Computational kernels are specified using lambda syntax in one of three data parallelism templates:
5. Higher-level languages for performance portability
In lines 20, 23, and 29 of Listing 12.23, we used the parallel_for template. The KOKKOS_LAMBDA macro replaces the [=] or [&] capture syntax. Kokkos takes care of its definition and does so in a much more readable form.
RAJA for a More Adaptable Performance Portability Layer
The RAJA performance portability layer aims to provide portability with minimal disruption to existing Lawrence Livermore National Laboratory source code. In many ways, it is simpler and easier to adopt than other comparable systems. RAJA can be built with the following options:
As shown in the following listing, RAJA also has good CMake support.
5. Higher-level languages for performance portability
The RAJA version of the streaming triad undergoes only a few changes, as shown in the listing below. RAJA also makes extensive use of lambdas to ensure portability across CPUs and GPUs.
5. Higher-level languages for performance portability
5. Higher-level languages for performance portability
The required changes for RAJA consist of inlining the RAJA header file at line 2 and replacing the computation loop with Raja::forall . It's clear that the RAJA developers have ensured a low entry barrier to improve performance portability. To run the RAJA benchmark, we included a script that builds and installs RAJA, as shown in the listing below. The script then proceeds to build the stream triad source code using RAJA and executes it.
5. Higher-level languages for performance portability
6. Exercises
1. Change the host memory allocation in the CUDA streaming triad example to use pinned memory (Listings 12.1–12.6). Did you see any performance improvement?
2. In the sum reduction example, try an array of 18,000 elements, all initially assigned an index value. Run the CUDA source code, then the version in SumReductionRevealed. You may want to adjust the amount of output.
3. Convert the CUDA reduction example to HIP by HIPifying it.
4. For the SYCL example in Listing 12.20, initialize arrays a and b on the GPU device.
5. Convert the two initialization loops in the RAJA example in Listing 12.24 to the Raja:forall syntax. Try running the CUDA example.