Graph Traversal on Tensor Cores: A BFS Framework for Modern GPUs

Unknown

Conceptio › Archive › arXiv CS

arXiv CSopen access

Graph Traversal on Tensor Cores: A BFS Framework for Modern GPUs

Unknown · 2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

cloud distributedcomputing parallelcomputing

distributed computing, parallel computing, cloud

Graph Traversal on Tensor Cores: A BFS Framework for Modern GPUs Deniz Elbek and Kamer Kaya

arXiv:2606.05081v1 [cs.DC] 3 Jun 2026

✦

Abstract—Modern GPUs have Tensor Cores (TCs) capable of extremely high-throughput matrix operations, yet graph algorithms remain difficult to accelerate because of their irregular and data-dependent execution patterns. This work presents B LEST, a TC–accelerated framework that reformulates Breadth-First Search (BFS) as a bit-level sparse matrix–vector computation while addressing the load imbalance, memory inefficiency, and synchronization overheads that limit prior approaches. B LEST introduces Binarized Virtual Slice Sets (BVSS), a graph representation that partitions work into balanced warp-level units and schedules only frontier-relevant regions of the graph. It further employs an optimized TC layout that maps neighbour checks onto binary MMA instructions without wasted outputs, reducing the number of required MMA calls by 8× compared with prior layouts. To mitigate atomic and cache bottlenecks, B LEST incorporates a lazy vertex-update scheme. We revisit the switching terminology for BFS and propose a mechanism that dynamically transitions from TCs to CUDA cores when it becomes more efficient. We also extend B LEST to multi-source BFS and closeness centrality workloads. Finally, we introduce a scalable graph reordering method that improves compression for scale-free-like graphs, while using RCM to improve locality for others. Across a broad set of real-world graphs, B LEST achieves average speedups of 22.0×, 7.7×, 8.1×, and 5.9× over GAP, Gunrock, GSWITCH, and BerryBees, respectively, establishing a new BFS baseline on GPUs. Thanks to its high performance, B LEST can compute the exact closeness centralities of 65.6M vertices in a social network with 3.6B edges in an hour using 100 H100 GPUs. Index Terms—BFS, GPUs, Tensor Cores, Sparse Matrix–Vector Multiplication, Multi-Source BFS, Direction Switching, Closeness Centrality.

I NTRODUCTION

ODERN GPUs have specialized matrix-multiplyaccumulate (MMA) units, i.e., Tensor Cores (TCs) that deliver high throughput for matrix operations. These units are designed for dense and highly regular workloads. Graph algorithms, however, live at the opposite end of the spectrum: they are driven by sparse, irregular, and data-dependent access patterns, and their performance is constrained not by peak arithmetic throughput but by memory behavior, synchronization, and load imbalance. Therefore, for graphs, the main question is not whether MMA units are fast, but under what conditions their dense computational model

1) Deniz Elbek ([email protected]) and Kamer Kaya ([email protected]) are with the Department of Computer Science and Engineering, Faculty of Engineering and Natural Sciences, Sabancı University, Istanbul, Turkey, and also with the Center of Excellence in Data Analytics, Sabancı University, Istanbul, Turkey. 2) Kamer Kaya is the corresponding author of this work.

can be exploited despite the irregularity inherent in stateof-the-art algorithms. In this work, our main focus will be on Breadth-First Search (BFS), however, we will also briefly discuss the use of TCs for other (iterative) graph algorithms. As one of the most important algorithms on graphs, BFS has applications across a broad range of domains— including network science [1], computer networks [2], recommendation systems [3], and compiler design [4]—where the BFS performance has a direct, end-to-end impact on the overall application throughput. Therefore, it is a popular research area in high-performance computing (HPC), and its performance has been optimized both algorithmically [5], [6], [7] and implementation-wise across many systems, including CPUs [6], [8], GPUs [9], [10], [11], and clusters [12], [13]. A single BFS can be represented as an iterative sparse matrix–sparse vector multiplication (SpMSpV) [14], where the matrix corresponds to the transpose of the bit-adjacency matrix of the underlying graph, and the vector is a bitfrontier vector with set bits identifying the nodes in the current frontier. Each SpMSpV produces an output frontier (bit) vector for the next level. The process continues until the frontier is all-zero, i.e., no reachable unvisited node remains. To harness the TCs, the multiplication must be recast as a specialized form of GEMM. Recent work [15] has introduced a methodology that successfully integrates the TCs into BFS pipelines and achieves speedups over the stateof-the-art across a diverse set of graphs. Nevertheless, several inefficiencies and extensions remain, which we address in this work. Beyond proposing the fastest BFS framework to date, we systematically analyze using TCs for graph algorithms. Our contributions are 5-fold: 1. We propose B LEST, which, to the best of our knowledge, is currently the fastest BFS implementation in the literature. • Thanks to its novel structure, Binarized Virtual Slice Sets (BVSS), it has near-perfect inter-warp load balance. • B LEST employs a TC multiplication layout that optimally maps the required edge operations onto these dense units, reducing the number of required MMA calls by a factor of 8× compared to state-of-the-art. • It introduces a lazy update scheme that (i) reduces the inherent cost of atomics and (ii) significantly improves cache locality across all levels of the memory hierarchy. • B LEST fuses the per-level BFS kernels into a single persistent kernel to eliminate host-side synchronization and kernel-launch overhead.

2. We reconsider switching in the TC era to specialize each BFS level execution, and enhance B LEST to perform a transition between TCs and CUDA cores when necessary. 3. We extend B LEST to multi-source BFS workloads, altering BVSS to support that transition as efficiently as possible. We evaluate these workloads on a representative multisource BFS application: closeness centrality. 4. We introduce J ACCARD W ITH W INDOWS, a novel yet simple reordering algorithm applicable to the largest realworld graphs. We use it for scale-free-like graphs to increase the effective compression ratio of BVSS. For other graphs, we use a traditional reordering algorithm, Reverse Cuthill–McKee (RCM) [16], to reduce L1 cache miss rates and excessive L2 cache line fetches. 5. We evaluate single-/multi-source B LEST on a broad set of real-world graphs and achieve average speedups of 21.96×, 7.70×, 8.12×, and 5.93× over GAP [17], Gunrock [10], GSWITCH [11], and BerryBees [15], respectively. The article is organized as follows. Section 2 provides background on BFS, TCs, and introduces the notation used throughout the paper. Section 3 introduces BVSS, our novel data structure for BFS on TCs. Section 4 presents the two reordering strategies employed in this work. Section 5 describes B LEST’s compute pipeline for single-source BFS and Section 6 extends it to multi-source BFS workloads. Both pipelines are evaluated extensively in Section 7. Related work is discussed in Section 8, and Section 9 concludes the article.

BACKGROUND AND N OTATION

Let G = (V, E) be a directed graph with n = |V| vertices and m = |E| edges, and let A ∈ {0, 1}n×n denote the transposed adjacency matrix, so that A[i][j] = 1 if and only if (j, i) ∈ E . Let x(k) ∈ {0, 1}n be the bit-frontier vector at BFS level k , (k) where xi = 1 if vertex i belongs to the frontier at level k . Wk Let v(k) = ℓ=0 x(ℓ) denote the visited vector after level k . The per-level BFS computation advances the frontier via x(k+1) = A · x(k) ∧ ¬v(k) ,   n−1 n−1 (1) _ (k) (k) =  A[i, j] ∧ xj  ∧ ¬vi  , j=0

i=0

where the multiplication operates over the Boolean semiring (∨, ∧); i.e., the dot product over (+, ×) is replaced by (∨, ∧), and already visited vertices are filtered from the result. The expression evaluates to 1 if and only if there exists at least (k) one incoming neighbour j of i such that A[i][j] = 1, xj = 1, (ℓ) and i has not been visited before, i.e., xi = 0 for all ℓ ≤ k . For a source s ∈ V , the recurrence in (1) is initialized with (0) (0) xs = 1, and xi = 0, ∀i ̸= s, and terminates at the first ∗ ∗ level k for which x(k ) = 0. The algorithm returns an array (k) level : V → Z≥0 , where level[i] = k such that xi = 1. (k) For an unreachable i, level[i] = ∞, so xi = 0 for all k . 2.1

GPU Tensor Core Architecture

GPUs comprise multiple compute units called Streaming Multiprocessors (SMs), each subdivided into SM SubPartitions (SMSPs). On recent NVIDIA GPUs, each SMSP

contains one TC, and a warp—a team of 32 threads executing the same instruction under the Single Instruction Multiple Thread (SIMT) paradigm—scheduled onto that partition can dispatch its work onto the partition’s private TC. TCs support multiple precisions and even structured sparsity configurations; however, their sparsity support imposes strict structural requirements not suitable for the irregular sparsity patterns of graphs. Still, TCs have found their way into a variety of domains that they were not originally designed for, including SpMV [18], [19], SpMM [20], [21], [22], [23], [24], [25], [26], [27], [28], [29], [30], [31], [32], [33], [34], triangle counting [35], quantized SpMM [36], fully-homomorphic encryption [37], epistatic detection [38], graph neural networks [39], [40], sparse deep learning [41], [42], [43], [44], sparse LLM inference [45], [46], [47], [48], stencil [49], reduction [50], and more recently BFS [15]. TCs can operate on the bit-level: given two k -bit vectors, a TC computes their product over the (popc, ∧) semiring, where popc returns the number of set bits produced by the bitwise AND of its two operands. Although Eq. (1) requires the Boolean semiring (∨, ∧), it can be faithfully emulated via the (popc, ∧) semiring supported by the hardware. Specifically, let ai ∈ {0, 1}n denote the i-th row of A; then h i (k+1) xi = 1 popc ai ∧ x(k) > 0 . (2) (k) Theemulation is exact: popc(ai ∧ x ) > 0 if and only if W (k) = 1. Here, the operands are not plain k -bit j ai [j] ∧ xj vectors, but m × k and k × n matrices, whose dimensions depend on the chosen precision configuration. Three binary MMA shapes are supported on recent NVIDIA architectures: m8n8k128, m16n8k128, and m16n8k256 [51]. In this work, we utilize the smallest available configuration, m8n8k128, to achieve the finest computation granularity.

2.2

Directions, Updates, and Switching in BFS

BFS is driven entirely by information transmission: the visited state of the vertices in the current frontier is passed to next frontier vertices. Following the terminology of Beamer et al. [6], in a level-synchronized BFS, this transmission can happen in two directions: (1) In top-down exploration, the vertices in the current frontier transfer their visited states to their unvisited (outgoing) neighbours. (2) On the other hand, in the bottom-up exploration, the unvisited vertices transfer the visited states from their (incoming) neighbour vertices in the current frontier. Hence, in top-down, the work queue is the current frontier, and in bottom-up, it is the set of unvisited vertices. The literature calls the former exploration direction push-based since the information is pushed from the current frontier to the next one. Similarly, the latter exploration is called pull-based. This coupling sounds natural since, in a top-down BFS, a thread naturally pushes updates to a shared state, whereas in a bottom-up one, it pulls the (incoming) neighbours’ states to its private state. In this work, we decouple the bottom-up mechanism from pull-based exploration, as well as top-down from push-based, yielding four distinct modes of a parallel BFS. Although inefficient on CUDA cores, pull-based update mechanics can be practically coupled with a top-down work queue, thereby removing the assumption that these terms are

Preprint submitted to IEEE Transactions on Parallel and Distributed Systems

4 x τ = start rowIds and masks

BVSS Data Structure Slice 0 will Sets be active

3 01010010

is active

virtualToReal 0 0 0

61 10010110

133 00000010

183 11110000

4 74 10001010 00011010

141 11001100

207 00001111

...

60 00011110

120 11111111

174 10000010

313 01100110

fra gB

gA fra

fra gB

gA fra

rowIds and masks (4 + 1) x τ = end

Virtual Slice 0 Sets

5 Slice 0

Q curr 4 Size

-1

Next Queue

Current Queue

Fcurr Qnext 0 Size

-1 -1

3 atomic add

0....0

1 0 ... 0

Slice 0

Fnext 00011000

n/8 - 1

11000000 .....

0 ... 0

3 01010010

0 ... 0

11000000

coming from the next frontier

t31 realPtrs 0 3 4 6

TC Compute

Processed by Warp ID: 0

n/8 - 1 .....

0 ... 0

if not visited

atomic or

pulling vertex: 3, result 1

