{'fraction_non_alphanumeric': 0.03365993467973072, 'fraction_numerical': 0.012219778266568908, 'mean_word_length': 4.565034619188922, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 12, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Accelerating Sparse Matrix-Matrix Multiplication with GPU Tensor Cores

Sparse general matrix-matrix multiplication (spGEMM) is an essential component in many scientific and data analytics applications. However, the sparsity pattern of the input matrices and the interaction of their patterns make spGEMM challenging. Modern GPUs include Tensor Core Units (TCUs), which specialize in dense matrix multiplication. Our aim is to re-purpose TCUs for sparse matrices. The key idea of our spGEMM algorithm, tSparse, is to multiply sparse rectangular blocks using the mixed precision mode of TCUs. tSparse partitions the input matrices into tiles and operates only on tiles which contain one or more elements. It creates a task list of the tiles, and performs matrix multiplication of these tiles using TCUs. To the best of our knowledge, this is the first time that TCUs are used in the context of spGEMM. We show that spGEMM, with our tiling approach, benefits from TCUs. Our approach significantly improves the performance of spGEMM in comparison to cuSPARSE, CUSP, RMerge2, Nsparse, AC-SpGEMM and spECK.


Introduction
Sparse general matrix-matrix multiplication (spGEMM), similar to its dense counterpart, performs the Matrix Multiplication (MM) of two sparse matrices. The main difference between sparse and dense MM is that we have to account for sparse matrices, which contain mostly zero elements. spGEMM is an important component in scientific and data analytics applications such as Graph analytics [1], Breadth-First-Search (BFS) [2], Algebraic Multigrid (AMG) [3], Schur complement [4], etc. Sparse matrices, which often contain millions of elements, require MM methods that do not waste computing resources on elements that are zero. The diverse structure and density of sparse matrices pose difficulties in regards to memory management and load balancing in parallel systems [5,6,7,8,9].
The re-emergence of deep learning motivated the creation of application specific integrated circuits (ASIC) that specialize in MM. MM is a core component of convolution [10] and these ASICs accelerate the calculation of the convolutional layers significantly in comparison to the normal, general purpose, processing cores. Such ASICs are Tensor Core Units (TCUs) from NVIDIA [11] and TPUs from Google [12]. We utilize the TCUs from NVIDIA to accelerate spGEMM for two reasons. First, accessibility. They are widely available as they are included in the new generation of GPUs from NVIDIA. Second, mixed precision. Mixed precision allows TCUs to mix 16-bit inputs and 32-bit multiplication and accumulation. Typically, in regards to deep learning 16-bits of precision (or less) is sufficient for training purposes and therefore tensor unit manufacturers opt only for 16-bit or lower precisions. Therefore, mixed precision of TCUs widens the application field to scientific problems which are more demanding w.r.t. precision.
Our work is motivated by three key observations. First, blocking sparse matrix storage formats [13], which group the elements of the matrix into rectangular tiles, are a good fit for TCUs which expect rectangular matrices as input. Second, TCUs are very efficient even when they are not fully occupied [14]. Third, even though TCUs support only low precision inputs, they can operate in mixed precision mode to perform operations that require higher precision, e.g., GEMM [15,16].
The key idea of our approach is to partition the input in tiles and multiply the tiles with TCUs. Tiles are sparse, but TCUs perform MM efficiently even when not fully occupied. Mixed precision mode is necessary in order to keep sufficient accuracy when multiplying large matrices.
Based on these observations, we propose a new GPU-based framework for spGEMM computation. Our novel methodology groups elements into tiles and uses the fast MM of TCUs to multiply the tiles. To the best of our knowledge, this is the first proposal of using TCUs in the context of spGEMM. Our TCU-based spGEMM methodology, which we name tSparse, has two advantages. First, it takes advantage of fast MM of TCUs. Second, by utilizing TCUs to process the MM, we leave the normal processing cores free for non-canonical workloads.
tSparse modifies Expand-Sort-Compress (ESC) methodology [3] of CUSP [17] to Sort-Expand and Compress (SEaC), i.e., tSparse brings both multiplication and accumulation after the sorting step. The benefits of this change are twofold. First, tSparse does not have to maintain in memory a large matrix for the intermediate products. Second, tSparse takes full advantage of Multiplication-Accumulation (MAC) operation of TCUs. To that end, we form a task list instead of calculating the intermediate products immediately. Our GPU kernels consume tiles from the task list and perform MM of the tiles using TCUs.
The rest of the paper is organized as follows. Section 2 gives a background on spGEMM. Section 3 presents a highlevel overview of tSparse, whereas Section 4 describes tSparse in-depth. Section 5 introduces our test configuration, which we use in Section 6 to evaluate the performance of tSparse. Section 7 presents related work. Finally, Section 8 concludes the paper.

Storage format
In sparse matrices, typically, the number of non-zero (NNZ) elements is much smaller than the number of zero elements. In order to save memory without degrading the performance of MM, we need an efficient way to store only non-zero (NZ) elements.

COO format
COO format stores each NZ element along with the coordinates the element would have in the dense representation of the matrix [22]. The COO format uses three arrays: for elements, for row indices and for column indices.

Bitmap format
TCUs simultaneously process multiple elements in rectangular structures. However, COO stores only single elements and has no concept of rectangular structures, therefore it is not sufficient by itself as a storage format for tSparse. For this reason, we use a bitmap-based block shaped storage format to store sparse matrices [13].
In our work, we use a bitmap format similar to [13] for three reasons: 1) it is simple and straightforward, 2) square tiles of fixed size fit well to TCUs, and 3) the performance of the format has been evaluated in [13]. The basic idea is to partition the matrices in a grid of 8×8 square blocks and work only on nonempty blocks. We refer to these blocks as tiles. The elements inside each tile have the same placement as they have in the dense representation of the matrix. Each element in the tile can be either zero or NZ. To keep track which elements are NZ we use a bitmap, a binary number of which each digit corresponds to one slot of the tile. If a slot contains a NZ we set the respective bit of the bitmap to "1", otherwise to "0". Then, tiles are stored in COO format. The difference with the standard COO format is that, instead of using single elements as values, now we use a tuplepair of two values: 1) an index to the element array. The element array holds the elements of the tile (elements of the same tile are in consecutive positions of the array), and 2) the bitmap of the tile. Fig. 1 shows how we convert a dense matrix to bitmap format. For simplicity we use 4×4 tiles. We represent the positions of NZ elements as "1"s in the bitmap. We store four values for each tile that has one or more NZ elements: 1) row index like in COO format, 2) column index like in COO format, 3) index into the element array, and 4) bitmap of the location of NZ elements in the tile.

