A Systematic Survey of General Sparse Matrix-Matrix Multiplication

In this section, we first give the symbolic representation commonly used in this paper, which is presented in Table 1. We use bold capital italic for matrices, lowercase bold italic for vectors, and lowercase italic for scalars.

There is no strict and formal definition of sparse matrix. The most popular one is given by Wilkinson: sparse matrix is any matrix with enough zeros that it pays to take advantage of them (Gilbert et al., 1992). Another common quantitative definition is given by Barbieri et al. (Filippone et al., 2017), that is, a matrix $\bm{A}$ is sparse if its number of non-zero entries (referred to as NNZ) is $O(n)$.

Storing matrices with a dense pattern often leads to a lot of useless calculations and redundant storage, because they usually have a few non-zero entries (referred to as non-zeros). The prevailing solution is to store each sparse matrix with a compression format. Table 2 summarizes the contributions using different formats.

COO, CSR, CSC, ELL, and DIA are five basic and popular compression formats. COO is the plainest format and stores the row index, column index, and value of each non-zero entry in three separate arrays. CSR is the most extensively used format in existing work. Instead of storing row indices, CSR stores the row pointers to the first non-zero entry per row. CSC replaces the column indices array of COO with column pointers. ELL compacts all non-zeros to the left side. DIA is specifically designed for diagonal sparse matrices. It stores non-zeros in each diagonal and the offset of each diagonal from the main diagonal.

In addition to the above five basic formats, some new sparse formats have been proposed over the past years. Buluç et al. (Buluç and Gilbert, 2008b)(Buluç and Gilbert, 2008a) propose double compressed sparse column (DCSC), an improved format based on CSC. It is designed for hypersparse matrix by removing all the repetitions in column pointers array. They also present DCSR format (Buluç and Gilbert, 2008b), which is a row-based dialect of DCSC. Bor$\check{s}$tnik et al. (Borstnik et al., 2014) design an efficient distributed sparse matrix multiplication algorithm using blocked CSR (BCSR) format. Park et al. (Park et al., 2022) propose huffman-coded non-zero indication (HNI) format, which is a bitmap-based data encoding. It uses non-zero indication bit-stream to replace row and column indices and encodes the stream with Huffman coding. Xie et al. (Xie and Liang, 2019) design compressed feature map (CFM) for efficiently storing sparse feature maps in convolution neural network (CNN). In SparTen proposed by Gondimalla et al. (Gondimalla et al., 2019), a sparse tensor is encoded into a two tuple of a bit-mask representation and a set of non-zeros. The bit-mask representation has 1’s for positions with non-zeros and 0’s otherwise. This bit-mask based compression idea is also used in (Kanellopoulos et al., 2019)(Pentecost et al., 2019)(Zachariadis et al., 2020)(Qin et al., 2020). Han et al. (Han et al., 2016)(Chen et al., 2017) use a CSR variation, run-length encoding (RLE), to store a sparse weight matrix in deep neural network (DNN). The new format stores a vector $\bm{v}$ and an equal-size vector $\bm{z}$ for each column in the weight matrix, where $\bm{v}$ saves the non-zero weights, and $\bm{z}$ stores the number of zeros before the corresponding entry in $\bm{v}$. Soltaniyeh et al. (Soltaniyeh et al., 2020) propose REAP intermediate Representation (RIP) to increase the throughput on FPGAs. It has three parts: shared feature, metadata and distinct features. To overcome the inefficient memory access of CSR, Srivastava et al. (Srivastava et al., 2020) propose a channel cyclic sparse row (C²SR) format, which assigns each matrix row to a fixed channel in a cyclic manner. In TileSpGEMM, Yu et al. (Niu et al., 2022) propose a sparse tile data structure which describes a sparse matrix using two levels of tile information.

Breadth-first search (BFS) is a key and fundamental subroutine in many graph analysis algorithms. The goal of the BFS is to traverse a graph from a given source vertex, which can be performed by a sparse matrix-vector multiplication (SpMV) between the adjacency matrix $\bm{A}$ of a graph $G=(V,E)$ and a sparse vector representing the source vertex (Zhang et al., 2011). Assume $n=|V|$, then the size of $\bm{A}$ is $n\times n$. Let $\bm{x}$ be a sparse vector with $\bm{x}_{i}=1$ and all other entries being zero, then the 1-hop vertexes from source vertex $i$, denoted as $\bm{v}_{i}^{1}$, can be derived by the SpMV operation: $\bm{v}_{i}^{1}=\bm{A}\times\bm{x}$. Repeating the operation from $\bm{v}_{i}^{1}$, we can receive the 2-hop vertexes from $i$. Finally, a complete BFS for the graph $G$ from vertex $i$ is yielded (Gilbert et al., 2006). In contrast, Multi-Source BFS (MS-BFS) runs multiple independent BFSs concurrently on the same graph from multiple source vertexes, which can be formulated as SpGEMM. A simple example for MS-BFS is presented in Figure 2. Let $\{1,2\}$ be two source vertexes, $\bm{A}$ is the adjacency matrix of the graph, and $\bm{X}=(\bm{x}^{1},\bm{x}^{2})$ is a rectangular matrix representing the source vertexes, where $\bm{x}^{1}$ and $\bm{x}^{2}$ are two sparse column vectors with $\bm{x}^{1}_{1}=1$ and $\bm{x}^{2}_{2}=1$ respectively and all other entries being zero. Then, the sparse matrix $\bm{B}^{1}$ representing the 1-hop vertexes (denoted as $\bm{v}^{1}$) from source vertexes $\{1,2\}$ is give by $\bm{B}^{1}=\bm{A}\times\bm{X}$. Repeat the multiplication of adjacency matrix and $\bm{B}^{1}$: $\bm{B}^{2}=\bm{A}\times\bm{B}^{1}$, we can derive the 2-hop vertexes from $\{1,2\}$. Finally, we get the results of BFS from vertices 1 and 2, which are $\{1,3,4,2,5,6\}$ and $\{2,3,4,1,5,6\}$ respectively.