Fig. 1: BVSS data structure and the flow of data reads, pull operations, and frontier updates: Slice set 2 (in green) is active, with Fcurr bits 11000000, i.e., the 1st and 2nd columns of slice set 2 are in the current frontier. Since B LEST’s queues operate over VSS indices, it retrieves the two VSSs corresponding to slice set 2, namely VSSs 4 and 5, directly from the current queue Qcurr (bottom left). For simplicity, the figure focuses on VSS 4 whose slices (rowIDs and masks) are read from memory and assigned Warp 0. The masks and the σ frontier bits of slice set 2 (accessed via virtualToReal from VSS 4 and then from Fcurr ) are fed to the TC as fragA and fragB, respectively. Pulls on the Boolean semiring are processed in two rounds, each processing half of the slices (colored purple and orange). For the slice with vertex/row 3 in the first half, the mask is 01010010, indicating that vertex 3 is an outgoing neighbour of the 2nd, 4th, and 7th vertices/columns of slice set 2. Since the 2nd vertex is active in the frontier, the popcount is nonzero, and an update is required (vertex 3 is unvisited). Since vertex 3 belongs to slice set 0 (covering vertices 0–7), its 4th bit is set in Fnext . Finally, B LEST locates all VSSs corresponding to slice set 0 via realPtrs and inserts their IDs (0, 1, and 2) into Qnext .

interchangeable; in its simplest form, a thread responsible for updating a private state can go over all the frontier queue vertices and check their neighbourhoods to see if its state variable needs to be updated. Since pull-based exploration admits a matrix multiplication formulation suitable for TCs, we adopt it; however, we maintain a work queue containing all the frontier vertices (and some more), which is common for the top-down exploration in the state-of-the-art. The literature shows that switching the direction, and hence the update mechanism, during an iterative graph algorithm pays off well not only for a BFS but also for others [52]. Adding the TCs to the list of available units yields another switching flexibility; an iteration can be processed not only by the CUDA cores but also by TCs. Furthermore, a hybrid processing is also possible, an avenue that has been explored for SpMM [30], [23], [32], [33]. Hence, as shown in this work, existing BFS frameworks that successfully apply excellent switching mechanisms, yet use only CUDA cores, such as GSWITCH [11], may not fully utilize the recent GPUs.

(k+1)

sets level[i] = k + 1 and xi = 1 if the mask-frontier multiplication produces a set bit, and i has not been visited at any prior level. The SotA is inefficient due to four factors: 1. Assigning one slice set per warp exposes an inter-warp load imbalance due to skewed degree distributions. 2. Slice sets are dispatched in a frontier-oblivious manner. The partition x(k) [sσ : (s + 1)σ] may be all zero, yet the assigned warp becomes active and must work. Even with early exits from warps due to all-zero frontiers, a warp may exit frequently while others do not. 3. When the information ratio in a mask, i.e., popc(mask) / σ , is low, the TCs perform redundant operations. (k+1) 4. The updates to level[i] and xi can cause high L1 cache miss rates and excessive L2 cache line fetches, respectively, due to scattered row IDs within a set. We address the first two problems in Sec. 3.1, where we propose our novel data structure. We address the third one in Sec. 4.1 and the last one in Secs. 4.2 and 5.2. 3.1

B LEST DATA STRUCTURE : BVSS

The Binarized Row Slice (BRS) structure [15] partitions A along its column dimension into sets of width σ , such that a nonzero A[i][j] falls in slice set s if sσ ≤ j < (s+1)σ . Row i is included in slice set s if it contains at least one nonzero in the column range [sσ, (s + 1)σ). Each such row is represented as a slice, a row ID and a σ -bit mask, where mask[j−sσ] is set for every nonzero A[i][j]. During computation, each slice set is assigned to a single warp, which assembles an intermediate matrix from the masks of that slice set and multiplies it against the corresponding partition x(k) [sσ : (s + 1)σ) of the current frontier vector. For a row i processed, the warp

Binarized Virtual Slice Sets (BVSS)

Let θ = 32/σ denote the number of slices assigned to each thread within a warp, and let τ = 32θ denote the total number of slices constituting one unit of warp work. Motivated by this, we partition each slice set into a variable number of virtual slice sets (VSS), each carrying at most τ slices. The upper-left portion of Figure 1 illustrates this partitioning: the three slice sets are split into three, one, and two VSSs, respectively, implying that the first slice set contains [2τ +1, 3τ ) slices, the third contains [τ +1, 2τ ) slices, and the second contains fewer than τ slices. Let Ns and Nv denote the total number of slice sets and VSSs, respectively. The mapping from slice sets to their virtual children is

Preprint submitted to IEEE Transactions on Parallel and Distributed Systems

maintained through an array realPtrs of size Ns + 1, where realPtrs[s + 1] − realPtrs[s] gives the number of VSSs emanating from slice set s. The inverse mapping is provided by virtualToReal of size Nv , where virtualToReal[v] identifies the parent slice set of VSS v . Throughout this work, we fix σ = 8, yielding θ = 4 and τ = 128. Consequently, each thread handles the multiplication of a 32-bit mask in collaboration with the rest of its warp, and is responsible for the update mechanism of 4 rows. VSSs containing fewer than τ slices—of which there are at most Ns such sets, at most one per slice set—are zero-padded with all-zero masks and arbitrary row IDs to simplify indexing: both row IDs and masks are accessed at position vτ +laneID·θ, where laneID is the thread index within the warp processing virtual slice set v . Consequently, each rowIds element packs four 4-byte rows, and each masks element holds one 32-bit mask. No VSSs are constructed for slice sets that contain no slices. We name this data structure Binarized VSSs (BVSS). Contrary to SotA [15], we maintain a work queue that contains only the VSSs that must be processed at the current BFS level, i.e., those whose parent slice set covers at least one column that belongs to the current frontier. In Fig. 1, the third slice set is active, and all of its child VSSs {4, 5} reside in the queue simultaneously. Taking VSS 4 as an example, the warp assigned to it resolves virtualToReal[4] = 2 and accesses the corresponding 8-bit frontier partition Fσcurr [2]1 . The multiplication of the τ = 128 slices assigned to this warp is performed in two rounds of TC multiplications; the exact layout of the matrix constructed for each round is detailed in Sec. 5.1. When a multiplication yields a set bit, and the pulling vertex has never been visited before, as visualized in the bottom-right portion of Fig. 1, Fσnext [0], containing the bit for the pulling vertex 3, is atomically updated. When a row within slice set s is discovered for the first time, all VSSs in the range [realPtrs[s], realPtrs[s + 1]) are enqueued into the next-level work queue. The algorithm terminates when the next-level queue is empty, i.e., when Fnext is all-zero. BVSS immediately resolves the first two problems stated in Sec. 3. Inter-warp load balance is resolved to near-perfect: any two warps selected from the grid can differ in the amount of work they perform by at most one virtual slice set, each containing a fixed and controlled number of slices. Such residual imbalance arises only when |Qcurr | mod P ̸= 0, where P is the number of warps in the grid, and is bounded by at most two TC multiplications across any two warps. The second deficiency, frontier-oblivious warp scheduling, is also resolved. Warps are dispatched only for VSSs whose parent slice set covers at least one active frontier column. Consequently, every item submitted to the work queue corresponds to a genuinely active region of the graph through which BFS is transmitting information at the current level. 3.2

Top-Down Exploration with Pull-Based Mechanics

The work-queue dynamics of the pipeline depicted in Fig. 1 align with a top-down exploration; the work is determined by the vertices in the current frontier (as in Fig. 2a). To be exact, B LEST processes the (outgoing) edges of the vertices 1. Throughout the paper, array subscripts indicate access granularity: for X ∈ {Fcurr , Fnext , Vcurr , Vnext }, X 1 [u], X σ [s], X 32 [w], and X κ [·] denote vertex, slice-set, 32-bit-word, and κ-bit accesses, respectively.

(a) All the outgoing edges of frontier vertices (pushed by two threads for j) Levels [0, k-1] Level k frontier (work queue) j .....

thread

.....

....

i’’ 01101000

00101100

01011101

xk 1 1 0 0 1 0 0 1

.....

01001111

....

thread

00100011

i’ 01101000

.....

thread

VSS

Unvisited vertices

01011101

(τ slices)

(b) All the incoming edges of unvisited vertices (pulled Levels by two threads [0, k-1] for i) Level k Frontier

window (σ vertices)

thread

window (σ vertices)

Frontier (work queue)

i x

k+1

Unvisited vertices (work queue)

edges of i

Outgoing edges of level k frontier vertices (and some more) offloaded to TC by the threads Filtered by the TC Reduced by the TC For each window, a single thread writes the information.

Fig. 2: The traditional a) top-down (push-based) and b) bottomup (pull-based) BFS-level processing. c) B LEST performs topdown exploration with a pull-based update mechanism.

that share a slice with at least one frontier vertex (as in Fig. 2c). However, it employs a pull-based approach; each thread is responsible for updating 4 (row) vertices. For each, the (incoming) neighbours determine whether the vertex resides in the frontier. In addition, since the (incoming) edges of a vertex i may be scattered to multiple slice sets, multiple threads can update the state private to i.

R EORDERING THE G RAPH

To reduce the total number of slices, i.e., to compress BVSS, and to resolve locality issues during updates, B LEST symmetrically reorders the graph G. We identify that graphs obeying power-law properties are highly compressible and introduce a scheme that increases the effective compression ratio by exploiting communities within slice sets. For graphs that are not scale-free-like, which are relatively less compressible, Reverse Cuthill–McKee (RCM) [16] is used to reduce the bandwidth of A, thereby better confining the updates. 4.1

Increasing Compression Ratio

In BVSS, a slice exists within a set when the row it represents has at least one incoming neighbour among the σ columns that the slice set covers. Consequently, the effective compression P ratio of a slicewith row ID i residing in slice σ−1 set s is j=0 A[i][sσ + j] /σ . To increase this ratio, one must construct slice sets such that the σ columns extensively overlap among their outgoing row neighbours. Jaccard similarity [53] is a popular approach to assess the overlaps between neighbourhoods. For two columns j and j ′ , the similarity is defined as

J (j, j ′ ) =

| nbrsA (j) ∩ nbrsA (j ′ ) | . | nbrsA (j) ∪ nbrsA (j ′ ) |

(3)

Although computing all pairwise Jaccard values may produce desirable slice sets for BVSS, its complexity is O(n2 δ),

Preprint submitted to IEEE Transactions on Parallel and Distributed Systems

