Performance Tuning for GPU-Embedded Systems: Machine-Learning-based and Analytical Model-driven Tuning Methodologies.

Regarding performance studies about GPU-embedded systems, previous works, such as those in [5, 6, 7], have provided valuable performance insights into the Jetson Tegra TK1. In [8], a performance characterization of both the Jetson TK1 and TX1 is conducted using roofline models, focusing on a matrix-multiplication application. However, the simplicity of the matrix multiplication fails to capture the complexities of real-world workloads. In [9, 10, 11], machine learning applications on the Jetson TX1, Nano, and AGX Xavier are respectively explored, offering significant performance insights but with a narrow focus on specific machine learning tasks. While some works have emphasized energy efficiency, comprehensive performance benchmarking has been lacking, except for [12], which also includes a performance variability analysis for the Jetson AGX Xavier. Notably, few studies have focused on performance, and the existing ones mostly rely on benchmarking approaches rather than developing performance models that can be applied to other embedded platforms.

In the field of desktop and server GPUs, there are examples of empirical tuning projects for the FFT such as [13, 14, 15, 16, 17, 18, 19]. GPU tridiagonal system solvers (TS) have been proposed in [20] and [21]. Nevertheless, these TS approaches do not exhibit significant optimization for modern architectures. More recent GPU TS can be found in [22, 23, 24, 25, 26, 27], while [28, 29, 4] introduce optimized GPU implementations for the scan primitive. Some tuning libraries, such as NVIDIA’s CUFFT [30] for FFT, CUSPARSE [31] for tridiagonal system solvers, and the CUB library [32] for the scan primitive, provide efficient implementations and are considered the state-of-the-art for these operations. In the realm of server GPU architectures, performance modeling has been explored in [4], while we extend their outlined principles to the context of GPU-embedded systems in this work, as detailed in the subsequent sections.

With respect to generic autotuning, where the goal is to tune any objective function rather than focusing on specific algorithms or platforms, we can find autotuners based on empirical searches such as OpenTuner [33] that use empirical approaches, being very slow. Some other generic auto-tuners, such as Kernel Tuner [34], implement various search optimization algorithms, but most of them have been demonstrated to not be more efficient than a random search [35]. Other research efforts have developed machine learning models [36, 37, 38, 39], but using random sampling indiscriminately to build the surrogate model. Instead, Bayesian optimization identifies the most informative candidates to train the surrogate model and achieves the same accuracy with fewer samples [40] [41]. Two generic auto-tuners, GPTune [42] and DeepHyper [43] use Bayesian optimization and are recognized as cutting-edge autotuners for HPC applications.

In [4], BPLG authors presented an analytical tuning methodology for small-medium parallel-prefix problems on CUDA desktop and server GPUs, while our work extends this methodology to tune embedded GPUs. When considering performance modeling for GPU-embedded systems, several factors come into play, such as the limited number of SMs, reduced memory bandwidth, and constrained power consumption. GPU-embedded systems necessitate a distinct approach to performance modeling.

We adhere to the identical set of performance premises and parameters established in [4]: (i) a reduction of high-latency communications; (ii) a high warp occupancy to hide latency, also taking into account a high block parallelism rate per SM; and (iii) a high granularity of work performed by threads (high instruction level parallelism), being aware of register consumption. The attainment of these premises can be accomplished by specifying certain values to the set of performance parameters exposed within the BPLG library: $S$, the number of elements stored in shared memory per threadblock; $L$, which is the number of threads per threadblock; and $P$ which is the number of elements processed by each thread. They are also related by $S=P\times L$ (except when shuffle operations can replace the communication pattern based on shared memory). These three parameters are represented by the tuple $(S,P,L)$ and their values depend on both the algorithm, the target architecture and $N$. Also, this methodology concentrates on batch execution to effectively utilize all available GPU resources, where each invocation to the library simultaneously executes $G$ batches of size $N=r^{n}$, $r$ power of two.

