<figure><img src="image_1.png" /></figure>

HAL
open science

# Latency hiding strategies for matrix multiplication on GPUs using pipelined multiword arithmetic

Alexis Aune, Oguz Kaya

<figure><img src="image_2.png" /></figure>

## To cite this version:

Alexis Aune, Oguz Kaya. Latency hiding strategies for matrix multiplication on GPUs using pipelined multiword arithmetic. 2024. hal-04619183

HAL Id: hal-04619183

https://hal.science/hal-04619183

Preprint submitted on 20 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
<figure><img src="image_1.png" /></figure>

HAL
open science

# Latency hiding strategies for matrix multiplication on GPUs using pipelined multiword arithmetic

Alexis Aune, Oguz Kaya

<figure><img src="image_2.png" /></figure>

## To cite this version:

Alexis Aune, Oguz Kaya. Latency hiding strategies for matrix multiplication on GPUs using pipelined multiword arithmetic. 2024. hal-04619183

HAL Id: hal-04619183

https://hal.science/hal-04619183v1

Preprint submitted on 20 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Latency hiding strategies for matrix multiplication on GPUs using pipelined multiword arithmetic

Alexis Aune and Oguz Kaya

Université Paris-Saclay

{oguz.kaya, alexis.aune}@universite-paris-saclay.fr

**Abstract.** Matrix multiplication is at the heart of many applications in AI, scientific computing, and quantum simulations, for which GPUs have been widely adopted due to the immense compute power and hardware acceleration they provide for such computations. In this work, we present a new technique for speeding up dense matrix-matrix multiplications (GEMM) on GPUs by decomposing the operation into several cascading lower-precision GEMMs and exploiting then-available overlap between data transfers and mixed-precision compute. The proposed implementation exploits the high performance of low-precision matrix multiply-accumulate (MMA) operations on modern GPUs and their multi-stream execution capabilities to overlap multiple low-precision GEMMs with data transfers. Our method respectively achieves a speedup of up to 39% compared to a standard transfer-compute-transfer and up to 16% compared to a 2 × 2 tiled GEMM implementation using tensor cores and leveraging stream parallelism to overlap transfers with the computation.

**Keywords:** GPU computing · matrix multiplication · numerical linear algebra · tensor cores · mixed-precision computing · pipelining

## 1 Introduction

The multiplication of two dense matrices $C = AB$ ($C \in \mathbb{R}^{m \times n}, A \in \mathbb{R}^{m \times k}, B \in \mathbb{R}^{k \times n}$) is a fundamental operation at the core of numerous applications from various disciplines ranging from the training and inference of deep neural networks [4] and quantum circuit simulations [5] to numerical solvers [13] and discrete Fourier transforms [14]. For this reason, there have been endless efforts both on the software and the hardware side to provide optimized GEMM kernels, including recently introduced MMA hardware in GPUs and other accelerators [7,10]. While developing such specific hardware as well as tailored compute kernels [9] are indispensable, providing new ways to achieve better GEMM performance in specific scenarios is also crucial to improve the performance of the applications encountering such cases.

In this work, we focus particularly on the scenario where the input matrices initially reside on the CPU main memory, and the result of the GEMM is to be potentially transferred back from the GPU memory to the CPU main memory
# Latency hiding strategies for matrix multiplication on GPUs using pipelined multiword arithmetic

Alexis Aune and Oguz Kaya

Université Paris-Saclay

{oguz.kaya, alexis.aune}@universite-paris-saclay.fr

**Abstract.** Matrix multiplication is at the heart of many applications in AI, scientific computing, and quantum simulations, for which GPUs have been widely adopted due to the immense compute power and hardware acceleration they provide for such computations. In this work, we present a new technique for speeding up dense matrix-matrix multiplications (GEMM) on GPUs by decomposing the operation into several cascading lower-precision GEMMs and exploiting then-available overlap between data transfers and mixed-precision compute. The proposed implementation exploits the high performance of low-precision matrix multiply-accumulate (MMA) operations on modern GPUs and their multi-stream execution capabilities to overlap multiple low-precision GEMMs with data transfers. Our method respectively achieves a speedup of up to 39% compared to a standard transfer-compute-transfer and up to 16% compared to a 2 × 2 tiled GEMM implementation using tensor cores and leveraging stream parallelism to overlap transfers with the computation.

**Keywords:** GPU computing · matrix multiplication · numerical linear algebra · tensor cores · mixed-precision computing · pipelining

## 1 Introduction

The multiplication of two dense matrices $C = AB$ ($C \in \mathbb{R}^{m \times n}, A \in \mathbb{R}^{m \times k}, B \in \mathbb{R}^{k \times n}$) is a fundamental operation at the core of numerous applications from various disciplines ranging from the training and inference of deep neural networks [4] and quantum circuit simulations [5] to numerical solvers [13] and discrete Fourier transforms [14]. For this reason, there have been endless efforts both on the software and the hardware side to provide optimized GEMM kernels, including recently introduced MMA hardware in GPUs and other accelerators [7,10]. While developing such specific hardware as well as tailored compute kernels [9] are indispensable, providing new ways to achieve better GEMM performance in specific scenarios is also crucial to improve the performance of the applications encountering such cases.

In this work, we focus particularly on the scenario where the input matrices initially reside on the CPU main memory, and the result of the GEMM is to be potentially transferred back from the GPU memory to the CPU main memory
<table><caption>**Table 1.** Peak performance of dense GEMM operations on Nvidia GPUs in TFlop/s. Bold values represent a tensor core operation whereas regular values denote the SIMT performance.</caption><thead><tr><th>GPU</th><th>FP16</th><th>TF32</th><th>FP32</th><th>FP64</th><th>FP ratio</th><th>Bandwidth (GiB/s)</th></tr></thead><tbody><tr><td>V100 (SXM)</td><td><b>125</b></td><td>-</td><td>15.7</td><td>7.8</td><td>1:8:16</td><td>897</td></tr><tr><td>A100 (SXM)</td><td><b>312</b></td><td><b>156</b></td><td>19.5</td><td><b>19.5</b></td><td>1:6:16</td><td>1550</td></tr><tr><td>H100 (SXM)</td><td><b>990</b></td><td><b>330</b></td><td>67</td><td><b>67</b></td><td>1:6:30:30</td><td>2039</td></tr></tbody></table>

as soon as it is available (e.g., in order to continue the work on the CPU, free the limited GPU memory, or communicate the result to other nodes in a distributed memory kernel, etc.) after executing the GEMM on the GPU. We introduce a new method, called multiword pipelined mixed-precision GEMM (PPGEMM), to accelerate a such scenario by overlapping host-to-device (H2D) and device-to-host (D2H) transfers with the GEMM kernel execution on the GPU in a particular way. In this regard, the PPGEMM could be seen as an alternative to a tiled GEMM using stream parallelism to carry out H2D transfers of small input tiles to the GPU rapidly, thereby enabling an early initiation of the GEMM kernel execution as well as the D2H transfer of output tiles. However, a potential drawback of this strategy, which we especially aim to avoid in PPGEMM, is the smaller matrix/tile size, reducing the arithmetic intensity of each GEMM execution. This can potentially slow down the overall computation unless the matrix dimensions are sufficiently large to saturate modern GPUs to achieve peak performance [15], which becomes increasingly arduous in many applications due to the dramatic increase in the compute power of modern GPUs with the MMA hardware, as presented in Table 1.

