PARALiA: A Performance Aware Runtime for Auto-tuning Linear Algebra on Heterogeneous Systems

Level-3 BLAS routines operate on matrices, and typically the problem is decomposed into data-parallel tasks, following some matrix decomposition scheme. On multi-GPU setups, this is necessary to exploit the multiple devices as parallel workers, similarly to multi-core or multi-node execution [11]. Decomposition and distribution schemes are important, as they determine the required amount of communication between workers. Early multi-GPU BLAS libraries, such as the CUDA-based cuBLASXt [7], and its LAPACK-compatible wrapper NVBLAS [8], use a simple round-robin scheme to distribute 2D tile-based sub-kernels to devices, resulting in unnecessary communication. State-of-the-art libraries such as BLASX [9] and XKBLAS [10] borrow the 2D block-cyclic decomposition and distribution from distributed computing, which achieves a good homogeneous distribution of communication on a virtual 2D-grid of workers.

In this work, we also employ the state-of-practice 2D block-cyclic decomposition/distribution. Figure 2 shows an example of the 2D block-cyclic distribution used for Level-3 BLAS matrix-matrix multiplication (GEMM). The available system devices are organized in a virtual 2D grid, which should be as square as possible - for example, a 2 × 2 grid for four devices, a 4 × 2 or 2 × 4 grid for eight devices, a 3 × 3 grid for nine devices, and so on.

The 2D block-cyclic decomposition in GEMM depicted in Figure 2 results in a favorable communication pattern for the read-only tiles of matrices A and B. Every GPU requires the same number of “row” tiles and “column” tiles from each of the input matrices. This is the basis for enabling communication optimization, as communication is the main bottleneck in multi-GPU BLAS performance on modern systems [9, 10, 12]. While different libraries have used different approaches for improving communication performance, the optimization targets can be roughly classified into communication overlap, avoidance, and routing.

Overlap : Overlap refers both to computation-communication overlap, as well as communication-communication overlap, when this happens between different devices. Computation-communication overlap is a common technique in GPU offload, both in single and multi-GPU setups and it has been extensively explored in the past [12, 13, 14, 15]. CUDA versions also frequently increase the overlap potential by enabling additional GPU copy-engines [16]. In this work, we also try to maximize overlap to improve the performance of Level-3 BLAS.

GPU Data caching : In a multi-GPU setup, the different GPUs have distinct memories. Because of the problem distribution, data that is necessary for computations need to be transferred from one device to another often. Although the effect of those transfers is mitigated by technologies like RMA, a common approach to reduce redundant communication is data caching/buffering on the GPU memory, as it also enables data reuse between subsequent subkernels. For example, in Figure 2, if GPU 0 does not cache data tiles, it needs to fetch \(C_{00}\), \(A_{00}\), \(A_{01}\) and \(C_{01}\) twice, resulting in 12 tile fetches instead of 8 (50% increase in communication volume), which can worsen depending on the problem size and the tiling size T.

BLAS libraries initially designed for CPUs [3, 4, 5, 6] use simple buffers in GPU memories to support some data caching, and cuBLASXt [7] follows the same design logic. Unfortunately, simple buffering is only sufficient for internal GPU caching, which fares well on a single GPU. Its performance degrades rapidly as the number of GPUs increases, since it neglects the very fast connections between discrete GPU memories modern interconnects offer [9]. In this work, we employ a GPU data caching scheme to reduce unnecessary communication and maximize data reuse.

Communication routing : Routing in a multi-GPU setup refers to selecting the fastest routes for transfers between GPUs. In the context of Level-3 BLAS, the multiple GPUs require different data tiles, which need to be transferred, and peer-to-peer transfers between GPUs (henceforth d2d) usually have considerably higher bandwidth than CPU-GPU transfers (henceforth h2d/d2h). To enable routing optimizations for multi-GPU BLAS, two components are necessary: 1) a cache consistency-like logic for the GPU buffers, to ensure that data tiles are always up to date and their Read/Write dependencies are respected, and 2) a way to distinguish the interconnect bandwidth levels in order to select the ‘closest’ fetch location when a read-only (RONLY) tile is available in multiple buffers. A solution like this is implemented in BLASX [9], which provides a hardware abstraction of the underlying interconnect in a hierarchical representation, designed for communication optimization in multi-GPU systems. The hierarchical representation abstracts the system as a cache hierarchy, with the CPU RAM being the main memory and the distinct GPU memories being levels of caches. This representation favors data reuse (last-level cache ‘hits’) and the faster GPU-GPU transfers (lower-level cache ‘hits’) to CPU-GPU transfers (cache ‘misses’), while a MESI-like protocol enables sharing RONLY blocks between GPUs. Still, BLASX performance is hindered by continuously writing back to the CPU and re-fetching output (WR) blocks. Moreover, recent multi-GPU clusters from NVIDIA are connected with NVLink lanes, which may offer more than one distinguishable bandwidth ‘levels’, compared to the simple distinction between ‘h2d’ and ‘d2d’ bandwidths. XKBLAS [10] overcomes the write-back bottleneck by performing lazy write-backs of output data. It additionally considers a heuristic ‘ranking’ of distinguishable bandwidth levels through information from the NVIDIA driver interface, and favoring transfers over device connections higher bandwidth. In our work, we also perform communication routing as a communication optimization, through a different hardware abstraction and caching scheme.