The application of the existing performance premises can be summarized in the following manner. The first premise looks for implementing coalescing patterns to achieve the maximum memory bandwidth. This is crucial in embedded GPUs where the memory bandwidth is reduced. Threads loading Node operators in a strided manner can help to achieve this. Regarding the second premise, GPU parallelism can be determined in terms of the number of threadblocks per SM (SM block parallelism), or by the number of warps per SM (SM warp parallelism). A trade-off must be sought between them depending on the problem features. The number of warps per SM, $W_{SM}$, is limited by the maximum number of warps per SM ($W^{max}$ architectural defined) and depends on the number of active threadblocks ($B_{a}$) being simultaneously executed per SM. Also, $B_{a}$ is bounded by the architecture design ($B^{max}$) and is constrained by the registers and shared memory available in the SM. With respect to the third premise, one thread is responsible for computing at least one Node operator from the parallel prefix circuit, where $P=max(fan\_{in},fan\_{out})$. It may be interesting to process more Node operators per thread, i.e. increase $P$, if extra register consumption does not affect the SM occupancy. This would reduce the number of required threads to compute the same amount of work. We could also increase the radix $r$ of the algorithm, when the communication pattern allows it. Each thread would continue working with a single Node operator, but its fan_in and fan_out values are increased, and consequently $P$ is also increased. This would change the definition of the algorithm owing to $N=r^{n}$, decreasing the number of steps taken.

The original methodology for desktop and server GPUs in [4] looks for maximizing both SM block and warp parallelism, and balancing them in the cases where is not possible to achieve their maximum values. Nevertheless, keeping higher warp occupancy is more important in CUDA embedded GPUs, compared to desktop or server GPUs due to several reasons. Higher warp occupancy reduces the number of idle resources and helps maximize the utilization of available hardware, resulting in more efficient computation with lower power consumption and reduced thermal output, leading to longer battery life. Additionally, many embedded applications, such as robotics, autonomous vehicles, and drones, have strict real-time requirements and low-latency demands. Higher warp occupancy helps ensure more consistent and predictable execution times, minimizing variations in thread execution and reducing overall latency. Prioritizing having a higher warp occupancy can yield greater benefits in GPU-embedded systems compared to having a high number of active thread blocks per SM due to high memory access latency. By keeping a higher warp occupancy, GPU performance can be maximized by effectively hiding memory latency and minimizing the impact of memory stalls, leading to improved overall throughput and resource utilization. Once a certain level of warp occupancy is achieved, hiding latency is not the main issue anymore. After that, focusing on increasing instruction level parallelism (ILP) can be particularly important for achieving optimal performance compared to server GPUs where resources are relatively more abundant. Thereby, increasing the radix, even when slightly reducing $B_{a}$, it reduces the number of taken steps (consequently synchronization points) and increases ILP.

Following this idea, this work proposes a performance guideline that we have found very accurate based on extensive empirical observations and experience with these devices. Specifically, our tuning guideline for GPU-embedded systems can be summarized as:

• •

  Choose $(S,L,P)$ values that achieve both $W_{SM}=W^{max}$ and $B_{a}=B^{max}$. If this is not possible:

• •

  Find the tuple that maximizes $B^{a}$, while the ratio $W_{SM}/W^{max}$ keeps between $[60\%-100\%]$. If this warp occupancy ratio cannot be reached:

• •

  Maximize the warp occupancy. If there are several possibilities, take the one that maximizes $P$.

• •

  In all previous steps, if the algorithm allows increasing its radix $r$, then select the configuration that increases $r$ even when reducing $B_{a}$.

Considering the GM20B GPU and proposed guideline, it is necessary to find the optimal values of the $(S,P,L)$ parameters for each algorithm. Figure 3 (a) summarizes the parallelism and the resource consumption achieved by different threadblock configurations.