Our method is based on a recently introduced *matrix multiword arithmetic* scheme [1] which expresses a high-precision matrix (e.g., $A$) as a sum of multiple low-precision matrices (e.g., $A_0 + A_1$), and computes the multiplication of two high-precision matrices through multiple low-precision-multiply and high-precision-accumulate operations (e.g., $C = A_0B_0 + A_1B_0 + A_0B_1$). In this paper, we aim to leverage the increased data transfer and compute granularity provided by such splitting with device-optimized GPU GEMM kernels that can effectively overlap computation and data transfers and analyze under which conditions such a scheme could yield better performances compared to a standard or tiled GEMM.

## 2 Related Work

Using multiple lower-precision GEMMs to achieve higher-precision is a strategy employed in many mixed-precision kernels [14,2,1,12], including precision levels that are not supported by the hardware [12]. The recent work [12] describes methods for performing a high-precision GEMM by performing multiple low-precision GEMMs, including double-double precision GEMMs (FP64x2) using ten FP64
<table><caption>**Table 1.** Peak performance of dense GEMM operations on Nvidia GPUs in TFlop/s. Bold values represent a tensor core operation whereas regular values denote the SIMT performance.</caption><thead><tr><th>GPU</th><th>FP16</th><th>TF32</th><th>FP32</th><th>FP64</th><th>FP ratio</th><th>Bandwidth (GiB/s)</th></tr></thead><tbody><tr><td>V100 (SXM)</td><td><b>125</b></td><td>-</td><td>15.7</td><td>7.8</td><td>1:8:16</td><td>897</td></tr><tr><td>A100 (SXM)</td><td><b>312</b></td><td><b>156</b></td><td>19.5</td><td><b>19.5</b></td><td>1:6:16</td><td>1550</td></tr><tr><td>H100 (SXM)</td><td><b>990</b></td><td><b>330</b></td><td>67</td><td><b>67</b></td><td>1:6:30:30</td><td>2039</td></tr></tbody></table>

as soon as it is available (e.g., in order to continue the work on the CPU, free the limited GPU memory, or communicate the result to other nodes in a distributed memory kernel, etc.) after executing the GEMM on the GPU. We introduce a new method, called multiword pipelined mixed-precision GEMM (PPGEMM), to accelerate a such scenario by overlapping host-to-device (H2D) and device-to-host (D2H) transfers with the GEMM kernel execution on the GPU in a particular way. In this regard, the PPGEMM could be seen as an alternative to a tiled GEMM using stream parallelism to carry out H2D transfers of small input tiles to the GPU rapidly, thereby enabling an early initiation of the GEMM kernel execution as well as the D2H transfer of output tiles. However, a potential drawback of this strategy, which we especially aim to avoid in PPGEMM, is the smaller matrix/tile size, reducing the arithmetic intensity of each GEMM execution. This can potentially slow down the overall computation unless the matrix dimensions are sufficiently large to saturate modern GPUs to achieve peak performance [15], which becomes increasingly arduous in many applications due to the dramatic increase in the compute power of modern GPUs with the MMA hardware, as presented in Table 1.

Our method is based on a recently introduced *matrix multiword arithmetic* scheme [1] which expresses a high-precision matrix (e.g., $A$) as a sum of multiple low-precision matrices (e.g., $A_0 + A_1$), and computes the multiplication of two high-precision matrices through multiple low-precision-multiply and high-precision-accumulate operations (e.g., $C = A_0B_0 + A_1B_0 + A_0B_1$). In this paper, we aim to leverage the increased data transfer and compute granularity provided by such splitting with device-optimized GPU GEMM kernels that can effectively overlap computation and data transfers and analyze under which conditions such a scheme could yield better performances compared to a standard or tiled GEMM.

## 2 Related Work

Using multiple lower-precision GEMMs to achieve higher-precision is a strategy employed in many mixed-precision kernels [14,2,1,12], including precision levels that are not supported by the hardware [12]. The recent work [12] describes methods for performing a high-precision GEMM by performing multiple low-precision GEMMs, including double-double precision GEMMs (FP64x2) using ten FP64
GEMMs, and provides a corresponding error analysis. In [1], the authors similarly investigate performing a high-precision GEMM by splitting the matrices into multiple low-precision words, which we later describe in Algorithm 1, then carrying out multiple mixed-precision GEMMs to obtain the result in high precision. They also provide a comprehensive error analysis, specifically with respect to the rounding algorithms used in Nvidia's FP16-multiply FP32-accumulate MMA units.

In this paper, we base our kernels on these same principles investigated in [12,1], but we instead explore the potential of this framework in the context of data transfer latency hiding and computation/communication overlap for a CPU-GPU GEMM execution. This involves optimized matrix splitting kernels on the CPU and GPU, overlapped CPU-GPU data transfers and mixed-precision GEMMs on the GPU leveraging the stream parallelism and tensor cores.

# 3 Multiword pipelined mixed-precision GEMM

The main objective of PPGEMM is to reduce the CPU-GPU transfer latency by increasing the data transfer granularity at the expense of potentially more operations performed in low-precision. This is in contrast to the tile parallelism that can similarly hide this latency, yet it sacrifices the arithmetic intensity instead due to multiple GEMMs performed on smaller tiles. As the low-precision compute is getting progressively less expensive on modern GPU architectures, as seen in Table 1, this low-precision compute overhead is expected to be negligible and overlapped with finer-grained data transfers in PPGEMM. This granularity is ensured by splitting the input matrices into a sum of multiple matrices of the same dimension but in lower-precision (e.g., $A = A_0 + A_1$ and $B = B_0 + B_1$) - hence of smaller memory size - which are called *words*. High-magnitude words ($A_0$ and $B_0$) can then be immediately transferred to the GPU. As a high-precision GEMM is performed upon their arrival ($A_0B_0$), lower-magnitude words can be transferred to the GPU ($A_1$), and the result can be updated with low-precision GEMMs ($A_1B_0$), which can typically be performed significantly faster than a high-precision GEMM on modern GPUs. This way, no sacrifice is made on the arithmetic intensity of the high-precision GEMM, and the small computational overhead caused by the extra low-precision GEMMs can be hidden by overlapping them with the subsequent word transfers ($B_1$) when computing the final result ($C = A_0B_0 + A_1B_0 + A_0B_1$).

The following subsection describes a general algorithm for splitting matrices at different precision levels. Next, we detail our PPGEMM kernel leveraging this matrix splitting strategy to achieve higher performance by better hiding the transfer latency.

## 3.1 Multiword splitting of a matrix

We define *splitting* as separating a high-precision floating-point scalar into multiple lower-precision ones whose sum yields the original value or a close approximation thereof. By extension, we define splitting a high-precision matrix as
GEMMs, and provides a corresponding error analysis. In [1], the authors similarly investigate performing a high-precision GEMM by splitting the matrices into multiple low-precision words, which we later describe in Algorithm 1, then carrying out multiple mixed-precision GEMMs to obtain the result in high precision. They also provide a comprehensive error analysis, specifically with respect to the rounding algorithms used in Nvidia's FP16-multiply FP32-accumulate MMA units.

In this paper, we base our kernels on these same principles investigated in [12,1], but we instead explore the potential of this framework in the context of data transfer latency hiding and computation/communication overlap for a CPU-GPU GEMM execution. This involves optimized matrix splitting kernels on the CPU and GPU, overlapped CPU-GPU data transfers and mixed-precision GEMMs on the GPU leveraging the stream parallelism and tensor cores.