where δ denotes the maximum column degree in A, making it infeasible for real-world graphs. A basic observation for B LEST is that only the σ columns within each slice set need to form a community; inter-slice-set column affinity has no impact on compression. Motivated by this, we restrict the computation within windows of vertices and propose a novel reordering, J ACCARD W ITH W INDOWS, given in Algorithm 1. The algorithm partitions A column-wise into nonoverlapping n/W windows of width W > σ where W mod σ = 0. Each window contains W/σ slice sets of size σ ; columns are greedily assigned to slice sets one at a time. For each slice set, a highest-degree unassigned column within the window is selected as a singleton seed j ∗ , whose outgoing neighbour set initializes R, with cost O(W ) (line 6). The inter[j] values used to form the slice sets are initialized with cost O(δW ) (line 9). The remaining σ−1 slots are filled by repeatedly selecting the column j ∈ Q—the set of unassigned columns within the window that share at least one row neighbour with R— that maximizes inter[j]/(|R| + degA (j) − inter[j]), where inter[j] is the number of common vertices within nbrsA (j) and R. When Q is empty, a highest-degree unassigned column is again chosen as a fallback (line 13). Upon each selection, R is extended with the newly covered rows, and inter[·] is incrementally updated via A⊤ for all column candidates in Q. Each slice set construction performs σ scans of Q of size at most W (lines 12–15). For the inter[·] updates, each selected column introduces δ new rows into R, and each such row contributes at most min(δ, W ) column candidate updates. Hence, the update cost is W n × × (δW + σ(W + δ min(δ, W )) (4) O W σ δW =O n× + W + δ min(δ, W ) (5) σ in total. With W ≪ n, the overall complexity is substantially below the O(n2 δ) of the naive implementation. We evaluate the window size W in Section 7 and how it affects both the execution time and the compression rate of the algorithm. 4.2

Boosting Cache Locality

In a GPU, each SM has a private L1 cache, whereas the L2 cache is shared among all SMs. All atomic stores are therefore coherent in L2, enabling concurrent memory updates to be visible across all SMs. When a TC multiplication produces a nonzero result for row/vertex i that has never been visited before, two updates indexed by row ID i occur: (i) level[i] = atomic k , and (ii) F1next [i] ← 1. The efficiency of these updates depends entirely on the distribution of row IDs within a VSS: to have high L1 and L2 cache hit rates, respectively, for level[i] updates and F1next [i] atomic stores, the row IDs constituting a VSS must be tightly clustered. For graphs that are not scale-free-like such as road networks, reordering G⊤ so that rows within a VSS are highly local is both feasible and effective. To quantify this locality, we introduce a metric called update divergence. For each of the θ columns of the VSS matrix shown in Fig. 1, we define the column divergence as the standard deviation of the row IDs with nonzero masks. A higher standard deviation indicates that the row IDs are more scattered and less clustered, whereas a lower value

Algorithm 1 J ACCARD W ITH W INDOWS Require: A ← G⊤ , A⊤ ← G, slice size σ , window size W Ensure: Inverse permutation π −1 1: for w = 0, . . . , ⌈n/W ⌉ − 1 in parallel do 2: ws ← wW ; we ← min(ws + W, n) 3: C←∅ 4: for s = 0, . . . , (we − ws )/σ − 1 do 5: ss ← ws + sσ ; se ← min(ss + σ, we ) 6: j ∗ ← arg maxj∈[ws , we )\C degA (j) 7: C ← C ∪ {j ∗ }; π −1 [j ∗ ] ← ss 8: R ← nbrsA (j ∗ ) 9: ∀ j ∈ [ws , we ) \ C : inter[j] ← |R ∩ nbrsA (j)| 10: Q ← {j ∈ [ws , we ) \ C : inter[j] > 0} 11: for ℓ = ss + 1, . . . , se − 1 do 12: if Q = ∅ then 13: j ∗ ← arg maxj∈[ws , we )\C degA (j) 14: else inter[j] 15: j ∗ ← arg maxj∈Q |R| + degA (j) − inter[j] 16: C ← C ∪ {j ∗ }; Q ← Q \ {j ∗ }; π −1 [j ∗ ] ← ℓ 17: for i ∈ nbrsA (j ∗ ) \ R do 18: R ← R ∪ {i} 19: for j ∈ nbrsA⊤ (i) ∩ ([ws , we ) \ C) do 20: if inter[j] = 0 then 21: Q ← Q ∪ {j} 22: inter[j] += 1 23: return π −1

reflects tighter locality. The set divergence of VSS s, denoted Udiv (s), is then the average column divergence over its nonempty columns. Finally, the update divergence is the mean Udiv = avg Udiv (s) computed over all Nv VSSs, where a s

lower Udiv indicates better locality. Reverse Cuthill-McKee (RCM) [16] assigns vertex IDs based on their visit order in a BFS-like traversal; vertices with the same parent are ordered by ascending degree, delaying frontier explosion and keeping connected vertices close. Table 1 reports Udiv for four graphs—two road networks (GAP-road and europe_osm), a planar, triangulation network (delaunay_n24), and one random network (rgg_24)—before and after RCM. For all, the dramatic reduction in Udiv confirms that RCM substantially tightens the clustering of row IDs within VSSs.

Graph GAP-road delaunay_n24

Udiv (unord.) 158,696 894,542

Udiv (ord.) 869 2,839

Graph europe_osm rgg_24

Udiv (unord.) 869,302 3,948

Udiv (ord.) 982 1,512

TABLE 1: Avg. update divergence Udiv before and after RCM. Unfortunately, RCM is not effective for scale-free-like graphs, which do not generally possess a recoverable low bandwidth. During preprocessing, we therefore first identify whether G is scale-free-like2 . If so, we apply J ACCARD W ITH W INDOWS (Section 4.1), which aims to maximize mask density within each slice set. Otherwise, we apply RCM to reorder G⊤ , which reduces Udiv and thereby improves 2. We classify a graph as scale-free-like if either its in- or out-degree distribution is heavy-tailed: the top 1% and 10% of vertices account for at least 5% and 40% of total degree, respectively, or a log-log degreehistogram fit for k ≥ 5 has slope −γ , R2 ≥ 0.70, and γ ∈ [1, 5].

Preprint submitted to IEEE Transactions on Parallel and Distributed Systems

10 14 18 22 26 30

11 1

27 3

31 3

64 96

34 66

10 11

10 10

....

19 23

13 14 15

12 12

13 13

14 14

18 19

16 16

18 18

19 19

....

σ=8

20 21

22 23

20 20 21 21

22 22

23 23

24 25 26 27

24 24 25 25 26 26

28 29 30 31

28 28 29 29 30 30 31

....

24 56 88

.... 31

63 95

VSS

35 ....

4 cols (each holding 2σ useful bits) σ σ

...

28 60 29 61 30 62 31 63

31 63 64 96 65

.... 95

fragC a) Tensor core computation (numbers are thread IDs)

8x8 popcount matrix 8 rows

2 3

21 25 29

Round 1

Round 2

fragB 3

16 20 24 28

3 0 1 2 0 1 0 0

σ σ 64 96 65

... 92

... -

Y Z

...

fragB Y Z

Z X

0 5 2 1 0 0 1 0

Y Z

... -

fragA

Shared fragB

Z X

fragA

Y Z

5 0 3 0 0 0 0 0 ... 2 0 0 0 0 0 4 0

Intermediate

fragC

b) BLEST VSS processing starting from slice-to-thread distribution (numbers are slice indices for VSS and fragA, frontier masks for fragB, pop-counts for fragC ).

Fig. 3: (a) The data layout for m8n8k128 on Tensor Cores. For fragA, fragB, and fragC, each box corresponds to a 32-bit word, and the number inside is the ID of the thread holding that word in its registers. (b) B LEST’s VSS processing pipeline starting from the 32 × 4 VSS matrix (numbers denote slice indices for the VSS and fragA, frontier masks for fragB, and popcounts for fragC). The pull operations for a VSS are handled in two rounds, each issuing one m8n8k128 instruction. In each round, two columns of the 32 × 4 VSS matrix are reshaped into a 32 × 2 intermediate matrix, then conceptually into the 8 × 8 fragA matrix (dashed lines), and finally packed into the 8 × 4 fragA matrix, where each 32-bit entry holds two σ -bit masks in its first 2σ bits. A single 4 × 8 fragB is shared across both rounds: letting α be the σ -bit frontier word and ϑ the all-zero σ -bit word, the 16-bit building blocks are X = αϑ, Y = ϑα, and Z = ϑϑ. All fragB words are zero except those held by threads 0, 4, 9, 13, 18, 22, 27, and 31, which carry XZ or Y Z to select the first or second packed slice mask, respectively. This layout is optimal: no fragC popcount is wasted, and every result is delivered directly to the thread that owns the corresponding slices, requiring no intra-warp communication.

cache efficiency during warp-level updates. Although Udiv measures per-virtual-slice-set locality, RCM also improves intra-virtual-slice-set locality: row IDs assigned to adjacent virtual slice sets are numerically proximate, so warps processing neighbouring VSSs also operate on clustered row IDs. Specifically, the level[i] update benefits when those warps are co-located on the same SM, since the relevant cache lines are then likely to reside in the shared L1 cache; the F1next [i] atomic update benefits unconditionally, as the L2 cache is shared across all SMs and clustered row IDs translate directly to higher L2 cache hit rates regardless of warp placement.

C OMPUTE M ECHANICS OF B LEST

As stated before, mapping the edge operations within a BFS onto TCs is non-trivial. Specifically, the matrix–vector primitive that underpins single-source BFS must be recast as a specialized form of matrix–matrix multiplication that these dense units are built for. Even a layout that efficiently achieves such a recasting is not sufficient to extract peak hardware performance from a BFS kernel, which is widely regarded as a memory-bound algorithm. To address both concerns, we first introduce our novel TC multiplication layout in Sec. 5.1, which optimally maps bit sparse-matrix– sparse-vector multiplication onto the dense units, and then in Sec. 5.2, we propose a lazy vertex-update scheme that substantially alleviates the inherent memory inefficiencies. 5.1

An Optimal Layout for TC Multiplication

Fig. 3 depicts the TC multiplication layout used in B LEST. In each of the two rounds per VSS, the warp processes 2/4 columns of the 32 × 4 VSS matrix. These two columns form a 32 × 2 intermediate matrix, conceptually reshaped into an 8 × 8 fragA matrix delimited by dashed lines in part (b). Since the m8n8k128 TC instruction expects fragA to have shape 8 × 4, two slices are packed into each fragA entry: for every row of the 8 × 8 fragA matrix, the two σ -bit

connectivity patterns belonging to the two selected columns are concatenated into the first 2σ bits of a 32-bit element, with the remaining 16 bits set to zero, reducing the 8 × 8 fragA matrix to the 8 × 4 fragA matrix. This packing is optimal: the 8 × 4 fragA matrix cumulatively holds data for exactly 64 slices, so multiplying against a single σ -bit frontier word fills all 64 entries of the 8 × 8 fragC matrix with useful popcounts, wasting no output. B LEST uses a single 4 × 8 fragB shared across both rounds, as shown in part (b). Let α denote the σ -bit frontier word and ϑ the all-zero σ -bit word. We define three 16-bit building blocks: X = αϑ, Y = ϑα, and Z = ϑϑ. Each row i ∈ {0, 1, 2, 3} of fragB contains exactly two nonzero 32-bit words, located at columns j = 2i and j = 2i + 1; all remaining entries are zero. For the eight nonzero entries, if j is even, the corresponding 32-bit word is XZ , and if j is odd, it is Y Z , as visualized in the fragB panel of part (b). Consequently, for even j , fragB selects the first packed slice via X , and for odd j it selects the second via Y , while the Z components mask out the unused halves. To avoid intra-warp communication for accumulation, we exploit the distribution of fragB words to warp threads, as in part (a) of Fig. 3. The 32-bit words of fragB are all zeros except those held by threads satisfying t mod 9 = 0, which hold XZ , and threads satisfying t mod 9 = 4, which hold Y Z . This distribution is unique: thread t owns the fragC entries at (i, j) and (i, j + 1), where i = ⌊t/4⌋ and j = 2 × (t mod 4), and t’s packed slices reside in fragA[i][j/2]. That fragA entry is multiplied by fragB[j/2][j] = XZ and fragB[j/2][j + 1] = Y Z , which select the first and second slice masks of t respectively, writing their outputs directly to fragC[i][j] and fragC[i][j + 1]—the very entries owned by t. Consequently, every thread retrieves its own slice–frontier dot-product results from its own fragC registers with no intra-warp communication while reducing the number of required MMA calls 8× compared to [15]. The basic B LEST algorithm based on our novel layout is summarized in Alg. 2.

Preprint submitted to IEEE Transactions on Parallel and Distributed Systems

Algorithm 2 B LEST Require: (1) BVSS data structure, (2) src: source vertex Ensure: (1) L: level array 1: L[v] ← ∞, ∀v ∈ V/{src}; L[src] ← 0; ▷ init 2: F1curr [v] ← 0, ∀v ∈ V/{src}; F1curr [src] ← 1; ▷ init 3: Fn ▷ init next [v] ← 0; 4: Qcurr ← [realPtrs[⌊src/σ⌋], realPtrs[⌊src/σ⌋ + 1]); ▷ init 5: Qnext ← ∅ 6: ℓ ← 0; ▷ current BFS level 7: cont ← true; ▷ non-empty state of the frontier 8: while cont do 9: ℓ ← ℓ + 1; 10: for w = warpID; w < |Qcurr |; w ← w + #warps do 11: vssin ← Qcurr [w]; ▷ virtual slice set ID 12: ssin ← virtualToReal[vssin ]; ▷ slice set ID 13: tile ← (vssin ≪ 5) + laneID; 14: (u0 , u1 , u2 , u3 ) ← rowIds[tile]; ▷ vectorized 128-bit read 15: mask ← masks[tile]; ▷ 4 × σ = 32 bit read 16: α ← Fσcurr [ssin ]; ▷ σ -bit frontier word 17: fragB ← 0; 18: if laneID mod 9 = 0 then fragB ← α; 19: if laneID mod 9 = 4 then fragB ← α; fragB ← fragB ≪ 8; 20: fragA ← low16(mask); ▷ the first two masks 21: fragC[0, 1] ← TC(fragA, fragB); ▷ 1st m8n8k128 22: fragA ← high16(mask); ▷ the last two masks 23: fragC[2, 3] ← TC(fragA, fragB); ▷ 2nd m8n8k128 24: for c ∈ {0, 1, 2, 3} do ▷ column in VSS matrix 25: if fragC[c] ̸= 0 then ▷ if dot-product is nonzero 26: u ← (u0 , u1 , u2 , u3 )[c]; ▷ row ID to update 27: ℓprev ← L[u]; 28: if ℓ < ℓprev then ▷ if u is not visited 29: L[u] ← ℓ; ▷ set u’s level 30: ssout ← ⌊u/σ⌋; ▷ u’s slice set index atomic

35: 36: 37: 38: 39: 40:

old ← Fσnext [ssout ] ∨ (1 ≪ (u mod σ)); if old = 0 then ▷ u’s set is seen for the first time [s, e) ← [realPtrs[ssout ], realPtrs[ssout + 1]); atomic Qnext ← Qnext ∪ [s, e); G RID S YNC(); ▷ level synchronization cont ← (|Qnext | > 0); ▷ check the frontier state swap(Fcurr , Fnext ); swap(Qcurr , Qnext ); ▷ swap arrays G RID S YNC(); |Qnext | ← 0; Fn ▷ clear next frontier data next ← 0; G RID S YNC();