We use two sets of parameters participating in the tuning search [39]: Input Parameters, $\mathcal{A}$, which characterize the input in terms of computational shapes or data layout of the application; and Performance Parameters, $\mathcal{B}$, the parameters to be optimized for the target architecture. Given an input $a_{j}\in\mathcal{A}$, our goal is to search for the optimal vector of performance parameters $x_{B}\in\mathcal{B}$ that meets $\operatorname*{arg\,min}_{x_{B}\in\mathcal{B}}f(X):X=(a_{j},x_{B})$. Here, $f(\cdot,\cdot)$ denotes the objective function that maps each configuration to a corresponding execution time. Calculating derivatives for $f$ is indoable [43], necessitating actual evaluations of different configurations. In this methodology, a Bayesian optimization search is employed, as it has been shown to be effective in exploring promising regions of the parameter space [49], thereby minimizing the number of required evaluations. The procedural workflow for Bayesian optimization is outlined as follows: First, a small set of configurations are randomly sampled from the search space and evaluated. The resulting data (configurations and execution times) are added to a dataset that trains the surrogate model. Based on the actual predictions of the surrogate model, the acquisition function optimizes a score, guiding the selection of the next configuration to be evaluated. The iterative process continues until a stopping criteria is met. To do this Bayesian-optimization search, we use the GPTune framework [42] that uses the Linear Coregionalization Model (LCM) [50] as surrogate model, and the Expected Improvement acquisition function [51].

As a summary of the process, the first step is to define the corresponding Input Parameters, $\mathcal{A}$, and Performance Parameters, $\mathcal{B}$, together with any constraint in order to define the search space. Then, this search space is explored with the searching workflow, which orchestrates the search by calling the Bayesian-optimization framework. The interactive search is run for each algorithm on the target platform until a user-defined stop criteria is met. For our BPLG parallel-prefix search, the Input and Performance parameters are described in Table I. The problem size describes the input problem we are solving, and its range of values depends on the BPLG implementation for each operation. Regarding the performance parameters, the range of taken values, which defines the search space, also depends on the parallel-prefix algorithm implementation being tuned, as it is seen in the next section. We also added shuffle as a binary performance parameter that chooses an alternative communication based on warp shuffles instead of shared memory, when available. With respect to the stopping criteria, we perform a sliding-window check and stop the search when no progress has happened in the last iterations, which saves extra evaluations once a good optimum has been found. In this search, we stop the search if no progress within the last 5 evaluations. Also, we set a high execution-time value for those executions with configurations that are invalid or are not finishing after 1 minute.

Both methodologies can be extended to larger problem sizes. When the problem size exceeds the capacity of shared memory ($N>S$), collaborative computation among multiple thread blocks becomes necessary. Two options are available in recent CUDA versions to synchronize these thread blocks: launching multiple kernels using a multi-kernel strategy or utilizing Cooperative Groups (available from CUDA 9.0 onwards). However, the number of thread blocks that can participate in Cooperative Groups is limited by the maximum number of resident thread blocks per Streaming Multiprocessor (SM). The resource consumption of BPLG kernels can significantly reduce the number of resident thread blocks per SM, resulting in fewer threadblocks than the number required to complete the collaborative computation. BPLG implements FFT for problem sizes exceeding shared memory capacity, and our work aims to tune this implementation as well.

In their work [52], the authors of BPLG presented an analytical model for tuning BPLG large-sized parallel-prefix operations on desktop and server GPUs, employing a multi-kernel strategy. Building on this idea, we extend our analytical heuristics for small to medium sizes on CUDA-embedded systems with the premise that The number of kernels $m$ needs to be minimized. If $N=r^{n}$ and $S=r^{s}$, then the number of required kernels $m$ can be calculated as $m=\lceil\frac{n}{s}\rceil$. Once the optimal $S$ is determined using this expression, our previous guideline for small to medium problems is applied to tune the $(S,P,L)$ and radix $r$ values for each resulting kernel. The inclusion of multiple kernels has introduced additional complexities in the tuning procedure, as the $(S,P,L)_{m}$ values are interdependent. The value of these parameters in one kernel impacts the workload of others.

It is noteworthy to mention that the ML-based methodology will treat the tuning process as a black box. We just have to extend the parameter space, but the underlying intricacies and interdependencies among kernels are transparent to this purpose.