# 3 Multiword pipelined mixed-precision GEMM

The main objective of PPGEMM is to reduce the CPU-GPU transfer latency by increasing the data transfer granularity at the expense of potentially more operations performed in low-precision. This is in contrast to the tile parallelism that can similarly hide this latency, yet it sacrifices the arithmetic intensity instead due to multiple GEMMs performed on smaller tiles. As the low-precision compute is getting progressively less expensive on modern GPU architectures, as seen in Table 1, this low-precision compute overhead is expected to be negligible and overlapped with finer-grained data transfers in PPGEMM. This granularity is ensured by splitting the input matrices into a sum of multiple matrices of the same dimension but in lower-precision (e.g., $A = A_0 + A_1$ and $B = B_0 + B_1$) - hence of smaller memory size - which are called *words*. High-magnitude words ($A_0$ and $B_0$) can then be immediately transferred to the GPU. As a high-precision GEMM is performed upon their arrival ($A_0B_0$), lower-magnitude words can be transferred to the GPU ($A_1$), and the result can be updated with low-precision GEMMs ($A_1B_0$), which can typically be performed significantly faster than a high-precision GEMM on modern GPUs. This way, no sacrifice is made on the arithmetic intensity of the high-precision GEMM, and the small computational overhead caused by the extra low-precision GEMMs can be hidden by overlapping them with the subsequent word transfers ($B_1$) when computing the final result ($C = A_0B_0 + A_1B_0 + A_0B_1$).

The following subsection describes a general algorithm for splitting matrices at different precision levels. Next, we detail our PPGEMM kernel leveraging this matrix splitting strategy to achieve higher performance by better hiding the transfer latency.

## 3.1 Multiword splitting of a matrix

We define *splitting* as separating a high-precision floating-point scalar into multiple lower-precision ones whose sum yields the original value or a close approximation thereof. By extension, we define splitting a high-precision matrix as
independently splitting every element composing it. Splitting a matrix results in multiple lower-precision matrices called *words* [1]. One can then use these words to perform more granular GEMMs at the word level, thus enabling the use of high-performance, lower-precision MMA hardware to achieve high precision.

The process of splitting a matrix is described in the SPLIT function provided in Algorithm 1, which produces a set of *n* lower-precision words $A_k$ and scaling factors $s_k$ ($0 \le k \le n$) from the input matrix A. Scaling factors are used to prevent the loss of precision caused by an overflow or underflow in low-precision words during the splitting process or subsequent multiword operations. While some algorithms exist to find effective scaling factors [3], we instead opted for a more straightforward approach by using the max norm of the pre-cast matrix at each step of the splitting process. A shifting factor $s_A$ can also be provided so that the algorithm splits $A/s_A$ instead, and reinjects this factor in every $s_k$. This can further assist in reducing the risk of overflows. Finally, for the PPGEMM kernel presented in this work, we only consider splitting an FP32 input matrix into two FP16.

<caption>**Algorithm 1:** SPLIT: Split a matrix into low-precision words</caption>

**Input:** A: Input matrix of precision $fl_A$
$fl_{0:n-1}$: Precision of each word $A_k$ for $0 \le k < n$
$s_A$: Shifting factor for A.
**Output:** $A_{0:n-1}, s_{0:n-1}|A \approx \sum_{k=0}^{n-1} s_k fl_A(A_k)$
1. $\{s_0, A_0\} \leftarrow fl_0(A/s_A)$
2. **for** $k \leftarrow 1$ **to** $n-1$ **do**
3. &nbsp;&nbsp;&nbsp;&nbsp;$$ \{s_k, A_k\} \leftarrow fl_k(A/s_A - \sum_{j=0}^{k-1} (fl_A(A_j)s_j)) $$
4. &nbsp;&nbsp;&nbsp;&nbsp;$$ s_k \leftarrow s_A s_k $$
5. **return** $\{A_{0:n-1}, s_{0:n-1}\}$

## 3.2 PPGEMM kernel

In Algorithm 2, we present our PPGEMM kernel, named 2f16h, that leverages multiword splitting of a matrix for latency hiding. We use the **upload** and **download** keywords to indicate an asynchronous data transfer from the CPU to the GPU or from the GPU to the CPU, respectively. An **on stream** block denotes the code executed on the GPU on a specific CUDA stream. All code outside these blocks runs on the CPU. Each **on stream** block is divided into three columns. The first column contains the events the stream awaits to continue its execution. The second column contains the code to be executed. The last column indicates events fired after the code is executed. The **sync** keyword expresses a synchronization barrier, requiring all previously launched GPU kernels to terminate before continuing the execution. We use HS (half-single) to denote an FP16-multiply FP32-accumulate matrix multiplication, for which we call the cuBLAS implementation [8].
independently splitting every element composing it. Splitting a matrix results in multiple lower-precision matrices called *words* [1]. One can then use these words to perform more granular GEMMs at the word level, thus enabling the use of high-performance, lower-precision MMA hardware to achieve high precision.

The process of splitting a matrix is described in the SPLIT function provided in Algorithm 1, which produces a set of *n* lower-precision words $A_k$ and scaling factors $s_k$ ($0 \le k \le n$) from the input matrix A. Scaling factors are used to prevent the loss of precision caused by an overflow or underflow in low-precision words during the splitting process or subsequent multiword operations. While some algorithms exist to find effective scaling factors [3], we instead opted for a more straightforward approach by using the max norm of the pre-cast matrix at each step of the splitting process. A shifting factor $s_A$ can also be provided so that the algorithm splits $A/s_A$ instead, and reinjects this factor in every $s_k$. This can further assist in reducing the risk of overflows. Finally, for the PPGEMM kernel presented in this work, we only consider splitting an FP32 input matrix into two FP16.

<caption>**Algorithm 1:** SPLIT: Split a matrix into low-precision words</caption>

**Input:** A: Input matrix of precision $fl_A$
$fl_{0:n-1}$: Precision of each word $A_k$ for $0 \le k < n$
$s_A$: Shifting factor for A.
**Output:** $A_{0:n-1}, s_{0:n-1}|A \approx \sum_{k=0}^{n-1} s_k fl_A(A_k)$
1. $\{s_0, A_0\} \leftarrow fl_0(A/s_A)$
2. **for** $k \leftarrow 1$ **to** $n-1$ **do**
3. &nbsp;&nbsp;&nbsp;&nbsp;$$ \{s_k, A_k\} \leftarrow fl_k(A/s_A - \sum_{j=0}^{k-1} (fl_A(A_j)s_j)) $$
4. &nbsp;&nbsp;&nbsp;&nbsp;$$ s_k \leftarrow s_A s_k $$
5. **return** $\{A_{0:n-1}, s_{0:n-1}\}$

## 3.2 PPGEMM kernel