Many applications run hundreds of BFSs over the same graph. Compared with running BFSs sequentially, running multiple BFSs concurrently in a single kernel allows us to share the computation between different BFSs without paying the synchronization cost (Gilbert et al., 2006). As one of the most expensive operations of MS-BFS, SpGEMM is worthy of more attention to be paid.

Clustering is one of the unsupervised learning methods for statistical data analysis, and MCL is one of the graph clustering algorithms proposed for biological data (Selvitopi et al., 2020b). Using MCL, the closely connected points are grouped into clusters by performing random walks on a graph based on Markov chains. Let $\bm{A}$ denote a probability matrix, and each entry in the matrix is the probability of each point reaching others. The sum of each column in this matrix is 1. There are three steps in the clustering iteration. In the first step of expansion, the probability matrix of reaching other points starting from any point after one step is given by $\bm{B}=\bm{A}\times\bm{A}$. This not only strengthens the connection between different areas, but also leads to the convergence of probabilities. In the second step, all entries smaller than a given threshold are pruned. Then, the third step of inflation is required to weaken the possibility of loosely connected points by computing the power of each element in the probability matrix. Next, the matrix is replaced with the new matrix and used for next iterations. After a number of iterations, the points in a graph are gradually clustered into groups.

Due to the high time and memory overhead of MCL, high performance MCL (HipMCL) algorithm is proposed for fast clustering of large-scale networks on distributed platforms (Selvitopi et al., 2020b).

Algebraic Multigrid (AMG) is a multi-grid method developed based on Geometric multigrid (GMG). AMG iteratively solves large and sparse linear system $\bm{A}\bm{x}=\bm{b}$ by automatically constructing a hierarchy of grids and inter-grid transfer operators (Falgout, 2006)(Briggs et al., 2000). Generally, AMG includes two processes: setup and solve. The setup phase constructs multiple components of multigrid algorithm. The solve phase executes multigrid cycling based on these components, and SpMV dominates this phase. SpGEMM is an important kernel in setup phase, whose critical steps are presented in Algorithm 1. It first constructs interpolation operator $\bm{P}_{l}$ based on input matrix $\bm{A}$ (line 3), then the restriction operator $\bm{R}_{l}$ is the transposition of $\bm{P}_{l}$ (line 4). Finally, the coarse-grid system $\bm{A}_{l+1}$ is constructed using Galerkin product (line 5), which is implemented with two SpGEMMs (Bell et al., 2012)(Ballard et al., 2016b). Generally, $\bm{P}_{l}$ is tall and skinny, and $\bm{R}_{l}$ is short and fat. Their non-zeros’ distribution is closely related to the used interpolation algorithms.

These SpGEMMs, taking more than $80\%$ of the total construction time, are the most expensive components. Moreover, the construction of operators (thus SpGEMM) is an expensive part in overall execution since it may occur at every time step (for transient problems) or even multiple times per time step (for non-linear problems), making it important to optimize SpGEMM (Elliott and Siefert, 2018).

SpGEMM is also one of the most important components for genome assembly (Selvitopi et al., 2020a)(Guidi et al., 2021a)(Guidi et al., 2021b), NoSQL database (Cohen et al., 2009)(Gadepally et al., 2015)(Hutchison et al., 2015), triangle counting (Cohen, 2009)(Wolf et al., 2015)(Wolf et al., 2017), graph contraction (Gilbert et al., 2008), graph coloring (Kaplan et al., 2006)(Deveci et al., 2016), the all pairs shortest path (Chan, 2007), sub-graph (Vassilevska et al., 2006)(Buluç and Madduri, 2011), cycle detection or counting (Censor-Hillel et al., 2015), and molecular dynamics (Weber et al., 2015)(Akbudak and Aykanat, 2014).

SpGEMM algorithms using precise prediction usually consist of two phases: symbolic and numeric phases. In the symbolic phase, the precise NNZ in each row/column of the output matrix is computed based on row and column indices of input sparse matrices, an example is shown in Figure 5. Real values of non-zeros are calculated in the numeric phase.

The implementations of SpGEMM in Kokkos Kernels (Rajamanickam et al., 2021), cuSPARSE (NVIDIA, 2021), MKL (Intel, 2021), and RMerge (Gremse et al., 2015) are typical representatives of this method. Besides, existing work (Demouth, 2012), (Nagasaka et al., 2017), (Demirci and Aykanat, 2020b) and (Akbudak and Aykanat, 2014) also exploit this approach. In order to speed up the symbolic phase, Deveci et al. (Deveci et al., 2017)(Deveci et al., 2018b) design a graph compression technique to compress the matrix $\bm{B}$ by packing its columns as bits. In (Gu et al., 2020), the authors estimate the memory requirement for $\bm{C}$ as well as the number of bins and allocate space for global bins in the symbolic phase. SpECK (Parger et al., 2020) uses size information collected in symbolic execution to guide the selection of accumulators in numeric phase.