5.2

Lazy Vertex Updates

31: 32: 33: 34:

When the graph bandwidth and update divergence are high and cannot be reduced by reordering, lines 27, 29, and 31 of Alg. 2 suffer from poor L1/L2 cache hit rates, respectively, as the row IDs updated per warp are highly scattered. Beyond locality, two additional factors degrade performance. First, lines 31 and 34 each issue an atomic operation in the hot path, stalling the warp until the operation completes in L2 cache and its return value is received. Second, lines 27/28 introduce a race condition: multiple threads may concurrently read the same u, pass the unvisited check, and proceed to line 31, inflating the number of atomic operations. The atomic at line 31 correctly prevents duplicate queue insertions, but at the cost of amplifying atomic contention. RCM partially mitigates the race at lines 27/28 for networks that are not scale-free-like by co-scheduling warps on the same SM over clustered, largely overlapping row ID ranges, increasing the likelihood that a concurrent level update is already visible in the shared L1 cache before a second thread checks. Furthermore, on the same networks, the per-level frontier is

small enough that the absolute number of atomics remains manageable. On scale-free-like networks, however, their cumulative effect is detrimental. To address these bottlenecks simultaneously, we introduce the lazy vertex update scheme. The GPU instruction set architecture (SASS) provides two instructions that ensure atomicity: ATOMG and REDG. Although both guarantee full atomicity, REDG is the asynchronous counterpart of ATOMG: the compiler emits REDG whenever it detects that the return value of an atomic operation is not used by the calling warp, allowing the warp to continue execution without waiting for the previous memory value to be returned—a capability that ATOMG does not offer. Motivated by this, the lazy vertex update scheme defers all necessary updates to the end of the current BFS level by replacing the synchronous atomic operations in Algorithm 2 with their asynchronous REDG counterparts, and resolves them fully only once the entire BFS level completes. The modified algorithm is provided in Algorithm 3. Alg. 3 introduces two additional arrays, Vcurr and Vnext , holding the cumulative visit status of all vertices up to and including the current and next BFS levels, respectively. The set difference Vnext \ Vcurr therefore identifies exactly the vertices joining the frontier at the next level. The per-level computation is split into two stages. Stage 1 (Lazy marking). The TC multiplication proceeds identically to Alg. 2, but for every pulling vertex u that 1 produces a nonzero dot-product, Vnext [u] is marked lazily via an asynchronous REDG atomic or, with no return value consumed. Crucially, not every such vertex genuinely requires an update: vertices already visited at prior levels are filtered out in Stage 2. Deferring this decision eliminates all synchronous ATOMG instructions from the hot path. Stage 2 (Frontier finalization). Once all warps complete Stage 1, a G RID S YNC separates the two stages, after which all threads sweep Vnext collectively in a fully coalesced manner. Each thread processes one 32-bit word per iteration, 32 32 computing diff = Vcurr [t] ⊕ Vnext [t] to isolate the bits corresponding to vertices that are genuinely new to the frontier. The σ -bit partition ssmask of diff for each of the 4 slice sets in the word is extracted and, if nonzero, iterated via the ffs instruction—which locates the lowest set bit in a single cycle—to identify each new frontier vertex u, assign L[u] ← ℓ, and enqueue the corresponding VSS range [realPtrs[ssout ], realPtrs[ssout + 1]). Since threads are 32 assigned to disjoint words, the F32 curr [t] and Vcurr [t] updates require no atomics, contrary to Alg. 2, where Fcurr updates require atomics with expensive cache costs. The level updates L[u] ← ℓ, on the other hand, exploit high spatial locality since u values within a word are consecutive. Furthermore, the VSS enqueuing is performed warp-atomically, making the only full atomic part of the algorithm immensely cheap, reducing the required number of atomic operations by a factor of 32× relative to only this part of Alg. 2. Although Stage 2 introduces a Θ(n) sweep per BFS level, we identify that on scale-free-like networks with only tens of BFS levels, it removes almost the entire atomic bottleneck and resolves cache problems, making the trade-off highly beneficial. 5.3

Switching Mechanism of B LEST

During a BFS, a frontier crowded with many VSSs may become a bottleneck even with the optimizations described

Preprint submitted to IEEE Transactions on Parallel and Distributed Systems

Algorithm 3 B LEST with Lazy Vertex Updates

Algorithm 4 B LEST Bottom-Up Kernel with CUDA Cores

Require: (1) BVSS data structure, (2) src: source vertex Ensure: (1) L: level array 1: L[v] ← ∞, ∀v ∈ V/{src}; L[src] ← 0; ▷ init 2: F1curr [v] ← 0, ∀v ∈ V/{src}; F1curr [src] ← 1; ▷ init 3: V1curr [v] ← 0, ∀v ∈ V/{src}; V1curr [src] ← 1; ▷ init 4: V1next [v] ← 0, ∀v ∈ V/{src}; V1next [src] ← 1; ▷ init 5: Qcurr ← [realPtrs[⌊src/σ⌋], realPtrs[⌊src/σ⌋ + 1]); ▷ init 6: ℓ ← 0; ▷ current BFS level 7: cont ← true; ▷ non-empty state of the frontier 8: while cont do 9: ℓ ← ℓ + 1; 10: Stage 1: Lazy marking 11: for w = warpID; w < |Qcurr |; w ← w + #warps do 12: ... 13: Lines 11–23 of Algorithm 2 14: ... 26: for c ∈ {0, 1, 2, 3} do ▷ column in VSS matrix 27: if fragC[c] ̸= 0 then ▷ if dot-product is nonzero 28: u ← (u0 , u1 , u2 , u3 )[c]; ▷ row ID to update

Require: CSR (Compressed Sparse Rows) of G⊤ 1: Stage 1: Lazy marking 2: for w = warpID; w < ⌈n/32⌉; w ← w + #warps do 32 3: unvisitedMask ← ¬ Vcurr [w] 4: if (unvisitedMask ≫ laneID) & 1 then 5: u ← 32w + laneID 6: if u < n then 7: for ptr = rowPtrs[u] to rowPtrs[u + 1] − 1 do 8: v ← colIds[ptr] 9: if F1curr [v] = 1 then

29: 30: 31: 32: 33: 34: 35: 36: 37: 38: 39: 40: 41: 42: 43: 44: 45: 46: 47: 48: 49: 50: 51: 52: 53: 54:

V1next [u] |Qnext | ← 0; G RID S YNC();

atomic

∨

▷ lazy mark ▷ prepare for Stage 2

Stage 2: Frontier finalization for t = threadID; t < ⌈n/32⌉; t ← t + #threads do next ← V32 ▷ vertices visited until here next [t]; diff ← V32 ▷ new frontiers in Stage 1 curr [t] ⊕ next; rssOffset ← 4t; ▷ 4 slice sets per word if diff ̸= 0 then ▷ next frontier is non-empty V32 curr [t] ← next; F32 ▷ set bits for frontier vertices curr [t] ← diff; for set ∈ {0, 1, 2, 3} do ▷ each 8-bit word in dif f ssmask ← (diff ≫ (8 · set)) & 0xFF; ▷ select 8-bits if ssmask ̸= 0 then ▷ change is in this set ssout ← rssOffset + set; while ssmask ̸= 0 do ▷ go over all bits b ← ffs(ssmask ) − 1; ▷ in {0, . . . , 7} u ← ssout · σ + b; ▷ new frontier vertex L[u] ← ℓ; ssmask ← ssmask & (ssmask − 1); ▷ unset u [s, e) ← [realPtrs[ssout ], realPtrs[ssout + 1]); warp-atomic

Qnext ← G RID S YNC(); cont ← (|Qnext | > 0); swap(Qcurr , Qnext ); G RID S YNC();

Qnext ∪ [s, e);

above. In such cases, bottom-up exploration of the next frontier with pull-based mechanics on CUDA cores may be preferable. Indeed, the graphs that may benefit from this are those for which lazy updates are useful. Conveniently, the 2nd stage, frontier finalization described above, traverses the entire Vcurr , and the number of unvisited vertices can be computed as the popcount of its bitwise complement. This value, #unvisited, is used by B LEST to decide whether to switch at the start of each level based on the condition

#unvisited < η · |Qcurr |,

(6)

where |Qcurr | is the current frontier queue size, and η is the switching constant, set to 10 based on our preliminary results. When the inequality (6) is met for a level, B LEST invokes Alg. 4 instead of the standard TC-based kernel. Finally, B LEST relies on G RID S YNC in Algs. 2, 3, and 4,

atomic

32 10: Vnext [w] ∨ (1 ≪ laneID) 11: break 12: |Qnext | ← 0; ▷ prepare for Stage 2 13: G RID S YNC(); 14: Stage 2: Frontier finalization (same as Algorithm 3)

which has been available via Cooperative Groups [54] and synchronizes the grid without returning control back to the host. This eliminates per-level kernel-launch overhead, whose performance impact is especially significant on graphs with thousands of BFS levels, such as road networks.

B LEST FOR M ULTI -S OURCE BFS W ORKLOADS

Although B LEST can trivially perform concurrent BFSs with distinct sources, a careful profiling [55] of Alg. 3 reveals that 39% of the total execution time is spent waiting on accesses to BVSS. With an efficient multi-source execution, it is possible to overlap these memory accesses across concurrent BFSs. 6.1

Multi-Source BFS

Let κ be the number of BFSs processed concurrently in a single kernel, and let ρ be ⌈n/σ⌉. In Alg. 3, the frontier and visited arrays encode a single BFS: bit u is set if vertex u is in the frontier or has been visited, respectively. The natural extension to κ concurrent BFS instances stacks these arrays in structure-of-arrays (SoA) format, so that each u is associated with a κ-bit element whose c-th bit encodes the frontier or visited status of u in the c-th BFS. Bitwise operations such κ κ as diff = Vcurr [u] ⊕ Vnext [u] then process all κ BFS instances simultaneously in a single instruction, as in [56]. This naive adoption, however, introduces a memory coalescing problem. In Alg. 3, each thread processes 4σ consecutive vertices (4 slice sets) and updates the corresponding frontier entries, with consecutive threads touching consecutive 4-byte words — an access pattern that is fully coalesced. With κ = 256, each κ-bit element occupies 32 bytes, so the stride between the elements touched by consecutive threads becomes 4σ κ8 bytes (1KB for σ = 8). We resolve this by (1) scheduling consecutive threads to consecutive slice sets (instead of words) and (2) introducing a bijective index remapping getVI(u, ρ) = (u mod σ) · ρ + ⌊u/σ⌋, which contiguously stores the vertices with the same u mod σ . For example, with n = 24 and ρ = 3, the layout becomes 0, 8, 16, 1, 9, 17, . . . , 7, 15, 23. Hence, the index u · κ + c is replaced by getVI(u, ρ) · κ + c. This allows consecutive threads to process consecutive slice sets, whose κ-bit elements are contiguous in memory, reducing the inter-thread stride to κ/8 bytes, i.e., 32 bytes for κ = 256, and restoring coalesced accesses with only reindexing overhead.

Preprint submitted to IEEE Transactions on Parallel and Distributed Systems

When a slice set is active in at least one BFS, all of its VSSs are enqueued into the frontier. However, being active only for a single BFS should not yield unnecessary lazy markings to be conducted for every other κ − 1 BFS instances. To address this, we introduce an array activeSets, where activeSets[s] is a κ-bit element whose c-th bit is set iff slice set s is active for the c-th BFS. During lazy marking, the warp iterates over the set bits of activeSets[s], and processes only the active BFSs. For efficiency, values from different BFSs are first accumulated into a thread-private buffer, and the resulting lazy atomic updates to Vnext are issued to L2 cache only after all active BFSs have been processed. Hence, neither the lazy marking nor the frontier finalization introduces additional overhead. The algorithm is given in Alg. 5. Without reindexing, Fσcurr [ssin ] is a contiguous σ -bit word storing the active/frontier status of the columns of slice set ssin . However, with it, these σ bits are no longer contiguously stored. To reconstruct them efficiently for the TC multiplication, the first σ lanes of the warp gather Fκcurr [laneID ·ρ+ssin ] (line 12 of Alg. 5). The remaining 32−σ lanes do not issue any memory access and initialize their buffer to zero (line 13). Each value fetched by the first σ lanes is a κ-bit word encoding the full multi-source BFS frontier status, with bit c of lane u’s word set if vertex ssin · σ + u is in the frontier for BFS c. In each iteration over an active BFS c (line 16), the warp performs a partial warp-level transpose via warp.ballot: each of the σ lanes extracts bit c from its local κ-bit buffer (line 19), contributing a single frontier bit to the collective ballot. Since lane u ∈ {0, . . . , σ − 1} holds the frontier status of vertex ssin · σ + u for all κ BFSs, its vote corresponds to whether vertex ssin · σ + u is in the frontier for BFS c. The result of warp.ballot, the σ -bit word α, assembling these σ votes, is broadcast to all lanes (line 20). Then, the TC multiplication runs as in the single-source BFS. The lazy vertex update scheme is a bottleneck for networks that are not scale-free-like; the per-level frontier is small and full swap in Stage 2 is expensive. To mitigate, we introduce dirty slice sets; a slice set is marked dirty during Stage 1 whenever any of its vertices is lazily marked in Vnext , signaling that it requires processing in Stage 2. Dirty slice sets are built during lazy marking and consumed during frontier finalization; if a slice set is not dirty, it is skipped entirely in Stage 2, since none of its vertices were marked and no new frontier entries can reside within it. The dirty and active slice sets are duals: the former is produced in Stage 1 to accelerate Stage 2, while the latter is vice versa. Alg. 5 is kept generic to use for kernels such as closeness centrality. Hence, the level updates which are not necessary in those are not included. If they were, the statement level[κ × b + idx] = ℓ, which sets u’s level for the bth BFS–for every set bit of dif f at index b–would be within the if at line 49.