Sparse matrix-matrix multiplication
The general matrix multiplication (GEMM) has the form: where A, B, C are the input matrices and D is the output. In spGEMM, similarly to dense matrices, to get one element of the output, we need to multiply the NZ elements of one row of A with the corresponding NZ elements of one column of B and then accumulate the intermediate products (i.e., calculate the inner product). The difference in spGEMM is that we multiply the corresponding elements only if the elements in the corresponding positions of the row of A and the column of B are both NZ and we accumulate only NZ products. The various ways to access the elements of A and B and perform MM are listed in [23]. The same multiplication method applies even if instead of elements we use tiles. In this case, the product of two corresponding tiles is their outer product (or equally MM). We make two important observations. First, (1) takes the form when accumulating the tiles, where C is both output and accumulator. We use only the form of (2) for the rest of this work. Second, the matrix Multiplication-Accumulation of small tiles is exactly what the TCUs were designed to do. Fig. 2 shows an example of the A * B sparse MM.

Challenges of spGEMM
Unlike other sparse matrix operations, in spGEMM both input and output are sparse. Therefore, it is very difficult to utilize the knowledge we infer from the sparsity structure of the input matrices to make arithmetic and memory optimizations. The reasons that make spGEMM more challenging than spMV (sparse matrix-vector multiplication) are three [5,6,7,8,9]: First, data access is highly irregular because it depends on the sparsity structure of both matrices and the interaction of both structures. MM is not trivial for two reasons: 1) it requires access of possibly distant memory locations to load the inputs, and 2) it requires inserting the result to the output with possibly irregular access patterns. Bitmap (hex) 0x829C 0x860 0x100 0x8230 0x8412 Bitmap storage format Second, it is computationally expensive to know the size of the output before the actual MM. According to Liu et al. [6], there are four methods to estimate how much memory we need to allocate for the output. First, the precise method, which makes an estimate that is very close to the actual size of the output, typically by doing a partial execution of the MM algorithm. Second, the upper bound method, which typically uses as upper bound the amount of intermediate products. Third, estimation using probability theory. Fourth, progressive memory increase, which allocates more memory if the previously allocated memory overflows. In all cases, at the end of spGEMM, we remove empty or unused entries from the allocated memory as necessary.
Third, load balancing. The sparsity structure of both matrices and the interaction of the structures determine the distribution of workload. NNZ elements in each row of the input or output may vary significantly, which makes it difficult to partition the workload among threads.

Expand-Sort-Compress and CUSP
CUSP [17] is a library that specializes in sparse matrix operations. It is open-source and easily accessible on GitHub. It is written in Thrust which makes it easy to read and port to other platforms. Therefore, it provides a good "boilerplate" to test our approach.
CUSP uses Expand-Sort-Compress ESC methodology. ESC performs the spGEMM in three steps [17,3]. First, Expand. ESC multiplies each NZ element a i, j of A with all NZ elements of row B( j, :) of B to get the intermediate products (no accumulation in this step) [23]. Second, Sort. ESC sorts the intermediate products of the previous step so that products that correspond to the same element of C are in consecutive positions. Third, Compress. ESC calculates each element of C by accumulating all respective products, which are in consecutive positions, thanks to the sorting step.