Cohen (Cohen, 1997) transforms the problem of estimating NNZ in the output matrix to the size estimation of reachability sets in a directed graph, and presents a Monte Carlo-based algorithm to estimate the size of reachability sets. The algorithm can be demonstrated using a hierarchical structure graph of matrix product, as shown in Figure 6. The algorithm starts by assigning a vector of same size, initialized with exponential-distribution random samples, to each node in $L_{1}$. The vector of each node in the higher layer is equal to the column-wise minimum of the vectors of its neighbors in the lower layer. Finally, NNZ in each row of $\bm{C}$ is estimated based on the vector of each node in $L_{3}$. For any tolerated relative error $\epsilon>0$, it can compute an $1\pm\epsilon$ approximation of NNZ in result matrix in time $\mathcal{O}(n/\epsilon^{2})$. Then, Amossen et al. (Amossen et al., 2014) improve this method to expected time $\mathcal{O}(n)$ for some particular $\epsilon$. Anh et al. (Anh et al., 2016) utilize a similar size estimation technology based on the rows sampling of $\bm{A}$ and $\bm{B}$. In (Cohen, 1997), the authors introduce an algorithm to determine the multiplication order of chain products with minimal number of calculation operations. Paper (Selvitopi et al., 2019) also uses this three-layer graph representation to estimate the size of result matrix.

The third method computes an upper-bound NNZ in the output matrix and allocates corresponding memory space. The most commonly used method is to count NNZ in the corresponding rows of $\bm{B}$ for each non-zero entry in $\bm{A}$. Taking the matrices presented in Figure 5 for example, the upper bound of NNZ in each row of the output matrix is stored in the array $U=\{u_{1},u_{2},u_{3},u_{4}\}$. The first row of $\bm{A}$ has two non-zeros, whose column indices are 2 and 4. Therefore, the upper-bound NNZ in the first row of $\bm{C}$ equals to the sum of NNZ in the second and forth rows of $\bm{B}$. The ESC algorithm (Bell et al., 2012)(Dalton et al., 2014) proposed by Bell et al. is a representative of this method. Nagasaka et al. (Nagasaka et al., 2019) also count a maximum of scalar non-zero multiplications per row of the result matrix. Then each thread allocates a hash table based on the maximum and reuses the hash table throughout the computation by re-initializing at the beginning of calculating each row.

The fourth method, also known as the progressive method, dynamically allocates memory as needed. It first allocates memory of proper size and then starts matrix multiplication. A larger memory block is re-allocated if the current memory is insufficient. The implementation of SpGEMM in Matlab (Gilbert et al., 1992) is a representative of this method. It first guesses the size, and then allocates a memory block that is larger by a constant factor (typically 1.5) than the current space if more space is required at some point.

Liu and Vinter (Liu and Vinter, 2014)(Liu and Vinter, 2015) propose a hybrid method which calculates the upper-bound NNZ for each row and groups all rows into multiple bins according to NNZ. The method allocates space of the upper-bound size for short rows and progressively allocates space for long rows. In TileSpGEMM proposed by Yu et al. (Niu et al., 2022), it first calls the symbolic implementation in NSparse (Nagasaka et al., 2017) to get the sparse tile structure of $\bm{C}$. Then, it uses binary search to find intersection sparse tiles from $\bm{A}$ and $\bm{B}$. Finally, bit mask operations are used to calculate the number of non-zeros of each tile in $\bm{C}$.

Of the four methods, precise method not only saves the memory usage, but also enables the sparse structure of $\bm{C}$ to be reused for different multiplies with the same structure of input matrices (Deveci et al., 2018b)(Gustavson, 1978). Moreover, it presents significant benefits in graph analytics, because most of them work only on the symbolic structure, no numeric phase (Wolf et al., 2017). However, the calculation of two phases means that it needs to iterate through the input matrices twice, leading to higher computation overhead than other methods. The accuracy and the speed of the second method depend on the probabilistic algorithm used, and additional memory allocation must be launched when the estimate fails. The upper-bound method is efficient and easy to implement, but it usually leads to memory over-allocation. The progressive method allocates memory dynamically as needed, and additional memory allocation must also be launched when the first allocation fails. In practice, the choice needs to be made according to the discussed problems.

Block partition of SpGEMM can be categorized into 1D, 2D and 3D algorithms based on how they partition the work among computing units (Azad et al., 2016)(Ballard et al., 2016a). We first introduce the workcube notation before diving into the details, as shown in Figure 8(b). The front, top, and side views of the workcube represent the two input matrices and the output matrix presented in Figure 8(a), respectively. The workcube can be divided into $3\times 3\times 3$ voxels. Each voxel inside represents the scalar multiplication of two non-zeros, which are mapped to cells of the voxel in the front and top views, and contributes to its projection in the side view (Ballard et al., 2013). Based on this projection, task partition of SpGEMM can be seen as the partition of workcube in different dimensions. Next, we introduce 1D, 2D and 3D partition based on this notation.

