MAGNUS: Generating Data Locality to Accelerate Sparse Matrix-Matrix Multiplication on CPUs

For a sparse matrix $X$, we define $n_{X}$, $m_{X}$, and $nnz_{X}$ as the number of rows, columns, and nonzero entries, respectively. The set $\mathcal{S}(X)$ denotes the column indices of all nonzero entries in $X$. Our notation also extends to individual rows. For example, for row $i$ of $X$, $nnz_{X_{i}}$ denotes the number of nonzero entries, and $\mathcal{S}(X_{i})$ denotes the set of column indices corresponding to the nonzero entries.

The general sparse matrix-matrix multiplication operation (SpGEMM) is defined as $C=AB$, where $A$, $B$, and $C$ are sparse matrices. Multithreaded implementations of SpGEMM are typically based on Gustavson’s row-by-row algorithm (gustavson, ):

i.e., for some row $i$ of $C$, the rows of $B$ are scaled by the nonzero values in row $i$ of $A$. These scaled rows are then summed together to give the final row of $C$. Since each row of $C$ is computed independently, multithreaded implementations typically partition rows among threads, which is the approach we take for MAGNUS.

There are two main ingredients for implementing Gustavson’s method: the matrix storage scheme and the algorithm that accumulates the scaled rows of $B$. Compressed sparse row (CSR) format is one of the most popular storage schemes and is especially useful for algorithms such as Gustavson’s method that traverse matrices row-wise. The CSR format requires three arrays: $C.col$, $C.val$, and $C.rowPtr$. Arrays $C.col$ and $C.val$ of size $nnz_{C}$ store the column indices and values, respectively, of the nonzero entries of $C$, and $C.rowPtr$, of size $n_{C}+1$, stores the starting positions of each row in $C.col$ and $C.val$.

Algorithm 1 shows the pseudocode for the numeric phase of Gustavson’s method. Variables in bold are global, meaning they are shared and visible across all threads. The scaled rows of $B$ are summed using a dense accumulator, defined as the combination of denseAccumBuff and bitMap. The column indices in row $i$ of $A$ are loaded as $idx=A.col[A.rowPtr[i]]$, and the rows of $B$ are loaded by reading $B.col$ from $B.rowPtr[idx]$ to $B.rowPtr[idx+1]$. Array $denseAccumBuff$ of size $m_{C}$ is updated for each column index of $B$, a companion bitmap stores the nonzero positions in $denseAccumBuff$, and $colBuff$ stores the running list of column indices in $C$. Besides $A$, $B$, and $C$, all variables are thread-local.

Algorithm 1 can be extended to other types of accumulators, e.g., hash map-based accumulators, where the dense accumulation array and bitmap are replaced by a hash map. In other Gustavson-based algorithms, such as expand-sort-compress (ESC) (ESC, ; pbSpGEMM, ), the intermediate product of $C$ is written to an array instead of directly updating the accumulator. This is shown in Algorithm 2, where the intermediate product is generated in the first loop by storing $B.col[k]$ and $A.val[j]\times B.val[k]$ in $colBuff$ and $valBuff$, respectively. The intermediate product is then sorted, duplicates are merged, and the result is written to $C$ (these steps take place in sortMergeWrite($colBuff$,$valBuff$)).

To compute $C.rowPtr$, which is an input to the numeric phase, an initial symbolic phase is required. The symbolic phase typically has the same high-level algorithm as the numeric phase but without performing the multiplication $A.val[j]\times B.val[k]$ and writing to $C$. For example, in Algorithm 1, the symbolic phase does not include the modifications to $denseAccumBuff$, $C.col$, and $C.val$. Instead, only the bitmap is updated along with a counter that outputs the exact number of nonzero entries for each row of $C$. Finally, $C.rowPtr$ is computed using a prefix sum on the counters. This type of symbolic phase is known as precise prediction, where the number of nonzero entries in $C$ is calculated exactly before the numeric phase.