Real number representation in digital computer systems
Computer systems have to store real numbers in bit representation. Floating point numbers are a common representation. The location of the decimal point and the number of bits determines the precision and range of the represented numbers. We denote the 32-bit representation as fp32, whereas the 16-bit as fp16. In contemporary systems, typically, we use floating point numbers as defined in IEEE 754 technical standard [24]. Usually, the fewer the number of bits, the faster the processing of the numbers is.

Tensor Core Units
NVIDIA, with the latest generations of Graphical Processing Units (GPUs) [25], brought Tensor Cores to the mainstream market. Nvidia TCUs are ASICs that have the purpose of accelerating MM. Therefore, our work on spGEMM has significant benefits by properly adapting spGEMM to TCUs.
TCUs mainly target deep learning, which is not very demanding precision-wise. Therefore, to accelerate MM, TCUs usually work with lower precision number representation (16bit or less). However, fp16 or lower precision is detrimental to the output because precision and range of fp16 numbers can be insufficient when dealing with physical problems. To rectify this problem NVIDIA provides mixed precision functionality. Mixed precision allows to mix numbers of different precision. The two defining characteristics of the mixed precision implementation of NVIDIA are as follows. First, although inputs A and B are in fp16 precision, their multiplication happens in full precision. Second, the product is stored as fp32 to accumulators C and D [25]. Fig. 3 shows the two characteristics. Markidis et al. and Haidar et al. [16,15] evaluate the performance and precision of GEMM and linear equation solving using the mixed precision mode of TCUs. They show that TCUs can be used in other physical problems, outside deep learning. We use mixed precision functionality to extend the applicability of our work to non deep learning workloads.

Overview of our technique
Our TCU-based sparse GEMM technique alters ESC methodology in order to efficiently work with tiles. To that end, tSparse creates a task list which delivers tiles to TCUs to accelerate the MM of the tiles. At this point it is important to make a distinction between the terms elements and tiles. Elements represent a single real or integer number, whereas tiles represent a group of elements (Section 2.1.2). This section gives a high level overview of the components of tSparse and how they relate to each other.
3.1. Creating the task list and allocating memory for tiles tSparse creates the task list in two steps. First, tSparse determines which tiles of A will be multiplied with which tiles of B in order to get each tile of C. Second, unlike the expansion phase of ESC, tSparse does not immediately multiply the corresponding tiles of A and B. Instead, tSparse creates a task list that holds only the locations of the corresponding tuples of A and B. The entries of the task list are sorted so that pairs that correspond to the same tile of C are in consecutive positions. With the help of the task list, tSparse also estimates the number of C tiles in the output and allocates memory accordingly.

Counting kernel
Dense matrices contain as many elements as the product of their dimensions. However, in sparse matrix-matrix multiplication, the NNZ of the output depend on the structure of the two input matrices and the interaction of these structures throughout the MM. Therefore, we know the exact size of the output only after the matrices are multiplied. However, in order to allocate memory for storing the result of the MM we first need to know the exact size to allocate because GPU kernels cannot easily reallocate memory during their execution. The purpose of the counting kernel is to make an estimation of the size of memory we need to allocate in order to store the element array of C, before calling the multiplication kernel. Non-blocking spGEMM approaches allocate memory only for elements. tSparse allocates memory for elements, in addition to memory for tiles. However, allocating memory for elements is not as simple as calling a few parallel primitives because tile multiplication is a sparse MM. The sparsity structure of both input tiles and the interaction of the structures define the memory allocation requirements.
In our implementation, we use the precise method (Section 2.3) to get an estimate of the count of elements the output has, i.e., the counting kernel is a partial implementation of the multiplication kernel. The counting kernel requires shorter execution time than the multiplication kernel because it neither loads nor stores any elements. Instead, this kernel uses the bitmap to put zeros and ones in the place of the elements and "simulates" the MM. We allocate memory equal to the NNZ of the output of the simulated MM.

Multiplication kernel
Once we know how much memory to allocate for tiles and elements we use the multiplication kernel to multiply and accumulate all pairs of tiles of A and B that correspond to each C tile.