1D Partition. 1D partition only divides the workcube in one of the three dimensions. As shown in Figure 9, each sub-figure represents a variant of 1D partition. A ”layer” of the workcube is assigned to one computing unit in all variations. Figure 9(a) divides the workcube by planes parallel to the top view, referred to as $T$-partition. Each computing unit is responsible for the multiplication of a set of rows in $\bm{A}$ and the entire matrix $\bm{B}$, producing the corresponding rows of $\bm{C}$. Similarly, Figures 9(b) and (c) divide the workcube by planes parallel to the front view and side view, respectively. We refer to these two partition variants as $F$-partition and $S$-partition.

$T$-partition. Liu et al. (Liu and Vinter, 2014) group rows of $\bm{C}$ to different bins according to their upper-bound NNZ to maintain a balanced load. Nagasaka et al. (Nagasaka et al., 2017) partition rows of $\bm{C}$ according to precise NNZ in each row in numeric phase. Kurt et al. (Kurt et al., 2017) propose to assign one thread block to a fixed number of rows in $\bm{A}$ and assign a same number of threads in the block to complete the row-wise multiplication of each row in $\bm{A}$ and corresponding entries in $\bm{B}$. Deveci et al. (Deveci et al., 2017) propose a hierarchical and parallel SpGEMM algorithm based on Kokkos library (Edwards et al., 2014). At the first level, one kokkos-team is assigned to calculate a set of rows of $\bm{C}$. Each kokkos-thread within the team is responsible to produce a subset of these rows at the second level. Multiple vector-lanes in one kokkos-thread are assigned to perform the row-wise multiplication of each non-zero entry in the subset of $\bm{A}$ and corresponding rows in $\bm{B}$ at the third level. They also explore $T$-partition but with different task assigning schemes for HPC architectures (Deveci et al., 2018b). Winter et al. (Winter et al., 2019) present a partition scheme assigning the same NNZ of $\bm{A}$ to each block while ignoring the row boundaries. Instead of splitting rows evenly (Elliott and Siefert, 2018), Li et al. (Li et al., 2019) split the matrix into multiple row blocks based on NNZ. To address the input or output reuse problem, Zhang et al. (Zhang et al., 2021) split matrix $\bm{A}$ into row fibers and dispatch them to processing elements (PEs) in a SpGEMM accelerator GAMMA, and each PE then performs a linear combination of row fibers of $\bm{B}$ to produce a row fiber of $\bm{C}$. Shivdikar et al. (Shivdikar, 2021) group multiple rows in a single window and assign one PIUMA (Programmable Integrated Unified Memory Architecture) block to a window. The size of a window depends on the scratchpad size.

$F$-partition. Azad et al. (Azad et al., 2015) assign processors to process columns of the upper triangular matrix in triangle counting. Lin et al. (Lin et al., 2013) design an architecture for SpGEMM on FPGAs, the computation is partitioned evenly to all PEs, and each PE is assigned to calculate multiple columns of matrix $\bm{C}$.

$S$-partition. Based on the outer-product multiplication, Deveci et al. (Deveci et al., 2018a) divide rows of $\bm{B}$ into blocks so that each block can be fitted into the HBM of KNL. This partition also induces a row-wise partition of $\bm{A}$. Zhang et al. (Zhang et al., 2020) optimize the outer-product SpGEMM by compacting non-zeros of each row in $\bm{A}$ to the left so that the data locality for both input and output matrices are jointly optimized. Gu et al. (Gu et al., 2020) develop an improved ESC SpGEMM based on outer product. They group the expanded triples into bins to saturate memory bandwidth.

Others. Ballard et al. (Ballard et al., 2016b) present a theoretical and detailed analysis for the communication cost of Galerkin triple product in the smoothed aggregation of AMG methods. They conclude that the row-by-row product ($T$-partition) is the best 1D method for the first two multiplications, and the outer product ($S$-partition) is the best 1D method for the third multiplication.

Most of SpGEMM can be implemented using 1D partition. However, when the NNZ in columns or rows of the matrix $\bm{A}$ or $\bm{B}$ increases, the communication overhead becomes substantially huge in distributed computing. This indicates that a ”layer” in 1D partition is required to be processed by more computing units. Therefore, the 2D algorithms with fine-grained partition are proposed.

2D Partition. 2D algorithms divide the workcube in two of the three dimensions, so there are also three variants. Figure 10(a) illustrates three partition variants that divide the workcube by planes parallel to two of the three views (top and front views, front and side views, top and side views), referred to as $TF$-partition, $FS$-partition and $TS$-partition hereafter. In $TF$-partition, each computing unit computes a block of $\bm{C}$. However, in $TS$-partition and $TF$-partition, each computing unit computes a block of intermediate results of $\bm{C}$.

$TF$-partition. SpSUMMA was first proposed by Buluç et al. (Buluç and Gilbert, 2008b)(Buluç and Gilbert, 2008a) based on dense SUMMA algorithm (van de Geijn and Watts, 1997). They also present a comparative SpCannon algorithm based on dense Cannon algorithm (Cannon, 1969). Both algorithms logically organize processors as a 2D grid, and map a block of $\bm{A}$ and $\bm{B}$ to each processor. In SpSUMMA, each processor broadcasts data to other processors in the same row/column of the grid, and also receives data broadcast by these processors. In SpCannon, processors exchange data using point-to-point communication. In both algorithms, each processor accesses a row block of $\bm{A}$ and column block of $\bm{B}$, and calculates a result block of $\bm{C}$, which conforms with $TF$-partition. Matrix partition similar to SUMMA is also used by Jin et al. (Jin and Ziavras, 2004). Patwary et al. (Patwary et al., 2015) partition $\bm{A}$ by row, while partition $\bm{B}$ by column when a certain condition is satisfied.

