GPU accelerated maximum cardinality matching algorithms for bipartite graphs

We design, implement, and evaluate GPU-based algorithms for the maximum cardinality matching problem in bipartite graphs. Such algorithms have a variety of applications in computer science, scientific computing, bioinformatics, and other areas. To the best of our knowledge, ours is the first study which focuses on GPU implementation of the maximum cardinality matching algorithms. We compare the proposed algorithms with serial and multicore implementations from the literature on a large set of real-life problems where in majority of the cases one of our GPU-accelerated algorithms is demonstrated to be faster than both the sequential and multicore implementations.


Introduction
Bipartite graph matching is one of the fundamental problems in graph theory and combinatorial optimization. The problem asks for a maximum set of vertex disjoint edges in a given bipartite graph. It has many applications in a variety of fields such as image processing [17], chemical structure analysis [15], and bioinformatics [2] (see also another two discussed by Burkard et al. [4,Section 3.8]). Our motivating application lies in solving sparse linear systems of equations, as algorithms for computing a maximum cardinality bipartite matching are run routinely in the related solvers. In this setting, bipartite matching algorithms are used to see if the associated coefficient matrix is reducible; if so, substantial savings in computational requirements can be achieved [7,Chapter 6].
Achieving good parallel performance on graph algorithms is challenging, because they are memory bounded, and there are poor localities of the mem-ory accesses. Moreover, because of the irregularity of the computations, it is difficult to exploit concurrency. Algorithms for the matching problem are no exception. There have been recent studies that aim to improve the performance of matching algorithms on multicore and manycore architectures. For example, Vasconcelos and Rosenhahn [18] propose a GPU implementation of an algorithm for the maximum weighted matching problem on bipartite graphs. Fagginger Auer and Bisseling [10] study an implementation of a greedy graph matching on GPU. Ç atalyürek et al. [5] propose different greedy graph matching algorithms for multicore architectures. Azad et al. [1] introduce several multicore implementations of maximum cardinality matching algorithms on bipartite graphs.
We propose GPU implementations of two maximum cardinality matching algorithms. We analyze their performance and employ further improvements. We thoroughly evaluate their performance with a rigorous set of experiments on many bipartite graphs from different applications. The experimental results conclude that one of the proposed GPU-based implementation is faster than its existing multicore counterparts.
The rest of this paper is organized as follows. The background material, some related work, and a summary of contributions are presented in Section 2. Section 3 describes the proposed GPU algorithms. The comparison of the proposed GPU-based implementations with the existing sequential and multicore implementations is given in Section 4. Section 5 concludes the paper.