The data placement hazard : While all the aforementioned optimizations target the full-offload scenario, where all data are initially located on the CPU, BLAS multi-GPU libraries support input data on any GPU. On modern systems, CPU-GPU (h2d/d2h) transfers are inherently much slower than peer-to-peer, GPU-GPU (d2d) transfers, therefore having part of the input data available on GPUs should lead to a less communication-bound problem.

However, as shown in Figure 1, the performance of BLASX and XKBLAS drops significantly in all but the full-offload scenarios and only PARALiA provides the expected increased performance. We demystify the cause for this counter-intuitive behavior in Figure 3, which shows the communication pattern, number of transfers and average achieved bandwidth for BLASX, XKBLAS and our proposed runtime, PARALiA. BLASX suffers from excessive costly write-backs to the CPU due to its caching policy, although these h2d/d2h could be completely avoided in the scenario where all data are initially on the GPUs [10]. Additionally, BLASX is bandwidth-agnostic, with transfers passing through a variety of bandwidth levels, since its hierarchical abstraction is only capable of recognizing the difference between h2d/d2h and d2d, and is not sufficient for modern interconnects. XKBLAS, on the other hand, avoids intermediate write-backs to the CPU due to its lazy write-back design and provides a much more balanced communication map in the full-offload scenario, with a clear preference for the higher bandwidths. This desirable behavior does not extend to the scenario where data initially populate GPUs 0, 1 and 2. On the contrary, the communication map becomes more dense around these locations and creates a communication bottleneck, with a lot of extremely low-bandwidth transfers. The cause of this meltdown lies in the core of current multi-GPU Level-3 BLAS optimization - the decomposition and distribution method. All usual distribution methods (round-robin, block-cyclic, 2D block-cyclic, etc.) target a homogeneous, balanced scenario, where data - and consequently the communication volume - is distributed between workers as equally as possible. In the full-offload scenario, all data must first be fetched from the host (\(id=8\) in the heatmaps), which in our testbed - and most modern HPC clusters - has the same bandwidth for all h2d/d2h connections and therefore favors a balanced distribution. On the other hand, in the second scenario, some GPUs are closer (higher d2d bandwidth) and some are further (lower d2d bandwidth) from the data, which results in transfers passing through a variety of connections, some of which are very slow. In our work, PARALiA, we take data placement into consideration and optimize communication, under the assumption of heterogeneity in the underlying interconnect, through an effective abstraction of the hardware, modeling, and autotuning. As shown in Figure 3, our work achieves higher average bandwidth compared to other techniques regardless of the initial data placement.

State-of-the-art multi-GPU libraries attempt to mitigate the aforementioned problem of imbalance between transfers, and also adapt to potential heterogeneous computational capabilities, by using task load-balancing through work-stealing [9, 10]. However, although work-stealing can improve load balancing, it disrupts the potential temporal locality offered by a good initial task decomposition and distribution, as is the 2D block-cyclic decomposition. Moreover, although work stealing may address the load imbalance coming from small performance differences between GPUs, it is not sufficient for heterogeneous workload distribution, as we will show in Section 4.4. Finally, the work stealing mechanisms integrated in existing approaches are performance-agnostic; they balance work between devices, without taking into account the computational capabilities of the devices or the communication properties of the problem. This design results in inefficient resource allocation and usage for BLAS execution, and is incompatible with future heterogeneous systems. This monolithic design does not lead to efficient execution since it can’t adapt resource allocation to BLAS problem characteristics and is incompatible with future heterogeneous systems.

In this work, we argue that a homogeneous distribution coupled with work-stealing is not able to effectively handle the built-in heterogeneity of modern HPC systems for the case of BLAS, and propose a model-based approach that adjusts the decomposition, communication, and task allocation to the characteristics of 1) each different system through offline micro-benchmarks and 2) each different BLAS problem and data placement during runtime. Our approach steps away from the typical homogeneous distribution coupled with work stealing, as we believe that this approach cannot simultaneously account for the interconnect layout, problem size, and data placement, and therefore cannot effectively handle the built-in heterogeneity of modern multi-GPU HPC systems for the case of BLAS.