$TS$-partition. Deveci et al. (Deveci et al., 2018a) partition $\bm{A}$ and $\bm{B}$ by row for GPU. When one partition of $\bm{A}$ cannot fit in GPU fast memory, column-wise partition is applied for row strip of $\bm{A}$.

2D algorithms perform reasonably well on a few hundred processes. However, as the number of processes increases, the communication cost becomes a bottleneck, and 3D algorithms were developed to reduce the communication cost.

3D Partition. 3D algorithms divide all three dimensions of the workcube. As shown in Figure 10(b), each computing unit owns a block of $\bm{A}$ and a block of $\bm{B}$ to compute an intermediate result block of the matrix $\bm{C}$.

Azad et al. (Azad et al., 2016) present a parallel implementation of the 3D SpGEMM algorithm. Hussain et al. (Hussain et al., 2021) use a similar 3D partition, while dividing each block of $\bm{B}$ into multiple batches to meet the size of available memory. To avoid excessive communications and logistic operations of 2D SpSUMMA, Weber et al. (Weber et al., 2015) propose a novel SpGEMM algorithm MPSM3 for the computation of the density matrix in electronic structure theory. All processes are organized as a 3D Cartesian topology, and each box in the workcube is assigned to a specific process on terms of the physical properties of the problem. Yu et al. (Niu et al., 2022) divide two input sparse matrices into a number of 16-by-16 tiles to utilize 8-bit unsigned char data type for local indices.

Traditional block-based matrix partitions rarely consider the workload (NNZ in each block) assigned to each processing element. For example, SpSUMMA presumes that non-zeros are uniformly and randomly distributed across the row/column, which may not hold for most sparse matrices.

Hypergraph Partition. The hypergraph representation for $\bm{C}=\bm{A}\times\bm{B}$, denoted as $\mathcal{H}=(\mathcal{V},\mathcal{N})=(\mathcal{V}^{\bm{AB}}\cup\mathcal{V}^{\bm{C}},\mathcal{N})$, is defined as a set of vertices $\mathcal{V}$ and a set of nets (hyperedges) $\mathcal{N}$. Here we take the outer-product SpGEMM as an example to introduce the hypergraph based partition. In $\mathcal{H}=(\mathcal{V}^{\bm{A}\bm{B}}\cup\mathcal{V}^{\bm{C}},\mathcal{N})$, each vertex $v_{i}$ in $\mathcal{V}^{\bm{AB}}$ denotes the outer product of the column $\bm{a}_{*i}$ of $\bm{A}$ with the row $\bm{b}_{i*}$ of $\bm{B}$. $\mathcal{V}^{\bm{C}}$ contains a vertex $v_{ij}$ for each non-zero entry $\bm{c}_{ij}$ in $\bm{C}$. $\mathcal{N}$ contains a net $n_{ij}$ for each non-zero entry $\bm{c}_{ij}$ in $\bm{C}$. Figure 11(a) shows an example of SpGEMM and Figure 11 (b) shows its hypergraph representation. The task of hypergraph partitioning is to divide a hypergraph into two or more roughly equal-sized parts such that a cost function on the nets connecting vertices in different parts is minimized (Çatalyürek et al., 2011).

Bipartite Graph Partition. A bipartite graph for $\bm{C}=\bm{A}\times\bm{B}$, denoted as $\mathcal{G}=(\mathcal{V}^{\bm{A}\bm{B}}\cup\mathcal{N}^{\bm{C}},\mathcal{E})$, is defined as two disjoint sets of vertices $\mathcal{V}^{\bm{A}\bm{B}}$ and $\mathcal{V}^{\bm{C}}$, and a set of edges $\mathcal{E}$. Also take outer-product SpGEMM as an example. The semantics of each vertex in $\mathcal{G}$ is the same with those in $\mathcal{H}$. The difference is that the dependence of $\bm{c}_{ij}$ is captured with edges, instead of a net. If an outer product represented by $v_{x}$ produces a partial result for $\bm{c}_{ij}$, an edge connecting vertices $v_{x}$ and $v_{ij}$ is added in $\mathcal{E}$. An example of a bipartite graph for the outer-product SpGEMM is presented in Figure 11(c).

Hypergraph-based partition was first applied in outer-product SpGEMM by Akbudak et al. (Akbudak and Aykanat, 2014). Besides, they also present two extended hypergraph model, which partitions the matrix $\bm{C}$ by row and column, respectively. Ballard et al. (Ballard et al., 2016a) present a fine-grained hypergraph model, in which each vertex is either a scalar multiplication or a non-zero entry. Based on this hypergraph model, Kurt et al. (Kurt et al., 2017) use an iterative method to construct hypergraph partitions to ensure that the data entries involved in each partition does not exceed the cache capacity. Akbudak et al. (Akbudak and Aykanat, 2017) use hypergraph model to improve outer-product and inner-product SpGEMM.