Firstly, it should be pointed out that each element is an equation composed of 4 simple-precision coefficients, and only WM allows modifying the radix definition. With respect to the analytical-based tuning, the optimal configuration which maximizes both block and warp parallelism, following Figure 3(a), is achieved with $L=64$ and fewer than 32 registers per thread (the 3rd row in Table). Starting with the CR pattern and taking into consideration additional variables and index calculation, $P$ must be equal to 2 in order to consume less than 32 registers per thread. When $N\leq 64$, each problem is solved with 32 threads and $P=2$ elements per thread. Thus, shuffle instructions can be employed within the warp instead of shared memory. Please observe that each threadblock is executed with $L=64$ threads, thus 2 batches are simultaneously solved. When $N=128$, there are two alternatives: either using $P=4$, 32 threads per problem (2 problems per threadblock) and shuffle communications; or $P=2$, with 64 threads per problem (one problem per threadblock) and the use of shared memory to perform the communications ($128\ el.\ \times 4\ coef.\ \times 4\ Bytes=2048\ Bytes$). In the first case, 63% of warp parallelism and 20 threadblocks are achieved (4th row in Figure 3(a)); whereas the second case obtains 100% of warp occupancy, but only 16 threadblocks (5th row in Table) and less bandwidth in communications. Following our performance recipe, the first alternative is chosen. For the remaining problem sizes, ($N>128$), shared memory is needed for storing the problem data; thus $S=N$. When $S$ occupies more than 3072 bytes, the case of $N>128$, it is not possible either the maximum warp parallelism or the maximum block parallelism, choosing the configurations that meet the tuning strategy requirements. Due to the high register consumption, $P=2$ is established for these problem sizes. The same tuning values are obtained for PCR. In the case of LF, each element is composed of two equations ($2\times 4\times 4$ bytes per element). As the number of registers must not exceed 32, considering registers employed for additional variables, $P$ must be strictly 2, but the same performance parameters are obtained despite of this restriction following Figure 3(a). Finally, WM allows modifying the radix of the algorithm due to its regular communication pattern. Hence, in the case of having several reasonable configuration alternatives, increasing $P$ must be prioritized. Specifically, when $N\leq 128$, $(S,P,L)=(0,4,64)$ is used, achieving 75% warp occupancy and 24 out of 32 threadblocks, performing shuffle communications within each warp. For the remaining sizes, $(S,P,L)=(N,4,N/4)$ is used.

Regarding the ML-based tuning, the methodology will search for $L$ values between 32 and 1024; $P$ values among $\{2,4,8\}$ (higher values are unrealistic due to high register pressure) and $S$ can take values up to $2048$ elements (32 KB per threadblock). Although it is possible to use up to $48KB$ per threadblock, that would not represent any $N=r^{n}$ size power of radix. The radix $r$ is set to 2 for CR, PCR and LF, since larger values are not allowed for these patterns, but the model does explore the set $r=\{2,4,8\}$ in WM. The $shuffle$ parameter is only allowed when $N/P\leq 32$, i.e. the computation of each problem can be performed within a warp. Also, the restrictions (!shuffle OR S==0) and (!shuffle AND S=PxL) have to be met, because the shuffle optimization releases the need of shared memory. Figure 4(a) shows the suggested configuration for the WM algorithm, together with the number of candidate evaluations required to train the surrogate model. Radix $r$ consistently aligns with the recommended $P$ value, and the radix $r$ is omitted in the figure. Also, shuffle optimizations were applied whenever $S=0$ in the figure. The higher problem sizes, the fewer configurations are valid, reducing the number of evaluated candidates in the search.

BPLG implements the scan primitive with shuffle instructions where 32 elements per problem are stored in shared memory [4]. Examining Table 3(a), the maximum warp and block parallelism is achieved with $L=64$ and fewer than 32 registers per thread. Considering auxiliary variables and index calculation, $P$ must be less than or equal to 4. Additionally, shared memory consumption per threadblock must be fewer than 2048 bytes; thus, $S\leq 512$ single-precision elements. As $L=64$ and $P=4$, the number of problems per threadblock can be expressed with $\frac{64}{N/4}=256/N$. Since each problem uses 32 shared memory elements, $S=32\times\frac{256}{N}=8192/N$. When $N>256$, $L$ has to be increased and only one problem can be solved per threadblock, thus $L=N/4$ and $S=32$. The KS pattern obtains the same values by following the same reasoning. Figure 3(c) summarizes the performance parameters for the LF and KS patterns.