Putting everything together
In order to perform the matrix multiplication of A with B, we 1) determine which products need to be accumulated for each tile of C, and 2) allocate memory for the tiles of C and the element array. Using the task list and the counting kernel we determine the memory allocation size for the tiles of C and the element array, respectively. Subsequently, the multiplication kernel has everything it needs to multiply A with B.
We expect three benefits. First, by placing both the multiplication and accumulation steps of MM in the multiplication kernel we can use TCUs for MM. By moving MM to TCUs, the computational heavy MM is no longer a bottleneck of the spGEMM algorithm. Second, the use of bitmap format reduces memory consumption because one row index and one column index represents up to 64 elements [13]. Third, by grouping elements to tiles, we reduce the amount of values we have to manipulate and therefore there are additional performance benefits (e.g., fewer values to sort during sorting phase of ESC).

Implementation details
This section describes the three main parts of our approach, tSparse, in detail: 1) the creation of the task list, 2) the estimation of how much memory to allocate for the elements of the output, and 3) the MM of the tiles. It also describes some smaller components of our implementation.

Creating the task list and allocating memory for tiles
The four main parts of our algorithm are: 1) finding corresponding tiles of A and B, 2) filtering of zero products, 3) creating a task list with entries that point to the tuples of A and B, and 4) estimating how much memory to allocate for the tuples of the output. In detail: First, similarly to CUSP, tSparse uses parallel primitives from Thrust library to find the correspondence among tiles of A and tiles of B. The main difference with CUSP is that instead of using single elements we use tiles. This is the part of our methodology that is the same between CUSP and tSparse. We do not change it because it is efficient and accounts for less than 15% of the total execution time of our spGEMM.
Second, the algorithm that finds the correspondence among tiles of A and B considers the whole tile as a single value, i.e., the result is the same regardless if the 8×8 tile contains one or 64 elements. Therefore, the correspondence algorithm finds a correspondence between tiles of A and B even if the elements inside the two tiles are not corresponding. The resulting tile of the MM of such tiles is a tile with only zero elements. tSparse has a routine that removes this type of corresponding tiles by applying boolean arithmetic on bitmaps. This routine is fast as it does not compute MMs. The culling significantly lightens the workload of the rest of our methodology (fewer tiles to sort, multiply etc.).
Third, tSparse creates the task list. An important consideration is that in ESC methodology the MM and accumulation are in different steps, i.e., MM is in the Expand step and accumulation in the Compress step. However, to fully take advantage of the combined MAC operation of TCUs we need to perform both MM and accumulation of tiles in the same step. In order to put MM and accumulation in the same step we sort the locations of the multiplicands, instead of sorting the intermediate products of the MM. By sorting the locations, we defer the MM step until after the Sorting step. Effectively, the locations of the multiplicands form a task list, of which each entry points to one tile of A and one tile of B. Merging the expansion and compression steps and moving them after the sorting step significantly reduces both memory allocations and costly data movement to/from global memory. The reason is that tSparse does not store the intermediate products. Therefore, we save memory by not allocating memory for the 64-bit bitmaps, an undefined number of up to 64 elements and the corresponding indices in the element array. We reduce data movements to/from global memory by not storing the intermediate products before the sorting step and loading them for a second time after the sorting step.
Fourth, our algorithm counts how many of the intermediate products correspond to the same tile of C using a segmented reduce parallel primitive on the sorted indices. We create an offset array from the prefix sum of the counted intermediate products. Our GPU kernels use the offset array for indexing purposes. The prefix sum also gives the total count of tuples in the output. We use this total count to allocate memory.
Sorting. Sorting a long task list is computationally demanding, therefore we need advanced optimization strategies. tSparse works with sorted matrices. We use this knowledge to sort the workload of each row of A separately. For this task, we use the segmented sort by Kaixi et al. [26], which employs a hybrid sorting scheme based on row length.

Counting kernel
The counting kernel works in four steps. First, it reads the bitmaps of A and B. Second, it creates tiles, wherein each element is set to "1" or "0" based on the corresponding position in the bitmap. Third, it multiplies and accumulates the tiles of A and B that correspond to each individual tile of C using TCUs. Fourth, we count how many elements of the resulting tile are NZ with the ballot instruction. We repeat for all tiles of C and accumulate the counts. The counting kernel returns an array of which each value holds an estimation of how many elements each tile of C has.
Each tile contains only "1"s or "0"s, therefore: 1) the estimation of memory requirements can be more than what is actually required because the counting kernel cannot account for possible numerical cancellation (after the actual multiplication there is a compaction stage that removes empty entries from the allocated element array), and 2) half precision is enough for executing the MM of the counting kernel because we work only with zero and NZ elements.