On modern CPUs, maximizing cache reuse is crucial to the performance of any application. In SpGEMM, the accumulator is the most frequently accessed, where the amount of reuse is determined by the sparsity pattern of $A$ and $B$. For optimal performance, the dense accumulator should be confined to the L2 cache, which is the highest level of private cache. This efficient cache utilization occurs naturally in specific matrix structures, such as banded matrices or matrices that yield a highly sparse $C$. However, for matrices that produce “large” intermediate products (where “large” refers to both the number of nonzero elements and a wide distribution of column index values), SpGEMM faces significant challenges. A prominent example is random power-law matrices that model social networks (rmat, ). For such matrices, the size of $denseAccumBuff$ often exceeds the capacity of the L2 cache and the large intermediate product results in frequent accesses to the entire $denseAccumBuff$ array. Consequently, $denseAccumBuff$ must frequently be evicted from and reloaded to the L2 cache, resulting in suboptimal performance. This breakdown in locality presents a substantial obstacle for current SpGEMM algorithms, as we will demonstrate through microbenchmarks and SpGEMM results.

Optimizations related to the accumulation step have been the primary focus of research on SpGEMM. For load balancing, a common approach is based on the observation that rows with different intermediate sizes require different accumulation strategies (liu1, ; liu2, ; ESC, ; pbSpGEMM, ; liu3, ; OpSparse, ; spECK, ; nagasaka3, ). In (ESC, ), $A$ is reordered by increasing intermediate size. Other approaches group rows based on the size of the intermediate product, where a different accumulator is used for different groups (liu1, ; liu2, ; liu3, ; OpSparse, ; spECK, ; nagasaka3, ; adaptLoadBalanceGPU, ). In (liu2, ), five different accumulation strategies are used, including priority queues and sorting algorithms. Developing new accumulators is another topic that has been widely studied, especially for modern multicore machines with vector processors (fevre, ; hierRowMerging, ; ASA, ; AVX512sort, ; nagasaka1, ; nagasaka2, ; gamma, ). A common approach is to optimize sorting algorithms (ESC, ; tricount3, ; liu1, ; AVX512sort, ; pbSpGEMM, ; registerAware, ) or data structures such as heaps (azad, ; nagasaka1, ; nagasaka2, ) and hash maps (nagasaka1, ; nagasaka2, ; cuSPARSE, ; registerAware, ; mcl, ; balancedHashing, ; kokkos, ). In Section 4.3, MAGNUS is compared with hash map and heap-based approaches, which are considered state-of-the-art (cseg, ; nagasaka1, ; nagasaka2, ; mkl, ).

Perhaps most relevant to MAGNUS are recent works on improving the cache behavior of accumulators, proposed in (partway, ) and improved in the CSeg method (cseg, ). The core concept is to partition the columns of $B$ into segments, where the number of segments is chosen so that the dense accumulator fits in the L2 cache. An additional high-level summary matrix is used to store the segmentation information. CSeg was shown to be overall faster than many state-of-the-art libraries mentioned above. For a more extensive overview of SpGEMM research, including distributed memory algorithms, see (spgemmSurvey, ).

MAGNUS uses a novel hierarchical algorithm to generate two levels of locality. The coarse-level algorithm reorders the intermediate product into discrete chunks and the fine-level algorithm further subdivides and accumulates the coarse-level chunks. The number of chunks used in both algorithms is based on optimal parameters that are computed using the input matrix properties and the target system specifications. These parameters are optimal in the sense that they minimize the storage cost of all frequently accessed arrays. The levels are generated using a set of basic operations, including histogramming and prefix summing. Combining the building blocks and the optimal parameters creates the locality required by the accumulators.

For sufficiently “small” matrices (as discussed later in the derivation of the MAGNUS parameters), the fine-level algorithm alone provides an adequate level of data locality. This algorithm is based on Gustavson’s method, similar to most SpGEMM algorithms. However, achieving scalability for “massive” matrices requires both fine- and coarse-level locality. Here, massive matrices are those in which the data structures required by the fine-level algorithm, including the accumulator, exceed the capacity of the L2 cache. The coarse-level algorithm employs an outer product-based approach, necessitating an additional pass over the intermediate product, which increases the total data volume. This additional cost means that using the standalone fine-level algorithm wherever possible is advantageous, which is why we reserve the coarse-level algorithm only for massive matrices. Since the intermediate product is generated in both approaches, MAGNUS can be classified as an ESC-type algorithm (ESC, ).

