Profiling on Nvidia GPUs
These notes document the steps we took to measure and improve the CanudaX_BSSNMoL thorn's performance on Nvidia GPUs. It is the hope that this guide can be useful starting point for developers of other CarpetX thorns to understand and improve the performance of their own thorns on GPUs.
Based on prior experience using the thorn, we should have an idea which routine takes the longest time in a simulation, and whose performance we want to improve. For example, the RHS evaluation is the most expensive part of the CanudaX_BSSNMoL thorn for a binary black hole inspiral simulation. This can be seen in the Cactus Timer output:
======================================================================
% Time/s Min/s Max/s Timer (gettimeofday)
======================================================================
100.0 28850.2 28681.7 28976.4 [0058] ODESolvers: ODESolvers_Solve_Subcycling in CCTK_EVOL
100.0 28849.2 28680.6 28975.3 CallFunction CCTK_EVOL: ODESolvers::ODESolvers_Solve_Subcycling
100.0 28848.1 28679.4 28974.1 ODESolvers::Solve
77.7 22413.2 22411.1 22416.3 EvolveRegrid
50.2 14470.8 14392.7 14637.1 Sync
24.1 6956.9 6765.0 7329.7 CallScheduleGroup
22.3 6441.8 5511.8 7181.9 ODESolvers::Solve::poststep
21.9 6321.6 6148.6 6674.4 ODESolvers::Solve::rhs
16.7 4827.3 4684.9 5101.3 [0059] CanudaX_BSSNMoL: CanudaX_BSSNMoL_calc_bssn_rhs in ODESolvers_RHS
16.7 4820.6 4677.7 5094.7 CallFunction CanudaX_BSSNMoL_RHSGroup: CanudaX_BSSNMoL::CanudaX_BSSNMoL_calc_bssn_rhs
In the above, the scheduled function CanudaX_BSSNMoL_calc_bssn_rhs is the most expensive part of ODESolvers_RHS group. It is a good target to profile and optimize, since we get better bang for our buck by speeding up this part of the thorn.
Note that when we offload the computation to the GPU, ideally we want to be able to saturate 100% of the GPU's computational throughput and only do work on the CPU when required, e.g. for input and output. If the required data for a GPU calculation is not available, then the GPU will remain idle until the data is available, and this data transfer between host (CPU) and device (GPU) has a very high latency.
With the above in mind, the first thing to look out for is excessive data transfers that happen outside of the beginning or the end of the simulation. In an evolution code using CarpetX on GPUs, most of the host/device memory transfer happens at the beginning and the end, since CarpetX enables AMReX's managed memory arena by default. This choice lets AMReX do all the memory transfers between host/device automatically, so that no explicit memory management is needed from the user side.
For GPU runs, we desire exactly one large continuous chunk of memory per GPU. To ensure this, we adjust max_grid_size_x,y,z in the parameter file such that each GPU has exactly one grid patch, and set max_tile_size_x,y,z to a very large number so that there is only a single tile. Oftentimes, if these parameters were not set appropriately, the GPU run could have very poor performance due to the overhead of launching kernels over many tiles. See the CarpetX thorn guide for a detailed explanation of these parameters.
Nvidia Nsight Systems measures the hardware use for the duration of the application, which can be a useful first step to understand the how different parts of the application make use of the hardware resources such as CPU usage, GPU usage, and memory bandwidth.
The timeline view presents hardware usage metrics over the duration of the application.
For a CarpetX simulation, it is easy to check which scheduled functions are being executed during the timeline by expanding the NVTX rows, and check the corresponding hardware usage for the duration of a particular function.
It may be useful to first look at the GPU's SM (streaming multiprocessor) throughput, SM warp occupancy, and memory bandwidth.
- If the average SM throughput is not high (say <60%), this might be a symptom that the GPU is waiting for data from the GPU
- If the SM throughput is high but the warp occupancy is low, then the GPU is not able schedule instructions to its full capability, possibly due to register pressure, waiting for global memory accesses on the GPUs, or branch divergence. In any case, this warrants a deeper look at the kernel level, and profiling the kernel using Nsight Compute can usually provide a good insight.
The stats system view presents execution statistics of different kernels or ranges of the program. For CarpetX simulations, it's probably most informative to start with the NVTX Start/End Range Summary. We can get an overview of the average, min/max, median, stddev of the duration of each NVTX range, which can let us know if any section of the application takes particularly long, or has a high variance in the execution time.
As seen in the Nsight Systems profiler output, CarpetX applications already have their scheduled functions annotated with NVTX. However, we may want to select the ranges to annotate in a more fine-grained manner, and the default annotated names can be long (e.g. CallFunction CanudaX_BSSNMoL_RHSGroup: CanudaX_BSSNMoL::CanudaX_BSSNMoL_calc_bssn_rhs!). Hence, we recommend adding custom NVTX ranges to mark out the start and end of a particular section of code you want to profile.
To set ranges in your application using NVTX annotation, a few source modifications are necessary in your source code:
Include the header in your C++ source code
#include <nvtx3/nvToolsExt.h>Mark the start and end of the range, which you can name.
const nvtxRangeId_t range = nvtxRangeStartA("RHS calculation");
[...]
nvtxRangeEnd(range);In the above example, we name this range "RHS calculation", and all of the code and GPU kernels executed in between the start and end are included in the range. For more advanced usage of NVTX annotation for Nsight Compute, you can find a list of examples with illustrations here
Finally, while not strictly necessary, it would be nice to wrap the above additions in an #ifdef guard so that the annotations are added only when compiled with nvcc, the CUDA compiler driver.
#ifdef __CUDACC__
// CUDA-specific additions here
#endifCompiling with different compilers, e.g. if you are compiling the code for CPU applications, will just result in the annotations being ignored.
Finally, we also suggest adding -generate-line-info to the CUDA compiler flags to your Optionlist. This allows you to resolve the source code in the Nsight Compute report and identify the lines in the code causing performance issues.
CUCCFLAGS = [...] -generate-line-info
Once we decide which routine to profile in detail, we should set up a small simulation which can be finished quickly and does not require a lot of memory. If the memory usage is too high, this can lead to the profiler hanging.
In our case, we set up our simulation such that we minimize other parts of the simulation which are not essential for the performance being evaluated:
- Use Cartesian Minkowski initial data, since it is quick to initialize and the RHS is evaluated independent of what the actual initial data is
- Set up a unigrid instead of multiple refinement levels
- Disable all analysis thorns and output
- End the simulation at iteration 1 With this setup, a simulation of 256^3 cells can finish in a few seconds on one GPU before doing any profiling, and finish in 10 minutes with the profiler attached.
Nsight Compute is a profiling tool that takes detailed hardware metrics at the kernel level, which can help the user understand performance bottlenecks and potential areas for optimization. For an introduction, we recommend checking out the recordings below:
- Short highlight of key features - Overview of functionalities
- Memory workload analysis - Tutorial explaining how optimizing memory use can impact performance
- Kernel Performance Analysis - Long tutorial session analyzing kernels with different performance issues
Let's say our CarpetX simulation is executed as follows
cactus_sim twopunctures.parand we want to profile our scheduled function in ODESolvers::Solve::rhs, which we had annotated in the NVTX range CanudaX_BSSNMoL_calc_bssn_rhs. Let's also say that there are two kernels launched in the scheduled function, each of which is launched via the loop_int_device function from the Loop thorn.
To launch the same simulation and attach Nsight Compute, we typically do something like the following:
ncu --set full --replay-mode kernel -s 2 -c 2 --kill yes \
--nvtx --nvtx-include "CanudaX_BSSNMoL_calc_bssn_rhs" \
cactus_sim twopunctures.parWhat the above command does is that it skips the first two kernel launches within the range, i.e. skipping the first ODESolvers::Solve::rhs instance. We make this skip because the first time the kernel is launched, there is extra work involving allocating and transferring data from CPU to GPU, so its performance is not representative of what happens during a simulation. This command then profiles the next two kernel launches, which happen during the second ODESolvers::Solve::rhs instance, and kills the application afterwards to end the task quickly.
The output will be shown in stdout, but one can also set the -o parameter to specify the output name, which will save it to a file ending in .ncu-rep. The file can be viewed in ncu-ui, the graphical version of Nsight Compute, and it contains much more information along with suggestions and explanations.
If you are profiling a parallel job over multiple MPI ranks on a cluster, you can target just one process using the following wrapper script
#!/bin/bash
# Use $PMI_RANK for MPICH, $OMPI_COMM_WORLD_RANK for openmpi, and $SLURM_PROCID with srun.
if [ $SLURM_PROCID -eq 0 ]; then
ncu -o %q{SLURM_JOB_NAME}.report.%q{SLURM_PROCID}.%q{SLURM_JOBID} \
--set basic --replay-mode kernel -s 1 -c 2 --kill yes \
--nvtx --nvtx-include "CanudaX_BSSNMoL_calc_bssn_rhs" \
"$@"
else
"$@"
fiThe script, which we save as ncu-wrap.sh, would then be run in a job script as follows:
srun bash ncu-wrap.sh bash -c "cactus_sim twopunctures.par"You may need to adjust the srun to mpirun, mpiexec, or something else depending on what your cluster supports.
Here we document the steps whereby we improved the performance of the RHS routine in CanudaX_BSSNMoL by profiling, making optimizations, and iterating over this process.
We started with the unoptimized RHS routine as our baseline result, which contained one giant kernel. The overall GPU utilization is poor because the compute and memory throughput are both very low. Usually a kernel would achieve either a high compute or memory throughput, and having neither hints at something deeper that prevents the GPU from efficiently doing the work.
In our case, our kernel has a very high GPU register usage that causes data to be spilled to global memory, which is off-the chip and has a very high latency. As a result, GPU threads can remain idle while waiting for data to be transferred back to the chip, which the profiler noted as "long scoreboard stalls".
| Duration | Compute throughput | Memory throughput | |
|---|---|---|---|
| 0 | 79.57 ms | 13.05 % | 46.97 % |
The first step we took was to separating the Ricci tensor calculation out as a standalone kernel, so we will profile two kernels at once in the following. Immediately, we see that the overall runtime improved from 79.57ms to 58.87ms.
We then went on several rounds of profiling and optimizing bits of the code by reducing temporary variables, making use of data locality, increasing cache hitrates. To identify hotspots in the source code, we checked for Live register use in the "Source" section of the profiler report. The optimizations where we saw the biggest improvements were from #6 (replacing division), #7 (Fuse derivative computations), and #8 (Refactor the Ricci calculation).
Left: Ricci tensor kernel Right: RHS calculation kernel reusing the Ricci tensor
| Duration | Compute throughput | Memory throughput | Duration | Compute throughput | Memory throughput | ||
|---|---|---|---|---|---|---|---|
| 1 | 27.96 ms | 18.00 % | 64.13 % | 30.91 ms | 23.30 % | 27.23 % | Split RHS kernel into Ricci and (rest of RHS) |
| 2 | 23.64 ms | 14.55 % | 61.82 % | 30.71 ms | 23.21 % | 27.37 % | Reduce number of temporary arrays |
| 3 | 22.93 ms | 14.99 % | 62.98 % | 30.88 ms | 22.41 % | 23.32 % | Reorder the Christoffel indices |
| 4 | 22.76 ms | 15.03 % | 63.11 % | 27.16 ms | 30.30 % | 18.56 % | Reorder calculations to move reusesd arrays closer in code |
| 5 | 23.53 ms | 14.82 % | 62.75 % | 23.17 ms | 30.37 % | 23.33 % | Compute finite difference first thing in kernels |
| 6 | 22.08 ms | 14.90 % | 65.71 % | 12.76 ms | 25.53 % | 45.16 % | Replace division with multiplication by reciprocals |
| 7 | 16.99 ms | 17.10 % | 57.01 % | 12.16 ms | 25.76 % | 37.62 % | Fuse computation of 1st and 2nd derivatives of the same grid functions |
| 8 | 6.64 ms | 33.40 % | 43.62 % | 12.06 ms | 25.69 % | 37.51 % | Refactor the Ricci calculations into smaller loops |
| 9 | 6.69 ms | 33.56 % | 43.84 % | 11.29 ms | 26.23 % | 40.44 % | Make use of tensor symmetries and not recompute the |
| 10 | 6.39 ms | 34.82 % | 46.20 % | 11.25 ms | 26.02 % | 46.20 % | Reduce the block size of kernel launch from 256->1280 |
At the end of the optimizations, we are still not able to move either of the kernels into the compute bound region (compute throughput >60%). Some of the reasons include the still-high register pressure, uncoalesced global memory loads coming from computing finite-difference stencils, and uncoalesced writes to grid functions to global memory. However, we ended up with a total execution time of 17.64 ms, a significant (~4x) speedup over the baseline of 79 ms!