Background and contributions
and a set of edges E such that for each edge, one of the endpoints is in V 1 and other is in V 2 . Since our motivation lies in the sparse matrix domain, we will refer to the vertices in the two classes as row and column vertices.
A matching M in a graph G is a subset of edges E where a vertex in V 1 ∪ V 2 is in at most one edge in M. Given a matching M, a vertex v is said to be matched by M if v is in an edge of M, otherwise v is called unmatched. The cardinality of a matching M, denoted by |M|, is the number of edges in M. A matching M is called maximum, if no other matching M with |M | > |M| exists. For a matching M, a path P in G is called an M-alternating if its edges are alternately in M and not in M. An M-alternating path P is called M -augmenting if the start and end vertices of P are both unmatched.
There are three main classes of algorithms for finding the maximum cardinality matchings in bipartite graphs. The first class of algorithms is based on augmenting paths (see a detailed summary by Duff et al. [8]). Pushrelabel-based algorithms form a second class [12]. A third class, pseudoflow algorithms, is based on a more recent work [13]. There are O( √ nτ ) algorithms in the first two classes (e.g., Hopcroft-Karp algorithm [14] and a variant of the push relabel algorithm [11]), where n is the number of vertices and τ is the number of edges in the given bipartite graph. This is asymptotically best bound for practical algorithms. Most of the other known algorithms in the first two classes and in the third class have the running time complexity of O(nτ ). These three classes of algorithms are described and compared in a recent study [16]. It has been demonstrated experimentally that the champions of the first two families are comparable in performance and better than that of the third family. Since we investigate GPU acceleration of augmenting-path-based algorithms, a brief description of them is given below (the reader is invited to two recent papers [8,16] and the original resources cited in those papers for other algorithms).
Algorithms based on augmenting paths follow the following common pattern. Given an initial matching M (possible empty), they search for an M-augmenting path P. If none exists, then the matching M is maximum by a theorem of Berge [3]. Otherwise, the augmenting path P is used to increase the cardinality of M by setting M = M ⊕ E(P) where E(P) is the edge set of the path P, and M ⊕ E(P) = (M ∪ E(P)) \ (M ∩ E(P)) is the symmetric difference. This inverts the membership in M for all edges of P. Since both the first and the last edge of P were unmatched in M, we have |M ⊕ E(P)| = |M| + 1. The augmenting-path-based algorithms differ in the way these augmenting paths are found and the associated augmentations are realized. They mainly use either breadth-first-search (BFS), or depth-firstsearch (DFS), or combination of these two techniques to locate and perform the augmenting paths.
Multicore counterparts of a number of augmenting-path based algorithms are proposed in a recent work [1]. The parallelization of these algorithms is achieved by using atomic operations at BFS and/or DFS steps of the algorithm. Although using atomic operations might not harm the performance on a multicore machine, they should be avoided in a GPU implementation because of very large number of concurrent thread executions.
As a reasonably efficient DFS is not feasible with GPUs, we accelerate two BFS-based algorithms, called HK [14] and HKDW [9]. HK has the best known worst-case running time complexity of O( √ nτ ) for a bipartite graph with n vertices and τ edges. HKDW is a variant of HK and incorporates techniques to improve the practical running time while having the same worst-case time complexity. Both of these algorithms use BFS to locate the shortest augmenting paths from unmatched columns, and then use DFSbased searches restricted to a certain part of the input graph to augment along a maximal set of disjoint augmenting paths. HKDW performs another set of DFS-based searches to augment using the remaining unmatched rows. As is clear, the DFS-based searches will be a big obstacle to achieve efficiency. In order to overcome this hurdle, we propose a scheme which al-ternates the edges of a number of augmenting paths with a parallel scheme that resembles to a breadth expansion in BFS. The proposed scheme offers a high degree of parallelism but does not guarantee a maximal set of augmentations, potentially increasing the worst case time complexity to O(nτ ) on a sequential machine. In other words, we trade theoretical worst case time complexity with a higher degree of parallelism to achieve better practical running time with a GPU.

Methods
We propose two algorithms for the GPU implementation of maximum cardinality matching. These algorithms use BFS to find augmenting paths, speculatively perform some of them, and fix any inconsistencies that can be resulting from speculative augmentations.
The overall structure of the first GPU-based algorithm is given in Algorithm 1, APsB. It largely follows the common structure of most of the existing sequential algorithms, and corresponds to HK. It performs a combined BFS starting from all unmatched columns to find unmatched rows, thus locating augmenting paths. Some of those augmentations are then realized using a function called Alternate (will be described later). The parallelism is exploited inside the InitBfsArray, BFS, Alternate, and if augmenting path f ound then The outer loop of Algorithm 1 iterates until no more augmenting paths are found, thereby guaranteeing a maximum matching. The inner loop is responsible from completing the breadth-first-search of the augmenting paths. A single iteration of this loop corresponds to a level of BFS. The inner loop iterates until all shortest augmenting paths are found. Then, the edges in these shortest augmenting paths are alternated inside Alternate function. Unlike the sequential HK algorithm, APsB does not find a maximal set of augmenting paths.
By removing the lines 9 and 10 of Algorithm 1, another matching algorithm is obtained. This method will continue with the BFSs until all possible unmatched rows are found; it can be therefore considered as the GPU implementation of the HKDW algorithm. This variant is called APFB.

Algorithm 2: BFS Kernel Function-1 (GPUBFS)
Data: bf s level, bf s array, cxadj, cadj, nc, rmatch, vertex inserted, augmenting path f ound 1 process cnt ← getP rocessCount(nc); 2 for i from 0 to process cnt − 1 do The GPU threads partition the columns vertices in a single dimension. Each thread with id tid is assigned a number of columns which is obtained via the following function: Once the number of columns are obtained, the threads traverse their first assigned column vertex. The indices of the columns assigned to a thread differ by tot thread num to allow coalesced global memory accesses. Threads traverse the neighboring row vertices of the current column, if its BFS level is equal to the current bf s level. If a thread encounters a matched row during the traversal, its matching column is retrieved. If the column is not traversed yet, its bf s level is marked on bf s array. On the other hand, when a thread encounters an unmatched row, an augmenting path is found. In this case, the match of the neighbor row is set to −2, and this information is used by Alternate later.