Figure 1 shows a simple workflow for MAGNUS, where two threads are used to multiply two $8\times 8$ matrices. The example shows how rows 1 and 2 assigned to thread $t_{0}$ are computed, where two chunks are used for both the coarse and fine levels. The outer product-based approach traverses the submatrix of $A$ (corresponding to rows 1 and 2) column by column, and the highlighted rows of $B$ are traversed row by row. This traversal generates the four coarse-level chunks. Each chunk is reordered again to get the eight fine-level chunks, where the accumulator is applied to get the final result. It is important to note that all coarse-level chunks are generated before executing the fine-level algorithm. However, the fine level is processed depth-first: for each coarse-level chunk, the fine level is generated and then immediately accumulated (similar to a Gustavson-based method) before proceeding to the next coarse-level chunk.

The key property of MAGNUS is that the range of column indices in each fine-level chunk is significantly smaller than $m_{C}$ (the number of columns of $C$). This allows the dense accumulator to fit in the L2 cache when $m_{C}$ exceeds the L2 cache capacity. In Figure 1, the per-chunk range of column indices is two, compared to $m_{C}=8$. For example, the column indices in the fourth fine-level chunk of either row fall within the range $[6,7]$. Although not shown in the figure, the column indices in each chunk are shifted into the chunk-local range $[0,1]$, which means that we only need a dense accumulator of length two. If we consider a theoretical system with an L2 cache capable of storing a dense accumulator with a maximum of two elements, each fine-level chunk can be accumulated with minimal L2 cache misses. In contrast, if we used Gustavson’s method with a conventional dense accumulator, L2 cache misses would occur frequently after loading the first two elements in the second row of $B$.

The high-level steps of MAGNUS are: (1) the setup phase, which includes calculating the optimal number of chunks (this will be discussed later in Section 3.5), computing the intermediate product size for each row, and categorizing each row; (2) the symbolic phase; and (3) the numeric phase. The setup phase is inexpensive compared to the symbolic and numeric phases: calculating the optimal number of chunks has constant time complexity, while the remaining setup steps involve a highly parallel single pass over the rows of $C$, with time complexity $O(n_{C}/t)$, where $t$ is the number of threads. Row categorization is necessary because not all rows require locality generation.

MAGNUS categorizes each row based on its structure and the system’s specifications:

1. (1)

   Sort: If the number of intermediate elements is less than the dense accumulation threshold, we can directly apply a sort-based accumulator, as in Algorithm 2. The dense accumulation threshold will be described later in this section.

2. (2)

   Dense accumulation: If the intermediate row length fits into the L2 cache, we can directly apply dense accumulation to the row, as in Algorithm 2. This is because the range of column indices does not exceed the size of the L2 cache. The intermediate row length is the difference between the minimum and maximum column index of the intermediate product.

3. (3)

   Fine level: If $s_{finelevel}<s_{L2}$, the fine-level algorithm can be applied, where $s_{finelevel}$ is the number of bytes required to store all necessary fine-level data structures (this is discussed later in this section).

4. (4)

   Coarse level: The coarse-level algorithm is applied to all remaining rows, where the fine-level algorithm is applied to each coarse-level chunk.

The accumulation parameters mentioned above will be discussed in Section 3.4, including a description of the sorting algorithm. The first two categories imply that some rows possess intrinsic locality and do not benefit from the locality-generation algorithms in MAGNUS. After the rows are categorized, the rows of $C$ are computed category-first to ensure that data specific to a particular category are cached for as long as possible. In our OpenMP implementation of MAGNUS, we use a parallel for loop with dynamic scheduling to traverse the rows in each category with a no wait clause. The no wait clause ensures that threads proceed to the next category without unnecessary synchronization.

For simplicity, we assume $m_{C}$ is a power of two in our descriptions of the algorithms in MAGNUS. We use precise prediction for the symbolic phase, but for brevity, we will only describe the numeric phase. See Section 2 for clarity on the differences between the symbolic and numeric phases.