Algorithm 5 B LEST (Generic) Multi-Source BFS Kernel

6.2

After accumulating the far values across all BFSs, the centrality can be computed via:

Closeness Centrality

To compute closeness centrality of all vertices, one can run an all-pairs shortest-path computation on G. When G is connected, undirected, and unweighted, this is equivalent to performing BFSs from all source vertices, which requires an efficient multi-source BFS kernel. Within the if statement at line 49, the required values can be captured as

far[idx] ← far[idx] + ℓ · popc(diff).

(7)

Require: (1) BVSS data structure, (2) srcs: κ source vertices {src0 , . . . , srcκ−1 } 1: init(srcs, Fcurr , Vcurr , Vnext , Qcurr , activeSets, dirtySets); 2: n̂ ← ⌈n/σ⌉ · σ ; ▷ σ -padded vertex count 3: while cont do 4: ℓ ← ℓ + 1; 5: Stage 1: Lazy marking 6: for w = warpID; w < |Qcurr |; w ← w + #warps do 7: ... 8: Lines 11–15 of Algorithm 2 9: ... 12: if laneID < σ then fv ← Fκcurr [laneID · ρ + ssin ]; 13: if laneID ≥ σ then fv ← 0; 14: B ← activeSets[ssin ]; ▷ κ-bit active BFS mask 15: M[0..3] ← 0; ▷ per-row register accumulation buffers 16: while B ̸= 0 do ▷ two TC calls per active BFS c 17: c ← ffs(B) − 1; ▷ in {0, . . . , κ − 1} 18: B ← B & (B − 1); ▷ unset c 19: inc ← (fv ≫ c) & 1; 20: α ← warp.ballot(inc); ▷ reconstruct frontier 21: ... 22: Lines 17–23 of Algorithm 2 23: ... 28: for q ∈ {0, 1, 2, 3} do ▷ column in VSS matrix 29: if fragC[q] ̸= 0 then 30: M[q] ← M[q] ∨ (1 ≪ c); ▷ accumulate 31: for q ∈ {0, 1, 2, 3} do ▷ flush register buffers to Vnext 32: if M[q] ̸= 0 then 33: u ← (u0 , u1 , u2 , u3 )[q]; 34: 35: 36: 37: 38: 39: 40: 41: 42: 43: 44: 45: 46: 47: 48: 49: 50: 51: 52: 53: 54: 55: 56: 57:

atomic

κ [getVI(u, ρ)] ∨ M[q]; Vnext dirtySets[u/σ] ← 1; |Qnext | ← 0; ▷ prepare for Stage 2 G RID S YNC(); Stage 2: Frontier Finalization for ub = threadID · σ; ub < n̂; ub ← ub + #threads · σ do ssout ← ub /σ if dirtySets[ssout ] = 0 then continue dirtySets[ssout ] ← 0; A ← 0; for u = ub , . . . , ub + σ − 1 do idx ← getVI(u, ρ); next ← Vκnext [idx]; diff ← Vκcurr [idx] ⊕ next; Fκcurr [idx] ← diff; if diff ̸= 0 then κ Vcurr [idx] ← next; A ← A ∨ diff; activeSets[ssout ] ← A; if A ̸= 0 then [s, e) ← [realPtrs[ssout ], realPtrs[ssout + 1]);

warp-atomic

← Qnext ∪ [s, e); ... Lines 51–54 of Algorithm 3 Qnext

cc[u] = (n − 1)/far[getVI(u, ρ)].

(8)

Since Alg. 5 processes κ concurrent BFS instances per launch, it is invoked ⌈n/κ⌉ times, each time from a distinct batch of κ sources, to perform all necessary distance computations. For a disconnected graph, a different normalization using component sizes can be applied to each centrality score.

Preprint submitted to IEEE Transactions on Parallel and Distributed Systems

E XPERIMENTS

B LEST3 is implemented in C++/CUDA and compiled with gcc 12.3.0 and CUDA 13.0. Experiments are conducted on two servers, Arch-1 and Arch-2. Arch-1 is equipped with a 64-core Intel Xeon Gold 6548Y+ (2.5 GHz) and 1 TB of host RAM, paired with an NVIDIA H200 GPU (141 GB HBM3e) on the device side. Arch-2 is equipped with two Intel Xeon Platinum 8460Y+ processors (40 cores each, 2.0 GHz) and 512 GB of host RAM, paired with NVIDIA H100 GPUs (64 GB HBM2e) on the device side. All experiments are conducted on Arch-1, with the exception of the Closeness Centrality experiment (Table 3), which is conducted on Arch-2. We use two benchmark suites: (1) the GAP Benchmark Suite [17], and (2) a custom benchmark suite containing large graphs from SuiteSparse [58], comprising all graphs 3. https://github.com/delbek/blest

TABLE 2: Performance (in ms.) of B LEST against SotA SS-BFS implementations. All runtimes are averaged over 64 random sources. All the speedups are normalized to [15]. Some baselines failed on large graphs due to unsupported graph size or out-ofmemory issues. These cases are marked with (err.) in the table.

57M 1.4B 1.9B 4.2B 4.2B

739.6 317.5 130.2 161.9 201.5 (err.) (err.) 63.7 327.1 (err.) 277.8 74.3 290.2 (err.) (err.) 164.9 468.8 (err.) (err.) 192.9

nlpkkt240 ✓ 27M 760M 234.3 25.1 66.4 33.3 8.9 3.7× uk-2005 ✓ 39M 936M 177.2 331.8 36.3 22.8 8.8 2.6× it-2004 ✓ 41M 1.1B 216.2 (err.) 97.1 21.4 10.1 2.1× europe_osm X 50M 108M 1573.6 942.1 417.2 481.9 137.5 3.5× com-Friends. ✓ 65M 3.6B 347.0 (err.) (err.) 135.6 11.5 11.8× Spiel._k600 X 72M 216M 520.9 40.5 231.2 181.3 19.4 9.3× webbase-2001 ✓ 118M 1.0B 292.0 72.8 74.7 29.2 12.1 2.4× kmer_V1r X 214M 465M 1493.7 60.0 244.6 317.9 45.3 7.0× mawi ✓ 226M 480M 272.9 882.9 107.1 457.3 55.5 8.3×

Speedups are given over the graphs the kernels could process.

Min Max Geomean

0.1× 0.1× 0.2× 1.7× 5.3× 4.3× 0.3× 0.8× 0.7×

1× 1× 1×

B LEST vs. [15]

|E|

X 23M ✓ 61M ✓ 50M ✓ 134M ✓ 134M

B LEST (this work)

|V|

GAP-road GAP-twitter GAP-web GAP-kron GAP-urand

GAP [17]

Graph

Gunrock [10] GSWITCH [11] BerryBees [15]

Times (in milliseconds) Scale-Free

A binary TC multiplication provides a single bit of information for two (binary) vectors: whether or not there exists an index at which both vectors simultaneously have a set bit. For graph algorithms in which information transmission can be conducted entirely on this single bit, TCs are well-suited and promise performance improvements. BFS, Multi-Source BFS, and Closeness Centrality—as well as the ones not covered in this work, such as Triangle Counting—satisfy this property. BFS focuses on the existence of two set bits at the same position of the frontier vector and the connectivity mask, i.e., whether the pulling vertex has an active frontier neighbour. Since this transmission involves only one bit, BFS maps onto TCs. Not every graph algorithm, however, depends on information that can be encoded in a single bit. Consider a hypothetical adaptation of Alg. 3 to Connected Components. Instead of L, one can maintain a comp array, where comp[u] = comp[v] if vertices u and v are connected. A trivial, synchronized implementation initially enqueues all VSSs in the frontier and sets comp[u] = u, u ∈ V . Whenever a TC produces a nonzero result at line 27 of Alg. 3, it signals a state change for u; comp[u] ← comp[v] for some column/vertex v of the VSS. However, although a state change signal is a one-bit value, one must also identify v to retrieve the actual information comp[v] to be transmitted into comp[u]. This necessitates an explicit bitwise traversal of the σ -bit slice to locate the column v , which makes the kernel inefficient. On the contrary, for BFS and closeness centrality, the information to be transmitted, L[v], i.e., the distance to the BFS source, is the same regardless of the VSS column v . A similar argument can also be made for kernels whose mechanics are more similar to BFS, such as Betweenness Centrality. In a Brandes-like implementation [57], the state change during the first stage implies npaths[u] += npaths[v] where npaths[u] is the number of shortest paths from the BFS source to u for all u ∈ V . To correctly perform this update, one needs to locate v . This is also true for BFS implementations that lack level synchronization, for which the active frontier queue can accommodate vertices belonging to different BFS levels. For all these graph kernels, the techniques used by B LEST need to be extended as they cannot be applied as efficiently as they currently are.

satisfying |V| ≥ 23M and |E| ≤ 232 − 1. When multiple graphs from the same matrix group meet these criteria, only the largest is retained. Graphs already included in the GAP suite are deduplicated, yielding 14 graphs in total. Table 2 presents the single-source BFS (SS-BFS) performance of B LEST against SotA implementations. B LEST decisively outperforms all baselines across all graphs, with average speedups of 22.0×, 7.7×, 8.1×, and 5.9× over GAP [17], Gunrock [10], GSWITCH [11], and BerryBees [15], respectively. Even [15]—the only other TC-based SS-BFS kernel and the fastest prior implementation in the literature— is on average 5.9× slower than B LEST, establishing B LEST as the new state-of-the-art in the SS-BFS literature.

GAP

Using Tensor Cores for Other Graph Algorithms

|V| ≥ 23M, |E| ≤ 232 − 1

6.3

41.3 3.9× 6.6 9.7× 15.7 4.7× 7.2 23.1× 13.4 14.4×

2.1× 23.1× 5.9×

For the Multi-Source BFS (MS-BFS), the most recent GPU-based implementation is not compatible with modern GPUs [59]. Although CPU-based, another well-known implementation [60] delivers good performance on 64 cores of Arch-1, yet it is 25.9× slower on average than B LEST for the graphs in Table 2. Given the current state of the MS-BFS literature, B LEST is all the more valuable: it is engineered for modern GPUs and delivers performance that is 2.7× faster on average than its own single-source variant. Since computing closeness centrality for all the vertices has O(nm) complexity, approximation algorithms [61], [62] that apply sampling, pivoting, pruning, or top-k approaches are used at the cost of accuracy. With B LEST, it is now possible to compute exact values for all vertices in real-world graphs. To demonstrate B LEST’s MS-BFS performance, we compute exact scores for com-Friendster from Table 2, which is, to the best of our knowledge, the largest graph for which exact centrality values have ever been computed. With 100 H100s, computing the closeness centrality of com-Friendster took 3,665 seconds, i.e., roughly 100 GPU-hours. For this experiment, sources are partitioned via MPI into disjoint batches across GPUs and the partial sums are reduced in a

Preprint submitted to IEEE Transactions on Parallel and Distributed Systems

final stage. We additionally compute the closeness centrality values for five other graphs—two social, one geometric, one random, and one road network—given in Table 3. TABLE 3: Closeness centrality performance of B LEST across six graphs. Execution times are reported in seconds. |V|

GPUs Graph

|E| Time Graph

|V|

|E| Time

5 com-LiveJ. 4.0M 69.4M 118 com-Orkut 3.1M 234.4M 131 10 delaun._n24 16.8M 100.7M 4,836 rgg_n_2_24 16.8M 265.1M 4,711 100 com-Friend. 65.6M 3.6B 3,665 GAP-road 23.9M 57.7M 1,701

7.1

Ablation Studies