Multiplication kernel
The multiplication kernel performs the actual multiplication and constructs the COO matrix of the output, i.e., it sets the row and column indices, the index and bitmap tuple and the elements of the element array. The multiplication kernel loads the actual elements from A and B and stores the result in the memory allocated by the counting kernel. We emphasize that we use only TCUs for the MM and not TCUs in addition to the standard CUDA cores.
There are two important considerations when multiplying the elements, which are real numbers. First, fp16 arithmetic has a very limited representation range of numbers , which we can easily exceed with multiplication. Therefore, we prefer the mixed precision functionality of TCUs. Second, unlike the counting kernel where we have only positive numbers, when accumulating real numbers, elements get canceled as a result of numerical cancellation. Many tiles may end up empty, something that the counting kernel, which acts on boolean elements, does not predict. For this reason, our multiplication kernel has the additional task of marking for removal tiles that are empty. Fig. 4 compares element loading between the counting and multiplication kernels. The counting kernel does not need to load any actual element. It just creates "1"s based on the bitmap. The multiplication kernel, on the other hand, loads the elements and it places them according to the bitmap.

16-bit Bitmap
Multiplication kernel Counting kernel

Other components
Arrangement of tiles in the TCUs. NVIDIA does not provide an Application Programming Interface (API) for using TCUs with small 8×8 tiles [11]. TCUs execute MM on 256 elements at a time. However, our tiles have a size of 8×8, which means that a large part of the TCU remains unused. Although a TCU does not have to be fully loaded in order to get performance benefits, we can fit two tiles in a single TCU. 16×16 is the only supported matrix configuration that can fit two 8×8 tiles [11]. To put two tiles in the same TCU two steps are necessary. First, we initialize the 16×16 matrix to zero. Second, the tiles must be placed in a diagonal of the 16×16 matrix (Fig. 5). If the tiles were not in the diagonal, but instead side-by-side (top-bottom), the same row (column) of the 16×16 matrix would have elements of two unrelated tiles, which would mix the inner products (of rows of A with columns of B) of the first tile with the inner products of the second tile. The API forces loading to/from TCUs through shared memory. We find the internal layout of the registers of the TCUs (fragments) and we access them directly instead. Loading through registers, which are faster than shared memory, we minimize data movement and increase the performance of MM. Load balancing. Both counting and multiplication kernels calculate an inner product of tiles (Section 2.2). The number of intermediate products required for each tile of C is different, because the number depends on the sparsity structure of A and B. Therefore the workload for the computation of each tile of C is different, leading to thread blocks with different amounts of work. Thread blocks with different execution times create imbalance among the SMs of the GPU. To tackle this issue, we assign each pair of tiles of C to a different thread block. When a thread block finishes and releases the resources, the scheduler of the SM schedules another thread block to take its place. Therefore, SMs take blocks according to their needs and stay fully occupied until completion of the MM.
Compaction of zeros. MM creates zero elements because of numerical cancellation. To store the output in a strictly sparse format we need to remove all zeros. Therefore, we need a way to detect zeros and remove them or, in other words, compact the arrays that hold the elements and the tuples. For the element array, which is just an array, we use a compaction parallel primitive from Thrust. For the array that holds the tuples, we first mark empty tiles in the multiplication kernel (see Section 4.3).

Putting everything together
Algorithm 1 summarizes tSparse. First, tSparse creates a task list (lines 1-8). Second, it estimates how much memory to allocate for tiles (lines 9-15). Third, it estimates how much memory to allocate for elements (lines [16][17] 16: element count ← CountingKernel(task list) 17: AllocateMemGPU(element count) 18: C tiles , C elements ← MultiplicationKernel(task list) 19: CompactZeroElements(C elements ) 20: CompactEmptyTiles(C tiles ) Figure 6 illustrates details of the multiplication kernel (Algorithm 1, line 18). In order to get a tile C of the output, we have to accumulate a varied number of products (MMs of tiles of A and B), e.g., C0 = A1 × B1 + A4 × B4 and C1 = A3 × B3 for the example of Fig. 2 (note the color code of tiles). Each TCU processes two tiles, e.g., TCU0 calculates C0 and C1. However, C1 has fewer addends. Therefore, when we load the second addend of C0, we use 0s in place of the second addend of C1. From the final result, we extract the bitmaps (using ballot) and elements of the output.
We compare tSparse with cuSPARSE from CUDA Toolkit [11], CUSP [17], RMerge2 [9], Nsparse [20], AC-SpGEMM [7], spECK [21]. In addition, to confirm the benefit of using TCUs, we create one implementation of tSparse without TCUs and we name it nonTCU. In nonTCU, we use the same method as [13] to multiply the tiles (Algorithm 2). We note that we include the time of memory allocations in the execution time of AC-SpGEMM. To evaluate the performance of tSparse, we perform the A * A MM, which has the benefit that both matrices have the same sparsity structure. The same structure makes it easier to make observations. In order to facilitate the presentation of the performance of tSparse, we select a subset of 16 matrices from SuiteSparse Matrix Collection [18] for our dataset. All selected matrices are square, as our A * A problem dictates. We select matrices which have elements in the fp16 range. Table 1 shows the characteristics of our dataset. We denote a matrix stored in bitmap storage format as C tiles . For nontiling (i.e., single element) approaches, we denote as C the total amount of intermediate products. Similarly, for our tiling approach, we denote as C tiles the total amount of intermediate products, where each product is a MM between tiles. NNZ(·) of a matrix denotes the NNZ values of the matrix. Finally, from the second and third columns of Table 1 derives the average row size of the input, which is defined as NNZ(A)/Dimensions, and we denote as RowA.
We divide our dataset into two parts. The first part (upper half of Table 1) consists of matrices that other works use [5,6,7,8,9,13,20,19,27,28]. The second part (lower half of Table 1) consists of matrices that we select after taking into consideration the two criteria that derive from the analysis in Section 6.1: 1) NNZ(A) > 300000, and 2) RowA > 42. .