In Algorithm 2, we present our PPGEMM kernel, named 2f16h, that leverages multiword splitting of a matrix for latency hiding. We use the **upload** and **download** keywords to indicate an asynchronous data transfer from the CPU to the GPU or from the GPU to the CPU, respectively. An **on stream** block denotes the code executed on the GPU on a specific CUDA stream. All code outside these blocks runs on the CPU. Each **on stream** block is divided into three columns. The first column contains the events the stream awaits to continue its execution. The second column contains the code to be executed. The last column indicates events fired after the code is executed. The **sync** keyword expresses a synchronization barrier, requiring all previously launched GPU kernels to terminate before continuing the execution. We use HS (half-single) to denote an FP16-multiply FP32-accumulate matrix multiplication, for which we call the cuBLAS implementation [8].
<table><caption>**Algorithm 2:** 2f16h PPGEMM</caption><thead><tr><th colspan="2">Input: A, B: Input matrices<br>Output: C: Output matrix such that C ≈ AB</th></tr></thead><tbody><tr><td>1</td><td>{A<sub>0:1</sub>, s<sub>0:1</sub><sup>A</sup>} ← SPLIT(A, {f16, f16}, ||A||<sub>∞</sub>)</td></tr><tr><td>2</td><td>{B<sub>0:1</sub>, s<sub>0:1</sub><sup>B</sup>} ← SPLIT(B, {f16, f16}, ||B||<sub>∞</sub>)</td></tr><tr><td>3</td><td>on stream S<sub>copy</sub> do</td></tr><tr><td>4</td><td>upload A<sub>0</sub></td><td>A<sub>0</sub></td></tr><tr><td>5</td><td>upload B<sub>0</sub></td><td>B<sub>0</sub></td></tr><tr><td>6</td><td>upload A<sub>1</sub></td><td>A<sub>1</sub></td></tr><tr><td>7</td><td>upload B<sub>1</sub></td><td>B<sub>1</sub></td></tr><tr><td>8</td><td>on stream S<sub>1</sub> do</td></tr><tr><td>9</td><td>A<sub>0</sub>, B<sub>0</sub> C<sub>t</sub> ← HS(s<sub>0</sub><sup>A</sup>s<sub>0</sub><sup>B</sup>A<sub>0</sub>B<sub>0</sub>)</td><td>C<sub>t</sub></td></tr><tr><td>10</td><td>{C<sub>0</sub>, C<sub>t</sub>}, {s<sub>0</sub><sup>C</sup>, s<sub>t</sub><sup>C</sup>} ← SPLIT(C<sub>t</sub>, {f16, f32}, s<sub>0</sub><sup>A</sup>s<sub>0</sub><sup>B</sup>)</td><td></td></tr><tr><td>11</td><td>download C<sub>0</sub></td><td></td></tr><tr><td>12</td><td>on stream S<sub>2</sub> do</td></tr><tr><td>13</td><td>C<sub>t</sub>, A<sub>1</sub> C<sub>t</sub> ← HS(s<sub>1</sub><sup>A</sup>s<sub>0</sub><sup>B</sup>A<sub>1</sub>B<sub>0</sub> + s<sub>t</sub><sup>C</sup>C<sub>t</sub>)</td><td></td></tr><tr><td>14</td><td>B<sub>1</sub> C<sub>t</sub> ← HS(s<sub>0</sub><sup>A</sup>s<sub>1</sub><sup>B</sup>A<sub>0</sub>B<sub>1</sub> + C<sub>t</sub>)</td><td></td></tr><tr><td>15</td><td>{C<sub>1</sub>, s<sub>1</sub><sup>C</sup>} ← SPLIT(C<sub>t</sub>, {f16}, s<sub>0</sub><sup>A</sup>s<sub>0</sub><sup>B</sup>2<sup>-11</sup>)</td><td></td></tr><tr><td>16</td><td>download C<sub>1</sub></td><td></td></tr><tr><td>17</td><td>sync</td><td></td></tr><tr><td>18</td><td>C ← f32(s<sub>0</sub><sup>C</sup>C<sub>0</sub> + s<sub>1</sub><sup>C</sup>C<sub>1</sub>)</td><td></td></tr></tbody></table>

The 2f16h kernel provided in Algorithm 2 multiplies two FP32 matrices A and B provided as input and computes $C = AB$ in FP32. First, in lines 1-2, both A and B are split into two FP16 words $A_0, A_1$ and $B_0, B_1$ together with their respective scaling coefficients $s_0^A, s_1^A$ and $s_0^B, s_1^B$ on the CPU. Next, in lines 3-7, all these words are scheduled for transfer to the CPU on the stream $S_{copy}$ in the order $A_0, B_0, A_1, B_1$ together with corresponding events $A_0, B_0, A_1, B_1$ so that the GEMM of high-magnitude words $A_0B_0$ can start as soon as possible. In lines 8-11, the stream $S_1$ waits for events $A_0$ and $B_0$ in parallel in order to launch FP16-multiply FP32-accumulate GEMM (i.e., $HS(A_0B_0)$), which is then split into an FP16 and a FP32 word $C_0$ and $C_t$, respectively. After this splitting, $C_0$ contains the high-magnitude word of the output matrix C, which can immediately be scheduled to transfer back to the CPU as the last operation on the same stream, and the event $C_t$ is fired. Here, it is important to keep $C_t$ in FP32 since the two remaining matrix multiplications will still need to be accumulated on it in FP32. Finally, in lines 12-16, the stream $S_2$ waits for events $C_t, A_1$, and $B_1$ to perform the remaining two multiplications $HS(A_0B_1)$ and $HS(A_1B_0)$, accumulate them on $C_t$. These refinement multiplications are expected to overlap with the transfer of $C_0$ to the CPU. We then utilize the SPLIT kernel to truncate $C_t$ to FP16 to obtain the low-magnitude word $C_1$ of the output matrix with the scaling factor $s_1^C$ and download it to the CPU.
<table><caption>**Algorithm 2: 2f16h PPGEMM**</caption><thead><tr><th></th><th>Input: A, B: Input matrices<br>Output: C: Output matrix such that C ≈ AB</th><th></th></tr></thead><tbody><tr><td>1</td><td>{A<sub>0:1</sub>, s<sub>0:1</sub><sup>A</sup>} ← SPLIT(A, {f16, f16}, ||A||<sub>∞</sub>)</td><td></td></tr><tr><td>2</td><td>{B<sub>0:1</sub>, s<sub>0:1</sub><sup>B</sup>} ← SPLIT(B, {f16, f16}, ||B||<sub>∞</sub>)</td><td></td></tr><tr><td>3</td><td>on stream S<sub>copy</sub> do</td><td></td></tr><tr><td>4</td><td>upload A<sub>0</sub></td><td>A<sub>0</sub></td></tr><tr><td>5</td><td>upload B<sub>0</sub></td><td>B<sub>0</sub></td></tr><tr><td>6</td><td>upload A<sub>1</sub></td><td>A<sub>1</sub></td></tr><tr><td>7</td><td>upload B<sub>1</sub></td><td>B<sub>1</sub></td></tr><tr><td>8</td><td>on stream S<sub>1</sub> do</td><td></td></tr><tr><td>9</td><td>A<sub>0</sub>, B<sub>0</sub> | C<sub>t</sub> ← HS(s<sub>0</sub><sup>A</sup>s<sub>0</sub><sup>B</sup>A<sub>0</sub>B<sub>0</sub>)</td><td></td></tr><tr><td>10</td><td>download C<sub>0</sub></td><td>C<sub>t</sub></td></tr><tr><td>11</td><td>download C<sub>1</sub></td><td></td></tr><tr><td>12</td><td>on stream S<sub>2</sub> do</td><td></td></tr><tr><td>13</td><td>C<sub>t</sub>, A<sub>1</sub> | C<sub>t</sub> ← HS(s<sub>1</sub><sup>A</sup>s<sub>0</sub><sup>B</sup>A<sub>1</sub>B<sub>0</sub> + s<sub>t</sub><sup>C</sup>C<sub>t</sub>)</td><td></td></tr><tr><td>14</td><td>B<sub>1</sub> | C<sub>t</sub> ← HS(s<sub>0</sub><sup>A</sup>s<sub>1</sub><sup>B</sup>A<sub>0</sub>B<sub>1</sub> + C<sub>t</sub>)</td><td></td></tr><tr><td>15</td><td>download C<sub>1</sub></td><td></td></tr><tr><td>16</td><td>download C<sub>0</sub></td><td></td></tr><tr><td>17</td><td>sync</td><td></td></tr><tr><td>18</td><td>C ← f32(s<sub>0</sub><sup>C</sup>C<sub>0</sub> + s<sub>1</sub><sup>C</sup>C<sub>1</sub>)</td><td></td></tr></tbody></table>