The autotuner is the backbone of PARALiA’s optimization. Its purpose is to improve (1) communication throughput and (2) workload distribution for arbitrary system/problem configurations. Due to the more generic nature of this problem, using a heuristic-based approach is bound to favor a subset of configurations, based on which the heuristics were designed, that being either specific system characteristics (e.g., number of CPUs/GPUs, inter-connectivity) or problem characteristics (e.g., data size, placement). For this reason, the autotuner uses a model-based approach instead, which looks at each configuration as a different problem, by combining its system and problem characteristics at runtime.

The autotuning algorithm that commences during each routine invocation is shown in detail in Figure 5. When a routine is invoked, the problem parameters are extracted from the routine. The autotuner loads pre-obtained transfer coefficients from the PARALiA database and uses them to construct an abstraction of the system characteristics called LinkMap. Then, the autotuner loops over candidate workload distributions, estimates their total performance and selects the best one. Each workload distribution consists of a) a list of \(dev_{num}\) devices (\(active\_dev_{ids}\)), which is a subset of the total system devices, b) a list of sub-kernel ratios (\(active\_dev_{ratio}\)) suggested for each device, and c) a transfer routing map optimized for this specific distribution. Regarding (a) and (b), since their combined search space is very large (\(active\_dev_{ratio}\) are float values), we instead decouple them by iterating on the possible device combinations \(active\_dev_{ids}\) (which are discrete) and selecting their \(active\_dev_{ratio}\) with a model-based method. Specifically, for each device combination we start with equal sub-kernel ratios, and iteratively adjust the ratios based on a performance prediction for each device (more in Section 3.3), until a \(active\_dev_{ratio}\) with similar performance per device (within 5%) is reached. Regarding (c), the autotuner adjusts and optimizes the LinkMap to each aforementioned scenario using its specific problem characteristics (more in Section 3.2). Finally, the best workload distribution is selected by using some metric-related aggregator (e.g., maximum for time, sum for energy, etc.) on the performance of each device obtained during the estimation of (b). We note that the autotuner also selects a tiling size T for tile decomposition (as depicted in Figure 4), but this process is disconnected from (a), (b) and (c) and performed based on CoCoPeLiA [12] due to its small impact in multi-GPU performance.

Since the hardware abstractions of previous libraries target homogeneous distributions in systems with similar device and interconnect capabilities, they are not suitable for any workload distribution. To mitigate this we assume the most generic system in an abstraction called LinkMap, capable of representing any system with arbitrary devices and connections between them. The LinkMap abstraction disconnects from the notion of “CPU” and “Main memory” and treats all parts of a system similarly; any candidate data location or available computational resource is categorized as a device and is connected via links with all other devices, which are responsible for data transfers between them. In the LinkMap representation each device is defined by a unique id (\(dev_{id}\)). While not common in current systems, different devices can share memory, in which case the transfer link time between them is always equal to zero. Additionally, this abstraction assumes a fully-connected virtual topology; even if an actual hardware connection does not exist between a device pair. Therefore, this creates a fully-connected graph, with where devices are the nodes and the links are the edges: the \(dev_{num}\) nodes are connected via a 2D grid \((dev_{num},dev_{num})\) of edges/links. The LinkMap representation is implemented in C++ as a class whose members and functions are shown in Table 1. It consists of five 2D matrices \(link_{\lbrace lat,bw,bw-shared,route,sl\rbrace }\) that hold its values and three functions that are used during auto-tuning to update them.

The LinkMap representation by itself does not contain any insights, it just represents the most general case. Its usefulness is its adaptability to any system and problem data placement, which happens during runtime and has three basic phases, implemented in the LinkMap functions. First, once per program during the first routine invocation, \(load\_link\_weights()\) loads the transfer coefficients \(link_{\lbrace lat,bw,sl\rbrace }\) from the database. This provides a basic System LinkMap containing empirically obtained estimations for the system in general. Then during the autotuning of any routine, \(estimate\_problem\_throughput()\) adjusts the LinkMap bandwidths (\(link_{bw-shared}\)) according to the current device/data configuration. Specifically, it assumes that all links that connect the \(dev_{num}\) devices (\(active\_dev_{ids}\)) to the \(data_{num}\) data locations (\(data_{locs}\)) perform transfers for the entire routine execution, and apply the slowdown of simultaneous usage (Equation (6)) to the bandwidth of each such link:

The final optimization phase is to calculate something similar to the shortest paths for this graph. In our case, we want to reroute transfers that would pass through links with low bandwidth to series of links of higher bandwidth. For example using three devices, if \(link_{bw-shared}(0\rightarrow 1) = 1 Gb/s\), \(link_{bw-shared}(0\rightarrow 2) = 3 Gb/s\) and \(link_{bw-shared}(2\rightarrow 1) = 4 Gb/s\), the shortest transfer route for \(0\rightarrow 1\) would be optimized to \(0 \rightarrow 2 \rightarrow 1\), since it would be faster to transfer data from device 0 to device 1 through device 2, instead of using their direct link. To avoid very long routes the re-routing algorithm (\(optimize\_problem\_routes()\)) uses a max_hops argument that limits the intermediate data locations (currently supported \(max\_hops=\lbrace 1,2\rbrace\)), and we use \(max\_hops=1\) in our evaluation. Performing these intermediate ‘hops’ during runtime has a very low overhead since PARALiA already holds tile buffers in all devices. Re-routing significantly improves performance since 1) bandwidth is increased for the otherwise slowest transfers, which are an important bottleneck, and 2) in level-3 BLAS, transferring a read-only data chunk with additional ‘hops’ (like device 2 in the example) also stores it to these devices for potential use.

As explained in Section 3.1, the ratio adjustment and the total performance aggregation in the autotuner use an estimation of the offload performance of each device (henceforth \(pred_{metric}(dev_{id})\)). Performance prediction in multi-GPU setups is considerably more complex than on a single GPU, as scheduling on multiple devices involves runtime decisions regarding data caching and simultaneous resource utilization that are not static or known beforehand. For this reason, we use a performance upper bound based on the full-overlap model [15], instead of using more advanced overlap models [12, 14, 15]. We note that, for simplicity, all equations presented below use time as the performance metric, but PARALiA supports more performance metrics that are later explained in detail. Table 2 summarizes the modeling notation used in this work.

First, we combine the full-overlap upper bound [15] with the PARALiA database to get a routine-specific, full-overlap prediction for each device’s total performance:where h2d stands for host-to-device and d2h for device-to-host transfers, and \(\sum _{\lbrace i,j\rbrace }^{is_{\lbrace R,W\rbrace }}\) are the subsets of the \(data_{num}\) matrices/vectors that are problem inputs and outputs, respectively. To adjust the model for multi-device offload, we need to replace ‘h2d’ and ‘d2h’ time with the transfer times of all links connecting \(data_{locs}\) to each device. To do this, first, we calculate the transfer time for each link (\(t_{link}\)) as a function of transferred bytes with:

by combining each link’s latency and bandwidth using the well-accepted latency/bandwidth model [12, 14, 15, 17, 18, 19]. Then, we assume the best-case scenario, where all input matrices/vectors are distributed equally between the \(dev_{num}\) devices by combining Equation (1) with Equation (2) to generalize for any initial data placement:

Equation (3) provides a more accurate prediction for the full-overlap performance of a routine, if multi-GPU execution does not involve additional transfers/data sharing between devices. This assumption works for level-1 and level-2 BLAS, but in level-3 BLAS decomposition each tile is reused by many sub-kernels and therefore transferred to multiple devices throughout a routine’s lifetime. Since PARALiA uses a 2D block-cyclic decomposition (\(DC_{row}\), \(DC_{col}\)) for level-3 BLAS, we consider this baseline scenario of 1) exchanging equal portions of RONLY bytes between all decomposition rows and columns and 2) no output data sharing. We estimate the proportional increase in transfer volume for each device as:

