Engineering A Workload-balanced Push-Relabel Algorithm for Massive Graphs on GPUs

Given a directed graph $G(V,E)$, where $V$ and $E$ are the vertex and weighted edge set, respectively, each edge $(u,v)$ has a weight, capacity $c(u,v)$, representing the maximum amount of flow that can be passed through this edge. The maximum flow problem tries to find the maximum flow value from the source vertex $s$ to the sink vertex $t$. Generally, algorithms solve the maximum flow problem operating on the residual graph $G_{f}(V,E_{f})$, which has the same vertex set as the given $G$ with residual edges $E_{f}$. The weight $c_{f}(u,v)$ of a residual edge $(u,v)$ represents how much-remaining flow can be passed from $u$ to $v$. Note that when there is a flow from $u$ to $v$, there will be a residual edge from $v$ to $u$, even if there is no $(v,u)$ in $E$.

The foundational algorithm solving the maximum flow problem is the Ford-Fulkerson algorithm [6], which finds augmenting paths from the source to the sink iteratively. It pushes the additional flow until no such path remains. Edmonds and Karp improved the efficiency of finding augmenting paths using breadth-first search (BFS). The Edmonds-Karp algorithm [5] reduces the runtime complexity of Ford-Fulkerson from $O(Ef)$ to $O(VE^{2})$, where $f$ is the maximum flow value. Dinic et al. [4] further improved upon the basic augmenting path approach by employing a level graph and blocking flows, which significantly reduced the number of augmentations needed; hence, reducing the complexity to $O(V^{2}E)$ and $O(V^{3/2}E)$ on general and unit-capacity graph respectively. Goldberg et al. [7] proposed the push-relabel algorithm (also known as the preflow-push algorithm), which focuses on local operations to adjust the flows. Some optimizations are based on it to improve the complexity to $O(V^{2}E)$.

The basic idea of the push-relabel algorithm is to generate as much flow as possible at the source, and gradually push it to the sink. It introduces the excess flow concept that allows a vertex to have multiple flows coming into it; namely, it allows a vertex to have more flow coming in than passing out during the push-relabel procedure, this vertex is called active. We denote $e(v)$ as the excess flow value on vertex $v$. The procedure of the push-relabel algorithm is to find active vertices and apply push and relabel operations on them until there are no remaining active vertices. To find the initial active vertices, the algorithm pushes flow from the source to all its neighbor vertices as much as possible, called preflow. A push operation then forces an active vertex to discharge its excess flow and pushes to its neighbors in $G_{f}$. To avoid the endless pushing, the push-relabel algorithm also introduces the height function $h$ on each vertex. At first, the source height is $|V|$, and the height of vertices except the source is 0. It forces an active vertex $u$ to push to the vertex $v$ iff $h(u)=h(v)+1$. If no neighbor vertex satisfies the constraint, the active vertex will relabel its height by finding the minimum height $h^{\prime}$ of its neighbor vertices and setting $h(u)\leftarrow h^{\prime}+1$. With the relabel operation, the height of an active vertex, which cannot push its excess flow at the end, will be increased. The active vertex will be deactivated once its height exceeds $|V|$. The whole procedure ends when there are no active vertices.

The push-relabel algorithm has not only the best theoretical and practical complexity, but it is well-suitable to parallelization due to its inherent structure and the local nature of its operations. Numerous works have developed the parallel version on the different hardware platforms, such as multiprocessor [2, 10], and graphics processing unit (GPU) [9, 24, 14]. Since there are different trade-offs and optimization used in different algorithms, we target the lock-free algorithm proposed by He et al. [9], which is the state-of-the-art algorithm operating on the GPU shown in 1. It parallelly checks whether a vertex $u$ is active by assigning a thread in a GPU (Line: 9). For convenience, we will use thread $u$ to represent the thread, which checks the vertex $u$. If the vertex is active, the thread $u$ will find its neighbor vertex $v$ whose height is minimum among other neighbor vertices (line:10 - 13). The thread then push flow from $u$ to $v$ (line:15 - 19) when $h(u)>h(v)$; otherwise, the thread will relabel the active vertex $u$. After cycles times of iteration, the algorithm applies global-relabeling heuristic, one of the push-relabel optimizations, to improve the practical performance [9]. The global relabeling updates the height of each vertex by performing a backward breadth-first search (BFS) from the sink to the source in the residual graph $G_{f}$. The height will be reassigned to the shortest distance from the sink. After global relabeling, the procedure subtracts Excess_total from the excess flow of those inactive vertices to guarantee the termination of the procedure (line: 6 in Algorithm 1). Since the lock-free algorithm assigns a thread to process a vertex, we also call this a thread-centric approach.