The 2f16h kernel provided in Algorithm 2 multiplies two FP32 matrices A and B provided as input and computes $C = AB$ in FP32. First, in lines 1-2, both A and B are split into two FP16 words $A_0, A_1$ and $B_0, B_1$ together with their respective scaling coefficients $s_0^A, s_1^A$ and $s_0^B, s_1^B$ on the CPU. Next, in lines 3-7, all these words are scheduled for transfer to the CPU on the stream $S_{copy}$ in the order $A_0, B_0, A_1, B_1$ together with corresponding events $A_0, B_0, A_1, B_1$ so that the GEMM of high-magnitude words $A_0B_0$ can start as soon as possible. In lines 8-11, the stream $S_1$ waits for events $A_0$ and $B_0$ in parallel in order to launch FP16-multiply FP32-accumulate GEMM (i.e., $HS(A_0B_0)$), which is then split into an FP16 and a FP32 word $C_0$ and $C_t$, respectively. After this splitting, $C_0$ contains the high-magnitude word of the output matrix C, which can immediately be scheduled to transfer back to the CPU as the last operation on the same stream, and the event $C_t$ is fired. Here, it is important to keep $C_t$ in FP32 since the two remaining matrix multiplications will still need to be accumulated on it in FP32. Finally, in lines 12-16, the stream $S_2$ waits for events $C_t, A_1$, and $B_1$ to perform the remaining two multiplications $HS(A_0B_1)$ and $HS(A_1B_0)$, accumulate them on $C_t$. These refinement multiplications are expected to overlap with the transfer of $C_0$ to the CPU. We then utilize the SPLIT kernel to truncate $C_t$ to FP16 to obtain the low-magnitude word $C_1$ of the output matrix with the scaling factor $s_1^C$ and download it to the CPU.
**2x2 tiled GEMM implementation using streams** To benchmark the PP-GEMM kernel, we use a tiled GEMM implementation which divides the matrices into $2 \times 2$ tiles of size $(N/2) \times (N/2)$ each to hide the transfer latency (e.g., $A_{00}, A_{01}, A_{10}, A_{11}$), uploads the input matrix tiles to the GPU in an optimal order, and performs tile multiplications as soon as they are uploaded where each tile of the output matrix gets two contributions (e.g., $C_{10} = A_{10}B_{00} + A_{11}B_{10}$) and is downloaded as soon as it is finalized. We employ multiple streams for uploads and downloads and to perform the GEMMs in parallel where possible to saturate the GPU. Stream parallelism could also be employed in the PPGEMM kernel to launch HS GEMMs in parallel. However, this would require multiple buffers (e.g., $C_t$) needing to be merged afterward, which resulted in an overall slowdown when tested. We then multiply the matrix tiles using either the SIMT cuBLAS SGEMM kernel, the 3TF32 kernel provided by Cutlass that uses three TF32 matrix multiplications, or a custom 3F16 kernel we developed that splits each FP32 tile into two FP16 words with scaling on the GPU, then performs three FP16-multiply FP32-accumulate GEMMs using cuBLAS.

# 4 Experiments

## 4.1 Setup

In our benchmarks, we use two GPUs from two widely used recent architectures: V100 (Volta [6]), and A100 (Ampere [7]). Note that the V100 tensor cores cannot perform TF32-based MMA operations, but both A100 and V100 provide FP16-FP32 MMA support. The A100 is installed on a server with two 20-core Intel Xeon 6346 CPUs clocked at 3.1Ghz with 36MB of total cache per socket and 1TiB of six-channel DDR4-3200 memory and uses PCIe 4.0. The V100 is installed on a server with two 16-core Intel Xeon 6230 CPUs clocked at 2.1Ghz with 27.5MB of total cache per socket and 768GiB of six-channel DDR4-2933 memory and uses PCIe 3.0. Since none of CPUs in these machines are equipped with AVX512F16 support, we use the Cuda Toolkit that provides software support for FP16 simulation in order to split matrices on the CPU on these machines. To accurately measure the CPU-side matrix splitting cost, we used another machine with two 56 cores Intel Xeon Platinum 8480+ with 105MB of total cache per socket and 512GiB of eight-channel DDR5-4400 memory, and used an optimized matrix splitting kernel we developed that leverages vectorization using AVX512 intrinsics and multicore parallelism using OpenMP.

We compare the PPGEMM implementation with the monolithic and a 2x2 tiled version of an FP32 tensor core-based GEMM. We report the results using 3TF32 and 3F16 kernels for monolithic and tiled multiplications on A100, and using only the 3F16 kernel on V100 due to the lack of TF32 support on this architecture. We perform our benchmarks on square matrices of sizes that are powers of two ranging from $2^6$ to $2^{14}$. We exclude smaller dimensions since using a GPU for very small matrices can hardly be beneficial, and the fixed overhead incurred by PPGEMM (due to multiple memory transfers, multiword splits, etc.) would render it even less attractive. Similarly, we do not measure our
**2x2 tiled GEMM implementation using streams** To benchmark the PP-GEMM kernel, we use a tiled GEMM implementation which divides the matrices into 2 × 2 tiles of size (N/2) × (N/2) each to hide the transfer latency (e.g., A₀₀, A₀₁, A₁₀, A₁₁), uploads the input matrix tiles to the GPU in an optimal order, and performs tile multiplications as soon as they are uploaded where each tile of the output matrix gets two contributions (e.g., C₁₀ = A₁₀B₀₀ + A₁₁B₁₀) and is downloaded as soon as it is finalized. We employ multiple streams for uploads and downloads and to perform the GEMMs in parallel where possible to saturate the GPU. Stream parallelism could also be employed in the PPGEMM kernel to launch HS GEMMs in parallel. However, this would require multiple buffers (e.g., Cₜ) needing to be merged afterward, which resulted in an overall slowdown when tested. We then multiply the matrix tiles using either the SIMT cuBLAS SGEMM kernel, the 3TF32 kernel provided by Cutlass that uses three TF32 matrix multiplications, or a custom 3F16 kernel we developed that splits each FP32 tile into two FP16 words with scaling on the GPU, then performs three FP16-multiply FP32-accumulate GEMMs using cuBLAS.

# 4 Experiments

## 4.1 Setup