Where \(RONLY_{num}\) is the number of matrix/vectors with \(is_R =1\) and \(is_W = 0\), and \(RONLY\_sum\_bytes\) is the sum of their corresponding bytes. This represents a lower bound of the added bytes due to multi-GPU BLAS3 data sharing for each device. We assume these bytes are equally distributed between devices, and use the average bandwidth of all links to estimate the additional transfer time:in which the extra communication in bytes for each device is multiplied by the inverse of its average receive bandwidth, which serves as an average estimate for the expected bandwidth of these transfers. We finally construct the full-overlap model used for the estimated performance of each GPU in a multi-GPU environment by adding the extra transfer time of Equation (4) to Equation (3):Performance metrics : All the aforementioned models return a time prediction for the execution on a single GPU. We use the maximum predicted execution time, \(total\_pred\_{t} = max(pred_t(dev_{id}))\) to evaluate different candidate workload distributions. PARALiA also supports utilization/energy-centric metrics, based on the total power consumption of all GPUs in the multi-GPU setup, \(total\_pred\_{W}\), which we combine with \(total\_pred\_{t}\) to further enhance the workload selection process. In this work, we use 1) performance (\(FLOPs=\frac{FLOP}{time}\)), 2) an inverse energy-delay product (\(EDP_i=FLOPs^2/W\)) and 3) an inverse power-delay product (\(PDP_i=FLOPs/W\)) for workload distribution (PARALiA select(\(\lbrace FLOPs,EDP_i,PDP_i\rbrace\), respectively). As evaluation metrics we use performance (FLOPs, in GFLOPs or TFLOPs) and energy efficiency (\(PDP_i\), in \(GFLOPs/W\)). To support rapid experimentation with additional metrics, we have simplified the addition and benchmarking of new metrics, requiring the modification of 3-4 lines of code only.

The PARALiA database stores the empirical measurements required to construct PARALiA’s system abstraction LinkMap, and to estimate performance in the autotuner. These measurements include transfer latencies and bandwidths, which are collected only once per system, and execution time/Watt measurements, collected for each routine. The database is automatically built by PARALiA at installation time, with a set of automated micro-benchmarks (we denote this process as the DB Builder), for all available devices and all connections between them.

Database Builder : The DB builder is an extension of the relevant component in CoCoPeLia [12], which performs single-device BLAS routine benchmarks for all system devices, and extends the set of system benchmarks to model transfers according to the requirements of the LinkMap representation. In PARALiA’s DB builder, we opt for ease of use, robustness and short benchmarking time. For ease of use, PARALiA provides a fully automated process for micro-benchmarking and for producing the empirical transfer models. Additionally, PARALiA is easily extensible to accommodate new backend BLAS library options, providing micro-benchmark template scripts, which can be easily modified with new routine invocations and any additional parameters. For robustness, for each benchmarked value, we repeat measurements until the 95% confidence interval of the mean falls within 5% of the reported mean value (taking at least 10 measurements to obtain these means). Finally, for short benchmarking times, we try to strike a balance between the number of measurements required for robustness, and their overall execution time. The DB builder benchmarks the computation time of the different backend BLAS routines. The minimum problem dimensions and steps are static and predefined for all benchmarks, while the maximum dimensions are calculated based on the available memory of the target device. The results of the DB builder are stored as a database and made available to the framework for all subsequent calls in the system.

Kernel lookup : The offload performance models used in the autotuner require an estimate for the routine execution time and average Watts per device. Using the same technique as in CoCoPeLia [12], we only collect measurements for the time/power of fine-grained chunks of specific, small tiling sizes, namely \(\lbrace t,W\rbrace _{exec}(routine, dev_{id}, T[,T[,T]])\), for which we then use value lookup in the database. The average GPU Watt consumption is obtained by sampling Watt values at regular intervals with the CUDA nvml-driver during each routine benchmark and averaging these. Micro-benchmarks are performed per routine and per device (\(dev_{num}\times routine_{num}\) times) and use separate BLAS backends depending on the target \(dev_{id}\). Devices with \(0 \le dev_{id} \lt cuda\_dev\_num\) are reserved for available CUDA devices and devices with \(dev\_{id} \gt =cuda\_dev\_num\) are reserved for available CPUs (usually one). The value lookup micro-benchmarks use cuBLAS for NVIDIA GPUs and OpenBLAS for CPUs, but are easily extensible with minimal adjustments to other \(dev_{id}\) ranges (e.g., for AMD devices) or for different BLAS implementations (e.g., a custom GPU implementation instead of cuBLAS) that follow the BLAS standard.

Transfer coefficients: To obtain \(link_{\lbrace lat,bw\rbrace }\), we follow the most widely used semi-empirical approach; we measure a set of transfer times and use them to fit the coefficients of basic linear models for transfer time. We obtain \(link_{lat}\) empirically as the average latency of multiple single-byte transfers. For \(link_{bw}\), we run benchmarks for square transfers with \(dtype=double\), for \(D1=D2=256 \xrightarrow []{step_{ad}=256} \sqrt {max\_device\_memory/2}\), and use least square regressions on the obtained samples in the manner of [20]. Then, we estimate the slowdown \(link_{sl}\) for simultaneous link usage (e.g., transfer overlap for any two links), assuming it imposes a constant throughput slowdown and does not affect latency. This slowdown is calculated with a single micro-benchmark for each link; first, for the link of interest, a large transfer (\(D1=D2=\sqrt {max\_device\_memory/2}\)) is tested isolated (\(t_{link1}\)), and then, it is tested overlapped (\(t_{link1-link2_{over}}\)) with multiple similar transfers on the other link, resulting in the slowdown:

Since the method of obtaining the sl is empirical, we assume a maximum slowdown of \(sl_{link1-link2} = 2.0\) (i.e., the effective bandwidth is halved), to avoid empirical errors spilling into the models. For all transfer experiments we use the PARALiA wrapped functions for transfers, which currently use the cudaMemcpy2DAsync routine in their back-end with pinned host memory, as required by these asynchronous calls.

The PARALiA preprocessor is responsible for the framework initialization and the transformation of problem data for BLAS execution in a multi-GPU system, which is broken down into the three basic operations described next.

Environment setup: This operation is performed when a BLAS routine is invoked for the first time or with a new set of parameters. It allocates buffers, initializes data structures, and performs backend-specific actions, like creating CUDA streams and events, to be used by the scheduler.

Tile decomposition: The bulk of preprocessing in BLAS libraries involves decomposing the problem data into smaller chunks, usually referred to as tiles. As most similar multi-GPU libraries [6, 7, 9, 10], in PARALiA we decompose vectors to 1D tiles and matrices to 2D square tiles, using a tiling size T provided by the autotuner. After data decomposition into tiles, PARALiA identifies all sub-kernels deriving from this decomposition, generating the relevant data/task dependencies.

Device initialization: This operation initializes the devices that will participate in BLAS execution (\(active\_dev_{ids}\)) and distributes sub-kernels to them, proportionally to their assigned problem ratios (\(active\_dev_{ratio}\)). PARALiA supports multiple sub-kernel distribution patterns (Sequential, Round-Robin, 1D-cyclic and 2D-cyclic) and by default uses Round-Robin for BLAS 1 & 2 and 2D-cyclic for BLAS3. Unlike scheduler-centric multi-GPU libraries [9, 10] that rely on dynamic load-balancing, PARALiA follows a static approach since load-balancing is based on effective performance estimation performed before scheduling.

The PARALiA scheduler manages data caching in distinct device memories, and data transfers between memories, invokes all sub-kernels on their assigned devices, and synchronizes their execution and results, as analyzed below.

Data caching: For this operation, PARALiA uses a Software_buffer C++ class, similar to BLASX and XKBLAS, which represents a buffer in each device with a distinct memory, and can store 1D and 2D tiles. Each Software_buffer holds a number of blocks depending on the problem size and per-device memory limitations, and employs a simple block-sequential write-back policy to swap tiles during sub-kernel execution, using a MESI-like protocol similar to BLASX [9]. This Software_buffer in each device is initialized the first time a program calls a PARALiA BLAS routine, and is updated/extended in subsequent calls to match their device and size requirements.

Sub-kernel invocation: The scheduler spawns \(active\_dev_{num}\) POSIX threads, that are responsible for invoking sub-kernels in their corresponding devices.

After a sub-kernel is invoked, it issues three sub-tasks: a) the requests for all its input tile dependencies (e.g., fetching tiles if they are not available in its device memory), b) the sub-kernel computations to be performed after dependencies are met, and c) potential data write-backs to the initial memory location for any tile it modified. The optimization of sub-kernel invocation order is important for multi-GPU BLAS scheduling, since it affects both task parallelism and the communication and data reuse pattern [9, 10]. We leave the sub-kernel order problem for future work because PARALiA focuses on model-assisted communication and workload distribution, not dynamic scheduling techniques.