The most significant difference between the original and lock-free push-relabel algorithms is that the lock-free one relaxes the push constraint of the height function $h(v)=h(u)+1$ to $h(v)>h(u)$ (line: 14). By finding the minimum-height neighbor and using atomic operations, the correctness of the lock-free algorithm has been proven in [10].

For the convenience of description, we use the CUDA (Compute Unified Device Architecture, proposed by Nvidia Corporation) terminologies. The modern GPU extends SIMD (Single-Instruction-Multiple-Data) to SIMT (Single-Instruction-Multiple-Thread), where a warp is the primary computing unit. A warp usually consists of 32 threads, which share a programming counter (PC) and execute in a lockstep manner; namely, all threads in a warp will execute and access the memory simultaneously; hence, the if-else condition inside a warp will cause warp divergence and serialize the execution of a warp, leading to tremendous overhead. On the top of warps, the thread block (TB) executes multiple warps concurrently to hide the latency of memory access on each warp. Note, all threads in a TB share two fast on-chip memory, shared memory and cache. The former memory is programmable, while the latter is automatic caching. In addition to this cached memory, the memory coalescing technique can also improve memory access efficiency. If threads inside warp access the memory continuously, the hardware will combine these accesses into one request to increase the memory bandwidth. Multiple TBs comprise an entire grid, which maps to a physical GPU hardware and equips up with global memory for data used in a GPU.

Despite the parallel lock-free algorithm mentioned in Section 2.2 has gained excellent performance improvement on GPU, the thread-centric approach still faces significant workload imbalance, especially when the underlying graph becomes larger and larger. Specifically, the prior approach fails to fully utilize the parallelism of the GPU. Before diving into the imbalance problem, we want to derive the cost model of the thread-centric approach on GPU to help us explain where the workload imbalance happens. The execution time model can be roughly represented by the following equation:

Figure 1 (b) shows the architecture for our workload-balanced push-relabel implementation. As mentioned in the Section 2, we aim to alleviate the two workload imbalances using the thread-centric approach. Fundamentally, we first assign all threads to scan all vertices to find the active ones and add the active vertices to the active vertex queue (AVQ). With the AVQ, each thread has an equal workload when finding active vertices. Besides, the active vertices can be coalesced in the AVQ; hence, we can assign a tile (a group of threads) to find the minimum-height neighbor vertex of an active vertex, which reduces the time of the searching time (Algorithm 1 line: 11-13) from $O(d(v))$ to $O(log_{2}{d(v)})$. Since we can process multiple vertices and use multiple threads in an active vertex simultaneously, we call this approach two-level parallelism. Furthermore, we designed two enhanced CSRs, reversed CSR and bidirectional CSR, to reduce significant time in scanning all neighbors or finding backward edges on the residual graph.

The prior approach used an adjacency matrix to represent a residual graph, so it takes $O(d(u))$ time to find all in- and out-neighbor vertices with a given vertex $u$. Using an adjacency matrix is convenient and less time-consuming in finding the minimum-height neighbor vertex, but it has $O(V^{2})$ memory complexity, which puts considerable pressure on the GPU. More and more people have embraced the compressed sparse representation (CSR) to reduce the memory consumption of a massive graph, especially the sparse graph. However, using a conventional CSR as the residual graph is inefficient for the push-relabel algorithm. As mentioned in Section 2.2, an active vertex must find the minimum-height neighbor $v$ in the residual graph. After finding the minimum-height neighbor, we need to find the backward edge of the $(v,u)$, so that we can decrease $flow(u,v)$ and increase $flow(v,u)$. If we directly put the backward edge below the forward edge shown in Figure 2 (b), it can access the backward edge in a constant time, but it takes $O(|E|)$ to find an active vertex’s neighbors. For instance, finding neighbors of vertex $2$ requires scanning all orange blocks to find $(2,0)$.

To alleviate the inefficiency of finding incoming neighbor vertices of a given vertex, we designed reversed CSR (RCSR) and bidirectional CSR (BCSR) for the residual graph, shown in Figure 2 (b), (c), respectively. The RCSR uses another CSR for backward edges. The flow_idx records the index of backward flow rather than the value. In this way, we can find all vertex 2’s neighbors by scanning the original CSR and reversed CSR (the orange part of Figure 2 (b)). Unfortunately, we found that accessing RCSR puts tremendous pressure on the memory bandwidth since neighbors of a vertex are stored in discontinuous addresses, causing uncoalesced memory access.