We evaluate the benefits of the proposed optimizations with an ablation study for both SS-BFS and MS-BFS variants. For SS-BFS, we construct five variants of B LEST, starting from the base B LEST (A) and incrementally adding one optimization at a time until the final version is reached. This ablation study is summarized in Table 4. B LEST (A) is the initial variant, incorporating only the BVSS data structure and kernel fusion, which achieves a 1.6× average speedup over [15] (col. 3). As optimizations are added one by one, the speedup increases, with the final version of B LEST achieving an average 5.9× speedup over [15] across all the graphs. The ablation study clearly demonstrates that each optimization may contribute, although the impact varies across graphs. For instance, reordering is effective on 8/14 graphs, neutral on 2/14, and degrades performance on 4/14. TABLE 4: SS-BFS ablation study of B LEST (in ms., averaged

Graph

Times (in milliseconds) Pseu. Berry B LEST Diam. Bees [15] (A) (AB) (ABC) (ABCD) (Full)

GAP

GAP-road GAP-twitter GAP-web GAP-kron GAP-urand

8.4K 16 64 8 7

161.9 62.8 49.3 41.3 63.7 66.2 54.9 41.3 74.3 19.1 17.0 16.7 164.9 177.7 154.4 74.3 192.9 236.0 202.6 201.3

N/A N/A 18.8 6.6 15.7 N/A 49.6 7.2 52.3 13.4

|V| ≥ 23M , |E| ≤ 232 − 1

over 64 random sources). (A): BVSS + kernel fusion. (AB): (A) + optimal TC layout. (ABC): (AB) + reordering. (ABCD): (ABC) + lazy vertex updates. (Full): (ABCD) + switching. N/A indicates the optimization is inapplicable. For each graph/column, the speedup w.r.t. [15] is computed using the latest non-N/A runtime among all variants up to and including that column.

nlpkkt240 242 uk-2005 210 it-2004 65 europe_osm 30.1K com-Friendster 26 Spielman_k600 601 webbase-2001 632 kmer_V1r 589 mawi 7

33.3 14.9 11.2 9.0 22.8 13.0 11.7 14.2 21.4 14.8 13.3 14.4 481.9 198.8 160.4 137.5 135.6 173.3 146.2 107.5 181.3 138.0 130.7 14.5 29.2 27.2 25.3 30.6 317.9 370.9 366.3 168.4 457.3 56.7 42.6 57.2

8.9 N/A 8.8 N/A 10.1 N/A N/A N/A 53.5 11.5 19.4 N/A 12.1 N/A 45.3 N/A 55.5 N/A

Min Speedup Max Speedup Geomean

1.0× 0.8× 0.9× 1.0× 1.0× 8.1× 10.7× 12.5× 1.0× 1.6× 1.9× 2.5×

2.1× 2.1× 9.3× 22.9× 3.9× 5.9×

Some optimizations are marked N/A for certain graphs, as they are designed for graphs meeting specific conditions. Direction switching is inapplicable to graphs whose per-level frontier never satisfies Eq. (6), and the corresponding (Full) column entry is marked N/A. Similarly, the lazy vertex update scheme (Sec. 5.2) is developed specifically for scalefree-like networks whose update divergence is high and cannot be reduced via reordering algorithms such as RCM

due to their inherently high bandwidth. For graphs that do not require the lazy scheme, Alg. 2 is used exclusively, and the corresponding ABCD entry in Tab. 4 is marked N/A. The dispatch decision is made based on Udiv : graphs with Udiv > 25,000 are dispatched to the lazy kernel (Alg. 3), while those with Udiv ≤ 25,000 are dispatched to Alg. 2. The difference between (A) and (AB) lies in the layout; in (A), the one in [15] is used, requiring 16 TC multiplications per 128 slices. (AB) replaces this by the layout proposed in Sec. 5.1, requiring only 2 TC multiplications. Both variants use TCs as the compute unit, though the latter utilizes them 8× more efficiently. Although the average speedup is only 1.2×, it is consistent over all the graphs. Although Tab. 4 shows that reordering degrades performance on 4/14 graphs, this does not imply reordering being harmful, but rather that the graphs favor natural orderings that reordering cannot improve upon. To validate this, we randomly reorder 6 graphs—including 4/14 where reordering is harmful and the 2/14 where it is neutral—and measure the performance of (ABC) against the randomly reordered variant (AB) in Tab. 5. Except for mawi and GAP-urand, reordering is effective, although its real benefit is limited due to good natural vertex orderings. Since the reordering requires a hyperparameter W to be set, we evaluate its sensitivity on GAP-web by varying W ∈ [23 , 218 ] and analyzing the compression ratio and BFS execution time in Figure 4. As Fig. 4 shows, increasing W has a positive impact on both the compression ratio and BFS performance, as expected. However, the relationship is clearly concavedown, signaling diminishing returns as W grows. Since larger W also increases the reordering cost, the choice of W requires balancing compression quality against preprocessing overhead. In all our experiments, we set W = 216 . TABLE 5: Effect of replacing the natural ordering by a random reordering on graphs where (ABC) degrades the performance over (AB). Runtimes are in milliseconds. Graph GAP-web webbase-2001 mawi

Rnd. Order B LEST B LEST (AB) (ABC) Graph 108.77 107.21 55.47

16.66 it-2004 30.63 uk-2005 57.24 GAP-urand

Rnd. Order B LEST B LEST (AB) (ABC) 75.82 14.42 64.07 14.16 201.83 201.30

Fig. 4: Effect of window size W on compression ratio and BFS runtime on GAP-web.

The impact of switching depends heavily on the perlevel frontier size, as expressed in (6). Consequently, graphs with a small diameter are the most promising candidates for switching to the bottom-up CUDA Core implementation described in Section 5.3. Figure 5 analyzes the switching

Preprint submitted to IEEE Transactions on Parallel and Distributed Systems

(a) com-Friendster

(b) mawi

(d) GAP-web

(e) GAP-kron

(f) GAP-urand

Fig. 5: Switching analysis across 6 graphs. Each subplot shows the execution time (ms) of Top-Down, Bottom-Up, B LEST, and Optimal for each BFS level. The Optimal curve is simply min(Top-Down, Bottom-Up) for each level, representing ideal switching. The misclassification rates and speedups indicate the fraction of levels at which B LEST is wrong and the speedup of Optimal over B LEST.

Finally, we evaluate the optimizations of Section 6 for the multi-source BFS extension of B LEST via the ablation study in Table 6. (Naive) serves as the baseline, executing the SS-BFS kernel independently over multiple sources. Three further

variants follow: (A) is Alg. 5 without reindexing or dynamic slice set tracking (activeSets/dirtySets); (AB) drops only the dynamic tracking; and (Full) is Alg. 5 in its entirety. TABLE 6: MS-BFS ablation study of B LEST (in sec., total of 256 BFSs over random sources). (Naive): multiple SS-BFSs. (A): Alg. 5 without reindexing or dynamic slice set tracking. (AB): Alg. 5 without dynamic slice set tracking. (Full): Alg. 5.

GAP

Graph

|V| ≥ 23M |E| ≤ 232 − 1

behavior of the 6 graphs from Tab. 4 with the lowest pseudodiameters, computed by finding 64 pseudo-peripheral vertices and executing BFS from each to identify the traversal with the greatest depth. The modes are: Top-Down (TC), Bottom-Up (CUDA Core), B LEST (the policy determined automatically via (6)), and Optimal (a hypothetical oracle selecting the correct mode at each level). For graphs where B LEST determines that switching is unnecessary,—by executing 3 BFS runs from random sources both with and without switching to assess its benefit—the B LEST curve always tracks the Top-Down curve. This is the case for mawi (Fig. 5b) and GAP-web (Fig. 5d), where the Top-Down curve aligns exactly with the Optimal curve, confirming that the preprocessing decision is correct. When switching is enabled, B LEST attempts to follow the Optimal curve by applying Eq. (6) at each level. For most, B LEST’s decisions are nearoptimal, e.g., on GAP-urand (Fig. 5f), the curves are nearly indistinguishable, with the notable exception of GAP-kron (Fig. 5e), where a misclassification causes B LEST to select the wrong mode at 2 out of 8 levels, resulting in a 1.69× higher execution time than Optimal. On GAP-twitter (Fig. 5c), the match is exact, demonstrating that η = 10 is an effective threshold on H100, although a graph-specific η may extract the remaining performance headroom. We note, however, that η = 10 assumes that processing one VSS with TCs costs the same as processing 10 vertices with CUDA Cores, a ratio calibrated specifically for the Hopper GPU. This constitutes a limitation of B LEST: the switching threshold is GPU-specific and may require retuning on other architectures.

Times (in seconds) B LEST B LEST (Naive) (A) (AB) (Full)

GAP-road GAP-twitter GAP-web

10.79 1.73 4.35

41.83 1.47 1.97

3.64 1.11 0.77

nlpkkt240 uk-2005 it-2004 europe_osm com-Friendster Spielman_k600 webbase-2001

2.30 10.61 10.12 2.76 2.37 1.84 2.90 2.10 1.78 33.30 146.14 100.51 3.00 4.87 4.20 4.97 32.64 30.16 5.17 9.30 5.66

1.68 0.70 0.64 9.05 2.75 2.14 1.56

Min Max Geomean

7.2

1.0× 1.0× 1.0×

49.18 1.80 2.39

0.15× 0.16× 1.09× 1.82× 2.21× 5.65× 0.52× 0.65× 2.69×

Overhead of B LEST

The preprocessing overhead of B LEST is threefold: (1) constructing the CSC (Compressed Sparse Columns) of A from the edge list, which is an I/O-bound operation; (2) reordering the CSC matrix via either Alg. 1 or RCM; and (3) constructing the BVSS data structure. The first overhead is shared by all state-of-the-art implementations listed in Table 2, whereas the latter two are specific to B LEST. Among the baselines, [15] also constructs a specialized data structure, BRS, whose construction cost is analogous to that of BVSS, as described

Preprint submitted to IEEE Transactions on Parallel and Distributed Systems

in Sec. 3. Table 7 reports the preprocessing costs for B LEST for all the graphs. Whenever Alg. 1 is selected, we set W = 216 . TABLE 7: Preprocessing overhead of B LEST (in sec.) across all

TABLE 8: Memory footprint of B LEST (in GB) across all experimental graphs for the Multi-BFS workload (κ = 256 sources). The CSR column is only applicable to graphs for which direction switching is enabled; it is marked N/A otherwise. Active CSR BVSS Frontier Queues Levels /Dirty Total

GAP

GAP-road GAP-twitter GAP-web GAP-kron GAP-urand

N/A 1.77 6.12 5.76 N/A 1.79 17.43 18.58 17.72 19.85

2.30 5.91 4.86 12.88 12.88

0.02 24.52 0.10 63.06 0.05 51.85 0.28 137.44 0.34 137.44

0.10 28.71 0.25 81.20 0.21 58.77 0.55 187.17 0.55 188.78

nlpkkt240 N/A 1.28 uk-2005 N/A 0.88 it-2004 N/A 1.06 europe_osm N/A 3.77 com-Friendster 14.71 15.93 Spielman_k600 N/A 5.34 webbase-2001 N/A 1.48 kmer_V1r N/A 15.84 mawi N/A 2.08

2.69 3.79 3.96 4.89 6.30 6.93 11.34 20.54 21.71

0.03 28.67 0.04 40.41 0.04 42.28 0.05 52.13 0.26 67.18 0.07 73.91 0.09 120.98 0.21 219.14 0.24 231.63

0.12 32.77 0.16 45.28 0.17 47.52 0.21 61.05 0.27 104.65 0.30 86.55 0.49 134.38 0.88 259.34 0.93 259.42

CSC Reord. BVSS

GAP-road R 9.3 9.5 5.3 GAP-twitter J 534.6 388.1 72.1 GAP-web J 340.8 871.8 20.5 GAP-kron J 1014.2 1804.4 159.6 GAP-urand J 1294.0 1942.9 181.3 nlpkkt240 J 79.9 28.9 8.4 uk-2005 J 85.8 306.3 10.4 it-2004 J 97.4 199.6 12.9 europe_osm R 11.8 19.0 22.1 com-Friend. J 720.6 750.6 149.2 Spiel._k600 R 42.9 23.5 15.7 webb.-2001 J 100.6 429.1 15.1 kmer_V1r R 70.5 164.6 74.6 mawi J 62.0 1312.7 81.5

We finally evaluate the memory footprint of B LEST in Table 8. As evident from the table, the structures required for BVSS and the computational mechanics of B LEST are remarkably compact, even for the largest real-world graphs. The size of level arrays in MS-BFS, however, grows proportionally with the number of concurrent BFS instances κ and becomes the dominant memory consumer. For 4 graphs, this results in an out-of-memory (OOM) error on the NVIDIA H200, which is why these graphs are excluded from the MSBFS experiments in Table 6. We note that this limitation is inherent to all GPU-based MS-BFS frameworks, and reducing κ straightforwardly can resolve this issue.