Algorithm 3: Alternate
Data: cmatch, rmatch, nc, nr, predecessor 1 process vcnt ← getP rocessCount(nr); 2 for i from 0 to process vcnt − 1 do  Algorithm 3 gives the description of the Alternate function. This kernel alternates the matching edges with unmatching edges of the augmenting paths found; some of those paths end up being augmenting ones and some are only partially alternated. Here, each thread is assigned a number of rows. Since rmatch of an unmatched row (that is also an endpoint of an augmenting path) has been set to −2 in the BFS kernel, only the threads whose row vertex's match is −2 start Alternate. Since there might be several augmenting paths for an unmatched column, race conditions while writing on cmatch and rmatch arrays are possible. Such a race condition might cause infinite loops (inner while loop) or inconsistencies, if care is not taken. We prevent these by checking the predecessor of a matched row (line-8). For example, in Fig. 1, two different augmenting paths that end with r 2 and r 3 are found for c 1 . If the thread of r 2 starts before the thread of r 3 in Alternate, the match of c 2 will be updated to r 2 (line-10). Then, r 3 's thread will read matched row of c 2 as r 2 (line-7). This would cause an infinite loop without the check at line-8. Inconsistencies may occur when the threads of r 2 and r 3 are in the same warp. In this case, the if-check will not hold for both threads, and their row vertices will be written on cmatch (line-10). Since only one thread will be successful at writing, this will cause an inconsistency. Such inconsistencies are fixed by This array holds the root (as the index of the column vertex) of an augmenting path, and this information is transferred down during BFS. Whenever an augmenting path is found, the entry in bf s array for the root of the augmenting path is set to L 0 − 2. This information is used at the beginning of BFS kernel. No more BFS traversals is done, if an augmenting path is found for the root of the traversed column vertex. Therefore, while the method increases the global memory accesses by introducing an extra array, it provides an early exit mechanism for BFS.
We further improve GPUBFS-WR by making use of the arrays root and bf s array. BFS kernels might find several rows to match with the same unmatched column, and set rmatch[·] to −2 for each. These cause Alternate to start from several rows that can be matched with the same unmatched column. Therefore, it may perform unnecessary alternations, until these augmenting paths intersect. Conflicts may occur at these intersection points (which are then resolved with FixMatching function). By choosing L 0 as 2, we can limit the range of the values that bf s array takes to positive numbers. Therefore, by setting the bf s array to −(neighbor row) at line 19 of Algorithm 4, we can provide more information to the Alternate function. With this, Alternate can determine the beginning and the end of an augmenting path, and it can alternate only among the correct augmenting paths. APsB-GPUBFS-WR (and Alternate function used together) includes these improvements. However, they are not included in APFB-GPUBFS-WR since they do not improve its performance.