In our benchmarks, we use two GPUs from two widely used recent architectures: V100 (Volta [6]), and A100 (Ampere [7]). Note that the V100 tensor cores cannot perform TF32-based MMA operations, but both A100 and V100 provide FP16-FP32 MMA support. The A100 is installed on a server with two 20-core Intel Xeon 6346 CPUs clocked at 3.1Ghz with 36MB of total cache per socket and 1TiB of six-channel DDR4-3200 memory and uses PCIe 4.0. The V100 is installed on a server with two 16-core Intel Xeon 6230 CPUs clocked at 2.1Ghz with 27.5MB of total cache per socket and 768GiB of six-channel DDR4-2933 memory and uses PCIe 3.0. Since none of CPUs in these machines are equipped with AVX512F16 support, we use the Cuda Toolkit that provides software support for FP16 simulation in order to split matrices on the CPU on these machines. To accurately measure the CPU-side matrix splitting cost, we used another machine with two 56 cores Intel Xeon Platinum 8480+ with 105MB of total cache per socket and 512GiB of eight-channel DDR5-4400 memory, and used an optimized matrix splitting kernel we developed that leverages vectorization using AVX512 intrinsics and multicore parallelism using OpenMP.

We compare the PPGEMM implementation with the monolithic and a 2x2 tiled version of an FP32 tensor core-based GEMM. We report the results using 3TF32 and 3F16 kernels for monolithic and tiled multiplications on A100, and using only the 3F16 kernel on V100 due to the lack of TF32 support on this architecture. We perform our benchmarks on square matrices of sizes that are powers of two ranging from 2⁶ to 2¹⁴. We exclude smaller dimensions since using a GPU for very small matrices can hardly be beneficial, and the fixed overhead incurred by PPGEMM (due to multiple memory transfers, multiword splits, etc.) would render it even less attractive. Similarly, we do not measure our
implementations' performance on matrices of dimensions past $2^{14}$. This is because the overall computation becomes entirely compute-bound in this scenario, eliminating the need for the PPGEMM kernel.

Every benchmark consists of a 10 seconds warmup phase, which allows the GPUs to reach a stable frequency dependent on the workload being tested, reducing the noise in the measurements. This is followed by a 20 seconds measurement phase. The data points shown in the figures are the median of all samples collected during the measurement phase to reduce the impact of outliers on the results, as spikes of latency are to be expected when measuring CPU and GPU workloads of extremely short duration. The CPUs used during the benchmark are locked to their base frequency to further reduce the possibility of outliers in the measurement. Further characteristics of these GPUs are provided in Table 1.

The performance of the various implementations shown in Figure 1 is measured in *effective flops/s*, which we define as $2n^3/t_{total}$, where *n* is the dimension of matrices *A*, *B*, and *C*, and *t*<sub>total</sub> is the total time spent on the GEMM operation including memory transfers but excluding some initialization steps such as memory pinning. We report the performance of the PPGEMM kernel with and without the CPU-side matrix splitting cost. The former can be pertinent for specific mixed-precision numerical algorithms that already perform such a split during their execution [2].

Each benchmark is performed on the same input matrices with uniformly random elements in the interval (-0.5, 0.5] and [0, 1] to keep the measurements consistent with [1]. This two intervals respectively represet the lowest and the highest rounding error accumulation.

## 4.2 Performance analysis of PPGEMM kernels

In Figure 1, we provide the relative performance measurements of the PPGEMM kernel as well as multiple cuBLAS- and CUTLASS-based GEMM implementations for varying matrix sizes. The results are normalized with the performance of a monolithic 3F16-based SGEMM, which can be performed on both A100 and V100 platforms. **2f16h** represents the PPGEMM kernel in Algorithm Algorithm 2 without the CPU-side matrix splitting cost. We also include two other curves that involve this cost for the case when the matrix is already in the cache or outside the cache. **SGEMM 3f16** and **SGEMM 3tf32** are the standard transfer-transfer-GEMM-transfer execution using the 3F16 and 3TF32 kernels, respectively. **SGEMM 3f16 2x2** and **SGEMM 3tf32 2x2** represent 2x2 tiled versions of these kernels.

The PPGEMM kernel achieves a speedup of up to 39% and 25% on A100 and V100, respectively, relative to **SGEMM 3f16** for small to medium matrix sizes ($n \le 2^{10}$), and is significantly faster than the tiled variants. The kernel still manages to achieve a speedup of up to 21% and 20% w.r.t the baseline on A100 and V100, respectively, when considering the splitting cost when the input matrices are inside the CPU cache. When the input matrices are outside the cache, the PPGEMM still provides a speedup over the tiled variants on the V100 platform. For large matrix sizes ($n \ge 2048$), the tiled variants catch
implementations' performance on matrices of dimensions past $2^{14}$. This is because the overall computation becomes entirely compute-bound in this scenario, eliminating the need for the PPGEMM kernel.

Every benchmark consists of a 10 seconds warmup phase, which allows the GPUs to reach a stable frequency dependent on the workload being tested, reducing the noise in the measurements. This is followed by a 20 seconds measurement phase. The data points shown in the figures are the median of all samples collected during the measurement phase to reduce the impact of outliers on the results, as spikes of latency are to be expected when measuring CPU and GPU workloads of extremely short duration. The CPUs used during the benchmark are locked to their base frequency to further reduce the possibility of outliers in the measurement. Further characteristics of these GPUs are provided in Table 1.

The performance of the various implementations shown in Figure 1 is measured in *effective flops/s*, which we define as $2n^3/t_{total}$, where *n* is the dimension of matrices *A*, *B*, and *C*, and *t*_{total} is the total time spent on the GEMM operation including memory transfers but excluding some initialization steps such as memory pinning. We report the performance of the PPGEMM kernel with and without the CPU-side matrix splitting cost. The former can be pertinent for specific mixed-precision numerical algorithms that already perform such a split during their execution [2].

Each benchmark is performed on the same input matrices with uniformly random elements in the interval (-0.5, 0.5] and [0, 1] to keep the measurements consistent with [1]. This two intervals respectively represet the lowest and the highest rounding error accumulation.

## 4.2 Performance analysis of PPGEMM kernels

In Figure 1, we provide the relative performance measurements of the PPGEMM kernel as well as multiple cuBLAS- and CUTLASS-based GEMM implementations for varying matrix sizes. The results are normalized with the performance of a monolithic 3F16-based SGEMM, which can be performed on both A100 and V100 platforms. **2f16h** represents the PPGEMM kernel in Algorithm Algorithm 2 without the CPU-side matrix splitting cost. We also include two other curves that involve this cost for the case when the matrix is already in the cache or outside the cache. **SGEMM 3f16** and **SGEMM 3tf32** are the standard transfer-transfer-GEMM-transfer execution using the 3F16 and 3TF32 kernels, respectively. **SGEMM 3f16 2x2** and **SGEMM 3tf32 2x2** represent 2x2 tiled versions of these kernels.

The PPGEMM kernel achieves a speedup of up to 39% and 25% on A100 and V100, respectively, relative to **SGEMM 3f16** for small to medium matrix sizes ($n \le 2^{10}$), and is significantly faster than the tiled variants. The kernel still manages to achieve a speedup of up to 21% and 20% w.r.t the baseline on A100 and V100, respectively, when considering the splitting cost when the input matrices are inside the CPU cache. When the input matrices are outside the cache, the PPGEMM still provides a speedup over the tiled variants on the V100 platform. For large matrix sizes ($n \ge 2048$), the tiled variants catch
<figure><img src="image_3.png" /><figcaption>Fig. 1. Performance comparison of all GEMM variants.</figcaption></figure>