The fine-level algorithm has the following steps for each row: histogram, prefix sum, reorder, and accumulate, where $n_{chunksFine}$ is the number of fine-level chunks. As in Gustavson’s method, each row (or coarse-level chunk in cases where the coarse-level algorithm is applied) is computed before moving on to the next row. This means that the intermediate product is generated only for a single row (or coarse-level chunk) at any given time, unlike in outer-product-based approaches.

Pseudocode for the fine-level algorithm is shown in Algorithm 3, where the input is the column indices and values of a single coarse-level chunk for row $i$ of $C$. For rows that only require fine-level locality, $A$ and $B$ are read directly as in Algorithm 1, i.e., the loop headers on lines 3 and 12 of Algorithm 3 are replaced with the nested loop headers on lines 4 and 6 of Algorithm 1, respectively. The notation $arr\leftarrow\alpha$ means that the entire array $arr$ is initialized to the value $\alpha$. The C-style notation $arr[i]$ means that element $i$ of $arr$ is accessed, and $\&arr[i]$ means that the array is accessed starting at element $i$.

The first step towards reordering the intermediate product is to compute $offsetsFine$ (the chunk offsets) using a histogram and a prefix sum operation. The array $offsetsFine$ stores the start and end locations of each fine-level chunk in $colFine$ and $valFine$. The buffers $colFine$ and $valFine$ are typical of any ESC-type algorithm where the intermediate product must be explicitly stored. In the case of MAGNUS, they store the reordered intermediate product. In the histogram step, the column indices are mapped to chunks as $col/chunkLenFine$, where $chunkLenFine=m_{C_{maxL2}}/n_{chunksFine}$ is the chunk length of the fine-level chunks. The chunk length is the local range of the column indices within a chunk, where column indices are shifted into the range $[0,chunkLenFine)$ as shown on line 15. The value of $m_{C_{maxL2}}$ is equal to $m_{C}$ if the fine-level algorithm is used alone, or equal to the coarse chunk length if the coarse-level algorithm is used ($m_{C_{maxL2}}$ is discussed in more detail in Section 3.5). To optimize the mapping, the division operation is replaced with a bitwise operation by restricting $n_{chunksFine}$ to a power of two, ensuring that $m_{C}/n_{chunksFine}$ is also a power of two. Therefore, the mapping becomes $col\texttt{>>}shiftFine$, where $shiftFine=\log_{2}(chunkLenFine)$ and >> is the bitwise right-shift operator.

After the histogram is computed, the chunk offsets are computed using a prefix sum (inclusive scan) of the histogram. To reorder the input (shown in the second loop), the histogramming phase is repeated, where the histogram is used to track the current number of elements in each chunk. Elements are reordered by writing the input coarse-level chunk at the position $offsetsFine[chunk]+countsFine[chunk]$ in $colFine$ and $valFine$. As mentioned previously, the column indices are shifted into the local range of each chunk as $col-chunk\times chunkLenFine$ to allow for cache-efficient accumulation.

Finally, a call to accum() for each chunk invokes either sort-based accumulation or dense accumulation, where the size of $denseAccumBuff$ (see Algorithm 1) is now reduced from $m_{C_{maxL2}}$ to $chunkLenFine$. After the accumulation step, the column indices are shifted back into the correct range before writing to $C$. The variable $rangeFine_{j}$ denotes the range of column indices of the fine-level chunk $j$ (i.e., $[j\times chunkLenFine,(j+1)\times chunkLenFine)$), and $rangeCoarse$ is the range of column indices of the input coarse-level chunk. Figure 2 shows the workflow of the fine-level algorithm in terms of the data structures from Algorithm 3, where the input is the first chunk from the example in Figure 1.

The fine-level algorithm requires two additional arrays: $countsFine$ and $offsetsFine$, both of size $n_{chunksFine}$. Alongside $denseAccumBuff$ and $bitMap$, the goal is to keep $countsFine$, $offsetsFine$, and the active cache lines of $colFine$ and $valFine$ in the L2 cache. The active cache lines must be considered since we are writing to $colFine$ and $valFine$ at noncontiguous positions. In our implementations, we use non-temporal streaming stores when writing to $colFine$ and $valFine$, which avoids polluting the L2 cache. Non-temporal stores are intrinsic functions used on Intel processors (e.g., _mm512_stream_si512()) that write to memory without evicting cached data, allowing us to retain the accumulator and fine-level data structures in the L2 cache while streaming the intermediate product. Section 3.5 shows how we choose the optimal number of chunks that minimizes the total storage cost of these cached arrays.