Results and Analysis
In this section, we find the criteria that define when tSparse is the recommended spGEMM approach using all qualified matrices from SuiteSparse. Then, we collect four types of measurements using the matrices in Table 1: 1) the speedup of tSparse, 2) the execution time breakdown, 3) the numerical precision, and 4) the memory consumption.

Finding the selection criteria
We measure the execution time of cuSPARSE, CUSP, RMerge2, Nsparse, AC-SpGEMM, spECK, nonTCU and tSparse on a collection of matrices which we create from SuiteSparse. For this collection, we keep matrices that meet the following four conditions: 1) the input matrix has more than 10000 NNZ elements, 2) the input is square with real numbers in the fp16 range, 3) the bitmap density of the input is more than one, and 4) matrices for which all approaches that participate in the comparison return correct results. We then define criteria based on the characteristics of the matrices. These criteria facilitate the selection of the most appropriate spGEMM approach among tSparse and the other approaches.
We find the criteria in two steps. First, we analyze how each approach works in order to identify which matrix characteristics affect the performance of said approach. Second, we adjust the criteria to our matrix collection manually. Alternatively, we could use a machine-learning method to find the criteria [29], which we leave for future work. We make six major observations.
First, tSparse generally outperforms cuSPARSE in the larger matrices of our collection (NNZ(A) > 200000). This happens probably because cuSPARSE does not have sufficient shared memory to store the hash tables and therefore global memory traffic increases. Nevertheless, the size of the matrix is not the only thing that decides the performance of cuSPARSE. To the best of our knowledge, the performance of cuSPARSE also depends on the sparsity structure of the matrix because the structure affects the number of hash conflicts.
Second, tSparse performs better than CUSP when the number of intermediate products of tiles is smaller than the number of intermediate products of elements of CUSP, i.e., NNZ(C tiles ) < NNZ(C). Otherwise, CUSP, which handles only single elements, is faster than tSparse which has an additional overhead for handling tiles.
Third, tSparse is faster than Nsparse when the average row size of A ( RowA) is greater than 42 and NNZ(A) > 100000. Nsparse follows a hash table approach (Section 7.2). Therefore, larger rows possibly create more hash conflicts and hash tables are difficult to keep in shared memory (instead of the slower global memory) with large rows. We also note that web graph   Table 1: Matrix characteristics. We list the size of the matrix (number of rows/columns), the number of non-zeros of: the input (NNZ(A)), the output (NNZ(C)), the intermediate matrix (NNZ(C), the number of tiles (NNZ(C tiles )), and the number of tiles of the intermediate matrix (NNZ(C tiles )), and the density of the tiles of the input matrix. The upper part corresponds to matrices that are commonly used in the literature, the bottom part to matrices we selected based on our criteria. matrices (webbase-1M and web-NotreDame) have high irregularity, i.e., although the average row size is small (3.1 and 4.6 respectively), there are rows with thousands of elements (maximum 4700 and 3445 respectively). Therefore they are also taxing for Nsparse.
Fourth, RMerge2 shows behavior similar to Nsparse. tSparse has better performance than RMerge2 when RowA > 42. To the best of our knowledge, RMerge2, which follows a hybrid approach (Section 7.3), uses faster kernels when the row size is smaller than the size of a CUDA warp [9]. We note that the low recall score is owed to an instability of RMerge2 (possibly caused by a driver issue) that, in a seemingly random manner, decreases the performance.
Fifth, AC-SpGEMM is an ESC-based approach, like CUSP. Therefore, a similar criterion applies to AC-SpGEMM. The only difference is that, as AC-SpGEMM has much better performance than CUSP, the criterion is much stricter, i.e., 9 * NNZ(C tiles ) < NNZ(C).
Sixth, spECK is a hash table approach and has similar behavior to Nsparse. spECK generally performs better than Nsparse, therefore we select a criterion similar to Nsparse's but stricter, i.e., RowA > 42 and NNZ(A) > 300000.
Global criterion. In summary, a range of matrix characteristics define the performance of each spGEMM approach. To help with selecting tSparse over other approaches, we define two criteria that generally work well to show when tSparse is the most appropriate approach: 1) NNZ(A) > 300000, and 2) RowA > 42. Table 2 summarizes the various criteria. We note that tSparse performs better on RTX 2070 than on Titan RTX, for the reasons we describe in Section 6.2. Therefore, we can relax the criterion for non-high-end GPUs as follows: RowA > 21. Finally, we note that the first criterion may limit the suitability of our approach for applications that reduce the  Table 2: The criteria that define when tSparse is faster in comparison to other approaches. This table evaluates the performance of each condition using the precision and recall metrics from the classification theory. We note that the collection size varies to account for the incorrect results of each approach.
6.1.1. Evaluating the performance of the criteria We evaluate the speedup of tSparse over the other spGEMM approaches after applying our global criterion to the matrix collection. Table 3 shows the speedup of tSparse over the other spGEMM approaches when applying the global criterion.
We make three major observations. First, tSparse is faster than all other approaches by 1.46× − 36.93×, which confirms that our global criterion is reliable. Second, the speedup of all approaches follows the same trends for both GPUs. The RTX 2070 GPU performs better w.r.t. execution time, with more relaxed criterion, as the CPU does not bottleneck the execution time as much as on Titan RTX (Sections 6.3, 6.2). Third, the Minimum columns of Table 3 show speedups < 1. This happens because of the wrong predictions of our proposed criterion (false positives).