Synchronization: The sub-tasks of each sub-kernel (i.e., fetch data, compute, writeback) are executed asynchronously and overlapped with sub-tasks from other sub-kernels (software pipelining) and on other devices (multi-GPU) using CUDA events to enforce data dependencies and CUDA streams to enable overlap. After all sub-kernels are invoked, the scheduler synchronizes all sub-tasks and returns the result and control to the user upon completion.

For the performance evaluation we use a single testbed: the “clx-ai” nodes of HLRS’ HPC cluster VULCAN, HPC cluster at HLRS [21]. System details are presented in Table 3, along with the interconnect bandwidths stored in the LinkMap for the 9 devices (8 GPUs + CPU which are illustrated in Figure 6). The interconnect utilizes a mix of NVlink-1 (24 GB/s) and NVlink-2 (48 GB/s) for inter-GPU connectivity and PCiE (12 GB/s) for all CPU-GPU communication. In addition, we note that CPUs share PCIe bandwidth in sets of two (e.g., GPU 0-1, 2-3, etc.). For time measurements we use wrapped timers based on clock_gettime, with device synchronization (cudaDevice Synchronize()) also included; both timer and synchronization overhead were less than 1% for all benchmarks. We perform 20 executions for large benchmarks and 100 for small ones, after 10 warm-up runs, and we report the median time/performance of these runs. The allocation time needed for CPU/GPU buffers is not modeled or included in the total time, and all matrices/vectors are initialized with random values before execution. We use pinned host memory to enable asynchronous CUDA calls and the caches/buffers are flushed between runs. The above configuration is consistent for all our experiments and all state-of-the-art libraries we include in this work.