up with PPGEMM as expected due to the operation becoming more and more
compute-bound. We conclude that using the PPGEMM can indeed improve
the performance of the GEMM operation in the matrix size range 128 × 128
to 2048 × 2048, depending on the GPU architecture, and can provide a faster
drop-in replacement for FP32 GEMM in this size interval. One can design a
poly-algorithm to decide when to use monolithic, tiled, or PPGEMM offline for
each architecture and matrix size.

## 4.3 Cost of multiword splitting on the CPU

In Figure 2, we provide the relative cost of the splitting operation relative to the
entire PPGEMM operation for each matrix size of interest in two scenarios: one
where the matrix was already loaded in the cache (if possible) and one where the
matrix was always outside the cache, which we emulated by flushing the cache
between two benchmark iterations using AVX intrinsics.

The splitting cost is significantly higher when the matrix is not in the CPU
cache. This is due to the lower bandwidth and the higher latency of the main
memory compared to the CPU cache, as splitting is a highly memory-bound com-
putation. This also means that the operation should scale with higher memory
bandwidth or cache size/bandwidth, which are becoming increasingly available
on modern platforms (i.e., high-bandwidth memory support, 3D cache, etc.).
We also observe that the relative cost increases significantly when operating on
very small matrices. When the matrix is already in the cache, the overhead of
splitting becomes reasonable, particularly for the size range in which PPGEMM
is advantageous (i.e., $2^8 \le n \le 2^{10}$).
<figure><img src="image_3.png" /><figcaption>Fig. 1. Performance comparison of all GEMM variants.</figcaption></figure>

up with PPGEMM as expected due to the operation becoming more and more
compute-bound. We conclude that using the PPGEMM can indeed improve
the performance of the GEMM operation in the matrix size range 128 × 128
to 2048 × 2048, depending on the GPU architecture, and can provide a faster
drop-in replacement for FP32 GEMM in this size interval. One can design a
poly-algorithm to decide when to use monolithic, tiled, or PPGEMM offline for
each architecture and matrix size.

## 4.3 Cost of multiword splitting on the CPU

In Figure 2, we provide the relative cost of the splitting operation relative to the
entire PPGEMM operation for each matrix size of interest in two scenarios: one
where the matrix was already loaded in the cache (if possible) and one where the
matrix was always outside the cache, which we emulated by flushing the cache
between two benchmark iterations using AVX intrinsics.

The splitting cost is significantly higher when the matrix is not in the CPU
cache. This is due to the lower bandwidth and the higher latency of the main
memory compared to the CPU cache, as splitting is a highly memory-bound com-
putation. This also means that the operation should scale with higher memory
bandwidth or cache size/bandwidth, which are becoming increasingly available
on modern platforms (i.e., high-bandwidth memory support, 3D cache, etc.).
We also observe that the relative cost increases significantly when operating on
very small matrices. When the matrix is already in the cache, the overhead of
splitting becomes reasonable, particularly for the size range in which PPGEMM
is advantageous (i.e., $2^8 \le n \le 2^{10}$).
<figure><img src="image_4.png" /><figcaption>Fig. 2. Splitting cost</figcaption></figure>

## 4.4 Numerical error of PPGEMM kernels

In this section, we provide experimental results on the A100 platform for the error obtained from the PPGEMM kernel and other variants for random matrices with elements in [-0.5, 0.5) and [-0.5, 0.5) intervals. Error analyses for multiword [1] or cascading GEMM [12] are pertinent and can apply to the PPGEMM kernel. In Table 2, we report the relative error of the resulting matrix C obtained from each kernel with respect to the reference result $C_{ref}$ obtained from the OpenBLAS DGEMM kernel [11], i.e., error = $||C - C_{ref}||/||C_{ref}||$ where $||.||$ represents the matrix Frobenius norm. We note that for the [-0.5, 0.5) and [0, 1) interval, the kernels provide comparable relative error and that those using the FP16-multiply FP32-accumulate kernels involve slightly more error, particularly for larger matrix sizes due to error accumulation.

# 5 Conclusion

In this paper, we presented a new kernel called PPGEMM that is capable of more effectively hiding CPU-GPU transfer latency by using a multiword splitting of the input/output matrices in lower precision and overlapping the data transfers with computations on the GPU. Our kernel shows a significant speedup when compared to traditional monolithic and tiled GEMM implementations for certain matrix sizes. Overall, our experiments demonstrated a speedup of up to 39% using PPGEMM compared to the best alternative.

Although this work focuses on hiding the CPU-GPU transfer latency, PPGEMM has a strong potential for distributed memory kernels as well for hiding
<figure><img src="image_4.png" /><figcaption>Fig. 2. Splitting cost</figcaption></figure>

## 4.4 Numerical error of PPGEMM kernels

In this section, we provide experimental results on the A100 platform for the error obtained from the PPGEMM kernel and other variants for random matrices with elements in [-0.5, 0.5) and [-0.5, 0.5) intervals. Error analyses for multiword [1] or cascading GEMM [12] are pertinent and can apply to the PPGEMM kernel. In Table 2, we report the relative error of the resulting matrix C obtained from each kernel with respect to the reference result $C_{ref}$ obtained from the OpenBLAS DGEMM kernel [11], i.e., error = $||C - C_{ref}||/||C_{ref}||$ where $||.||$ represents the matrix Frobenius norm. We note that for the [-0.5, 0.5) and [0, 1) interval, the kernels provide comparable relative error and that those using the FP16-multiply FP32-accumulate kernels involve slightly more error, particularly for larger matrix sizes due to error accumulation.

# 5 Conclusion

In this paper, we presented a new kernel called PPGEMM that is capable of more effectively hiding CPU-GPU transfer latency by using a multiword splitting of the input/output matrices in lower precision and overlapping the data transfers with computations on the GPU. Our kernel shows a significant speedup when compared to traditional monolithic and tiled GEMM implementations for certain matrix sizes. Overall, our experiments demonstrated a speedup of up to 39% using PPGEMM compared to the best alternative.

Although this work focuses on hiding the CPU-GPU transfer latency, PPGEMM has a strong potential for distributed memory kernels as well for hiding
<table><caption>**Table 2.** PPGEMM kernel precision comparison on A100</caption><thead><tr><th rowspan="2">Kernel</th><th rowspan="2">Matrix Type</th><th colspan="4">Matrix Size</th></tr><tr><th>9</th><th>10</th><th>11</th><th>12</th></tr></thead><tbody><tr><td rowspan="2">2f16h</td><td>[-0.5,0.5)</td><td>9.05e-07</td><td>9.06e-07</td><td>1.73e-06</td><td>6.65e-06</td></tr><tr><td>[0,1)</td><td>5.20e-06</td><td>5.20e-06</td><td>1.09e-05</td><td>4.52e-05</td></tr><tr><td rowspan="2">sgemm 3f16</td><td>[-0.5,0.5)</td><td>9.47e-07</td><td>9.65e-07</td><td>1.77e-06</td><td>6.65e-06</td></tr><tr><td>[0,1)</td><td>5.20e-06</td><td>5.20e-06</td><td>1.09e-05</td><td>4.52e-05</td></tr><tr><td rowspan="2">sgemm 3f16 2x2</td><td>[-0.5,0.5)</td><td>5.64e-07</td><td>9.65e-07</td><td>9.68e-07</td><td>1.77e-06</td></tr><tr><td>[0,1)</td><td>2.41e-06</td><td>5.20e-06</td><td>5.20e-06</td><td>1.09e-05</td></tr><tr><td rowspan="2">sgemm 3tf32</td><td>[-0.5,0.5)</td><td>3.88e-07</td><td>4.38e-07</td><td>4.68e-07</td><td>5.15e-07</td></tr><tr><td>[0,1)</td><td>4.01e-07</td><td>4.16e-07</td><td>4.17e-07</td><td>4.35e-07</td></tr><tr><td rowspan="2">sgemm 3tf32 2x2</td><td>[-0.5,0.5)</td><td>3.80e-07</td><td>4.27e-07</td><td>4.47e-07</td><td>4.74e-07</td></tr><tr><td>[0,1)</td><td>3.89e-07</td><td>4.03e-07</td><td>3.95e-07</td><td>3.95e-07</td></tr></tbody></table>