Speedup
To show the benefits of our approach that uses TCUs for MM, we find the speedup over cuSPARSE, CUSP, RMerge2, Nsparse, AC-SpGEMM, spECK and nonTCU for the 16 matrices of our dataset. Figure 7 and Fig. 8 show the speedup of tSparse over the seven approaches when calculating A * A on the Titan RTX GPU. Figure 7 corresponds to the first part of our dataset (randomly selected matrices), whereas Fig. 8 corresponds to the second part of our dataset (matrices selected based on criteria). Figures 9 and 10 present the speedup on RTX 2070. GeoMean indicates the geometric average of the speedup for the matrices of the respective figure.
We make five major observations. First, tSparse that uses TCUs performs faster than our nonTCU implementation, by an average of 1.68×. This speedup may seem low considering that TCUs in mixed precision promise 4× more flops than normal fp32 operations [25]. There are two reasons that keep it low: 1) our counting and multiplication kernels, without TCUs, occupy about 50% of the total execution time of spGEMM, so according to Amdahl's law we do not expect more than 2× speedup of the total execution time, and 2) tSparse is memory bound, rather than arithmetic bound, because the sorting step results in a task list with entries that point to non-continuous memory locations. Fourth, although the spGEMM approaches generally maintain their relative ranking, the speedup of tSparse over the other approaches is greater on the RTX 2070 system. The reason is that, although the GPU parts become significantly faster on Titan RTX, the CPU tasks of tSparse (e.g., memory allocation/deallocation) take the same time on both systems. tSparse scales slightly worse than the other approaches on the higherend GPU due to the CPU bound tasks.
Fifth, tSparse shows better performance with denser tiles. Grouping NZ elements in tiles reduces the number of values to work on. According to [5], the main cost of ESC is sorting. Consequently, fewer values equals to less time sorting. Generally, tSparse has good performance for bitmap density greater than five.

Execution time analysis
In this section, we present the relative execution time of the main parts of tSparse, i.e., Task      We note that in the execution time we include all allocations, including the allocation of the output. We also include the time for conversion from fp32 to fp16 (if the input is stored in fp32).

Numerical precision
TCUs accept 16-bit numbers as input, whereas TCUs perform MM and accumulation in 32-bit precision (mixed precision -Section 2.6). In this section, we show how tSparse fares w.r.t. precision in comparison to a full 32-bit approach (we use CUSP as base). Figure 12 illustrates the precision of spGEMM as Symmetric Mean Absolute Percentage Error (SMAPE -Equation (3)).
We make two observations. First, when the input matrices contain patterns of "1" and "0" the SMAPE is 0%. Second, when the inputs are real numbers the SMAPE is on average 0.02%. Exception are the cases with inputs close to the limits of the fp16 range, where the round-off error is bigger (nemeth21 -2.49%, patents main -0.66%).

