Despite drawing great attention recently [1]- [11], accelerating graph analytics on GPUs remains challenging in that: (i) graphs are often irregular in connectivity, causing warp divergence, memory non-coalescing, and low arithmetic intensity; (ii) graphs are often large, many of which cannot directly fit into the caches, shared memory, or even the DRAM of GPUs. To address these challenges, the matrix-centric approach has been proposed and increasingly employed by GPU-based graph frameworks [3], [4], [12], [13]. Unlike traditional graphcentric frameworks that go through relevant nodes or edges iteratively, this approach employs sparse storage formats-such as compressed sparse column (CSC) or compressed sparse row (CSR)-to represent the adjacency matrix of a graph and then uses highly optimized linear algebra kernels (sparse matrix-vector and matrix-matrix multiplications, SpMV and SpGEMM) on the sparse formats for computation.
Although GraphBLAS has shown improved performance over traditional approaches in handling large graphs, there is still considerable potential to tap into. This work aims to unlock the potential by exploiting bit-level representations and optimizations on GPUs. The underlying observation is that for a large class of graphs (homogeneous graphs), a single bit is sufficient for indicating the adjacency relation of two vertices in a graph (1: adjacent, 0: not adjacent). Given that many graph algorithms center around computations upon adjacency matrices, using a bit-level representation can potentially reduce storage usage and improve computation efficiency.
Bit representations (bitmaps, bitvectors) have been used in vertex-based graph frameworks [14]- [16] for representing frontiers (i.e., active nodes); bit-level optimizations have, however, not yet been systematically explored in matrixbased graph frameworks. Existing GraphBLAS [17], [18] frameworks typically build on existing linear algebra libraries, which offer no bit-level representations of matrices or bit-level implementations of linear algebra functions.
This work answers three key research questions. RQ-1: What storage format should be used for a binary adjacency matrix?
Unlike the Boolean data structures (frontiers) in graphcentric frameworks, for matrix-based graph frameworks, the binary data structure in focus is the entire adjacency matrix, which sits at the center of the linear algebra operations in GraphBLAS. Systematic studies are needed for manipulating it; simply representing it as a bitmap cannot tap into the full potential of space savings by accommodating many unnecessary zeros; that also causes difficulties for the graph operations to leverage the highly tuned matrix-based libraries.
Based on the properties of adjacency matrices and various tradeoffs, we design a storage format, namely Bit-Block Compressed Sparse Row (B2SR). B2SR is inspired by the Block Compressed Sparse Row (BSR) format [19]. It takes a twolevel structure: The upper level is similar to BSR's upper level, using a sparse format to represent the locations of non-zero blocks (or called submatrices); the lower level differs from BSR in that it represents each non-zero block as a dense bit matrix-that is, each element in the block becomes one bit in the representation. The representation allows it to efficiently harvest the hardware computation capability at the low level and at the same time leverage the (regional) sparsity at the high level. Balancing the space savings and the indexing overhead lays the foundation for tapping into the potential of sparse binary matrices.
RQ-2: How to efficiently compute on the new representation?
We propose several new algorithms to implement the critical linear algebra kernels (SpMV and SpGEMM) for manipulating sparse matrices represented in B2SR. The design carefully tailors the kernel implementations around the efficient low-level bit-manipulation intrinsics on GPUs. It exploits the new optimizations and hardware-specific capability brought by the latest GPUs. These kernels lay the foundation for efficient manipulations of sparse binary matrices for graph analytics.
RQ-3: What are the performance implications? We evaluate the B2SR-based SpMV and SpGEMM on 521 binary matrices and five graph algorithms. The benefits are significant. B2SR provides up to 32⇥ space savings. On two generations of GPUs (NVIDIA's Pascal and Volta), we observe 40⇥ and 6555⇥ maximum speedups over the stateof-the-art sparse linear algebra libraries cuSPARSE [20] and GraphBLAST [4]. On graph algorithms, it offers up to 433⇥ acceleration on Breadth-first-search (BFS), 55⇥ on Single-Source-Shortest-Path (SSSP), 28⇥ on PageRank (PR), 69⇥ on Connected Component (CC) algorithms, and 52⇥ on Triangle Counting (TC) algorithm over GraphBLAST [4], a GPU graph processing framework with state-of-the-art performance. As is well known [21]- [23], no sparse format fits all matrices. So despite the effectiveness of B2SR, there are matrices that fit other sparse formats better. We provide a brief discussion and a simple sampling approach to assisting users in applying B2SR.
Graph programming frameworks are based on either graphcentric abstraction (vertex-or edge-centric) [1], [24]- [28] or matrix abstraction. For the performance advantages and direct leverage of advances in high-performance linear algebra libraries, matrix abstraction has received increasing interest in recent years. GraphBLAS [17], [18] is the mathematical core of matrix-based graph frameworks. It models graph traversal as operations on semi-rings. Frameworks that implement the standard include nvGraph [12], cuGraph [29], SuiteSparse [13], GraphBLAS template library (GBTL) [3], GraphBLAST [4], and so on. Among them, GraphBLAST [4] represents state of the art, achieving high performance on GPU by exploiting input and output sparsity [30] and enhanced load balance by exploiting the memory access patterns of sparse matrix multiplication.
These frameworks are mainly built under the same line of unified graph construct-using CSC or CSR to establish floating-point element space and perform matrix operations with underlying linear algebra libraries. They have not exploited bit-level optimizations. Even though in the sparse linear algebra libraries, bitmaps or bitvectors may be used to index the non-zero elements in a sparse matrix [31]- [33], those libraries fundamentally assume that the sparse matrices are general rather than binary matrices. They hence leave an immense performance potential untapped (as our comparison in Section VI shows).
In graph-centric frameworks (e.g., GraphMat [14], Graphlt [15], SpbLA [16]), there are some Boolean data structures (e.g., frontiers or active nodes), which are sometimes represented in bitmaps or bitvectors. Some works [34]- [36] exploit binary encoding or compression to achieve storage reduction and algorithm acceleration. For matrix-centric graph frameworks, the binary data structure in focus is the entire adjacency matrix, which sits at the center of the linear algebra operations in GraphBLAS-simply representing it as a bitmap leaves not only lots of potential 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 nRows = 16 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 nCols = 16 nTileRow = 2 0 1 2 1.00 1.00 1.00 1.00 1.00 1.00 1.00
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 for space savings but also causes difficulties for the graph operations to leverage the highly tuned matrix-based libraries. Therefore, systematic studies are needed for manipulating it. Bit-manipulation primitives have been used in enhancing the performance of deep neural networks (DNNs) on GPUs [37]- [39], ASICs, and FPGAs [40]- [42]. They are about dense binary operations, so they cannot efficiently handle sparse ops as those in graph processing. The irregularity of sparse matrices and accesses makes the issue more complex.
This section presents B2SR, the format we have designed for representing an adjacency matrix. Figure 1 illustrates the design. The principle we followed is that the representation should reduce the space cost as much as possible and at the same time facilitate the acceleration of the core graph operations. Drawing on the inspiration of BSR, we create the two-level representation of B2SR. The top-level takes advantage of well-proven effective sparse formats (CSR or CSC) to represent non-empty blocks. The bottom-level treats each non-empty block as a dense block and packs its elements into a bit representation. The combination of sparse formats and bit representation minimizes space usage, while the blocklevel dense bit format preserves low-level regularity, making efficient computation possible. We will explain the format in detail.
Since adjacency matrices are all square, it is natural to have the number of tile rows (nTileRow) set as nRows+tileDim 1 tileDim , where tileDim is the dimension (or bit-width) of the tile (e.g., 4,8,16,32). The number of non-empty bit-tiles can then be inferred from the nonzeros' coordinates of the sparse matrices. In our implementation, we use cuSPARSE's cusparseXcsr2bsrNnz() API to obtain the number of non-empty tiles from the CSR format. We utilize indexing arrays to record the coordination of the non-empty bytealigned bit-tiles. The proposed format comprises three arrays:
• Tile row indices (TileRowPtr): an integer array with the size of the number of tile-rows. It records the bittiles' row indices. It is an accumulated array with the i-th element recording the sum of total non-empty bittiles counting from the first tile row to the (i-1)-th row. Therefore, T ileRowP tr[i+1] T ileRowP tr[i] suggests the number of non-empty bit-tiles in the i-th tile-row. Table I: Binarized packing format. • Tile column indices (TileColInd): an integer array with the size equal to the number of non-empty bit-tiles. This array is for recording the tile column indices in the tile coordination system. • Bit-tiles storage (BitTiles): a bit-packing type (unsigned char, unsigned short, unsigned int, or unsigned long long int) array with a size equal to the tileDim⇥numofTiles (number of non-empty bit-tiles). It stores the binarized non-empty bit-tiles' layout in the order of their tile column indices. The proposed format has several merits: (1) It allows simpler transpose of the sparse matrix. By transforming the TileRowPtr and TileColInd from CSR to CSC, the sparse matrix is transposed. (BitTiles do not require transpose since we default the mxv() and mxm() algorithm to access the content of a tile always in row-by-row order.) We use cuSPARSE's cusparseScsr2csc() API to enable this function in our implementation. (2) The format has a better data accessing locality for SpGEMM and SpMV when computing in tile-row by tile-row order; Since the storage format is similar to CSR, we can use existing optimization on CSR-based algorithms to spearhead the linear algebra kernels in use. (3) In BitTiles, the format provides storage-saving compared to CSR. The proposed format carries extra zeros than CSR and COO, which store only the nonzeros. Nevertheless, with the binarized packing yielding up to 32⇥ space-saving per tile, most sparse matrices can still benefit from storage compression, especially when configured in small tile size.
Figure 2 presents the column-major and row-major packing of a bit-tile. We adopt column-major packing as default when converting CSR to B2SR. The transpose of B2SR tiles can be achieved by storing the additional row-major layouts. To figure out the proper byte-addressable data types to carry bits in proximity, we explore the packing granularity from 8-bit to 32-bit (reference Table I). Additionally, we use half of the space in an unsigned char to allow 4-bit (nibble) packing, which further reduces half in one dimension to carry unnecessary zero when the matrix is extremely sparse. The packing result yields four variants of B2SR format: B2SR-4 (for 4⇥4 tile size), B2SR-8 (for 8⇥8 tile size), B2SR-16 (for 16⇥16 tile size), and B2SR-32 (for 32⇥32 tile size). Generally, the space-saving depends on how the adjacency matrix is initially being stored. State-of-the-art GPU graph frameworks mostly use float to carry the elements, so our bit-packing can generally provide up to 32⇥ storage savings per square tile; for frameworks that use double to carry the elements, the bit-packing results in up to 64⇥ savings in space. This indicates an ability to store 32⇥ or 64⇥ larger graphs using the same amount of space. This also brings a potential reduction in required data accessing bandwidth during computation to enable higher throughput.
Bit-packing overhead To transform CSR to B2SR, we parallelize each tile-row's encoding procedure for the large graph. For the 4, 8, 16, or 32 continuous elements in CSR's RowPtr, we use cusparseXcsr2bsrNnz() and cusparseScsr2bsr() to obtain tile-row index and fullprecision tiles along a tile-row. Next, column-major or rowmajor bit-packing kernels are applied to process each tile's encoding. The routine's overall cost is about 3 to 34 ms. In practical applications, a graph can be reused by many users; even within one execution, a graph is often used repeatedly (e.g., for iterative processing). So despite format conversion may be needed, such a one-time cost can be greatly amortized in these conditions.
While the proposed method provides storage compression and performance gains in BLAS operators, it is evident that not all graph matrices are suitable for converting to Bit-GraphBLAS binarized format. For example, graph matrices with randomly distributed connections (nonzeros) or relatively dense patterns may not necessarily benefit from this format. When each bit tile does not capture enough nonzeros, we can have many empty bit-rows in tiles after bit packing. Converting the matrices from CSR to B2SR is not ideal since it expands the total amounts of storage. In addition, it adds additional per-thread workloads (OPS) compared to CSR's BLAS kernels in terms of operators. We observe that the tile size selection trade-off can be considerably different in various matrix patterns. Figure 3a shows that when the tile size equals 4⇥4, there are less than 30% non-empty tiles; when the tile size equals 32⇥32, non-empty tiles reach more than 80% for some matrices. The reason is that although increasing tile size can decrease the number of non-empty tiles, the reduction is often less than 4⇥ (the times of per tile size increment), causing the ratio of non-empty tiles to increase ultimately. We find that despite the increment, the total B2SR byte size can sometimes decrease because of the reduction in the number of tiles and indexing arrays. For instance, in matrix mycielskian12, we have CSR storage at 3.12 (MB), B2SR-4 at 675.70 (KB), B2SR-8 at 361.46 (KB), B2SR-16 at 358.89 (KB), and B2SR-32 at 429.89 (KB). The total byte size of the format does not monotonically increase as tile size increases. From the other angle, Figure 3b presents the reduction in the average occupancy of nonzeros in the non-empty tiles. The percentage of actual nonzeros in tiles can drop from 20% to less than 5% as the tile dimension differs. If the bit-tile is too large, the computation may waste processing too many tiles; if the bit-tile is too small, the indexing array may carry more unit workloads. Therefore, to fully utilize the double benefits (compression and computation efficiency) of Bit-GraphBLAS, graph users can first identify the potential benefits offline through our sampling profile method. In this way, when tackling extensive graphs or cumbersome, repetitive graph computations, users can experience worthy time and labor cost saving with only one-time format conversion to B2SR. The sampling works as in Algorithm 1. The user first specifies the number of rows to sample as N for N N umOf Rows (sampling more rows can accurately capture matrix characteristics but induces more significant overheads). We then have N random indices as the random index set S ✓ [0, 1, 2, ..., N]. Algorithm 1 Sampling Profile Scheme 1: for i in S do 2:
The sampling result provides a rough estimation of the compression rate of Bit-GraphBLAS on B2SR-4 to B2SR-32. Users can select the affordable compression rate and tile size configuration.
This subsection discusses the code patterns and scheme designs that leverage the proposed storage format. Our imple- We briefly introduce the GPU bit operation intrinsics that will be used in this work, with respect to the BSTC bitblock abstraction [37]: (1) __popc(): The population count function is for efficient bit-accumulation across a singe bit-row. CUDA supports population count along a 32-bit unsigned int. Paired with the logical AND operation, it can perform bit-dotproduct for two 32-bit bit-rows. ( 2) __shfl_sync(): This intrinsic is for exchanging the bit-row across the lanes of a warp. For BMM, it facilitates a faster bit-dot-product between a bit-row and multiple bit-columns. ( 3) __ballot_sync(): This warp-vote intrinsic returns a 32-bit unsigned integer whose N-th bit indicates a predicate setting by the N-th thread of a warp (assuming all threads are active). This is essentially equal to transposing a bit-column to a bit-row. Since bits are indexed from right to left in the bit-row packing, the function is equivalent to a 90 clockwise transposition to a bit-row. ( 4) __brev(): This intrinsic is used in bit-packing. Together with __ballot_sync(), this function rotates a bit-column 90 anti-clockwise into a bit-row.
The kernels of matrix-centric graph computing are matrixvector and matrix-matrix computations. Table II and III list the core schemes we have implemented. They correspond to the different scenarios of the inputs and outputs (binary "1-bit" or full precision "32-bit"). We use two of the schemes to explain our implementations:
Listing 1. As an example scheme of Binarized Sparse Matrix Multiply Vector (BMV), Listing 1 shows the code of bmv_bin_bin_full() in 32⇥32 tile size. The function demonstrates the bit multiplication between the binarized bittiles and the binarized vector. The result is a full-precision vector. The computation is as follows:
Before computation, the sparse matrix is packed into the hierarchical storage format -the vector is binarized into the column-major order with 32 consecutive elements compacted as an unsigned int. This allows the bit-columns to be fetched according to the same indexing system and enables fast bitdot-product with each bit-row in the bit-tiles. The A variable is a bit-row in a tile, and the B variable is the binarized vector. The output array C is a vector in full precision. In each warp of the thread block, the number of bit-tiles to be computed is indicated by the TileRowptr. In each iteration, the bit-tile and the 32-binarized vector perform bit-matrixvector-product using bit-wise AND and __popc(). Each lane Listing 2. As an example of Binarized Sparse Matrix Multiply Binarized Sparse Matrix (BMM), Listing 2 shows the code of bmm_bin_bin_sum() for B2SR-32. The A and B variables are the two bit-vector in tiles of the input sparse matrices. The output C is a single variable in full precision, summing up the nonzeros (1s) of the resulting bit matrix. In each warp of the thread block, the number of bittiles in A's tile row is indicated by the TileRowptr of A. In the outer "for" loop, the TileColInd of each bit-tile is used to retrieve the corresponding tile rows and bit-tiles of these tile-rows. The inner "for" loop performs the bit-tile-matrixmatrix-multiplication, where __shfl_sync() is used to retrieve the B's bit-vectors in each lane. The temporary result of each bit-vector in B is accumulated in separate registers for avoiding race conditions. Finally, the content in the 32 registers is summed up. The sum is atomically added to a single destination element in C.
Allowing efficient full-precision vector load In BMV and BMM, the bit-tiles in a tile-row are handled in a warp of 32 threads, following the warp-consolidation model [43]. By default, each of the thread blocks contains only one warp. So up to 64 thread blocks can be freely scheduled by a single SM scheduler. This works fine when kernels are computing in a binarized vector or matrix. However, in bmv_bin_full_full(), when engaging a full-precision vector as the multiplier, the one-warp-per-thread-block design impedes the flexibility to preload the common vector portions for the neighboring tile rows into shared memory. Thus, we implement the kernel scheme with 32 warps processing consecutive 32 tile-rows in a thread block. This further enhances the spatial locality of the computation workload. Figure 4 shows the thread mapping for B2SR-32, B2SR-16, B2SR-8, B2SR-4. The execution latency is generally shorter when 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 … 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 … 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 cooperatively loading 128 bytes (equal to the cache line size of GPUs in our experiment) of binarized tiles from global memory. At the same time, the corresponding number of fullprecision subvectors have to be loaded for multiplication. We set the thread block to contain 1024 threads to load the vectors into shared memory before multiplication. Bit-vectors on tilerow with the same column index can share the preload vectors from shared memory.
V. GRAPH ALGORITHMS This section explains how graph programs can be implemented upon the core operations.
Matrix-centric graph computing models graph traversals as operations on semirings [44]. As shown in Table IV our implementation can support the key semiring domain operations when performing vxm(), mxv(), mxm(). After the adjacency matrix is in B2SR, it remains binary throughout all operations. The vectors representing the frontier nodes are all in dense format. They can be either binarized for binary semiring or full-precision (float, unsigned, bool, etc.) for the non-binary domains to support a variety of graph algorithms. We also implement efficient masking schemes for both BMV and BMM: bmv_bin_bin_bin_masked(), bmv_bin_bin_full_masked(), bmv_bin_full_ full_masked(), and bmm_bin_bin_sum_masked().
We next use two graph algorithms to illustrate how to utilize these kernel backends when writing graph algorithms.
Breadth-First-Search (BFS) Breadth-first-search uses boolean semiring. In each iteration, the vxm() performs onedegree edge traversal to all the connected vertices. A mask of visited vertices is applied at the end to filter out the visited results. We introduce bmv_bin_bin_bin_masked() to enable this. GraphBLAST uses early exit to eliminate the masked element operations in their masked vxm(). Yet, a similar strategy does not apply to our case. In our implementation, the consecutive rows in a tile row are operated in the same warp. Early exit causes a performance penalty because of warp divergence. Therefore, in the bmv_bin_bin_bin_masked() kernel, the bitmask is applied right before the output store, having bit-wise AND with the negation of visited vertex vector (indicates unvisited vertices).
Single-Source Shortest-Path (SSSP) We implement the algorithm with delta-stepping SSSP [8] as in Graph-BLAST. SSSP utilizes tropical min-plus semiring. The intermediate vectors are reduced by minimum operation. bmv_bin_full_full() maintains the multiplier vector in full-precision, allowing it to carry minimum distance values. To realize the relaxation, we set an extra condition in the bmv_bin_full_full() such that the 0s in the adjacency matrix are identified as infinite (1), indicating unreachable. Only 0s along the diagonal are treated as actual zeros, which we omit their self-connectivity. Within a warp, a thread reduces all non-zero full-precision values along the multiplier vector by Min() (for B2SR-32). In B2SR-4, B2SR-8, and B2SR-16, since we use more than one thread to process the values along the multiplier vector, atomicMin() is applied to avoid race conditions.
PageRank (PR) PR uses arithmetic semiring. In each iteration, the page rank vector is multiplied by the column stochastic adjacency matrix. The column stochastic adjacency matrix is the adjacency matrix with each out-vertex connectivity divided by the vertex's out-degree. Since the page rank vector is in full-precision, we use bmv_bin_full_full() with an auxiliary vector v_out_degree to accommodate each vertex's out-degree. For each 1 on the matrix, the corresponding value on the page rank vector is divided by its out-degree on v_out_degree. Eventually, the intermediate vector is summed up with add operation to the output, indicating the weighted sum. Likewise, B2SR- process the workload along the vector cooperatively.
Connected Component (CC) We follow the CC implementation in GraphBLAST, which is based on the FastSV linear-algebraic connected component alogrithm [45], [46]. Similar to SSSP, CC uses tropical min-plus semiring. We adopt bmv_bin_full_full() since the frontier vector should be in full-precision. The mxv() is achieved by reducing the nonzero full-precision values along the intemediate vector using Min() and atomicMin().
Triangle Counting (TC) We implement TC as in Graph-BLAST, following Azas and Buluc's [47] and Wolf's [48]. The TC uses arithmetic semiring. It is achieved by multiplying the lower triangle of the adjacency matrix (L) with the transpose of itself (L T ) and then applying (L) as the mask to generate the output matrix. Ultimately, the non-zeros is summed up into one full-precision value; therefore, we fuse the reduction sum kernel with mxm() and directly perform atomicAdd() to global sum once a bitmap subroutine finishes. Since both input matrices and mask matrix can be sufficiently represented in binary format, so we realize the kernel through bmm_bin_bin_sum_masked().
Dataset We use all 521 binary square matrices in the SuiteSparse Matrix Collection [13]. The set of matrices contains the number of rows and columns ranging from 2 to 214,005,017 and the number of nonzeros from 2 to 11,588,725,964. To better summarize the similarity between the matrices with higher or lower performance in the evaluation, we further classify the matrices into six categories based on their patterns in Table V. GPU Environments We evaluate the proposed format and computation core functions on two NVIDIA GPU architectures, including Pascal and Volta. With the compute capability 6.0 and 7.0 configured, respectively, how the proposed format adapts to each hardware-specific variance is worth seeing. Table VI shows the configured SM information and memory hierarchy size of the two GPU architectures in the evaluation. We use CUDA version 10.0 across all our evaluations.
Algorithm Parameters The optimization setup of Graph-BLAST's algorithms is based on the default running script provided in their GitHub repository [49]. BFS is with earlyexit, structure-only, and operand reuse enabled. PR is limited to a maximum iteration of 10. The alpha parameter is set to 0.85, and pdfilon is set to 1e 9. The runtime for all algorithms is measured by the average of 5 runs.
B2SR brings significant storage savings for large sparse matrices. We show the compression ratio of the 521 binary graph matrices with respect to the default 32-bit floating-point CSR. The compression ratio thus is defined as: B2SR size CSR size . A lower value indicates a better compression rate. The compression ratio depends mainly on the nonzero distribution of the binary matrices. Figure 5a shows the compression ratio on the x-axis and the histogram recording the number of matrices using the four B2SR formats on the y-axis. In Figure 5b, the y-axis shows the number of matrices that belong to their: (1) optimal size (colored in blue): the least storage size required among the four B2SR formats of a matrix. (2) compressed size (colored in green): the B2SR format can provide a compression ratio < 100% for a matrix. For optimal, 162 matrices appear at B2SR-4, 291 matrices at B2SR-8, 26 matrices at B2SR-16, and 12 matrices at B2SR-32. For compressed, 491 matrices can have a compression ratio < 100% on B2SR-4, 421 on B2SR-8, 329 on B2SR-16, and 263 on B2SR-32.
There are multiple factors that contributed to the significant speedups B2SR has achieved. In addition to the gains by using native bit-level intrinsics such as __popc(), we have observed that more performance is from the reduced memory transactions and enhanced data locality. For example, for the matrix mycielskian8, by using B2SR, the number of global memory load transactions reduces by 4⇥ from 6630 to 1826, while the L1 cache hit-rate increases by 24% from 65.63% to 81.83%. We also observed different sweet areas for different B2SR tile sizes-such as, the smaller tile sizes (e.g., 4, 8) draw better L1 hit rates while the larger tile sizes (e.g., 32, 64) favor coalesced memory access-and also the impact from the profiles of individual graphs. In the next two sections, we describe the performance evaluation of the linear algebra kernels (BMV and BMM) and the five graph algorithms, respectively.
In this subsection, we evaluate the basic arithmetic cores BMV and BMM in terms of different schemes. We evaluate the speedups of the kernels over cuS-PARSE's SpMV (cusparseScsrmv()) and SpGEMM (cusparseScsrgemm()) with CSR in 32-bit floating-point nonzero storage. We compare the performance of each kernel scheme on B2SR-4, B2SR-8, B2SR-16, and B2SR-32. In Figure 6 and 7, the y-axis is the average speedups (of 5 runs) over cuSPARSE and the x-axis is the nonzero density of the matrices which is defined as # of nonzeros # of elements . A higher nonzero density (more to the right in the figure) implies a denser matrix while a lower one (more to the left) implies a sparser matrix.
BMV In BMV, we implement three schemes as the basic SpMV kernels (see Table II). We evaluate the performance of bmv_bin_bin_bin(), bmv_bin_bin_full(), bmv_bin_full_full() and compare it with cuSPARSE's cusparseScsrmv(). bmv_bin_bin_bin()'s performance is shown in Figure 6a and7a. Although the arithmetic capability of this scheme is limited to only binary operations, it allows a minimal memory footprint by keeping all value storage in binarized format. On GTX1080, bmv_bin_bin_bin() achieves an average speedup of 2.41⇥ in B2SR-4, 2.74⇥ in B2SR-8, 2.91⇥ in B2SR-16, and 2.85⇥ in B2SR-32. The max speedup over baseline is 40.47⇥ that appears at matrix ins2 with B2SR-32. On Titan V, bmv_bin_bin_bin() achieves an average speedup of 2.04⇥ in B2SR-4, 2.17⇥ in B2SR-8, 2.18⇥ in B2SR-16, and 2.26⇥ in B2SR-32. The max speedup over baseline is 25.16⇥ that appears at matrix ins2 with B2SR-32.
In bmv_bin_bin_full(), the vector input is binarized with the same tile dimension of the binary adjacency matrix. Compared to bmv_bin_full_full(), it requires less vector load bandwidth per matrix-vector multiplication. The inner product of each bit-row is done by bit-wise AND and population count of the resulting bit-row. In the 521 binary matrices, the runtime speedup of bmv_bin_bin_full() over cuSPARSE's full-precision CSR SpMV is shown in Figure 6b and 7b. On GTX1080, bmv_bin_bin_full() achieves an average speedup of 2.06⇥ in B2SR-4, 2.36⇥ in B2SR-8, 2.22⇥ in B2SR-16, and 2.97⇥ in B2SR-32. The For bmv_bin_full_full(), the performance is present in Figure 6c and7c. Unlike bmv_bin_bin_full(), the average performance gain decreases when enlarged the tile size. On GTX1080, it achieves an average speedup of 2.06⇥ in B2SR-4, 1.92⇥ in B2SR-8, 1.43⇥ in B2SR-16, and III). Figure 6d and 7d presents the performance of the BMM kernel compared to cuSPARSE's cusparseScsrgemm(). On GTX1080, it achieves an average speedup of 33.96⇥ in B2SR-4, 27.84⇥ in B2SR-8, 17.81⇥ in B2SR-16, and 10.22⇥ in B2SR-32. The maximum speedup is 6555⇥ that happens at matrix ins2 with B2SR-4. On Titan V, the performance gain is moderate compared to GTX1080. We accounts the reason for cuSPARSE's APIs have better performance gain on Volta than Pascal, while our implementation perform similar or evenly slightly poor on Volta than on Pascal. bmm_bin_bin_sum() on Titan V achieves an average speedup of 5.34⇥ in B2SR-4, 3.65⇥ in B2SR-8, 9.03⇥ in B2SR-16, and 12.25⇥ in B2SR-32. Interestingly, the most significant speedup is 3243⇥ also happens at matrix ins2 but with B2SR-32 instead of B2SR-4 as on GTX1080.
We compare the implemented B2SR-based Bit-GraphBLAS algorithms with GraphBLAST [4], a state-of-the-art GPUbased GraphBLAS framework. GraphBLAST switches between sparse and dense vector/matrix computation depending on sparsity degree across algorithm interactions with optimized CUDA kernels. For the GraphBLAST configuration, we use the same environment (e.g., CUDA Runtime 9.1) as indicated on the Github repository [49]. The graph algorithms are implemented following the GraphBLAS convention. For iteration-based algorithms, such as BFS, SSSP, PR, and CC, each iteration contains a frontier vector that performs a matrixvector multiplication with the adjacency matrix and several element-wise scalar operations. They are used to update the frontier vector for indicating neighbor aggregation in each iteration, through the mathematical semi-ring operation. The number of iterations depends on when the algorithm is converged at runtime.
Since matrix-vector multiplication is the major performance concern per iteration (>80% of the workload), in Bit-GraphBLAS, our major goal is to improve the efficiency of SpMV through our B2SR based BMV kernel. Table VII and VIII list the algorithm and kernel execution latency in ms for the proposed Bit-GraphBLAS with respect to GraphBLAST. As can be seen, under both conditions, Bit-GraphBLAS achieves considerable speedups through B2SR and the strong bit computation capability of modern GPUs.
For Bit-GraphBLAS, in the 521 binary matrices dataset, the patterns with better performance fall into three main categories: diagonal, block, and stripe (reference the classification in Table V). The performance of matrices from the three categories are shown in Table VII and VIII. In the subset of matrices, delaunay n14, se, and debr are stripe patterns; Erdos02, mycielskian9, EX3, net25, and myciel-skian10 are block patterns; ash292, netz4504 dual, minnesota, jagmesh6, uk, whitaker3 dual, rajat07, 3dtube are diagonal patterns. In the algorithm evaluation, BFS relies on the kernel bmv_bin_bin_bin_masked() with boolean semiring. On Pascal, diagonal pattern matrices can achieve up to 433⇥ acceleration in the whole algorithm and 1414⇥ in kernel; On Volta, it achieves up to 349⇥ speedup to GraphBLAST in algorithm and 1344⇥ in kernel. Matrices in stripe and block generally perform a moderate speedup from 2⇥ to 7⇥ on both GPU architectures. Other three SpMV-based algorithms (SSSP, PR, and CC) are fulfilled by bmv_bin_full_full() with relaxation and extension (details are described in Section V). They mainly achieve a acceleration over GraphBLAST with no more than 20⇥ algorithm-wise and 40⇥ kernel-wise.
In Table IX, we demonstrate the performance improvement of TC algorithm, which is essentially a one-time execution of the bmm_bin_sum_masked() kernel. On both Pascal and Volta, we still have diagonal patterns with the best performance. It can achieve up to 52⇥ on Pascal and 27⇥ on Volta GPUs.
It is noteworthy that the same matrix can often find lower GraphBLAST runtime on Volta than on Pascal. A similar effect can be found in both kernel and algorithm evaluations. For example, the 3dtube runs 151.89 (ms) on Pascal but only 79.49 (ms) on Volta. Nevertheless, Bit-GraphBLAS can have a larger or similar runtime on Volta than on Pascal GPUs. We have a 0.04 (ms) increase in runtime for the 3dtube case. We attribute the effect to that Volta has updated the warp execution model and eliminated implicit warp synchronous. This poses a little performance slowdown of binary intrinsics like __shfl_sync() and __ballot_sync() compared to the non-synchronizing __shfl() and __ballot() in Pascal GPUs.
Limitations of Bit-GraphBLAS: As Bit-GraphBLAS relies on the bit-operation and bit-data, it only directly applies to homogeneous graphs (i.e., the adjacency matrix is a binary matrix). Throughout the SuiteSparse data collection, we have observed that ⇠20% of graphs are homogeneous that can directly benefit from Bit-GraphBLAS. These graphs cover a wide range of domains including mathematics, power-grid, physics, electronics, material science, economics, thermal, fluid dynamics, etc. Additionally, as the weights for many heterogeneous graphs can be expressed by integers or fixedpoints (e.g., through normalization), similar to the recent effort decomposing a quantized-neural-network into several concurrent binary-neural-networks for acceleration [39], Bit-GraphBLAS can also be extended to support heterogeneous graphs with short bit-width. We set this as future work.
Platform Portability: Although in the evaluation we showcase Bit-GraphBLAS on NVIDIA GPUs due to hardware availability, the bit intrinsics that Bit-GraphBLAS relies on, such as popc(), shfl(), ballot(), and brev() are also available (despite using different names) in other GPUs and CPUs (like AMD's GPU and X86 CPUs). Therefore, we did not see significant challenges in supporting Bit-GraphBLAS on alternative hardware platforms (e.g., AMD GPUs through HIPIFY [50]).
We present Bit-GraphBLAS, a linear algebra-based graph framework that utilizes a Bit-Block Compressed Sparse Row (B2SR) format and bit manipulation primitives on GPUs to enable dense bit-operations on bit tiles within large sparse adjacency matrices. We explore different tile size configurations from 4⇥4 to 32⇥32, and suitable bit-packing types. We implement BMV and BMM schemes to support parse kernel operations and demonstrate their effectiveness on graph algorithms by reducing memory footprint. The result shows significant performance gain over full-CSR-based GPU graph frameworks. In sum, the novel storage format and algorithms compress the graphs' storage and accelerate the linear algebra kernels SpMV and SpGEMM through finer-grained bit-wise parallelism.