With respect to the $(S,P,L)$ exploration in the ML-based tuning, the search space is similar to the one seen with tridiagonal solvers. However, the radix of these algorithms cannot be explored, while $P$ can take any value from $P=\{2,4,8,16,32\}$. Due to the given BPLG implementations, $shuffle$ is always set to one, and $S=L/32$. Therefore, the ML-based tuning search only explores values for $P$ and $L$ parameters, and when raising $N$, the amount of possible configurations is highly reduced due to restrictions, leading to a minimal number of candidate evaluations. Figure 4(b) displays the suggested configurations for the scan LF algorithm.

The FFT communication pattern allows extending its radix. Thus, increasing $P$ reduces the number of computing steps, internal communications and synchronizations. Referring to Figure 3(a) as a point of reference once again, $L=64$ should be chosen. Using $P=4$ exceeds the number of 32 registers, but the benefits of having $P=4$ are high in terms of reducing computing steps, as mentioned in Section V-A. For $N\leq 256$, there are several problems being solved in each threadblock, as many as $\frac{64\times 4}{N}$. As $S=P\times L$ and each element is composed of two floats, $S=256$ elements, achieving 75% of warp parallelism. If $S=4096$, it implies $4096\times 2\times 4bytes=32KB$ of shared memory per threadblock and leads to only one active threadblock per SM. To avoid it, BPLG first exchanges the real part, performing the computations, and then exchanging the imaginary part, reducing the shared memory to 16KB per threadblock. For the remaining sizes, $L=N/4$ and $S=N$, as Figure 3(d) shows.

Regarding the ML-based search, radix $r$ can take any value between $r=\{2,4,8,16\}$, and consequently $P=\{2,4,8,16\}$. Also, $shuffle$ is always set to 0, as BPLG does not support shuffle communications for the FFT implementations. Considering $8$ bytes per element and the explained shared-memory multiplexing technique, $S\leq 8192$ elements. Figure 4(c) shows the achieved configurations for the FFT algorithm with the ML-based search. It always coincides the suggested radix $r$ with the recommended $P$ value, resulting in the omission of the radix $r$ in the Figure.

We consider problem instances within the range of $4096<N\leq 8388608$ elements to represent this type of problems. Starting with the analytical-based tuning, $S$ must be equal to 2048 (the maximum number of elements allowed) to minimize the number of required kernels with the multi-kernel strategy. As $2048=P\times L$, if we apply the given performance guideline, then the tuple $(2048,8,512)$ with radix-8 is chosen as maximizes the number of active warps (75%) per SM over other solutions for this equation. When $N\geq 524288$, three kernels have to be launched.

For small and medium problem sizes ($N\leq 4096$), the ML-based search could be deemed overkill, as a exhaustive search with few evaluations could have sufficed. However, when the search space grows, the ML-based approach fits better within the study work, offering a coherent and efficient approach to tackle the search. Figure 4(d) shows the required number of candidate evaluations to train the model. When three kernels must be tuned, the exhaustive search have to explore hundreds of valid combinations, while the ML-based methodology finds a good optimum with few evaluations.

In the context of tridiagonal systems, the data performance is quantified in MRows/s, utilizing the formula: $N\cdot b\cdot 10^{-6}/t$, where each equation operates on four coefficients. Table II presents the achieved performance for the CR, PCR, LF, and WM implementations in the BPLG library using the two proposed methodologies. Furthermore, the $\Phi$ metric is computed with the optimal configuration obtained through an exhaustive search. The results demonstrate that both methodologies yield comparable performance, closely approaching the optimal configuration identified by the exhaustive search. The analytical model developed for these computation patterns exhibits excellent performance, while the ML-based search also proves highly effective, as few combinations are valid with high problem sizes requiring minimal exploration to identify the global optimal solution. It is worth noting that slight variations in execution times among GPU runs account for most of the observed differences. To mitigate variability, we conducted 100 executions for each configuration, thereby reducing the impact of such fluctuations.