Device selection in homogeneous systems is meaningful for communication-bound problems but can be even more beneficial in heterogeneous systems, where devices can have different computation capabilities. While heterogeneous multi-device systems are not common nowadays, computational heterogeneity will probably be more commonplace in the future. Heterogeneous-like execution scenarios can also appear in current homogeneous multi-device systems, by applying different power management policies or sharing devices between users/processes. For this reason, we include a proof-of-concept application of our approach to an artificial heterogeneous system, which we emulate by loading the GPUs of the Vulcan clx-ai cluster with different computation workloads running in other processes. We configure these workloads empirically to represent GPUs with lower performance, resulting in the following double FP peak adjustments: \(GPU_{\lbrace 0,1,4,6\rbrace } = 3.5\) Tflops, \(GPU_{\lbrace 2,3\rbrace } = 5\) Tflops and \(GPU_{\lbrace 5,7\rbrace } = 7\) Tflops (original peak). We also do not adjust the power consumption of each device, resulting in different energy efficiency for each device category (e.g., \(GPU_{\lbrace 5,7\rbrace }\) are more energy efficient than \(GPU_{\lbrace 2,3\rbrace ,}\) etc.). This leads to a total system DP peak of 38 Tflops (vs 56 Tflops for the original system), and a \(PDP_i\) peak of 17.5 (vs 26 in the original system). We also note that a homogeneous-distribution DP peak (without load-balancing) is \(3.5 \cdot 8 = 28\) Tflops for future reference. Additionally, we limit the dataset to less than half the problems, by doubling the data size iteration steps and the minimum size for a total of 250 problems, and exclude cublasXT due to extreme benchmark times and having no load-balancing mechanism whatsoever. The results for this heterogeneous-emulated system are shown in Table 5 and Figure 9.

Figure 9 contains the aggregated performance (left) and energy efficiency-\(PDP_i\) (right) box-plots for the entirety of the aforementioned heterogeneous dataset of 250 problems (HS - 90, HL - 160), which now leans more towards computation-bound problems since the peak performance has lowered considerably while the interconnect is the same. BLASX, XKBLAS and PARALiA comm_opt follow similar performance and \(PDP_i\) patterns with the homogeneous system, albeit at lower peak as expected. Communication optimization is still an issue for BLASX and XKBLAS resulting in superior performance for PARALiA comm_opt, but all three approaches are limited by the aforementioned homogeneous-distribution peak around 28 Tflops. On the other hand, PARALiA select(\(EDP_i\)) manages to provide a better workload distribution, which results in better performance by approaching the peak of 38 Tflops for large problems, and has far superior \(PDP_i\) both due to better performance and due to awareness of the different characteristics of each emulated device. Summing up, PARALiA select(\(EDP_i\)) vastly outperforms previous approaches both in terms of performance and energy efficiency in the heterogeneous system as well, now further boosted by a better workload distribution in different devices.

Table 5 shows the achieved mean performance in a similar layout with Section 4.3 for the small (HS), large (HL) and total (HT) dataset also displaying results for PARALiA select(PERF). PARALiA select(\(PDP_i\)) is omitted due to always selecting from \(GPU_{\lbrace 5,7\rbrace }\) as expected without adding any additional insights. Since the performance is already visualized in Figure 9 we use Table 5 for problem size-related insights and mean comparison, making the following observations. First, PARALiA versions still outperform previous approaches in all subsets, but with a smaller performance difference added to the baseline PARALiA comm_opt from communication optimization. Specifically, PARALiA select(PERF) has a (geo)mean performance improvement in (S, L, T) of (2.2×, 1.7×, 1.8×) over BLASX and (2.0×, 1.9×, 1.9×) over XKBLAS, and PARALiA select(\(EDP_i\)) has (1.8×, 1.6×, 1.6×) over BLASX and (2.0×, 1.7×, 1.7×) over XKBLAS. This is expected, since by reducing peak performance and removing smaller problems from the total dataset, the impact of communication optimization is limited since many problems now become compute-bound. Second, the impact of workload selection increases performance and energy efficiency both for PARALiA select(PERF) and PARALiA select(\(EDP_i\)) in respect to PARALiA comm_opt, since devices with lower computational power are used for a smaller part of the problem or omitted by the auto-tuning runtime. This is more visible in the large (HL) problems which are compute-bound, since in small (HS) problems the communication optimization is still more important. The difference between PARALiA select(PERF) and PARALiA select(\(EDP_i\)) becomes more evident when comparing them with PARALiA comm_opt; PARALiA select(PERF) offers 1.12× mean performance and 1.3× mean \(PDP_i\) improvement while PARALiA select(\(EDP_i\)) offers 1.06× performance and 1.5× \(PDP_i\). Summing up, both PARALiA select(PERF) and PARALiA select(\(EDP_i\)) benefit from workload selection, offering different insights and balance between performance and energy in the heterogeneous system, outlining the increased importance of additional metrics for such future systems.

Since BLAS is an important part of a plethora of scientific code and solvers, from the early days of computing it played a vital role in scientific code optimization. Dongarra [1] firstly defined the BLAS standard as a set of “black box” routines that should follow a specific input/output layout and should be optimized by vendors and library providers transparently to the user, without any additional performance tuning from the side of the user. This black-box approach required performance engineering effort to reimplement or tune BLAS routines for every new generation of hardware. To offset this, Whaley et al. [11] implemented ATLAS, the most well-known BLAS autotuning framework, which uses their automated empirical optimization of software (AEOS) technique. ATLAS explores many possible routine implementations, tuning itself to each new system empirically by testing them and timing them, in order to improve cache utilization, increase parallelism and load balance sub-problems. In a similar notion, Van Zee and Vsan de Geijn [23] defined a set of a few kernels that can be used to optimize all BLAS 2 and BLAS 3 operations using empirical knowledge, and later Low et al. [24] proved that an analytic approach is sufficient for replacing the empirical part with auto-tuning. These techniques are not used directly in our work, but they form the groundwork for most BLAS-related autotuning and model-based optimization.