Experiments
The running time of the proposed implementations are compared against the sequential HK and PFP implementations [8], and against the multicore  is used to initialize all tested algorithms. We compare the running time of the matching algorithms after this common initialization. Two different main algorithms APFB and APsB can use two different BFS kernel functions (GPUBFS and GPUBFS-WR). Moreover, each of these algorithms can have two versions (i) CT: uses a constant number of threads with fixed number of grid and block size (256 × 256) and assigns multiple vertices to each thread; (ii) MT: tries to assign one vertex to each thread. The number of threads used in the second version is chosen as MT = min(nc, #threads) where nc is the number of columns, and #threads is the maximum number of threads of the architecture. Therefore, we have eight GPU-based algorithms.
The algorithms are run on bipartite graphs corresponding to 70 different matrices from variety of classes at UFL matrix collection [6]. We also permuted the matrices randomly by rows and columns and included them as a second set (labeled RCP). These permutations usually render the problems harder for the augmenting-path-based algorithms [8]. For both sets, we report the performance for a smaller subset which contains those matrices in which at least one of the sequential algorithms took more than one second. We call these sets O S1 (28 matrices) and RCP S1 (50 matrices). We also have another two subsets called O Hardest20 and RCP Hardest20 that contain the set of 20 matrices on which the sequential algorithms required the longest running time.
First, we compare the performance of the proposed GPU algorithms.   Table 1 shows the geometric mean of the running time on different sets. As we see from the table, using constant number of threads (CT) always increases the performance of an algorithm, since it increases the granularity of the work performed by each thread. GPUBFS-WR is always faster than GPUBFS. This is because of the unnecessary BFS traversals in the GPUBFS algorithm. GPUBFS cannot determine whether an augmenting path has already been found for an unmatched column, therefore it will continue to explore. This unnecessary BFS traversals not only increase the BFS time, but also reduce the likelihood of finding an augmenting path for other unmatched columns. Moreover, the Alternate scheme turns out to be more suitable for APFB than APsB, in which case it can augment along more paths (there is a larger set of possibilities). For example, Figs. 2(a) and 2(b) show the number of BFS iterations and the number of BFS levels in each iteration for, respectively, Hamrle3 and Delanuay n23 graphs. As clearly seen from both of the figures, APFB variants converges in smaller number of iterations than APsB variants; and for most of the graphs, the total number of BFS kernel calls are less for APFB (as in Fig. 2(a)). However, for a small subset of the graphs, although the augmenting path exploration of APsB converges in larger number of iterations, the numbers of the BFS levels in each iterations are much less than APFB (as in Fig. 2(b)). Unlike the general case, APsB outperforms APFB in such cases. Since APFB using GPUBFS-WR and CT almost always obtains the best performance, we only compare the performance of this algorithm with other implementations in the following. Figures 3(a) and 3(b) give the log-scaled speedup profiles of the best GPU and multicore algorithms on the original and permuted graphs. The speedups are calculated with respect to the fastest of the sequential algorithms PFP and HK (on the original graphs HK was faster; on the permuted ones PFP was faster). A point (x, y) in the plots corresponds to the probability of obtaining at least 2 x speedup is y. As the plots show, the GPU algorithm has the best overall speedup. It is faster than the sequential HK  algorithm for 86% of the original graphs, while it is faster than PFP on 76% of the permuted graphs. P-DBFS obtains the best performance among the multicore algorithms. However, its performance degrades on permuted graphs. Although P-PFP is more robust than P-DBFS to permutations, its overall performance is inferior to that of P-DBFS. P-HK is outperformed by the other algorithms in both sets.  show the performance profiles of the GPU and multicore algorithms. A point (x, y) in this plot means that with y probability, the algorithm obtains a performance that is at most x times worse than the best running time. The plots clearly show the separation among the GPU algorithm and the multicore ones, especially for original graphs and for x ≤ 7 for the permuted ones, thus marking GPU as the fastest in most cases. In particular, GPU algorithm obtains the best performance in 61% of the original graphs, while this ratio increases to 74% for the permuted ones. Figure 5 gives the overall speedups. The proposed GPU algorithm obtains average speedup values of at least 3.61 and 3.54 on, respectively, original and permuted graphs. The speedups increase for the hardest instances, where the GPU algorithm achieves 3.96 and 9.29 speedup, respectively, on original and permuted graphs.  on the complete set of graphs can be found at http://bmi.osu.edu/hpc/ software/matchmaker2/maxCardMatch.html). As seen from this table, except six instances among the original graphs and another two among the permuted graphs, the GPU algorithm is faster than the best sequential algorithm. It is also faster than the multicore ones in all, except five original graphs.

Concluding remarks
We proposed a parallel BFS based GPU implementation of maximum cardinality matching algorithm for bipartite graphs. We presented experiments on various datasets, and compared the performance of the proposed GPU implementation against sequential and multicore algorithms. The experiments showed that the proposed GPU implementations are faster than the existing parallel multicore implementations. The speedups achieved with respect to well-known sequential implementations varied from 0.03 to 629.19, averaging 9.29 on a set of 20 hardest problems with respect to the fastest sequential algorithm. A GPU is a restricted memory device. Although, an out-of-core or distributed-memory type algorithm is amenable when the graph does not fit into the device, a direct implementation of these algorithms will surely not be efficient. We plan to investigate the techniques to obtain good matching performance for extreme-scale bipartite graphs on GPUs.