Figure 5 illustrates the performance of the BPLG proposals in comparison to the CUSPARSE library. Notably, the jagged outline observed for BPLG-WM stems from the CUDA implementation. When the number of elements in a thread workload (radix $r$) does not evenly divide the total number of elements ($N$), a mixed-radix technique is employed. This technique involves performing the computation in multiple steps, with a lower $r$ value for the first step and the given $r$ value for the remaining steps. Consequently, additional synchronizations are required, and the limited shared memory availability restricts the number of simultaneous threadblocks per SM. On the other hand, CUSPARSE exhibits subpar performance, which is not exclusive to the Jetson platform. Since CUDA 7.0, the performance of the CUSPARSE library for solving tridiagonal systems has experienced notable degradation.

The performance evaluation of the Scan primitive is conducted in terms of MData/s, using the formula $N\cdot b\cdot 10^{-6}/t$. Table II demonstrates that both methodologies suggest configurations that align with the global optimum. Figure 6 provides an analysis of the performance of tuning BPLG implementations in comparison to the CUB library. Notably, the BPLG implementations exhibit similar performance characteristics, maintaining a constant amount of shared memory regardless of the problem size $N$. This is accomplished through the utilization of shuffle instructions for interthread communications, ensuring consistent performance across different problem sizes.

The default implementation of the CUB library does not inherently support multi-batch execution. However, it can be adapted by applying a segmented-scan transformation technique [54]. Nonetheless, invoking the library multiple times yielded the best performance in our experiments. Therefore, in the interest of fairness, the experimental results were obtained by invoking the CUB library $b$ times, resulting in an average performance of $14.97$ MData/s. When compared to the BPLG implementations, BPLG outperforms CUB for problem sizes smaller than $N=2048$. This behavior aligns with expectations, as the number of library invocations corresponds to the number of executed batches. With a greater number of smaller problem sizes, there are more library invocations, incurring additional overhead. Conversely, for larger problem sizes, the number of batches decreases, leading to improved performance. Despite the peak of performance for large problem sizes, BPLG implementations obtain an average speed-up of $1.22x$ over the CUB library.

The performance evaluation of the complex FFT is quantified in GFlops/s using the well-established formula: $5N\cdot log_{2}(N)\cdot b\cdot 10^{-9}/t$, where $N$ represents the input size, $b$ denotes the number of processed batches, and $t$ is the time measured in seconds. The obtained performance results for each tuning methodology are summarized in Table II, showcasing competitive configurations suggested by both methodologies. Figure 7 provides a comprehensive comparison between the BPLG implementation and the cuFFT library in terms of performance. It is evident that the tuning BPLG version achieves highly comparable performance to cuFFT for problem sizes up to $N=2048$. However, for larger problem sizes, the shared memory limitation becomes the primary bottleneck for the BPLG version. Notably, while this Maxwell architecture offers an increased shared memory capacity of 64 KB per SM, each threadblock can only utilize up to 48 KB. Consequently, a performance drop is observed beyond $N>2048$. It is important to note that despite this decline, the cuFFT library maintains comparable performance on average with the BPLG-FFT implementation.

In the context of small to medium problem sizes, the analytical model has consistently demonstrated comparable or superior performance compared to the ML-based search. This finding emphasizes the pragmatic advantage of utilizing the analytical model, as it eliminates the need for costly evaluations of numerous candidates to train a surrogate model. The ML-based search may not have appeared as a favorable option due to the relatively small search spaces, making exhaustive search or analytical modeling more practical. However, as depicted in Table II, a notable shift in behavior emerges for larger problem sizes. This discrepancy can be attributed to the increased complexity of the optimization problem, where up to three kernels require tuning, leading to a considerably larger search space. Consequently, the ML-based search becomes a more viable and beneficial approach. Figure 8 shows the performance of the tuned BPLG implementation compared to cuFFT, where comparable performance is achieved. The highest difference in performance happens at $N=8192$, for which cuFFT launches only one kernel while BPLG launches two. This evidences the overhead associated to increase the number of kernels in a computation. The growing trend in performance along sizes observed for small and medium problem sizes is reserved for large problem sizes, as shared memory becomes a limiting factor due to its reduced capacity and more frequent high-latency global-memory accesses.