Graph

|V| ≥ 23M , |E| ≤ 232 − 1

CSC Reord. BVSS Graph

Ord.

Graph

Ord.

graphs. The columns report: the ordering algorithm used, the cost of CSC construction from the edge list, reordering via J AC CARD W ITH W INDOWS (J) or RCM (R), and BVSS construction.

R ELATED W ORK

The idea of using linear-algebraic primitives in graph algorithms transformed graph analytics [63], [14], [64]. This became prominent with the availability of MMA units in modern accelerators. BerryBees [15] is an early effort to realize BFS using TCs via a bitmap-based data structure. In contrast to B LEST, it uses both SpMV and SpMSpV as a primitive and switches from the latter to the former whenever the frontier is sufficiently dense4 . Both variants of their implementation, however, suffer from major bottlenecks as explained in this work. BerryBees’ SpMV applies a frontier-oblivious slice-set distribution, which in turn severely exacerbates load imbalance across the warps. When the frontier becomes sparse, BerryBees switches to SpMSpV. BerryBees tries to solve the load-balancing problem in the SpMSpV by traversing the frontier and redistributing the slice sets based on the number of slices. This per-level, frontier-dependent redistribution limits the benefits of SpMSpV and highlights the importance of the proposed BVSS data structure, which achieves near-perfect load balance by construction, without requiring the algorithmic layer to be aware of or to react to such dynamic load variations. Gunrock [10] is a GPU library designed to balance performance, scalability, and programmability by expressing the graph algorithms through a data-centric “frontier” abstraction and a bulk-synchronous loop. It applies GPU-specific optimizations such as load-balanced workload mapping, idempotence to reduce atomics under concurrent discoveries, and push vs. pull direction switching. Hence, it specializes 4. This mechanism does not switch between push and pull. Instead, it improves the work efficiency by skipping pulls for vertices whose incoming neighbours contain no frontier vertices.

the execution flow with respect to the graph at hand. The literature, and also this work, follows the message that there is no single “best” GPU BFS frontier expansion policy that does not vary w.r.t. the graph structure. GSWITCH [11] dynamically reconfigures (with low overhead) the execution by selecting among optimization patterns, most critically push/pull switching, plus complementary choices such as load-balancing and active-set representation. It uses runtime/workload features to decide when to switch and which kernels to run, searching a BFS variant space in each iteration.

C ONCLUSION AND F UTURE W ORK

In this paper, we presented B LEST, an ultra-efficient framework for executing BFS, Multi-Source BFS, and Closeness Centrality on modern GPUs at a speed unprecedented in the literature. By systematically engineering every stage of the pipeline and developing novel approaches tailored specifically to the unique characteristics of BFS, we establish B LEST as the fastest BFS implementation in the era of Tensor Cores, achieving a 5.9× geomean speedup over the prior state-of-the-art across 14 diverse real-world graphs. Although B LEST demonstrates strong performance across the vast majority of experimental configurations, the direction-switching mechanism occasionally misclassifies the optimal execution mode on certain graphs, leaving extractable performance unrealized. Developing a more robust, potentially graphadaptive switching criterion that minimizes misclassification without incurring prohibitive preprocessing cost constitutes a natural direction for future work.

ACKNOWLEDGMENTS The numerical calculations reported in this paper were fully performed at 1) TUBITAK ULAKBIM, High Performance and Grid Computing Center (TRUBA resources) and 2) the EuroHPC Joint Undertaking (EuroHPC JU) supercomputer MareNostrum 5, hosted by the Barcelona Supercomputing Center (BSC). Access to MareNostrum 5 was provided through a national access call coordinated by the Scientific and Technological Research Council of Turkey (TÜBİTAK). We are grateful to TÜBİTAK and EuroHPC JU for providing access to these resources and supporting this research.

Preprint submitted to IEEE Transactions on Parallel and Distributed Systems

R EFERENCES [1]

F. W. Takes and W. A. Kosters, “Determining the diameter of small world networks,” in Proc. 20th ACM Int. Conf. on Information and Knowledge Management (CIKM’11). New York, NY, USA: ACM, 2011, pp. 1191–1196. [2] V. K. Akram and O. Dagdeviren, “Breadth-first search-based singlephase algorithms for bridge detection in wireless sensor networks,” Sensors, vol. 13, no. 7, pp. 8786–8813, 2013. [3] V. Pandey and P. Bonde, “BFSMpR: A BFS graph based recommendation system using MapReduce,” Int. J. on Recent and Innovation Trends in Comp. and Communication, vol. 5, no. 6, pp. 445–449, 2017. [Online]. Available: https://ijritcc.org/index.php/ ijritcc/article/view/794 [4] J. Barnat, L. Brim, and J. Chaloupka, “Parallel breadth-first search LTL model checking,” in Proc. 18th IEEE Int. Conf. on Automated Software Eng. (ASE’03). Montreal, Canada: IEEE, 2003, pp. 106–115. [5] C. E. Leiserson and T. B. Schardl, “A work-efficient parallel breadthfirst search algorithm (or how to cope with the nondeterminism of reducers),” in Proc. 22nd ACM Symp. on Parallelism in Algorithms and Archs. (SPAA’10). Santorini, Greece: ACM, 2010, pp. 303–314. [6] S. Beamer, K. Asanović, and D. Patterson, “Direction-optimizing breadth-first search,” in Proc. 2012 ACM/IEEE Int. Conference for High Perf. Comp., Netw., Storage and Analysis (SC’12). Salt Lake City, UT, USA: IEEE Computer Society, 2012, pp. 1–10. [7] L. Dhulipala, G. E. Blelloch, and J. Shun, “Theoretically efficient parallel graph algorithms can be fast and scalable,” ACM Trans. Parallel Comp., vol. 8, no. 1, pp. 4:1–4:70, 2021. [8] J. Shun and G. E. Blelloch, “Ligra: A lightweight graph processing framework for shared memory,” in Proc. 18th ACM SIGPLAN Symp. on Principles and Practice of Parallel Programming (PPoPP’13). Shenzhen, China: ACM, 2013, pp. 135–146. [9] D. Merrill, M. Garland, and A. Grimshaw, “High-performance and scalable GPU graph traversal,” ACM Trans. Parallel Comp., vol. 1, no. 2, pp. 14:1–14:30, 2015. [10] Y. Wang, Y. Pan, A. Davidson, Y. Wu, C. Yang, L. Wang, M. Osama, C. Yuan, W. Liu, A. T. Riffel, and J. D. Owens, “Gunrock: GPU graph analytics,” ACM Trans. Par. Comp., vol. 4, no. 1, pp. 3:1–3:49, 2017. [11] K. Meng, J. Li, G. Tan, and N. Sun, “A pattern based algorithmic autotuner for graph processing on GPUs,” in Proc. 24th ACM SIGPLAN Symp. on Principles and Practice of Parallel Programming (PPoPP’19). Washington, DC, USA: ACM, 2019, pp. 201–213. [12] A. Buluç and K. Madduri, “Parallel breadth-first search on dist. memory systems,” in Proc. 2011 ACM/IEEE Int. Conference for High Perf. Comp., Netw., Storage and Analysis (SC’11). New York, NY, USA: ACM, 2011, pp. 65:1–65:12. [13] F. Checconi and F. Petrini, “Traversing trillions of edges in real time: Graph exploration on large-scale parallel machines,” in Proc. 28th IEEE Int. Parallel and Dist. Processing Symp. (IPDPS’14). Phoenix, AZ, USA: IEEE, 2014, pp. 425–434. [14] J. Kepner, P. Aaltonen, D. Bader, A. Buluç, F. Franchetti, J. Gilbert, D. Hutchison, M. Kumar, A. Lumsdaine, H. Meyerhenke, S. McMillan, C. Yang, J. D. Owens, M. Zalewski, T. Mattson, and J. Moreira, “Mathematical foundations of the GraphBLAS,” in Proc. 2016 IEEE High Perf. Extreme Comp. Conference (HPEC’16). Waltham, MA, USA: IEEE, 2016, pp. 1–9. [15] Y. Niu and M. Casas, “BerryBees: Breadth first search by bit-tensorcores,” in Proc. 30th ACM SIGPLAN Annual Symp. on Principles and Practice of Parallel Programming (PPoPP’25). New York, NY, USA: ACM, 2025, pp. 339–354. [16] W. Liu and A. H. Sherman, “Comparative analysis of the Cuthill–McKee and the reverse Cuthill–McKee ordering algorithms for sparse matrices,” SIAM Journal on Numerical Analysis, vol. 13, no. 2, pp. 198–213, 1976. [Online]. Available: https://doi.org/10.1137/0713020 [17] S. Beamer, K. Asanović, and D. Patterson, “The GAP benchmark suite,” CoRR, vol. abs/1508.03619, 2015. [Online]. Available: https://arxiv.org/abs/1508.03619 [18] Y. Chen and J. X. Yu, “Bitmap-based sparse matrix-vector multiplication with tensor cores,” in Proc. 53rd Int. Conf. on Parallel Processing, ser. ICPP ’24. Sweden: ACM, Aug. 2024, pp. 1135–1144. [19] Y. Lu and W. Liu, “DASP: Specific dense matrix multiplyaccumulate units accelerated general sparse matrix-vector multiplication,” in Proc. Int. Conf. for High Perf. Comp., Netw., Storage and Analysis, ser. SC ’23. Denver, CO, USA: ACM, Nov. 2023.

[20] R. Fan, W. Wang, and X. Chu, “DTC-SpMM: Bridging the gap in accelerating general sparse matrix multiplication with tensor cores,” in Proc. 29th ACM Int. Conf. on Architectural Support for Programming Languages and Operating Systems, Volume 3, ser. ASPLOS ’24. La Jolla, CA, USA: ACM, Apr. 2024, pp. 253–267. [21] P. Okanovic, G. Kwasniewski, P. S. Labini, M. Besta, F. Vella, and T. Hoefler, “High performance unstructured SpMM computation using tensor cores,” in Proc. Int. Conf. for High Perf. Comp., Netw., Storage and Analysis, ser. SC ’24. Atlanta, GA, USA: IEEE, Nov. 2024, pp. 154:1–154:14. [22] H. Zhao, S. Li, J. Wang, C. Zhou, J. Wang, Z. Xin, S. Li, Z. Liang, Z. Pan, F. Liu, Y. Zeng, Y. Wang, and X. Chi, “Acc-SpMM: Accelerating general-purpose sparse matrix-matrix multiplication with GPU tensor cores,” in Proc. 30th ACM SIGPLAN Annual Symp. on Principles and Practice of Parallel Programming, ser. PPoPP ’25. Las Vegas, NV, USA: ACM, Mar. 2025, pp. 326–338. [23] Y. Dong, W. Jiang, X. Shen, H. Guo, Z. Shao, and H. Jin, “BRPSpMM: Block-row partition based sparse matrix multiplication with tensor and CUDA cores,” in Proc. 39th Int. Parallel and Dist. Proc. Symp. (IPDPS ’25). Milan, Italy: IEEE, May 2025, pp. 901–912. [24] L. Xiang, O. Asudeh, G. Sabin, A. Sukumaran-Rajam, and P. Sadayappan, “cuTeSpMM: Accelerating sparse-dense matrix multiplication using GPU tensor cores,” CoRR, vol. abs/2504.06443, 2025. [25] J. Shi, S. Li, Y. Xu, R. Fu, X. Wang, and T. Wu, “FlashSparse: Minimizing computation redundancy for fast sparse matrix multiplications on tensor cores,” in Proc. 30th ACM SIGPLAN Annual Symp. on Principles and Practice of Parallel Programming, ser. PPoPP ’25. Las Vegas, NV, USA: ACM, Mar. 2025, pp. 312–325. [26] Y. Chen, J. Xie, S. Teng, W. Zeng, and J. X. Yu, “Groot: Graph-centric row reordering with tree for sparse matrix multiplications on tensor cores,” in Proc. 20th European Conf. on Comp. Sys., ser. EuroSys ’25. Rotterdam, Netherlands: ACM, Apr. 2025, pp. 803–817. [27] G. Li, Y. Liu, W. Luo, R. Fu, X. Wang, and T. Wu, “FastSpMM: Leveraging tensor cores for sparse matrix multiplication,” in Proc. 22nd ACM Int. Conf. on Comp. Frontiers, ser. CF ’25. ACM, 2025, pp. 195–204. [28] Y. Xia, W. Wang, D. Yang, X. Zhou, and D. Cheng, “Voltrix: Sparse matrix-matrix multiplication on tensor cores withs asynchronous and balanced kernel optimization,” in 2025 USENIX Annual Tech. Conf., ser. USENIX ATC ’25. Boston, MA, USA: USENIX Assoc., Jul. 2025, pp. 699–714. [29] K. Zhang, X. Liu, H. Yang, T. Feng, X. Yang, Y. Liu, Z. Luan, and D. Qian, “Jigsaw: Accelerating SpMM with vector sparsity on sparse tensor core,” in Proc. 53rd Int. Conf. on Parallel Processing, ser. ICPP ’24. Gotland, Sweden: ACM, Aug. 2024, pp. 1124–1134. [30] Q. Wang, Y. Wang, Y. Luo, R. Luo, and P. Tang, “HR-SpMM: Adaptive row partitioning and hybrid kernel design for sparse matrix multiplication,” in Proc. 39th ACM Int. Conf. on Supercomputing, ser. ICS ’25. Salt Lake City, UT, USA: ACM, Jun. 2025, pp. 161–172. [31] C. Ma, D. Wu, Z. Deng, J. Chen, X. Huang, J. Meng, W. Zhu, B. Wang, A. C. Zhou, P. Chen, M. Deng, Y. Wei, S. Feng, and Y. Pan, “NM-SpMM: Accelerating matrix multiplication using N:M sparsity with GPGPU,” in Proc. 39th IEEE Int. Parallel and Dist. Processing Symp., ser. IPDPS ’25. Milan, Italy: IEEE, Jun. 2025, pp. 926–937. [32] Z. Li, X. Ke, Y. Zhu, Y. Gao, and Y. Tu, “HC-SpMM: Accelerating sparse matrix-matrix multiplication for graphs with hybrid GPU cores,” in Proc. 2025 IEEE 41st Int. Conf. on Data Eng. (ICDE). Hong Kong, China: IEEE, May 2025. [Online]. Available: https://ieeexplore.ieee.org/document/11112857 [33] A. Li, J. Sun, H. Li, W. Ji, and G. Sun, “RSH-SpMM: A row-structured hybrid kernel for sparse matrix-matrix multiplication on GPUs,” February 2026. [Online]. Available: https://arxiv.org/abs/2603.08734 [34] Y. Dong, Z. Shen, W. Jiang, Z. Liu, Y. Xu, B. He, R. Zheng, and H. Jin, “Bridging the gap between unstructured SpMM and structured sparse tensor cores,” in Proc. Int. Conf. for High Perf. Comp., Netw., Storage and Analysis (SC ’25). New York, NY, USA: ACM, 2025, pp. 645–660. [Online]. Available: https://dl.acm.org/doi/10.1145/3712285.3759849 [35] Y. Chen and J. X. Yu, “ToT: Triangle counting on tensor cores,” IEEE Trans. Parallel and Dist. Systems, vol. 36, no. 12, pp. 2679–2692, 2025. [36] S. Li, K. Osawa, and T. Hoefler, “Efficient quantized sparse matrix operations on tensor cores,” in Proc. Int. Conf. for High Perf. Comp., Netw., Storage and Analysis, ser. SC ’22. Dallas, TX, USA: IEEE, Nov. 2022, pp. 1–15.