Akbudak et al. (Akbudak et al., 2018) propose and compare multiple computational partitioning models based on hypergraph and bipartite graph partitioning, with the aim of reducing message volume. Further, three communication hypergraph partitioning models for three SpGEMM formulations are proposed to reduce the latency cost. The bipartite graph model for row-by-row SpGEMM was later used by Demirci et al. (Demirci and Aykanat, 2020b). Instead of associating two weights with each vertex, they propose a three-constrained partitioning in which each vertex is associated with three weights. Selvitopi et al. (Selvitopi et al., 2019) apply bipartite graph and hypergraph models for simultaneous scheduling of the map and reduce tasks for MapReduce jobs. Demirci et al. (Demirci and Aykanat, 2020a) propose hypergraph models for 2D and 3D partitioning to improve the performance of 2D and 3D SpGEMM.

Each partition algorithm usually has its own design considerations and advantages. As the most commonly used sparse formats, CSR and CSC usually store sparse matrices by rows and columns, respectively. It is compatible with the 1D partition. The workload of 1D partition is closely related to the non-zeros distributions in the divided rows or columns, and the non-zeros distributions in the other sparse matrix. Taking $T$-partition for example, we assume that each thread or a group of threads are responsible for the multiplication of one row of $\bm{A}$ and entire $\bm{B}$. If $\bm{B}$ is a regular sparse matrix and has uniform non-zeros distribution, then the NNZ in each row of $\bm{A}$ determines the workload of each computing unit. Therefore, the longest row of $\bm{A}$ dominates the performance. The problem of load imbalance is more serious and complicated if $\bm{B}$ is irregular. 2D and 3D partitions are fine-grained and produce a more balanced workload than the 1D partition. Communication cost is an important metric for the distributed system. In the 1D partition, each work process accesses an entire input matrix ($\bm{B}$ for $T$-partition, $\bm{A}$ for $F$-partition) or output matrix ($\bm{C}$ for $S$-partition), which requires a huge communication overhead. Differently, 2D and 3D partitions just access some rows/columns, or partial non-zeros. The graph partitioning algorithms build computation and communication graph models to ensure load balance and minimum communication cost. For SpGEMM with an irregular computation pattern, graph partitioning algorithms achieve a more balanced load and lower communication cost than block partitioning algorithms, but at the cost of higher partitioning overhead.

Dense accumulators use dense vectors to cache the intermediate results of current ”active” columns or rows in the result matrix. The most popular dense accumulator, was first proposed by Gustavson (Gustavson, 1978), usually has three vectors. The first one stores real values. The second one is used as a tag array to check if a column index has been inserted before. The third one stores column indexes. An example of the dense accumulator is presented in Figure 13.

Gustavson’s algorithm allows random access to a single entry in a specified time, and it is widely used and improved by many successive researches (Gilbert et al., 1992)(Patwary et al., 2015)(Siegel et al., 2010)(Lee et al., 2020). For example, the SPA accumulator in Matlab (Gilbert et al., 1992) uses a dense vector of the same size as the column size of $\bm{C}$, a dense vector with true/false ”occupied” flags, and an unordered list of the indices whose ”occupied” flags are true. Elliott et al. (Elliott and Siefert, 2018) note that prolongation matrices in the AMG method are usually tall and skinny, so they use a local dense lookup array for each thread that supports fast memory space allocation for result matrix by page-aligned and size-tracked techniques.

In sparse accumulator, intermediate results for a row are stored in compact data structures. There are two types of method for merging: list and hash-based methods.

List-based accumulating usually sorts the intermediate results by column indexes. According to the sorting algorithms utilized, we classify researches as follows. Radix sort based methods allocate the entries to be sorted to some ”buckets”, so as to achieve the goal of sorting. Bell et al. (Bell et al., 2012) propose the ESC algorithm, which sorts many scalar products in the second phase. It has been improved in (Liu and Vinter, 2015) and (Li et al., 2019). Dalton et al. (Dalton et al., 2015) use the B40C radix sort algorithm (Merrill and Grimshaw, 2011) which allows specifications in the number and location of the sorting bits to accelerate the sort operation. Gu et al. (Gu et al., 2020) use an in-place radix sort to group keys by a single byte sharing the same valid byte position. In (Liu et al., 2019), a sparse accumulator, which sorts data in the GPU registers, is proposed. Merge sort based methods first divide the sequence to be sorted into several subsequences. Then each subsequence is ordered, and all subsequences are combined into an overall ordered sequence. In (Gremse et al., 2015), the sparse rows of $\bm{B}$ selected and weighted by corresponding non-zeros of $\bm{A}$, are merged in a hierarchical way similar to merge sort. Besides, Liu et al. (Liu et al., 2019) propose a GPU register-based merge algorithm, which demonstrates significant speedup over its original implementation. Heap sort based accumulator is proposed by Azad et al. (Azad et al., 2016). The intermediate result is represented as a list of triples, and each triple includes row index, column index and value of each non-zero entry. A $k$-way merge on $k$ lists of triples is performed by maintaining a heap of size $k$, which stores the current minimum entry in each triple list. Then a multi-way merge routine finds the minimum triple from the heap and merges it into the current result. Nagasaka et al. (Nagasaka et al., 2019) implement a heap-sort based sparse accumulator on multicore architectures.