As the columns of $C$ increase, the storage of the fine-level data structures eventually exceeds the size of the L2 cache. An initial coarse level must be generated for such matrices, providing the first level of locality. To generate the coarse level, we use a modified outer product-based algorithm, where the intermediate product is generated and reordered for all rows that require coarse-level locality before any accumulation occurs. The reordered intermediate product is organized into discrete coarse-level chunks that can be handled independently by the fine-level algorithm.

The outer product-based approach is used because the coarse-level algorithm computes only the intermediate product without performing any accumulation. Therefore, maximizing the reuse of the input matrices rather than the accumulator is beneficial, which is a well-known property of outer product-based SpGEMM algorithms (pbSpGEMM, ; gamma, ). The combination of our reordering algorithm with the outer product formulation makes the coarse-level algorithm similar to propagation blocking-based algorithms (pbSpGEMM, ).

Conventional outer product algorithms typically generate the intermediate product by multiplying the columns of $A$, stored in CSC format, with the rows of $B$, stored in CSR format. The coarse-level algorithm in MAGNUS follows the same scheme, but on the subset of the rows categorized as coarse-level rows. Therefore, a CSC version of the submatrix $\hat{A}$ of $A$ is constructed, where $\hat{A}$ only includes these coarse-level rows. Each thread performs the following steps on its list of coarse-level rows, which are stored in the array $coarseRowsC$:

1. (1)

   For all $i\in coarseRowsC$, generate $coarseRowsB$, i.e., the unique set of rows in $B$ required to perform the outer product. This is done by iterating through the nonzero entries in $\hat{A}$ and setting a bitmap, where $coarseRowsB$ is the list of set bits.

2. (2)

   Construct the thread-local CSC submatrix $\hat{A}^{CSC}$ using the well-known approach for the conversion of CSR to CSC: a histogramming stage computes the number of nonzero entries per row, a prefix sum computes $\hat{A}^{CSC}.colPtr$, and $\hat{A}^{CSC}.colPtr$ is then used to construct the arrays $\hat{A}^{CSC}.row$ and $\hat{A}^{CSC}.val$. This step is inexpensive since the time complexity of outer product-based SpGEMM algorithms is dominated by the generation of the intermediate product in the following step. Specifically, the time complexity of constructing $\hat{A}$ is $O(nnz_{\hat{A}})$, while generating the intermediate product requires $O(\sum_{(i,j)\in\mathcal{S}(\hat{A})}nnz_{B_{j}})$ time, which is best-case $O(nnz_{\hat{A}})$ (i.e., each row in $B$ contains only one nonzero value) and worst-case $O(nnz_{\hat{A}}m_{B})$ (i.e., each row in $B$ is dense).

3. (3)

   Perform a histogram step by reading $\hat{A}^{CSC}$ and the rows of $B$ corresponding to $coarseRowsB$. A prefix sum of the resulting histogram generates the coarse-level offsets.

4. (4)

   Generate the coarse level by again reading $\hat{A}^{CSC}$, reading the rows of $B$ corresponding to $coarseRowsB$, and writing the result in $colCoarse$ and $valCoarse$ at the positions determined by the offsets. As in the fine-level algorithm, we use non-temporal streaming stores for the reordering.

5. (5)

   For each coarse-level chunk, apply the fine-level algorithm.

Steps 3-5 are shown in Algorithm 4. The same set of basic building blocks is used as in the fine-level algorithm (histogramming, prefix summing, and reordering). The key difference is that the chunks for all rows are generated at the same time. The prefix sum is modified to compute the chunk offsets within a row based on the offsets from the previous rows. Just as in the fine-level algorithm, we generate the reordered intermediate product using non-temporal streaming stores. In the final loop, the fine-level algorithm is applied to each coarse-level chunk for each row. Figure 3 shows the data structures for the example in Figure 1.