Preprint submitted to IEEE Transactions on Parallel and Distributed Systems

[37] D. Jiao, X. Deng, Z. Wang, S. Fan, Y. Chen, D. Meng, R. Hou, and M. Zhang, “Neo: Towards efficient fully homomorphic encryption acceleration using tensor cores,” in Proc. 52nd Int. Symp. on Comp. Arch., ser. ISCA ’25. Tokyo, Japan: ACM, Jun. 2025, pp. 107–121. [38] R. Nobre, M. Graça, L. Sousa, and A. Ilic, “EPIClear: Exploiting domain-specific features for epistasis detection acceleration on tensor cores,” in Proc. 39th ACM Int. Conf. on Supercomputing, ser. ICS ’25. Salt Lake City, UT, USA: ACM, Jun. 2025, pp. 293–307. [39] Y. Wang, B. Feng, Z. Wang, G. Huang, and Y. Ding, “TC-GNN: Bridging sparse GNN computation and dense tensor cores on GPUs,” in 2023 USENIX Annual Tech. Conf., ser. ATC ’23. Boston, MA, USA: USENIX Assoc., Jul. 2023, pp. 149–164. [40] J.-A. Chen, H.-H. Sung, R. Zhang, A. Li, and X. Shen, “Accelerating GNNs on GPU sparse tensor cores through N:M sparsity-oriented graph reordering,” in Proc. 30th ACM SIGPLAN Annual Symp. on Principles and Practice of Parallel Programming, ser. PPoPP ’25. Las Vegas, NV, USA: ACM, Mar. 2025, pp. 16–28. [41] Y. Han, I. Kim, J. Kim, and G. E. Moon, “Tensor core-adapted sparse matrix multiplication for accelerating sparse deep neural networks,” Electronics, vol. 13, no. 20, p. 3981, 2024. [42] R. L. Castro, A. Ivanov, D. Andrade, T. Ben-Nun, B. B. Fraguela, and T. Hoefler, “VENOM: A vectorized N:M format for unleashing the power of sparse tensor cores,” in Proc. Int. Conf. for High Perf. Comp., Netw., Storage and Analysis, ser. SC ’23. Denver, CO, USA: ACM, Nov. 2023. [43] Z. Ye, R. Lai, J. Shao, T. Chen, and L. Ceze, “SparseTIR: Composable abstractions for sparse compilation in deep learning,” in Proc. 28th ACM Int. Conf. on Architectural Support for Programming Languages and Operating Systems (ASPLOS ’23), Volume 3. Vancouver, BC, Canada: ACM, March 2023, pp. 660–678. [Online]. Available: https://dl.acm.org/doi/10.1145/3582016.3582047 [44] N. Zheng, B. Lin, Q. Zhang, L. Ma, Y. Yang, F. Yang, Y. Wang, M. Yang, and L. Zhou, “SparTA: Deep-learning model sparsity via Tensor-with-Sparsity-Attribute,” in 16th USENIX Symp. on Operating Systems Design and Implementation (OSDI ’22). Carlsbad, CA, USA: USENIX Assoc., July 2022, pp. 213–232. [45] D. Joo, H. Hosseini, R. Hadidi, and B. Asgari, “Coruscant: Codesigning GPU kernel and sparse tensor core to advocate unstructured sparsity in efficient LLM inference,” in Proc. 58th IEEE/ACM Int. Symp. on Microarchitecture, ser. MICRO ’25. Seoul, South Korea: ACM, Oct. 2025. [46] Y. Wang, X. Guo, J. Xiao, D. Chen, and G. Tan, “GeneralSparse: Bridging the gap in SpMM for pruned large language model inference on GPUs,” in 2025 USENIX Annual Tech. Conf., ser. USENIX ATC ’25. Boston, MA, USA: USENIX Assoc., Jul. 2025, pp. 417–432. [47] H. Xia, Z. Zheng, Y. Li, D. Zhuang, Z. Zhou, X. Qiu, Y. Li, W. Lin, and S. L. Song, “Flash-LLM: Enabling cost-effective and highly-efficient large generative model inference with unstructured sparsity,” Proc. VLDB Endowment, vol. 17, no. 2, pp. 211–224, 2023. [Online]. Available: https://dl.acm.org/doi/10.14778/3626292.3626303 [48] R. Fan, X. Yu, P. Dong, Z. Li, G. Gong, Q. Wang, W. Wang, and X. Chu, “SpInfer: Leveraging low-level sparsity for efficient large language model inference on GPUs,” in Proc. 20th European Conf. on Computer Systems (EuroSys ’25). Rotterdam, Netherlands: ACM, March 2025. [Online]. Available: https://dl.acm.org/doi/10.1145/3689031.3717481 [49] Y. Chen, K. Li, Y. Wang, D. Bai, L. Wang, L. Ma, L. Yuan, Y. Zhang, T. Cao, and M. Yang, “ConvStencil: Transform stencil computation to matrix multiplication on tensor cores,” in Proc. 29th ACM SIGPLAN Annual Symp. on Principles and Practice of Parallel Programming, ser. PPoPP ’24. Edinburgh, United Kingdom: ACM, Mar. 2024, pp. 333–347. [50] A. Dakkak, C. Li, J. Xiong, I. Gelado, and W. mei W. Hwu, “Accelerating reduction and scan using tensor core units,” in Proc. ACM Int. Conf. on Supercomputing, ser. ICS ’19. Phoenix, AZ, USA: ACM, Jun. 2019, pp. 46–57. [51] Parallel Thread Execution ISA, NVIDIA Corporation, 2025, PTX ISA Version 9.1, online documentation. [Online]. Available: https://docs.nvidia.com/cuda/parallel-thread-execution/ [52] M. Besta, M. Podstawski, L. Groner, E. Solomonik, and T. Hoefler, “To push or to pull: On reducing communication and synchronization in graph computations,” in Proc. 26th Int. Symp. on High-Performance Parallel and Dist. Comp., ser. HPDC ’17.

New York, NY, USA: ACM, 2017, p. 93–104. [Online]. Available: https://doi.org/10.1145/3078597.3078616 [53] P. Jaccard, “The distribution of the flora in the alpine zone. I,” New Phytologist, vol. 11, no. 2, pp. 37–50, 1912. [Online]. Available: https://doi.org/10.1111/j.1469-8137.1912.tb05611.x [54] M. Harris. (2017, may) Cooperative groups: Flexible CUDA thread programming. NVIDIA Dev. Blog, acc.: 2025-12. [Online]. Available: https://developer.nvidia.com/blog/cooperative-groups/ [55] NVIDIA Corporation, “NVIDIA Nsight Compute,” https:// developer.nvidia.com/nsight-compute, 2024, version 2024.3, accessed: 2026. [56] A. E. Sarıyüce, E. Saule, K. Kaya, and Ü. V. Çatalyürek, “Regularizing graph centrality computations,” Journal of Parallel and Distributed Computing, vol. 76, pp. 106–119, 2015. [Online]. Available: https://doi.org/10.1016/j.jpdc.2014.07.006 [57] U. Brandes, “A faster algorithm for betweenness centrality*,” The Journal of Math. Soc., vol. 25, no. 2, pp. 163–177, 2001. [Online]. Available: https://doi.org/10.1080/0022250X.2001.9990249 [58] T. A. Davis and Y. Hu, “The University of Florida sparse matrix collection,” ACM Trans. Mathematical Software, vol. 38, no. 1, pp. 1:1–1:25, 2011. [Online]. Available: https://doi.org/10.1145/2049662.2049663 [59] H. Liu, H. H. Huang, and Y. Hu, “iBFS: Concurrent breadth-first search on GPUs,” in Proceedings of the 2016 International Conference on Management of Data, ser. SIGMOD ’16. San Francisco, California, USA: ACM, 2016, pp. 403–416. [Online]. Available: https://doi.org/10.1145/2882903.2882959 [60] M. Then, M. Kaufmann, F. Chirigati, T.-A. Hoang-Vu, K. Pham, A. Kemper, T. Neumann, and H. T. Vo, “The more the merrier: Efficient multi-source graph traversal,” Proceedings of the VLDB Endowment, vol. 8, no. 4, pp. 449–460, 2014. [Online]. Available: https://www.vldb.org/pvldb/vol8/p449-then.pdf [61] E. Cohen, D. Delling, T. Pajor, and R. F. Werneck, “Computing classic closeness centrality, at scale,” in Proc. 2nd ACM Conf. on Online Social Networks (COSN 2014). ACM, 2014, pp. 37–50. [62] E. Bergamini, M. Borassi, P. Crescenzi, A. Marino, and H. Meyerhenke, “Computing top-k closeness centrality faster in unweighted graphs,” ACM Transactions on Knowledge Discovery from Data, vol. 13, no. 5, pp. 53:1–53:40, 2019. [63] A. Buluç and J. R. Gilbert, “The combinatorial BLAS: Design, implementation, and app.” The Int. J. of High Perf. Comp. App., vol. 25, no. 4, pp. 496–509, 2011. [Online]. Available: https://doi.org/10.1177/1094342011403516 [64] C. Yang, A. Buluç, and J. D. Owens, “Graphblast: A high-perf. linear algebra-based graph framework on the GPU,” ACM Trans. on Mathematical Software, vol. 48, no. 1, pp. 1:1–1:51, 2022. [Online]. Available: https://doi.org/10.1145/3466795

Deniz Elbek is an undergraduate student in Computer Science and Engineering at Sabanci University. His areas of research focus on High Performance Computing, Parallel Algorithms, and Massively Parallel Computer Architecture.

Kamer Kaya is an Associate Professor of Computer Science at Sabanci University. His research interests are Parallel Algorithms, High Performance Computing, and Graph and Sparse Matrix Algorithms.

Preprint submitted to IEEE Transactions on Parallel and Distributed Systems

Record · ID 259426 · SHA-256 dff909c11fc8ce88
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.