Hash-based. Hash-based sparse accumulator first allocates a memory space based on the upper bound estimation as the hash table and uses the column indexes of the intermediate results as the key. Then the hash table is required to be shrunk to a dense state. Finally, we sort the values of each row of the result matrix according to their column indexes to obtain the final result matrix compressed with sparse format (Liu et al., 2019). An example is shown in Figure 14. SpGEMM implementation in cuSPARSE (Demouth, 2012)(NVIDIA, 2021) is a representative of this accumulator. Deveci et al. (Deveci et al., 2013) design a hashmap-based, two-level, sparse accumulator that supports parallel insertions and merges from multiple vector lanes. Weber et al. (Weber et al., 2015) use the hash table to store each row of the matrix $\bm{C}$ based on the distributed BCSR format. In (Anh et al., 2016), the row and column indices of an intermediate result, normally stored in two 32-bit numbers, are first compacted into a single 32-bit value and then inserted into hash table as key. It is helpful to accelerate the later sorting operation. Nagasaka et al. (Nagasaka et al., 2017) utilize vector registers in KNL and multi-core architectures for hash probing to accelerate hash-based SpGEMM. Their hash function is defined as the reminder after division of the column index multiplied by a constant number by hash table size. Liu et al. (Liu et al., 2019) optimize the data allocations of each thread in a hash-based accumulator utilizing GPU registers. Selvitopi et al. (Selvitopi et al., 2020b) compare the heap and hash sort in distributed HipMCL, and conclude that the latter method performs better.

Deveci et al. develop KKMEM (Deveci et al., 2017) and KKSpGEMM (Deveci et al., 2018b) algorithms, which support two hash-based accumulators and one dense accumulator. Which accumulator to use depends on the features of the discussed matrix. KKTri-Cilk algorithm (Yasar et al., 2018) uses a similar accumulator to that in the KKSpGEMM algorithm. Meanwhile, it is optimized in the Cilk implementation.

Sparse or dense accumulators are selected according to the distribution of non-zeros in different parts of a matrix. Liu et al. (Liu and Vinter, 2014)(Liu and Vinter, 2015) propose a hybrid parallel result accumulating method. Three sort-based accumulators are selected according to the upper-bound NNZ in each row. In the spECK algorithm (Parger et al., 2020), a local balancer decides which accumulator to use for each block to achieve the best expected performance from the following three accumulators: direct reference, hashing, or dense accumulation. The hash function of its hash-based accumulator multiplies the element index with a prime number and then divides the result by hash table size. TileSpGEMM (Niu et al., 2022) uses the sparse accumulator working on a sparse tile and the dense accumulator working on a dense tile. Whether a tile is dense or sparse is determined by a preseted threshold.

Compared with sparse accumulators, dense accumulator requires a large memory space, especially in high concurrent scenarios. They usually require to maintain a private dense accumulator for each thread or thread group, which greatly reduces the scalability of the algorithm. There are three cases in which dense accumulators are preferred: (1) in hybrid accumulators to process the rows with a large number of non-zeros per row; (2) in SpGEMM algorithms on CPU whose large cache and modest computing cores are preferred by dense accumulator; (3) in SpGEMM whose result matrix has small column size. In contrast, the sparse accumulators have lower memory requirement and are more suitable for SpGEMM whose result matrix has a large number of columns. At the same time, sparse accumulators are also the primary choice for GPU-based SpGEMM algorithm, because GPU has limited shared memory. Among sparse accumulators, hash-based accumulators have low memory space requirement, with the cost of handling collisions. For three list-based sparse accumulators, the performance of the used sorting algorithm determines their performance. Heap and merge sort has the same time complexity, while the latter presents a higher space complexity. Radix sort is preferred when processing large-scale array.

Patwary et al. (Patwary et al., 2015) propose to divide $\bm{A}$ and $\bm{B}$ by row and column, respectively, aiming to alleviate the problem of high L2 cache misses caused by dense accumulator. Elliott et al. (Elliott and Siefert, 2018) present a single-phase OpenMP variant of the Gustavson algorithm. They use page-aligned allocations and track used size and capacity separately to accelerate the allocation of thread-local memories. Gu et al. (Gu et al., 2020) use the propagation block technique to improve the memory bandwidth utilization of OP SpGEMM. They first save the multiple tuples generated by the outer product into multiple partially sorted bins, and then sort and merge each bin independently. Chen et al. (Chen et al., 2019)(Chen et al., 2020) design multiple SpGEMM kernels based on different sparse formats. They all utilize coalesced DMA transmission instead of discrete memory access to achieve fast data loading from main memory. Furthermore, they reserve a minimum partition of $\bm{B}$ in local data memory to improve data reuse and avoid redundant data loading.

In (Patwary et al., 2015), matrices are divided into small partitions and dynamic load balancing is used over the partitions. They find that better results are achieved when the total number of partitions is $6$-$10$ times the number of threads. Chen et al. (Chen et al., 2019)(Chen et al., 2020) present a three-level partitioning scheme for Sunway architecture. At the first level, each core group in the SW26010 processor performs the multiplication of a sub-$\bm{A}$ and $\bm{B}$. At the second level, each CPE core in a core group multiplies the sub-$\bm{A}$ by a sub-$\bm{B}$. At the third level, to fit the 64 KB local data memory in each CPE core, sub-$\bm{A}$ and sub-$\bm{B}$ are further partitioned into several sets. To achieve a balanced load, they divide $\bm{A}$ and $\bm{B}$ according to the computational loads rather than rows.