The coarse-level algorithm utilizes a set of data structures similar to that of the fine-level algorithm, but without an accumulator. The absence of a cached accumulator enables the coarse-level algorithm to generate more chunks than the fine-level algorithm. However, this comes at the cost of increased data volume due to the additional pass over the intermediate product. Despite this tradeoff, for a sufficiently large matrix, the increased data volume proves more efficient than the frequent cache misses incurred by the fine-level algorithm, as we will demonstrate in our microbenchmark results.

Similarly to other outer product-based approaches, the memory requirement for the coarse-level algorithm is higher than the fine-level algorithm since the intermediate product for multiple rows must be stored. In the worst case, the memory requirement is proportional to the sum of the outer products of all rows, which may exceed the memory capacity of certain memory-constrained systems. Therefore, we use a batching approach, where we collect rows of $C$ into $coarseRowsC$ until one of two conditions is met: (1) the memory limit of our system is reached, or (2) the storage requirement for generating the coarse-level chunks ($countsCoarse$, $offsetsCoarse$, and the small buffers for the non-temporal streaming stores) exceeds the L2 cache size. If either condition is met, the batch of rows in $coarseRowsC$ is computed via steps 2-5. This batching process is repeated until all rows requiring coarse-level locality have been computed.

MAGNUS is accumulator-agnostic, allowing for portability and flexibility: only the storage requirements of the desired accumulators are needed to compute the optimal MAGNUS parameters (the optimal parameters are derived in Section 3.5. For portability, accumulators optimized for specific architectures can be selected without changing the locality-generation algorithms. This is evident in Algorithm 3, where the accum() function requires only an input chunk. For flexibility, MAGNUS allows for a hybrid approach in which each chunk chooses an accumulator based on the chunk characteristics. In this paper, we consider two accumulators: AVX-512 vectorized bitonic sorting (AVX512sort, ) and classical dense accumulation. For chunks with a small number of elements, sorting is performed. Otherwise, dense accumulation is used. When visiting each chunk, a threshold is used to choose the accumulator. This threshold is based on how the sorting algorithm is implemented: quicksort partitions the array, and then hard-coded, vectorized bitonic sorting networks sort the partitions. We found experimentally that dense accumulation is faster than sorting unless the sort size is small enough to bypass the quicksort algorithm and directly use bitonic sorting. For more details, see Section 4.1.1 and (AVX512sort, ).

The optimal number of chunks is chosen based on the following input parameters, which are readily available for end users: $s_{cacheLine}$, $s_{L2}$, and $m_{C}$. Parameters $s_{cacheLine}$ and $s_{L2}$ are the cache line and L2 cache sizes, respectively, which are retrieved by querying the underlying system, e.g., by using standard Linux commands. Parameter $m_{C}$ is the number of columns of $C$, which is already included in the CSR data structure of $B$ since $m_{C}=m_{B}$.

As explained in Section 3.2, the goal is to retain in the L2 cache certain data structures needed by the fine-level algorithm. For simplicity, assume $m_{C}$ is a power of two ($m_{C}$ is ceiled to the nearest power of two otherwise). Choosing the optimal number of fine-level chunks corresponds to selecting the value of $n_{chunksFine}$ that minimizes the convex function

where $s_{fineLevel}$ is the storage requirement for the L2-cached fine-level data structures. The first term is the storage requirement of the dense accumulator, where the number of elements in the underlying dense accumulator array is $m_{C}/n_{chunksFine}$. For the numeric phase, $s_{denseAccum}=s_{val}+1$, since we need $denseAccumBuff$ and $bitMap$ (see Algorithm 1). For the symbolic phase, $s_{denseAccum}=1$ since we only need $bitMap$. The second term is the storage requirement for reordering, where the storage cost per fine-level chunk is $s_{chunkFine}=s_{histoType}+s_{prefixSumType}+2s_{cacheLine}$. The parameters $s_{histoType}$ and $s_{prefixSumType}$ denote the size of the histogram and prefix sum array data types, respectively, which are both four bytes. The term $2s_{cacheLine}$ accounts for the storage of active cache lines when noncontiguously writing to $colFine$ and $valFine$ during the reordering phase. The value of $n_{chunksFine}$ that minimizes Equation 2 is