Memory requirements
In this section we present the peak memory consumption of all approaches. Fig. 13 presents the memory consumption for the 16 matrices of our dataset.
We make four major observations. First, hash table approaches require the smallest memory area as they do not have to store huge intermediate matrices. Second, ESC approaches require large amounts of memory. Third, tSparse, which "packs" elements into tiles, typically requires more memory than hash approaches and less than ESC. The amount of memory is proportional to the bitmap density. When the average bitmap density is very low, like in patents main and wiki-vote, tSparse requires a lot of memory due to the additional overhead for storing tiles. Fourth, RMerge2 performs well as it usually requires only 5 bytes for each row of the left-hand side.  We examine the performance of our approach in A * A and compare to cuSPARSE, CUSP, RMerge2, Nsparse, AC-SpGEMM and spECK over the 16 sparse matrices of our dataset. We draw two important conclusions. First, tSparse generally outperforms the other approaches for the larger matrices of our dataset when we have dense tiles and enough work per matrix row. Second, TCUs play an important role in the performance of our approach, speeding up nonTCU 1.68×. tSparse outperforms the other state-of-the-art approaches 1.53× to 31.77× when our selected criteria are met, without significant loss in accuracy. Therefore tSparse is a suitable alternative for spGEMM.

Related work
To our knowledge, this is the first paper to use TCUs in the context of spGEMM. With this approach, we group NZ into tiles. Unlike previous methods, which use normal fp32 cores, we use tensor cores to multiply the tiles. We show that, with this approach, we can increase the performance of spGEMM. In this section we summarize other spGEMM implementations.

Expansion Sorting Compression
There exists a substantial body of work on improving ESC methodology that we describe in Section 2.4. Dalton et al. [5] improve on ESC by optimizing sorting, the most time consuming step, and by localizing processing to shared memory. Kunchum et al. [28], in HybridSparse, implement variants of ESC. Winter et al. [7], in AC-SpGEMM, perform ESC locally in shared memory. They use dynamic scheduling of iterations of ESC to keep data longer in shared memory. Thus, they reduce global memory traffic and the cost of sorting a huge intermediate matrix in global memory.
tSparse moves expansion and compression steps after sorting. This way we reduce memory allocation and movements and utilize MAC of TCUs.

Hash tables
Hash tables can mitigate the cost of sorting and storing huge intermediate matrices. Demouth et al. [19] present one of the first implementations of spGEMM with hash tables. cuSPARSE [11] is based on the work of Demouth et al. Their approach has two drawbacks. First, there is imbalance between threads because different threads of the warp might have to insert a different number of values in the hash table. Second, shared memory space is limited, which results in frequent data movement to global memory. Anh et al. [27], in BalancedHash, and Nagasaka et al. [20], in Nsparse, reduce the consumption of shared memory by partitioning the rows of the input or output, respectively. Nsparse improves on BalancedHash in two ways: 1) by using hash tables of variable size in shared memory, less shared memory is required and more thread blocks can run, and 2) by using fewer auxiliary matrices, it keeps memory traffic low and reduces memory storage requirements. Deveci et al. [30] use two-level hash tables. They adapt the hash tables to the number of threads in order to create a method that is portable to many platforms. Parger et al. [21], in spECK find a trade-off between analysis cost for load balancing and expected gain.
In tSparse, 1) we have less values to sort because elements are grouped into tiles, and 2) we avoid storing the intermediate products by using the task list to directly accumulate them.

Hybrid
Hybrid methods select among multiple methods/kernels during spGEMM depending on the workload. Dalton et al. [5] improve on ESC methodology by changing the thread granularity of the sorting method based on the size of rows of C. Liu et al.
[6], in bhSparse, change the sorting-merging method depending on the size of rows of C. Kunchum et al. [28], in HybridSparse, choose the spGEMM method based on the workload of each row of A. They select among variants of ESC and their own method (a GPU implementation of scatter vectors). Gremse et al. [8,9], in RMerge create different kernels for different row sizes of A. Hybrid methods achieve good load balance thanks to their adaptability to the workload.
We can divide tSparse in two parts. The part that performs ESC on tiles to form the task list and the part that uses our kernels to perform MAC operations. For the first part, we use Thrust which achieves good load balance. For the second part, we let the GPU hardware scheduler manage the workload. This is possible because we employ many thread blocks, which the GPU schedules to SMs based on the available resources of each SM.

Conclusion
In this work, we utilize TCUs to increase the performance of spGEMM. To that end, we modify the ESC method and we create a task list of MMs of tiles. The key advantages of our approach, tSparse, are two. First, tiles reduce the number of values that the computationally demanding parts of ESC have to act on. Second, the task list sends the tiles to TCUs, which not only perform MM faster than normal computation cores, but also leave the normal cores free for different workloads.
The results confirm that TCUs increase the performance of MM and the combination of our tiling approach with TCUs provides significant benefits to spGEMM. TCUs increase the performance of tSparse by 68% in comparison to our nonTCU implementation. Our approach is, on average, 1.53× to 31.77× faster than cuSPARSE, CUSP, RMerge2, Nsparse, AC-SpGEMM and spECK when NNZ(A) > 300000 and RowA > 42. We conclude that our methodology improves the performance of spGEMM by making efficient use of tiles and TCUs. The source code of our approach is available at https://github.com/oresths/tSparse.