Patwary et al. (Patwary et al., 2015) store each column partition of $\bm{B}$ in an individual CSR format. This requires to change the data structure of $\bm{B}$ from CSR to blocked CSR in advance, which brings a significant format conversion overhead. Therefore, they use a simple upper-bound method to estimate the proportion of non-zeros per row that exceeds the size of the L2 cache, and format conversion occurs only when this proportion is higher than 30%.

In (Dalton et al., 2015), Dalton et al. propose an improved ESC method, which stores the non-zeros of $\bm{A}$ in shared memory to avoid repeated loading from global memory. In (Winter et al., 2019), Winter et al. propose an adaptive chunk-based GPU SpGEMM approach (AC-SpGEMM), which improves ESC algorithm by reducing intermediate results of $\bm{C}$ within shared memory. Gremse et al. (Gremse et al., 2015) also use shared memory to merge intermediate results. The main idea is to reduce the overhead of global memory accesses by merging rows using sub-warps. Then in (Nagasaka et al., 2017), Nagasaka et al. propose a fast SpGEMM algorithm that only requires a small amount of memory and achieves high performance. Since the memory subsystems of NVIDIA${}^{\text{\textregistered}}$ GPU differ from generation to generation significantly in both bandwidth and latency, two different data placement methods are proposed in (Deveci et al., 2018a). One is to keep the partitions of $\bm{A}$ and $\bm{C}$ in fast memory and stream the partitions of $\bm{B}$ to shared memory. The other is the opposite. Which method to choose often depends on the features of the input matrices. In (Liu and Vinter, 2014)(Liu and Vinter, 2015), memory pre-allocation for the result matrix is organized using a hybrid approach that saves a large amount of global memory space.

In (Demouth, 2012), one warp is assigned to one row of $\bm{A}$, and each thread in the warp is responsible for the multiplication of one non-zero entry in the row and the corresponding row of $\bm{B}$. To address the load imbalance between thread blocks, Lee et al. (Lee et al., 2020) propose Block Reorganizer based on the OP SpGEMM. The Block Reorganizer first gets the calculation load of all column-row pairs of input matrices, then combines, splits and assigns these pairs to thread blocks. Kurt et al. (Kurt et al., 2017) also use one thread block to process multiple rows of $\bm{A}$ simultaneously. Yu et al. (Niu et al., 2022) store input matrices and output matrix as multiple non-empty tiles and assign one warp to process one sparse tile of $\bm{C}$ in both symbolic and numeric phases on GPU, which greatly alleviates the load imbalance caused by irregular non-zero distributions of matrix rows.

To achieve a better load balancing, some work uses thread block to process the same NNZ. In (Winter et al., 2019), Winter et al. split the non-zeros of $\bm{A}$ uniformly and assign each thread block to process the same NNZ. However, the workload for each block varies based on the intermediate products generated due to the irregular non-zero distribution in the rows of $\bm{B}$. They present a fine-grained load balancing strategy within thread block. In addition, some researches assign tasks by intermediate results. For example, the authors of (Liu and Vinter, 2015)(Liu and Vinter, 2014) group the rows of the result matrix to multiple bins according to upper-bound NNZ estimation of each row, and assign different computing units to each bin to have a better load balancing.

In recent years, some authors have suggested two-level load balancing. SpECK (Parger et al., 2020) takes advance of both global load balance which splits the work into blocks, and local load balance which decides the number of threads assigned for each row of $\bm{B}$. Xia et al. (Xia et al., 2021) also adopt this load balance strategy in their own GPU implementation of SpGEMM. Also, SpGEMM in Kokkos Kernels (Rajamanickam et al., 2021) picks the best task assignment method according to non-zeros distribution of input matrices. Specifically, for the matrices with fewer non-zeros per row, it assigns one thread to a row. For the matrices with more non-zeros, multiple threads are assigned to one row, and vector parallelism is used. In (Dalton et al., 2015), SpGEMM is formulated as a layered graph model, and the expansion phase of ESC method is considered as the BFS in the levels of the layered mode. Then, they parallel the expansion phase over the multiplication of each non-zero entry of $\bm{A}$ and the corresponding row of $\bm{B}$ at the granularity of thread, warp, or thread block, according to the length of the row referenced from $\bm{B}$.

The emerging tensor core unit (TCU) in GPU attracts the attention of researchers. In (Zachariadis et al., 2020), Zachariadis et al. utilize TCUs to improve the performance of SpGEMM. They partition the input matrices into tiles and operate only on tiles which contain one or more non-zero entries.

Liu et al. (Liu and Vinter, 2014)(Liu and Vinter, 2015) present efficient parallel merging methods for rows with a large number of intermediate results. It consists of four steps: binary search and duplicate entries reducing, prefix-sum scan, non-duplicate entries copy, and merging of two sorted sequence in one continuous memory space. Dalton et al. (Dalton et al., 2015) propose an optimization scheme for the sorting process of ESC method. They replace global memory-based sorting operations over a large number of intermediate results with multiple parallel shared memory-based sorting operations over a small number of ones. When the shared memory space is insufficient, the global memory-based sorting is launched. Liu et al. (Liu et al., 2019) utilize GPU registers to optimize three typical sparse accumulators: sort, merge, and hash.