the node-to-node communication latency. One can also extend this idea to develop PPGEMM kernels for FP64 target precision; we leave such perspectives to a future work.
<table><caption>**Table 2.** PPGEMM kernel precision comparison on A100</caption><thead><tr><th rowspan="2">Kernel</th><th rowspan="2">Matrix Type</th><th colspan="4">Matrix Size</th></tr><tr><th>9</th><th>10</th><th>11</th><th>12</th></tr></thead><tbody><tr><td rowspan="2">2f16h</td><td>[-0.5,0.5)</td><td>9.05e-07</td><td>9.06e-07</td><td>1.73e-06</td><td>6.65e-06</td></tr><tr><td>[0,1)</td><td>5.20e-06</td><td>5.20e-06</td><td>1.09e-05</td><td>4.52e-05</td></tr><tr><td rowspan="2">sgemm 3f16</td><td>[-0.5,0.5)</td><td>9.47e-07</td><td>9.65e-07</td><td>1.77e-06</td><td>6.65e-06</td></tr><tr><td>[0,1)</td><td>5.20e-06</td><td>5.20e-06</td><td>1.09e-05</td><td>4.52e-05</td></tr><tr><td rowspan="2">sgemm 3f16 2x2</td><td>[-0.5,0.5)</td><td>5.64e-07</td><td>9.65e-07</td><td>9.68e-07</td><td>1.77e-06</td></tr><tr><td>[0,1)</td><td>2.41e-06</td><td>5.20e-06</td><td>5.20e-06</td><td>1.09e-05</td></tr><tr><td rowspan="2">sgemm 3tf32</td><td>[-0.5,0.5)</td><td>3.88e-07</td><td>4.38e-07</td><td>4.68e-07</td><td>5.15e-07</td></tr><tr><td>[0,1)</td><td>4.01e-07</td><td>4.16e-07</td><td>4.17e-07</td><td>4.35e-07</td></tr><tr><td rowspan="2">sgemm 3tf32 2x2</td><td>[-0.5,0.5)</td><td>3.80e-07</td><td>4.27e-07</td><td>4.47e-07</td><td>4.74e-07</td></tr><tr><td>[0,1)</td><td>3.89e-07</td><td>4.03e-07</td><td>3.95e-07</td><td>3.95e-07</td></tr></tbody></table>

the node-to-node communication latency. One can also extend this idea to develop PPGEMM kernels for FP64 target precision; we leave such perspectives to a future work.
References

1. Massimiliano Fasi, Nicholas J. Higham, Florent Lopez, Theo Mary, and Mantas Mikaitis. Matrix multiplication in multiword arithmetic: Error analysis and application to gpu tensor cores. *SIAM Journal on Scientific Computing*, 45(1):C1-C19, 2023.

2. Nicholas J. Higham and Theo Mary. Mixed precision algorithms in numerical linear algebra. *Acta Numerica*, 31:347-414, 2022.

3. Nicholas J. Higham, Srikara Pranesh, and Mawussi Zounon. Squeezing a matrix into half precision, with an application to solving linear systems. *SIAM Journal on Scientific Computing*, 41(4):A2536-A2551, 2019.

4. Bengio Yoshua LeCun Yann and Geoffrey Hinton. Deep learning. *Nature*, 2015.

5. Igor L. Markov and Yaoyun Shi. Simulating quantum computation by contracting tensor networks. *SIAM Journal on Computing*, 38(3):963-981, 2008.

6. NVIDIA. Nvidia tesla v100 gpu architecture. 2017.

7. NVIDIA. Nvidia a100 tensor core gpu architecture. 2020.

8. NVIDIA. cublas - cuda basic linear algebra subroutine library. 2023.

9. NVIDIA. Cutlass - cuda templates for linear algebra subroutines. 2023.

10. NVIDIA. Nvidia h100 tensor core gpu architecture. 2023.

11. OpenMathLib. Openblas - an optimized blas library.

12. Devangi N. Parikh, Robert A. van de Geijn, and Greg M. Henry. Cascading gemm: High precision from low precision. 2023.

13. Grégoire Pichon, Eric Darve, Mathieu Faverge, Pierre Ramet, and Jean Roman. Sparse supernodal solver using block low-rank compression: Design, performance and analysis. *Journal of Computational Science*, 27:255-270, 2018.

14. Louis Pisha and Łukasz Ligowski. Accelerating non-power-of-2 size fourier transforms with gpu tensor cores. pages 507-516, 2021.

15. Samuel Williams, Andrew Waterman, and David Patterson. Roofline: an insightful visual performance model for multicore architectures. *Communications of the ACM*, 52(4):65-76, 2009.
References

1. Massimiliano Fasi, Nicholas J. Higham, Florent Lopez, Theo Mary, and Mantas Mikaitis. Matrix multiplication in multiword arithmetic: Error analysis and application to gpu tensor cores. *SIAM Journal on Scientific Computing*, 45(1):C1-C19, 2023.

2. Nicholas J. Higham and Theo Mary. Mixed precision algorithms in numerical linear algebra. *Acta Numerica*, 31:347-414, 2022.

3. Nicholas J. Higham, Srikara Pranesh, and Mawussi Zounon. Squeezing a matrix into half precision, with an application to solving linear systems. *SIAM Journal on Scientific Computing*, 41(4):A2536-A2551, 2019.

4. Bengio Yoshua LeCun Yann and Geoffrey Hinton. Deep learning. *Nature*, 2015.

5. Igor L. Markov and Yaoyun Shi. Simulating quantum computation by contracting tensor networks. *SIAM Journal on Computing*, 38(3):963-981, 2008.

6. NVIDIA. Nvidia tesla v100 gpu architecture. 2017.

7. NVIDIA. Nvidia a100 tensor core gpu architecture. 2020.

8. NVIDIA. cublas - cuda basic linear algebra subroutine library. 2023.

9. NVIDIA. Cutlass - cuda templates for linear algebra subroutines. 2023.

10. NVIDIA. Nvidia h100 tensor core gpu architecture. 2023.

11. OpenMathLib. Openblas - an optimized blas library.

12. Devangi N. Parikh, Robert A. van de Geijn, and Greg M. Henry. Cascading gemm: High precision from low precision. 2023.

13. Grégoire Pichon, Eric Darve, Mathieu Faverge, Pierre Ramet, and Jean Roman. Sparse supernodal solver using block low-rank compression: Design, performance and analysis. *Journal of Computational Science*, 27:255-270, 2018.

14. Louis Pisha and Łukasz Ligowski. Accelerating non-power-of-2 size fourier transforms with gpu tensor cores. pages 507-516, 2021.

15. Samuel Williams, Andrew Waterman, and David Patterson. Roofline: an insightful visual performance model for multicore architectures. *Communications of the ACM*, 52(4):65-76, 2009.