Finite difference method – Laplacian part 4

Performance across different hardware and sizes

In addition to the MI250X GPU, the following AMD GPUs are considered:

1. RX 6900 XT

2. RX 7900 XTX

3. MI50

4. MI100

5. MI210

6. MI250

We will perform a scaling study for Kernel 5 starting with a problem size nx,ny,nz = 256, 256, 256 across all of the above GPUs. Each dimension is increased in multiples of 256 until we reach a global problem size of nx,ny,nz = 1024, 1024, 1024. Recall that the Figure of Merit (FOM) is defined as the effective memory bandwidth:

theoretical_fetch_size = ((nx * ny * nz - 8 - 4 * (nx - 2) - 4 * (ny - 2) - 4 * (nz - 2) ) * sizeof(double); 
theoretical_write_size = ((nx - 2) * (ny - 2) * (nz - 2)) * sizeof(double); 
effective_memory_bandwidth = (theoretical_fetch_size + theoretical_write_size) / average_kernel_execution_time;

This is the average rate of data transfer considering the least amount of data that needs to move through the entire memory subsystem. The objective here is to obtain a FOM as close to the achievable device HBM bandwidth as possible. Through prior experimentation of rerunning Kernel 5 across various thread block sizes, launch bounds, and tiling factors, we found the following parameters provided the highest FOM for each individual GPU on a 512 \times 512 \times 512 grid of points:

From here on out, experiments on both the MI250 and MI250X GPUs shall be conducted on a single GCD. The figure below depicts the FOM[1] across the various AMD GPUs and problem sizes.

Figure 1: Performance of Kernel 5 across various AMD GPUs and problem sizes

Note that the RX 6900 XT GPU does not have enough memory to execute the largest example nx,ny,nz = 1024, 1924, 1024 so it is not shown in the figure above. Performance across all sizes on the RX 6900 XT and RX 7900 XTX GPUs appear relatively consistent, whereas all Instinct™ GPUs exhibit degraded performance when the problem size crosses a certain threshold. In the next sections we will investigate the source of the performance degredation on Instinct™ GPUs.

Root cause of performance degradation

Let us refocus our attention back on the MI250X GPU. To understand the cause of the performance drop, we need to collect rocprof stats. We will first collect the FETCH_SIZE and compare it against the theoretical_fetch_size for each combination of nx,ny,nz. Ideally, the FETCH_SIZE should be equivalent to the theoretical amount as the elements loaded by each thread could be reused up to six times by neighboring threads. It should be noted that the MI100, MI210, MI250, and MI250X GPUs all have 8 MB of L2 cache shared among the compute units (CUs) per GCD, so the L2CacheHit we collect will be useful in determining how well cached data is reused. The figure below illustrates our findings across the same spectrum of problem sizes:

Figure 2: Fetch efficiency and L2 cache hit rates of Kernel 5 on a single MI250X GCD across various problem sizes

There is a clear correlation between the FOM, fetch efficiency, and L2 cache hit rate. Recall from the Laplacian Part 1 post that the baseline HIP kernel exhibited a 50% fetch efficiency while maintaining a relatively high L2 cache hit rate of 65%. This led us to believe that there was a problem with the kernel design. In this last experiment, both the fetch efficiency and L2 cache hit rate leveled around 33%, suggesting that the limiting factor for performance is the size of the L2 cache.

Another interesting observation from the above figure is that the dips in performance occur only when either nx or ny increases – increases in nz plays virtually no role in the performance. As seen from the thread block configuration and grid indexing in Kernel 5 above, all elements from a single xy-plane will be fetched from global memory before fetching elements from the next xy-plane. Stencil computations per thread have z-direction strides of -nx*ny, 0, and nx*ny so three xy-planes of data are needed. If the three xy-planes are too large to simultaneously fit int the L2 cache, then each thread will cycle the cache to fetch extra xy planes from global memory – for example, when only a single xy-plane fits (or partially fits) the L2 cache, the thread needs to cycle through two extra xy planes from global memory hence fetch efficiency dropping to as low as 33%.

Let us conduct a second experiment to further illustrate this point. We shall examine various values of nx and ny to identify the optimal z-direction stride before performance begins to degrade. From a starting z-stride of ~1 MB (equivalent to nx,ny = 512,256) with nx fixed to three different sizes, we incrementally increase ny until the xy-plane exceeds 8 MB (nx*ny*sizeof(double)). See the figure below:

Figure 3: Performance impact from an increasing nx*ny stride for various fixed values of nx on a single MI250X GCD

Regardless of the size of nx, performance degrades when a single xy-plane exceeds ~2.5 MB in memory. It turns out 2.5 MB is small enough to fit three consecutive xy planes in an 8 MB L2 cache, which explains why the FOMs are higher in value. However, the larger the stride the more elements start dropping out of of cache thereby lowering the L2CacheHit rate and subsequently the FOM. Performance levels out after nx*ny exceeds 8 MB because at this point only one plane fits (or partially fits) the cache – this is where the FETCH_SIZE is three times larger than the theoretical estimate.

Addressing the cache size limitation

So how do we address the issue of the cache sizes for these larger problems? This is especially important on hardware like the MI50 where the LLC is only 4 MB of L2 cache. Even if nz is small, large values of nx and ny, the L2 cache could become a memory traffic bottleneck. Let us consider when nx,ny,nz = 1024,1024,1024 for the next several experiments on an MI250X GCD. Below are two possible strategies for addressing L2 cache bottleneck:

dim3 grid(
    (nx - 1) / block.x + 1, 
    (ny - 1) / (block.y * m) + 1, 
    (nz - 1) / block.z + 1);

dim3 grid(
    (ny - 1) / (block.y * m) + 1, 
    (nz - 1) / block.z + 1, 
    (nx - 1) / block.x + 1);

int i = threadIdx.x + blockIdx.x * blockDim.x;
int j = m*(threadIdx.y + blockIdx.y * blockDim.y);
int k = threadIdx.z + blockIdx.z * blockDim.z;

int i = threadIdx.x + blockIdx.z * blockDim.x;
int j = m*(threadIdx.y + blockIdx.x * blockDim.y);
int k = threadIdx.z + blockIdx.y * blockDim.z;

 
dim3 grid(
    (nx - 1) / block.x + 1, 
    (ny - 1) / (block.y * m) + 1 , 
    (nz - 1) / block.z + 1);

...

    
laplacian_kernel<<<grid, block>>>(d_f, d_u, nx, ny, nz, 
    invhx2, invhy2, invhz2, invhxyz2);

int bny = (ny - 1) / 4 + 1;
dim3 grid(
    (nx - 1) / block.x + 1, 
    (bny - 1) / (block.y * m) + 1 , 
    (nz - 1) / block.z + 1);

...

for (int i = 0; i < 4; i++)
  laplacian_kernel<<<grid, block>>>(d_f, d_u, nx, ny, nz, 
      invhx2, invhy2, invhz2, invhxyz2, bny*i);

__global__ void laplacian_kernel(T * f, 
    const T * u, int nx, int ny, int nz, 
    T invhx2, T invhy2, T invhz2, T invhxyz2) {

...

int j = m*(threadIdx.y + blockIdx.y * blockDim.y);

__global__ void laplacian_kernel(T * f, 
    const T * u, int nx, int ny, int nz, 
    T invhx2, T invhy2, T invhz2, T invhxyz2, int sy) {

...

int j = sy + m*(threadIdx.y + blockIdx.y * blockDim.y);

Summary and Recap

In this part, we first examined the performance of Kernel 5 across various AMD GPUs and problem sizes and noted that there appeared to be sharp drop off in performance when the grid crosses a certain size threshold. It was determine through experimentation that the size of the LLC was the limiting factor of performance for problems with large xy-planes. Two different workarounds were presented to circumvent the issue regarding the cache size, both of which required modifications to only a few lines of code.

Throughout the course of the finite difference method Laplacian series, we started with a simple implementation of the HIP kernel and incrementally applied several different optimizations to dramatically improve the performance of the kernel. We shall now quickly recap all the different optimizations undertaken in this work:

1. Loop tiling (Part 2): If a kernel is known to load multiple reusable values, consider adding loop tiling aka loop unrolling. Several steps are required to refactor the code, but this ultimately reduces the number of thread blocks launched and increase the number of stencils computed per thread. The FETCH_SIZE could be greatly reduced with this approach.

2. Reordered access patterns (Part 2): Frequently accessing elements of a device by going both forward and backward in the memory address may prematurely evict reusable data from cache. The kernel should reorder all load instructions in the kernel to access data in a single direction i.e., accessing memory by ascending address. This optimization is needed in order to properly leverage the loop tiling optimization.

3. Launch bounds (Part 3): If the register usage of a kernel is extremely high, applying launch bounds may enable the compiler to make more appropriate decision in register allocation, thus improving the register pressure. The default is 1024 threads, but we strongly recommend setting this to the number of threads your kernel will be using.

4. Nontemporal loads (Part 3): Elements of the output array f are not shared by any other threads or kernels, so the builtin intrinsics allow the finite difference stencil computation to be written back to device memory with an increased priority for eviction from cache, thus freeing up more cache for the input array u.

5. Grid block configuration (Part 4): A quick way to circumvent problems that spill over the L2 cache can be to configure the thread block or grid indexing configuration. This lets smaller segments of the xy-plane fill the cache first allowing some threads to perform computations without having to fetch from global memory multiple times.

6. Splitting into subdomains (Part 4): If the entire problem is bottlenecked by the LLC size, consider splitting the problem into several subdomains and launching the HIP kernels to process each subdomain sequentially. This method is robust and can be tailored as needed for any LLC or problem size.

The finite difference stencil for the Laplacian on a structured grid is relatively straightforward to implement in HIP. In a future lab notes post, we will consider other partial differential equations like that of the acoustic wave equation which will require higher order stencils than the ones shown here. This most likely requires optimization strategies not outlined in these last four posts. Nonetheless, we believe all the optimization strategies outlined in this Laplacian series should be leveraged to get the best performance out of stencil kernels.

Accompanying code examples

If you have any questions or comments, please reach out to us on GitHub Discussions