Regarding GPU BLAS implementations, cuBLAS [2] was the first library to provide high performance on a single GPU, expecting, however, the data to be available in the GPU. cuBLAS autotunes the CUDA kernel block sizes internally, based on pre-performed offline empirical testing for each NVIDIA GPU. When the data reside outside the GPU memory, the execution of a BLAS routine requires data transfers as well. Gregg and Hazelwood [25] were the first to highlight the problem of data transfer overheads, arguing against the trend to exclude transfer overheads in the scientific reporting of the performance of GPU applications, and proposed a taxonomy for data transfers and their impact on offload performance. Numerous works model CPU-GPU transfers, using variants of the linear latency-bandwidth model for PCIe transfers [17, 26, 27, 28] and machine learning [27] to provide user insights, without, however, specifically targeting BLAS, or integrating them in an autotuner. The inclusion of transfers in GPU performance prediction improved accuracy, but whenever there was communication/computation overlap, simplistic models still failed to predict the actual performance. As highlighted by Hoefler et al. [13], modeling overlap areas is a crucial step in performance modeling. Towards this direction, Gómez-Luna et al. [14] explored the use of CUDA streams for 3-way concurrency, but they consider the stream creation time as the only overlap overhead. Werkhoven et al. [15] enhanced this work by offering multiple performance models for communication/computation overlap for various common offload scenarios (RMA, 2-way, 3-way), introduced stream transfer overlap latency, and provided methods to obtain the optimal number of CUDA streams for a given problem. In a similar notion, Liu et al. [19] offered a mathematical framework for software pipelining on GPUs using non-equal tiles, which focused on partitioning, scheduling, and granularity. All these models offer high accuracy, however, their modeling approach does not capture all problem characteristics present in BLAS. Moreover, they were never used in practice for autotuning, and therefore never faced the practical problems of transfer prediction in real systems, where empirical bandwidths defer from their theoretical counterparts depending on many parameters like the type of src/dest memory (normal, unified, pinned) or the underlying sharing of interconnect resources [29, 30, 31]. In our previous work CoCoPeLia [12], which aimed to create an autotuner that would selected good tiling sizes for domain decomposition in single-GPU BLAS, we had to overcome both of these problems to enable performance estimation. To overcome the first problem, we modified these models to account for common BLAS characteristics, like data reuse and bidirectional transfers. To overcome the second problem, we created an automated pipeline which fueled these models with offline empirical tests and integrated these into an end-to-end BLAS framework, which provides state-of-the-art performance for single-GPU BLAS.

Several research works focus on CPU-GPU hybrid BLAS execution, a topic relevant to our work, because of the problem of heterogeneous splitting, namely the splitting of a BLAS problem into subproblems, based on performance estimations for the execution on the CPU and the GPU, integrated in an end-to-end solution. Luk et al. [32] first proposed such an end-to-end solution for model-based CPU-GPU splitting for a subset of BLAS routines. Their solution employs an initial benchmark run, through which their backend BLAS functions obtain linear models for CPU and GPU performance, store them in a database (file) and utilize them in subsequent runs, to split a problem depending on its predicted CPU and GPU performance. Tomov et al. [4] envisioned a heterogeneous LAPACK and made use of BLAS-level parallelism, where the program is represented as a collection of BLAS-based tasks and dependencies. To this end, they use task graphs to represent BLAS task dependencies and use the CPU for small kernels, the execution of which is inefficient on GPUs, and the GPU for everything else, overlapping CPU and GPU work and transfers. Humphrey et al. [33, 34, 35] developed CULA, a CUDA framework which enables GPU-CPU simultaneous execution, with each one running predefined routine parts most fit to its paradigm. Bernabé et al. [18] employed CPU-GPU hybridization auto-tuning based on micro-benchmarks run at the time of the library installation, which were used to feed polynomial model coefficients. Finally, Ma et al. [36] explored energy efficiency, splitting the problem during execution until a balance between CPU-GPU execution time is reached, and then scaling frequencies to limit energy consumption. All the above approaches do not take into account transfers and overlap in their modeling/assumptions, favoring the GPU even in cases that the problem is transfer bound if its computational capability is larger than the CPU, and do not target multi-GPU systems. This leaves a wide gap between multi-GPU libraries and hybrid/heterogeneous BLAS research. In our work, PARALiA, we bridge this gap with high-quality subproblem scheduling and communication optimization on multiple GPUs, together with performance awareness and device selection, to ensure resource utilization for any type of heterogeneity.