We aggregated incoming and outgoing neighbors to achieve better locality and proposed BCSR shown in Figure 2 (d). Since the neighbor is aggregated, we cannot access the backward edge in a constant amount of time. For instance, the edge $(4,2)$ belongs to vertex $4$’s neighbor; Hence, we need to do an additional search to find the backward edge of $(2,4)$. If we sort the column list in ascending order by vertex ID, we can reduce the searching time from $O(d(v))$ to $O(log_{2}{d(v)})$ with binary search, where $d(v)$ is the degree of a given vertex $v$. We evaluate the performance of both RCSR and BCSR in Section 4.

Algorithm 2 shows the two-level parallelism algorithm for an iteration in Algorithm.1 (line: 9-21). In each iteration, we first use the atomic_add() operation to add the active vertex to the AVQ. Since we put a global synchronization in Line:5 (grid_sync()), we can re-organize the thread assignment for the search of minimum-height neighbor. Besides, we can early break the while loop (Algorithm 1, line: 8) when there is no active vertex remaining in the AVQ, which avoids the redundant iteration. We use a warp as a tile for an active vertex to parallelize the search of a minimum-height neighbor vertex. Since the neighbor of an active vertex stored in the CSR is continuous, using a warp can accelerate the process and reduce the memory request with inherent memory coalescing mentioned in Section 2.3. We use the parallel reduction proposed in [8] to find the vertex with minimum height. We chose the Kernel 7 implementation in [8] since it has the best bandwidth speedup among the seven kernels. Because the finding minimum operation has a very low arithmetic intensity, the peak bandwidth grain can benefit our implementation the most. After finding the minimum-height neighbor, we will use the delegated thread in a warp, whose localIdx is 0, to execute the push or relabel operations.

• •

  TC+RCSR. The thread-centric approach with the RCSR.

• •

  TC+BCSR. The thread-centric approach with the BCSR.

• •

  VC+RCSR. Our vertex-centric approach with the RCSR.

• •

  VC+BCSR. Our vertex-centric approach with the BCSR.

We measured the kernel execution time of different algorithms. Table 1 and Table 2 show the performance of the algorithms in maximum flow and bipartite matching tasks, respectively. At first glance, our vertex-centric approach can improve the execution of both RCSR and BCSR in the two tasks. It gains an average 2.49x and 7.31x execution time speedup in maximum flow tasks on RCSR and BCSR, respectively. In bipartite matching, our VC approach achieves 2.29x and 1.89x execution time speedup on RCSR and BCSR. We observed two critical points in these experiments. First, The VC method is suitable for graphs with a high standard deviation of degree, and these graphs should be a manageable size, or the overhead of synchronization can offset the workload balance benefits of VC. For example, the B0-B2 graphs are too small, so the speedup is marginal. The performance degradation of graphs R0-R2, S0, S1, and B12 comes from the characteristics of the graph itself. The Amazon0302 (R0) is the co-purchased product network on Amazon.com. If a product $v$ is frequently co-purchased with product $u$, the graph contains a directed edge from $i$ to $j$. In this network, almost all nodes are within the same SCC (Strongly Connected Component), and the degrees of these nodes are very close to each other; hence, the workload of this network is naturally balanced. The R1 and R2 are two real-world road networks, so their maximum degree is less than 10. Since we use a tile (usually more than 32 threads per tile) to find neighbors of an active vertex, the active node with a small degree will idle most of the threads within a tile and decrease the utilization of the GPU. On the contrary, if there are some nodes whose degrees are enormous, our vertex-centric approach can effectively alleviate the imbalanced workload, such as R5, B7, and B8.

The second observation is that the BCSR outperforms the RCSR in the maximum flow problem. Nevertheless, there is a slight degradation in the bipartite matching problem. This result recalls the notion in Section 3 that the RCSR has constant accessing time of backward edges, while the BCSR has a better locality than the RCSR. In the bipartite matching task, since the average degree of a graph is higher than that in the maximum flow task, the RCSR performs well because of its fast access to backward edges without searching. On the other hand, the better locality of BCSR can reduce the memory request number by memory coalescing, leading to better performance.

To prove that the performance improvement of our vertex-centric approach comes from the workload balance during execution, we used the method proposed in [1] and measured each warp’s execution time. By assigning the delegated (first) thread within a warp for recording the execution timestamp in the kernel function, we can draw the workload distribution for the TC and VC approach across 13 bipartite graphs in Figure 3. Note that both configurations operate on RCSR. Based on the result, we can make two observations.

First, even though the chart does not directly indicate a shorter execution time for VC (after being normalized by the mean), the vertex-centric approach does indeed reduce the standard deviation of execution times across warps, thus achieving a more even distribution of work.

Second, in smaller graphs, even if we equalize the uneven distribution of work among warps, the overall performance may still decrease due to excessive synchronization. The B0, B1, and B2 are the cases.