Behrouz Maham
email: bmaham@ut.ac.ir

Walid Saad

Mérouane Debbah
email: merouane.debbah@supelec.fr

Zhu Han

Interference Analysis for Spatial Reused Cooperative Multihop Wireless Networks

We consider a decode-and-forward based wireless multihop network with a single source node, a single destination node, and N intermediate nodes. To increase the spectral efficiency and energy efficiency of the system, we propose a cooperative multihop communication with spatial reuse, in which interference is treated as noise. The performance of spatial-reused space-time coded cooperative multihop network is analyzed over Rayleigh fading channels. More specifically, the exact closedform expression for the outage probability at the nth receiving node is derived when there are multiple interferences over noni.i.d. Rayleigh fading channels. Moreover, in high SNR scenario, closed-form asymptotic formulas for the outage probability are derived, from which, we show that the full-spatial diversity is still achievable given interferences from the transmission of concurrent packets. In addition, we propose a simple power control scheme which is only dependent on the statistical knowledge of channels. Finally, the analytic results were confirmed by simulations. It is shown by simulations that the spatialreused multihop transmission outperforms the interference-free multihop transmission in terms of energy efficiency in low and medium SNR scenarios.

I. INTRODUCTION

Cooperative multihop wireless systems have been considered as the promising technique to extend coverage area and reduce power consumption [START_REF] Khandani | Cooperative routing in static wireless networks[END_REF], [START_REF] Maham | Energy-efficient spacetime coded cooperation in outage-restricted multihop wireless networks[END_REF]. In [START_REF] Maham | Energy-efficient spacetime coded cooperation in outage-restricted multihop wireless networks[END_REF], three cooperative multihop transmission protocols were proposed that compromise between spectral and energy efficiencies. To further increase of multiplexing gain and energy efficiency, in this paper, we consider a cooperative multihop transmission with interference due to the simultaneous transmission of multiple packets. The idea of multihop transmission with spatial reuse is proposed in [START_REF] Sikora | Bandwidth-and power-efficient routing in linear wireless networks[END_REF]. To facilitate concurrent transmission of several packets in the network, the available bandwidth is reused among transmitters, with a minimum division of K nodes between simultaneously transmitting nodes. Therefore, we have to deal with a type of co-channel interference (CCI).

The performance analysis of multihop transmission in Rayleigh fading channels under CCI were recently studied in the number of literature such as [START_REF] Suraweera | Performance analysis of two hop amplify-and-forward systems with interference at the relay[END_REF]- [START_REF] Ikki | Multihop wireless relaying systems in the presence of cochannel interferences: Performance analysis and design optimization[END_REF]. However, to the best of our knowledge, this is the first work that investigate the performance analysis of multihop networks with multiple interferences over non-i.i.d. Rayleigh fading. This is of primarily importance for the study of interference due to spatial reused cooperative multihop transmission.

In this paper, we study the performance analysis of the decode-and-forward based cooperative multihop transmission with interference due to the concurrent transmission of multiple data. We show that the capacity of the cooperative multihop transmission can be improved by using spatial reuse scheme. The achievable rate of the multihop transmission can be increased up to ⌊ N +1 K ⌋ times, where K is the minimum separation of concurrently transmitting nodes in a network with N relays, in expense of performance degradation. Moreover, we derive a closed-form expression for the outage probability of the cooperative multihop system in presence of interferences due to the spatial reuse over Rayleigh fading channels. The simplicity of the calculated expression can give insights on performance of the system and ways to optimize the system. In addition, the asymptotic formulas for different signal-tonoise ratio (SNR) and interference-to-noise ratio (INR) conditions are derived. Furthermore, we formulate the problem of minimizing the transmit power for an outage-restricted equal power multihop network under the the assumption of no instantaneous CSI knowledge at the transmitters.

II. SYSTEM MODEL AND PROTOCOL DESCRIPTION

Consider a wireless communication network in which the source s intend to transmit its data to the destination d with the help of N cascaded intermediate nodes. Due to the broadcast nature of the wireless channel, some intermediate relays can overhear and retransmit the received packets. The channel between any two nodes in the network is assumed to be a Rayleigh fading. Similar to [START_REF] Khandani | Cooperative routing in static wireless networks[END_REF], each transmission could be either a broadcast transmission where one node transmits the signal that is heard by multiple receivers, or a cooperative transmission where multiple nodes concurrently transmit the signal to a single receiving node. Here, we employ the cooperation protocol proposed in [START_REF] Maham | Energy-efficient spacetime coded cooperative routing in multihop wireless networks[END_REF] where consists of N + 1 transmission phases. We assume there is no CSI knowledge at transmitters and only statistical CSI is available at the transmitters. Thus, distributed space-time coded transmissions like codes proposed in [START_REF] Maham | Distributed GABBA spacetime codes in amplify-and-forward relay networks[END_REF] is a feasible method to be employed for the cooperative transmission.

The transmission protocols have two major phases of noncooperative and cooperative stages. Depending on the requirements, non-cooperative phase may contain one are multiple steps. The next phases employ space-time cooperated transmission. As an example, Fig. 1 depicts a protocol employing distributed quasi-orthogonal space time code (D-QOST) with M = 4, where M is the number of cooperating nodes. The detailed description of cooperative multihop protocols, i.e. Broadcast-then-Cooperate, Multihop-then-Cooperate, and Full-Cooperation, is studied in [START_REF] Maham | Efficient cooperative protocols for general outage-limited multihop wireless networks[END_REF]. To be consistent through the paper, from now on, we consider the Multihop-then-Cooperate protocol illustrated in Fig. 2. However, the proposed procedure can be easily modified using two other protocols. Assuming the usage of full-rate distributed space-time codes, the number of cooperating nodes is equal to the transmitting packets. Hence, the source node intends to transmit M packets to the destination. The signals transmitted by the source terminal during the mth time slot of Phase 1 is denoted as s m (t), m = 1, . . . , M where t is the time index and is indicated as a group of M packets transmitted at a given time, and E{s m (t)} = 0 and E{|s m (t)| 2 } = 1 for m = 1, . . . , M . In Phase 1, the source transmits the information, and the signal received at the ith node in the first M time slots is given by

y i,m (t) = √ P 0,1 h 0,i s m (t) + v i,m , for m = 1, . . . , M, (1)
where P 0,1 is the average transmit power of the source symbol in the first phase, by assuming that the data symbols s m (t)'s have zero mean and unit variance, and v i,m denotes complex zero-mean white Gaussian noise with variance N 0 . The link coefficients from the jth node to the ith node h j,i , j = 0, 1, . . . , N , i = 1, 2 . . . , N + 1, are complex Gaussian random variables with zero-mean and variances σ 2 j,i , where the (N + 1)th node is the destination d. We assume coherence times of the channels are such that channel coefficients h j,i are not varying during M consecutive time slots. Note that σ 2 j,i = (d 0 /d j,i ) δ is the path-loss coefficient, where d j,i is the distance between nodes j and i, δ is the path-loss exponent, which typically lies in the range of 2 ≤ δ ≤ 6, and d 0 depends on the operating frequency.

For the first M -1 phases, we repeat the non-cooperative transmission described above for the nth relay, 1 ≤ n < M -1 to retransmit the source data. After initial stages, the cooperative transmission is utilized for routing the source packets to the destination. In contrast to [START_REF] Khandani | Cooperative routing in static wireless networks[END_REF] where transmitters are able to modify their phases, here, the instantaneous CSI is not known at the transmitter nodes. This assumption is realistic for most wireless systems. Hence, space-time coded cooperation is the appropriate choice to achieve the spatial diversity gain. In Phase n, M ≤ n, the previous M nodes transmit their signals concurrently toward the next node using an appropriate distributed space-time code, with a simple power allocation.

To facilitate simultaneous transmission of several packets in the network, the available bandwidth is reused among transmitters, with a minimum separation of K nodes between concurrently transmitting nodes. Fig. 4 shows a spatial-reused cooperative multihop network with M = 2. Since in Phase n ≥ M , M nodes transmit the intended data to the nth node, and by assuming half-duplex transmission, the minimum value of the spatially-reused factor K is M + 1, and thus, M + 1 ≤ K ≤ N + 1. For the message detection, the nth node consider all received signals not coming from the M previous nodes as Gaussian interference. In this work, we treat the interference as additive Gaussian noise. In the presence of inter-network interference from the spatial-reused nodes, for n ≥ M the received signal at the n-th receiving node can be represented as

y n,m = M ∑ m=1 P n-m,n h n-m,n x n,m (t) + ∑ u∈Un M ∑ m=1 P u-m,u h u-m,n x u,m (t -τ ) + v n,m , (2)
where x n,m (t) is the zero-mean space-time coded signal, normalized as E{|x n,m (t)| 2 } = 1 during the whole packet transmission, and P n-m,n , m = 1, . . . , M , is the average transmit power of node nm during the mth time slot of Phase n. In (2), U n denotes the set of nodes transmitting simultaneously with Nodes nm, m = 1, . . . , M , due to spatial reuse, i.e.,

U n = {u ∈ {1, 2, . . . , N + 1}|u ̸ = n and K divides n -u} , (3)
and τ = n-u K . The objective of the system In this system, we are interested in the reliable delivery of messages at a rate of R bits/second/Hertz by consuming the minimum total transmit power.

III. PERFORMANCE ANALYSIS OF SPATIAL REUSED COOPERATIVE MULTIHOP TRANSMISSION

In the following, the outage probability ρ out n Pr{r n < R} of the nth receiving node at the nth hop in spatial-reused system is derived, which describes the probability that the transmit rate R is larger than the supported rate r n . This probability depends on the fixed transmission parameters and the channel condition within the hops.

In the cooperative multihop transmission with spatial reuse factor of K, from (1) and ( 2), the instantaneous achieved rate at the nth hop becomes

r n = 1 K log ( 1 + ∑ M m=1 P n-m,n |h n-m,n | 2 N 0 W + ∑ u∈Un ∑ M m=1 P u-m,u |h u-m,n | 2
) ,

(4) where P n-i,n = 0 for n = 1, . . . , M -1, i = 2, . . . , M .

A. Outage Probability

Now, we calculate an exact closed-form expression for the outage probability at the nth receiving node in presence of interference from the multiple-antenna secondary BS. By defining γ th 2 RK -1, the outage probability of noncooperative transmission can be represented as

ρ out n = Pr { ∑ M m=1 P n-m,n |h n-m,n | 2 N 0 W + ∑ u∈Un ∑ M m=1 P u-m,u |h u-m,n | 2 < γ th } .
(5) Thus, the receiver can reliably decode the source data whenever r n ≥ R. For decoding the message correctly, the outage probability must be less than a desired end-to-end outage probability ρ max .

Lemma 1: Considering a set of independent exponential random variables X = {X 1 , . . . , X M } with mean of σ 2 xm , m = 1, . . . , M , the cumulative distribution function (CDF) of the summation of independent-not-identical exponentially distributed random variables, i.e., X = ∑ M m=1 X m is given by

Pr {X < x} = M ∑ m=1 α m ( 1 -e -x σ 2 xm ) (6)
where

α m = M ∏ j=1 j̸ =m σ 2 xm σ 2 xm -σ 2 xj . (7)
Proof: The proof is given in [START_REF] Amari | Closed-form expressions for distribution of sum of exponential random variables[END_REF]. Using the inductive reasoning, the following lemma can be obtained:

Lemma 2: For α m defined in [START_REF] Ikki | Multihop wireless relaying systems in the presence of cochannel interferences: Performance analysis and design optimization[END_REF], the following properties hold:

M ∑ m=1 α m = 1, (8)
M ∑ m=1 α m σ 2k xm = 0, for k = 1, . . . , M -1, (9)
M ∑ m=1 α m σ 2M xm = (-1) M +1 ∏ M m=1 σ 2 xm . ( 10
)
Proposition 1: Given finite sets of independent random variables X = {X 1 , . . . , X M } and Y = {Y 1 , . . . , Y Q } with non-identical exponential distribution and mean of σ 2

xm , m = 1, . . . , M , and σ 2 yq , q = 1, . . . , Q, respectively, the CDF of

SINR = ∑ M m=1 X m 1 +
∑ Q q=1 Y q can be calculated as

P {SINR < γ} = 1 - M ∑ m=1 α m e -γ σ 2 xm Q ∏ q=1 ( σ 2 yq σ 2 xm γ + 1 ) -1 . ( 11
)
Proof: The proof is given in Appendix I. From Proposition 1 and by defining

X m = Pn-m,n|hn-m,n| 2 N0W , m = 1, . . . , M , Y q = Pu-m,u|hu-m,n| 2 N0W
, q = 1, . . . , Q, and Q = |U n |M where |U n | denotes the cardinality of the set U n , the outage probability in (5) can be written as

ρ out n = P {SINR < γ th } = 1 - M ∑ m=1 α m,n e - γ th N 0 W P n-m,n σ 2 n-m,n × ∏ u∈Un M ∏ i=1 ( P u-i,u σ 2 u-i,n P n-m,n σ 2 n-m,n γ th + 1 ) -1 , (12)
where

α m,n = M ∏ j=1 j̸ =m P n-m,n σ 2 n-m,n P n-m,n σ 2 n-m,n -P n-j,n σ 2 n-j,n . (13)
The outage probability ρ n at the nth receiver is affected by all previous n nodes. An upper-bound expression for the outage probability at the destination, i.e., at the (N + 1)th hop can be found as [10, Eq. (29)]

ρ out ≤ 1 - N ∏ ν=0 ( 1 -ρ out N -ν+1 ) ΩM (ν) , (14)
where

Ω M (ν) = 1, 0 ≤ ν < M , and Ω M (ν) = ∑ M i=1 Ω M (ν -i). Note that for the case of M = 2, {Ω 2 (ν)} is a Fibonacci sequence, i.e., Ω 2 (ν) = Ω 2 (ν -1) + Ω 2 (ν -2).
In addition, for the extreme case of M = N + 1, we have Ω M +1 (ν) = 2 ν-1 . In addition, when M = 1, i.e., in the non-cooperative multihop transmission scenario, we have Ω M (ν) = 1, for ν = 0, . . . , M . It is important to note that assuming the equality in (14) implies that the outage at the destination happens even if one intermediate node experience an error. This guarantees that by using the power control strategies proposed in the next section, the outage probability QoS at the destination is satisfied. To get an insight into the relationship between the end-to-end outage probability of ρ des ρ N +1 and ρ 0 , we have

ρ des = 1 - N ∏ ν=0 (1 -ρ 0 ) Ωm(ν) = 1 -(1 -ρ 0 ) ∑ N ν=0 Ωm(ν)
.

Thus, the target outage probability at each hop ρ 0 can be represented in terms of the desired probability of error at the destination ρ des .

Furthermore, assuming ρ out ≪ 1, the outage probability at the destination in (14) can be approximated as follows:

ρ out ≈ N ∑ ν=0 Ω(ν) ρ out N -ν+1 . (16

N0W

, the outage probability at the nth receiving node can be stated as

ρ out n ≈ γ M th ∏ M m=1 SNR n,m ∑ u∈Un M ∑ m=1 α ′ m,u,n M ∑ i=0 INR i n,u,m (M -i)! .
(17) where α ′ m,u,n is defined as

α ′ m,u,n = ∏ i∈Un M ∏ j=1 (i,j)̸ =(u,m) P u-m,u σ 2 u-m,n P u-m,u σ 2 u-m,n -P i-j,i σ 2 i-j,n . (18)
Proof: The proof is given in Appendix II. From Proposition 2, and by using the definition of diversity order

G d = lim SNR→∞ -log(ρ out n ) log(SNR) [12, Eq. (1.19)],
we have the following corollary:

Corollary 1: In a spatial-reused multihop network with cooperation order of M , even in existence of inter-network interference due to spatial reuse, we can achieve the full diversity order of M .

Corollary 2: In interference-free conditions and high SNR regime, the outage probability in (17) can be modified as

ρ out n ≈ γ M th M ! ∏ M m=1 SNR n,m . (19)
From Proposition 1 and by using the facts that e -x ≈ 1x and 1 1+x ≈ 1x, for x ≪ 1, we have Corollary 3: In high SNR regime, for a non-cooperative spatial-reused multihop transmission, i.e., when M = 1, the outage probability can be approximated as

ρ out n ≈ γ th SNR n,1 ( 1 + ∑ u∈Un INR n,u,1 ) . ( 20
)
If the interference due to spatial reuse is strong, the following corollary can be obtained from Proposition 2:

Corollary 4: In high SNR and high interference scenario, i.e., when SNR n,m ≫ 1 and INR n,u,m ≫ 1, the system becomes interference-limited due to an error floor, and the outage probability can be approximated as

ρ out n ≈ 1 - M ∑ m=1 M ∏ j=1 j̸ =m SNR n,m SNR n,m -SNR n,j × ∏ u∈Un M ∏ i=1 ( INR n,u,i SNR n,m γ th + 1 ) -1 . (21)
IV. POWER ALLOCATION FOR MULTIHOP TRANSMISSION WITH SPATIAL-REUSE INTERFERENCE In the previous section, the impact of spatial-reuse interference on the performance of cooperative multihop system is studied. Here, we introduce a power control strategy to improve the system performance. In this section, we derive the required power for the multihop transmission scheme discussed in Section II in order to achieve a certain rate R with a given an outage probability QoS. Finding the optimum value of the transmit powers is complicated due to the complexity of outage probability Pr{r n < R} derived in ( 12) and ( 17). However, by assuming an equal power in every nodes, in the following, a suboptimal power allocation strategy is proposed.

In the case of interference-free transmission, as stated in [10, Theorem 1], the cooperative transmit power coefficients should be equal in each transmission phase. To get a more accurate result, we further assume equal transmission power in all phases to achieve a target outage probability QoS. Thus, assuming the equal transmit power, i.e., P n-m,n = P u-m,u = P 0 , for m = 1, . . . , M , n = 1, . . . , N + 1, and u ∈ U n , we have

ρ out n = 1 -e - γ th N 0 W Pn σ 2 n-1,n ∏ u∈Un ∏ i∈Mu ( σ 2 u-i,n σ 2 n-1,n γ th + 1 ) -1 , (22)
for n = 1, . . . , M -1, where M u = {1} if u < M , and

M u = {1, 2, . . . , M }, if u ≥ M .
For n = M, . . . , N + 1, the outage probability can be rewritten as

ρ out n = 1 - M ∑ m=1 A m,n e - γ th N 0 W Pn σ 2 n-m,n , (23)
where

A m,n = M ∏ j=1 j̸ =m σ 2 n-m,n σ 2 n-m,n -σ 2 n-j,n ∏ u∈Un ∏ i∈Mu ( σ 2 u-i,n σ 2 n-m,n γ th + 1
) -1 .

(24) Since ρ out n is an decreasing function of the power coefficient P 0 for P 0 ≥ 0, to find the minimum value of the problem in P 0 , the constraint ρ out n ≤ ρ n is turned into the equality. Thus, the positive root of ρ out nρ n = 0 should be calculated. Hence, from (22), for n = 1, . . . , M -1, we have

P n = -γ th N 0 W σ -2 n-1,n ln(1 -ρ n ) + ∑ u∈Un ∑ i∈Mu ln ( σ 2 u-i,n σ 2 n-1,n γ th + 1
) .

(25) For n = M, . . . , N +1, and for a given initial value, P n can be calculated from (23) using the following recursive equation:

P (t+1) n = -γ th N 0 W σ -2 n-1,n ln [ 1-ρn A1,n - ∑ i∈Mu-{1} Am,n A1,n e - γ th P (t) n σ 2 n-m,n ] (26) where P
(t) n
is the updated version of the power coefficient in the t-th iteration. Since ρ out n is a decreasing function of P n , to guarantee that ρ out n ≤ ρ n where ρ n is a target outage probability per hop, we have

P 0 = max {P * n } . (27)
where P * n is the solution of recursive equation in (26). Assuming a fixed per-hop outage target of ρ n = ρ 0 and using (15), we can represent P 0 in terms of ene-to-end outage probability of ρ des by replacing ρ n with ρ 0

= 1 -[1 -ρ des ] 1 ∑ N ν=0 Ωm(ν) in (26).
Proposition 3: In the spatial-reused multihop transmission, the minimum allowed target outage requirement at the destination is given by

ρ des ≥ 1-   M -1 ∏ n=1 ∏ u∈Un ∏ i∈Mu ( σ 2 u-i,n σ 2 n-1,n γ th + 1 ) -1   × N +1 ∏ n=M ( M ∑ m=1 A m,n ) ΩM (N -n+1) . ( 28
)
Proof: The minimum amount of permissable target outage per hop can be obtained by putting P n → ∞ in ( 22) and ( 23) to get

ρ n ≥ 1 - ∏ u∈Un ∏ i∈Mu ( σ 2 u-i,n σ 2 n-1,n γ th + 1 ) -1 , (29)
for n = M, . . . , N + 1, and

ρ n ≥ 1 - M ∑ m=1 A m,n , for n = M, . . . , N + 1. (30)
Combining ( 14), (29), and (30), the minimum feasible outage probability QoS at the destination is obtained as (28).

In addition, for a given desired outage probability ρ des at the destination, one can find the minimum spatial-reused factor, i.e., nodes distance K, using Proposition 3. Moreover, it can be observed from ( 8), (28), and (30) that when there is no interference, we have ρ des ≥ 0, and thus, there is no limitation in choosing ρ des .

For the case of non-cooperative multihop transmission, i.e., when M = 1, the closed-form solution for P 0 is obtained in the following proposition: Proposition 4: Assuming the equal power transmission from all nodes, the minimum transmit power P * 0 per node to achieve a per-hop outage probability of ρ n in a noncooperative spatial-reused multihop system over Rayleigh fading channels can be expressed as

P * 0 = max n      γ th N 0 W σ -2 n-1,n ρ n -γ th ∑ u∈Un σ 2 u-1,n σ 2 n-m,n      (31)
Proof: From the approximation given in Corollary 3, which is actually an upper-bound, and by the fact that P n = P 0 , for n = 1, . . . , N + 1, we have

ρ out n ≤ γ th N 0 W P n σ 2 n-m,n + γ th ∑ u∈Un σ 2 u-1,n σ 2 n-m,n . (32)
Then, combining ( 27) and (32), the result in (31) can be yielded. V. NUMERICAL ANALYSIS In this section, numerical results are provided to analyze the performance of the the proposed spatial-reused cooperative multihop scheme. A regular line topology is considered where nodes are located at unit distance from each other on a straight line. The optimal non-cooperative transmission in this network is to send the signal to the next closest node in the direction of the destination. Assume that rate R is 1 2 , bandwidth W is normalized to 1, and the path-loss exponent is assumed to be 3.

In Fig. 3, we compare the outage probability curves of the spatial-reused multihop transmission with respect to the interference-free multihop scenario. The depicted curves are outage probabilities at the last transmission phase, i.e., ρ N +1 , and the non-cooperative (M = 1) and cooperative (M = 2) cases with different spatial reuse factors (K = 3, 4) were compared with interference-free case. As it can be seen, in low and medium SNR regimes, spatial-reused multihop transmission outperforms the interference-free case. For instance, for the cooperative transmission case, when the outage probability of 10 -1 is required at each step, using the spatial-reused scheme with K = 3, around 5 dB saving in transmit power is achievable compared to the interference-free case. However, in high SNR conditions, one can observe that the interferencefree transmission performs better than spatial-reused schemes. Since higher K means lower concurrent transmissions, and thus, lower interferences, as K increases, curves get closer to the interference-free case. Furthermore, Fig. 3 confirms the correctness of our analytical results derived in [START_REF] Jafarkhani | Space-Time Coding Theory and Practice[END_REF], since the curves are exactly match the simulations results.

APPENDIX I PROOF OF PROPOSITION 1

The PDF of Y q is given as p q (y q ) = e . Moreover, the CDF of X = ∑ M m=1 X m is calculated in Lemma 1. By marginalizing over the independent random variables Y q , the CDF of SINR can be calculated as

P {SINR < γ} = ∫ ∞ 0;Q-fold Pr { X < γ + γ Q ∑ q=1 y q } Q ∏ q=1 p q (y q ) dy q = M ∑ m=1 α m ∫ ∞ 0;Q-fold   1 -e - γ(1+ ∑ Q q=1 yq ) σ 2 xm   Q ∏ q=1 p q (y q ) dy q = 1 - M ∑ m=1 α m e -γ σ 2 xm Q ∏ q=1 ∫ ∞ 0 e -yq ( γ σ 2 xm + 1 σ 2 yq ) dy q σ 2 yq , (33)
where in the third equality, we used the first property of Lemma 2 in [START_REF] Maham | Energy-efficient spacetime coded cooperative routing in multihop wireless networks[END_REF]. Thus, the closed-form solution for integral in (33) is obtained as [START_REF] Amari | Closed-form expressions for distribution of sum of exponential random variables[END_REF].

APPENDIX II PROOF OF PROPOSITION 2

We express the CDF of X = ∑ M m=1 X m in Lemma 1 in terms of its Taylor series as

Pr {X < x} = M ∑ m=1 α m ∞ ∑ k=1 -1 k! ( -x σ 2 xm ) k . (34)
In addition, Y = ∑ Q q=1 Y q , where Y q is defined in Lemma 1, has a distribution similar to X with different parameters, and its PDF can be represented as

where Ψ k is defined as

Ψ k = ∫ ∞ 0 ( -
where in the third equality, the binomial series expansion of (1 + y) k is used. Combining (36) in (38) and the closed-form solution for the integral is obtained. Then, by using the second property of Lemma 2, i.e., [START_REF] Maham | Distributed GABBA spacetime codes in amplify-and-forward relay networks[END_REF], the first M -1 terms in (36) becomes zero, and the outage probability is simplified as

Pr {SINR < γ} = ∞ ∑ k=M Ψ k M ∑ m=1 α m σ 2k xm , (39)
Finally, by the fact that σ

Fig. 1 .

 1 Fig. 1. Multihop-then-Cooperate Protocol: (a) Initial phases; (b) Subsequent phases using a distributed space-time code with M = 4.

Fig. 2 .

 2 Fig. 2. Spatial-reused cooperative multihop network with M = 2.

) B. Asymptotic Analysis Proposition 2 :N0W≫ 1 ,

 21 In high SNR conditions, i.e., when SNR n,m Pn-m,nσ 2 n-m,n and medium or low interference scenario due to spatial-reuse where interference terms is defined as INR n,u,m Pu-m,uσ 2 u-m,n

Fig. 3 .

 3 Fig. 3. The outage probability curves versus the transmit SNR in a wireless multihop network with N = 5 for non-cooperative (M = 1) and cooperative (M = 2) cases, different spatial reuse factors, and R = 1 2 bits/sec/Hz.

.0

  By marginalizing over the random variable Y and using (34), the integral in (33) can be rewritten asPr {SINR < γ} = ∫ ∞ Pr {X < γ(1 + y)} p y (y) dy =

1 ) 0 (- 1 )

 101 k+1 γ k (1 + y) k k! p y (y) dy.(37)Now, by replacing p y (y) from (35) in (37), we haveΨ k = ∫ ∞ k+1 γ k (1 + y)

  2 xm ≫ 1 is equivalent to SNR n,m ≫ 1, we can ignore higher order terms, and thus, we have Pr {SINR < γ} ≈ Ψ M

	∑ M m=1 probability can be further simplified as Pr {SINR < γ} ≈ αm σ 2M . Using (10), the outage xm (-1) M +1 Ψ M ∏ M m=1 σ 2

xm

Hence, the result in (

17

) is obtained.


Mehmet Deveci
email: mdeveci@bmi.osu.edu

Kamer Kaya

Bora Uc ¸ar
email: bora.ucar@ens-lyon.fr

Ümit V C ¸atalyürek

A Push-Relabel-Based Maximum Cardinality Bipartite Matching Algorithm on GPUs

Keywords: Push-relabel, GPU, maximum cardinality matchings, bipartite graphs

We design, develop, and evaluate an atomic-and lock-free GPU implementation of the push-relabel algorithm in the context of finding maximum cardinality matchings in bipartite graphs. The problem has applications on computer science, scientific computing, bioinformatics, and other areas. Although the GPU parallelization of the push-relabel technique has been investigated in the context of flow algorithms, to the best of our knowledge, ours is the first study which focuses on the maximum cardinality matching. We compare the proposed algorithms with serial, multicore, and manycore bipartite graph matching implementations from the literature on a large set of real-life problems. Our experiments show that the proposed pushrelabel-based GPU algorithm is faster than the existing parallel and sequential implementations.

I. INTRODUCTION

Bipartite graph matching is a classical problem in graph theory and combinatorial optimization. Given a bipartite graph, we want to find a set of vertex disjoint edges with the maximum cardinality. The problem is of interest in a variety of fields such as bioinformatics [?], scheduling [?, Section 3.8], and image processing [?]. Maximum cardinality bipartite matching is also employed routinely in sparse linear solvers to see if the associated coefficient matrix is reducible; if so, substantial savings in computational requirements can be achieved [?,Chapter 6].

Different approaches for computing maximum matchings in bipartite graphs exist in the literature. Duff et al. [?] discuss the design, analysis and sequential implementation of a class of algorithms which are based on augmenting paths. Push-relabel-based algorithms [?] form the second class. Their design and careful sequential implementation have been recently investigated. A rigorous set of experiments showed that the push-relabel approach is very promising and performs as good as the augmenting-path-based approach [?]. The implementation of the third, more recent, approach based on pseudoflow algorithms [?], have been described by Chandran and Hochbaum [?].

Algorithms dealing with graph problems are usually memory bounded, hence it is relatively hard to obtain a good parallel performance. Moreover, because of the computation irregularities, it is difficult to exploit concurrency. The matching problem is 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 [?] propose a GPU implementation of an algorithm for the maximum weighted matching problem on bipartite graphs. Fagginger Auer and Bisseling [?] study an implementation of a greedy graph matching on GPU. C ¸atalyürek et al. [?] propose different greedy graph matching algorithms for multicore architectures. Azad et al. [?] introduce several multicore implementations of the maximum cardinality matching algorithms on bipartite graphs. In a recent study [?], we investigated the performance of augmentingpath-based approach on GPU architectures.

In this work, we investigate the GPU parallelization of the push-relabel-based maximum cardinality bipartite matching algorithm on GPUs. We develop an atomic-and lock-free implementation to obtain a good performance. Our approach employs the global and gap relabeling heuristics which are also implemented in the GPU. We thoroughly evaluate the performance of the algorithm with a rigorous set of experiments on many bipartite graphs from different applications. The experimental results conclude that the proposed GPUbased implementation is faster than the existing sequential, multicore, and manycore approaches.

The rest of this paper is organized as follows. The background material, some related work, and a summary of contributions are presented in Section ??. Section ?? describes the proposed GPU algorithm and the techniques we develop to make it faster. The comparison of the proposed algorithm with its existing sequential, multicore, and manycore counterparts is given in Section ??. Section ?? concludes the paper.

II. BACKGROUND AND RELATED WORK

In a bipartite graph G = (V 1 ∪ V 2 , E) with two disjoint vertex sets V 1 and V 2 , each edge in E has one endpoint in V 1 and one in V 2 . Let m and n be the number of vertices in V 1 and V 2 respectively. And let Γ(v) = {u : {u, v} ∈ E} represent the neighborhood of a vertex v ∈ V 1 ∪V 2 . To better distinguish the vertices in V 1 and V 2 in the text, we will employ the matrix notation and use V R and V C for V 1 and V 2 , respectively, and refer the vertices in V R as the set of rows and the vertices in V C as the set of columns.

A subset of the edges, M, is called a

matching if no two edges in M share an endpoint. If v ∈ V is an endpoint in M we say v is matched. Otherwise, it is unmatched. A matching M is called maximal, if there exists no other matching M such that M ⊃ M. A maximal matching M is called maximum if |M| ≥ |M | for all possible matchings M where |M| is the cardinality of M. If |M| = |V R | = |V C |, M is called perfect matching.
The deficiency of a matching M is the difference between the cardinality of a maximum matching and |M|.

A. Maximum cardinality bipartite matching algorithms

A path in G is M-alternating if its edges alternate between those in matching M and those not in M. An M-alternating path P is called M-augmenting if the start and end vertices of P are both unmatched. Augmenting path-based algorithms are based on the following theorem.

Theorem 1 ([?]): Let G be a graph (bipartite or not) and let M be a matching in G. Then M is of maximum cardinality if and only if there is no M-augmenting path in G.

Given a possibly empty matching M, augmenting-pathbased algorithms searches for an M-augmenting path P. If none exists then the matching M is maximum by Theorem ??. Otherwise, P is used to increase the cardinality of M by setting M = M⊕E(P) where E(P) is the edge set of a path P, and M ⊕E(P) = (M∪E(P))\(M∩E(P)) is the symmetric difference. Since both the first and the last edge of P were unmatched in M, we have |M ⊕ E(P)| = |M| + 1. The way in which augmenting paths are searched constitutes the main difference between the algorithms based on augmenting paths, both in theory and in practice. These algorithms mainly use graph traversal techniques such as the depth-first search (DFS) and breadth-first search (BFS). We refer the reader to [?] for a recent survey on these algorithms.

B. The push-relabel algorithm for bipartite matching

Push-relabel algorithms search and augment simultaneously. They do not explicitly construct augmenting paths but they follow them with a clever heuristic. That is, they repeatedly augment a prefix of a hypothetical augmenting path

P 2 = (v, u, w) in G where w ∈ V C is matched to u ∈ V R , and v ∈ V C is an unmatched column.
Augmentations are performed by unmatching w and matching v to u. If w has an unmatched neighbor, the suffix of an augmenting path has been found, allowing the augmentation of |M|. These hypothetical augmentation operations are performed until no further suffixes can be found. The algorithms are guided by assigning a label to every vertex which provides a lower bound on the distance to the nearest unmatched row.

Goldberg and Tarjan's original push-relabel algorithm was designed for the maximum flow problem [?]. The bipartite matching problem is a special case of maximum flow. Hence, the algorithm can be simplified to be adapted to the bipartite matching problem. In fact, it is known to be one of the fastest algorithms for bipartite matching [?]. In the following, we portray this algorithm in a ready-to-implement pseudocode form as shown in Algorithm ??. This algorithm will be referred to as PR throughout the paper. Let ψ(v) denote a lower bound on the length of the shortest alternating-path from the vertex v to an unmatched row. Hence, ψ can be used to find the direction of the closest unmatched row for each vertex. If v is an unmatched column, this path is also an augmenting path. During the initialization of PR, the algorithm sets ψ(v) = 1 for all v ∈ V C and ψ(v) = 0 for all v ∈ V R . Throughout the algorithm, the unmatched columns are called active.

As long as an active column v exists, the algorithm performs the push operation. Γ(v) is searched for a row u ∈ Γ(v) with the minimum ψ(u). Note that ψ(v) -1 is the infimum for the value of ψ(u). This property, ψ(u) ≥ ψ(v) -1 for u with the minimum ψ in the neighborhood of v, holds after the initialization (ψ(u) = 0 and ψ(v) = 1 for all u ∈ V R and v ∈ V C ) and is maintained throughout the algorithm as the neighborhood invariant. Hence, when an edge {v, u} with ψ(v) = ψ(u) + 1 is found the search stops.

Let u be the row with minimum ψ in Γ(v). If u is unmatched, it can be matched to v, and the cardinality of M can be increased by one (single push). Note that an unmatched row vertex u has ψ(u) = 0, and it will always have the minimum ψ value. If u is matched to a column w, a double push operation is performed which removes {w, u} from M, adds {v, u} to M, and makes w active. Hence, once a row is matched, it never becomes unmatched again.

After the push operation, PR relabels v with ψ(v) ← ψ(u)+ 1. By definition, ψ(u) is a lower bound on the alternatingpath distance from u to a closest unmatched row. Since the path between v and its closest unmatched row must contain some u ∈ Γ(v), and ψ(u) has the minimum ψ among all the neighbors of v, ψ(u) + 1 is a lower bound on the length of a path between v and its closest unmatched row.

After relabeling v, PR relabels u by ψ(u) ← ψ(u) + 2. For a single push, ψ(u) becomes 2. For a double push, any alternating path from u to a closest unmatched row now contains v. Such a path always uses the edge (u, v), and hence, ψ(u) must be at least ψ(v) + 1 which is actually ψ(u) + 2 due to relabeling of v. Clearly, this new value maintains the neighborhood invariant ψ(u) ≥ ψ(v) -1.

The maximum length of any augmenting path in G is at most min(2m, 2n) -1. Since ψ is a lower bound on the length of a path to an unmatched row, when the minimum ψ(u) among all u ∈ Γ(v) becomes m + n, the vertex v cannot be on an augmenting path. In this case, no push operations are realized and v is considered unmatchable, i.e., inactive. The push and relabel operations are repeated until there is no active column left, either because all vertices have been matched or marked as inactive. Using Theorem ??, it is easy to show that in this case M is a maximum matching. The time complexity of the algorithm is O((n + m)τ ) [?].

A push-relabel implementation does not need to store a ψ(u) for a row u since the value is 0 when u is unmatched and ψ(w) + 1 when (u, v) ∈ M. However, such an approach increases the amount of jumps and arithmetic operations. Thus, we store a full ψ array for both the rows and the columns.

C. Global relabeling

As various studies showed, PR's performance can be improved by periodically setting all labels to exact distances. This heuristic is called global relabeling and realized via a BFS starting from all unmatched rows as shown in Algorithm ??. During the traversal, the label of each vertex v is set to the minimum distance from v to any unmatched row. And each vertex w not visited by the BFS is assigned a label ψ(w) = m+n, hence it is removed from further consideration.

Algorithm 2 GR: Global Relabeling

Input: A bipartite graph G = (VR ∪ VC , E) and a matching M in G Output: An accurate distance labeling ψ w.r.t. M 1: Q ← u for all unmatched u ∈ VR 2: Set ψ(v) = m + n for all v ∈ VC 3: Set ψ(u) = m + n for all matched u ∈ VR 4: while Q not empty do 5: u ← POP u from Q 6: for all v ∈ Γ(u) do 7: if ψ(v) = m + n then 8: ψ(v) ← ψ(u) + 1 9: if {v, w} ∈ M then 10: ψ(w) ← ψ(v) + 1 11: PUSH w to Q 12: return ψ
In practice, a counter is incremented every time the value of ψ(v) is changed in Line ?? of Algorithm ?? in order to keep track of the number of pushes. The usual approach is to call GR when this counter reaches a predetermined threshold. A threshold of n was suggested as the standard frequency of global relabels [?]. The performance of the push-relabel approach and its several sequential implementations for the bipartite cardinality matching problem have been well studied and various improvements have been proposed [?], [?], [?]. For example, it has been shown that when the active columns are visited with a FIFO order, the push-relabel approach obtains a better performance. The most recent implementation and a comparison with the augmenting-path-based algorithms are given by Kaya et al. [?]. The authors showed that PR is slightly better than the augmenting-path-based algorithms. In our experiments, we use their implementation for comparison purposes.

D. Parallel algorithms for bipartite cardinality matchings

Efficient multicore counterparts of a number of augmentingpath based algorithms have been proposed in a recent study [?]. The parallelization of the multicore algorithms is achieved by using atomic operations at BFS and/or DFS steps of the algorithm. Among these algorithms, P-DBFS, which employs vertex disjoint BFSs to find the augmenting paths, obtained the best performance in our experiments. Its relative performance was also good in the original experiments [?]. Thus, in this paper, we compare the performance of our GPUbased push relabel algorithm with P-DBFS.

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. In recent work, we designed and developed atomic-free GPU implementations, G-HK and G-HKDW [?], corresponding to two augmenting-path-based algorithms, HK [?] and HKDW [?]. HK has the best known worst-case running time complexity of O(τ √ n + m) for a bipartite graph with τ 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 DFS-based 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. In this paper, we compare the performance of the proposed algorithm with that of G-HKDW which performed better than G-HK [?].

GPU-based push-relabel implementations for the maximum flow problem exist in the literature [?], [?], [?]. Vineet and Narayanan used the push-relabel approach to find graph cuts in GPUs [?]. Their implementation does not have the global relabeling heuristic which is later added by Huseein et al. [?]. However, the implementation does not have a relabeling operation except the global one. He and Hong proposed a GPU-based implementation of the push-relabel approach with global relabeling on hybrid CPU-GPU architectures for the maximum flow problem [?]. To integrate both push and relabel operations into a CUDA kernel, they used the atomic fetchand-add function supported by CUDA-based GPUs.

To the best of our knowledge, ours is the first GPU implementation of the push-relabel approach for the maximum bipartite cardinality matching problem. The algorithm does not need locks or atomic instructions to perform push-relabel operations concurrently. Instead, the algorithm allows matching inconsistencies throughout its execution, which are resolved at the end. Furthermore, it is implemented fully on the GPU with CUDA including a global relabeling phase.

III. PUSH-RELABEL-BASED BIPARTITE CARDINALITY MATCHING ON GPUS

A high-level structure of our GPU-based algorithm is given in Algorithm ??. The main while loop at line ?? iterates until all active columns are consumed by the GPU kernel executions at line ?? where in a single execution of G-PR-KRNL, an active column is processed by a single GPU thread. In our implementation, we use an array ψ(•) to store the labels of vertices in V R ∪ V C and another array µ(•) to store their matching vertices. During the course of the algorithm, if two

vertices u ∈ V R and v ∈ V C are matched µ(u) = v and µ(v) = u. If a row u ∈ V R is unmatched, µ(u) = -1. If a column v ∈ V C is inactive, either µ(v) = -2, i.e.,
v cannot reach to an unmatched row with an augmenting path, or µ(v) = u and µ(u) = v, i.e., v is matched with a u ∈ V R . The algorithm maintains these as its matching invariants throughout the execution. Note that a column vertex v can have µ(v) > -1 even though it is unmatched.

Algorithm 3 G-PR: GPU-based Push-Relabel Matching Input: G = (VR ∪ VC , E), µ 1: loop ← 0 2: actExists ← true 3: iterGR ← 0 4: while actExists do 5: if loop = iterGR then 6: maxLevel, ψ ← G-GR(G, µ) 7:
iterGR ← GETITERGR(maxLevel, loop, iterGR)

8:
actExists ← false 9:

µ, ψ, actExists ← G-PR-KRNL(G, µ, ψ)

10:

loop ← loop + 1 11: µ ← FIXMATCHING(µ) 12: return µ

A. Global relabel operation on GPUs

The proposed algorithm G-PR periodically uses a globalrelabeling function to set exact ψ values. The pseudocode of this function is given in Algorithm ??. The function first calls a GPU kernel INITRELABEL. For each unmatched row u, the initialization kernel sets ψ(u) ← 0. For all other vertices, ψ values are set to m + n. After the initialization, a global relabeling is performed by multiple G-GR-KRNL calls where each call processes a level of a BFS starting from all unmatched rows. The details of the kernel is given in Algorithm ??. The row vertices are distributed to the threads and if a row u has ψ(u) = cLevel all of the unvisited columns v ∈ Γ(u) are visited. For each such column v, ψ(v) is set to cLevel + 1 and its matched row's ψ value is set to cLevel + 2. Although the same ψ(•) value can be set by multiple threads, since all threads use the same value for the same location there is no concurrency issues for G-GR-KRNL. ψ, uAdded ← G-GR-KRNL (G, ψ, µ, cLevel)

7: cLevel ← cLevel + 2 8: maxLevel ← cLevel Algorithm 5 G-GR-KRNL Input: G = (VR ∪ VC , E), ψ, µ, cLevel 1: for all row vertex u ∈ VR in parallel do 2: if ψ(u) = cLevel then 3: for all v ∈ Γ(u) do 4: if ψ(v) = m + n then 5: ψ(v) ← cLevel + 1 6: if µ(v) > -1 and µ(µ(v)) = v then 7: ψ(µ(v)) ← cLevel + 2 8: uAdded ← true
Following the practice, in our experiments, we used an initial matching M obtained by a standard greedy matching heuristic [?]. When the initial matching is empty, we start with ψ(u) = 0 for each u ∈ V R , and

ψ(v) = 1 for each v ∈ V C .
Even with a nonempty matching, the algorithm runs correctly after this initialization. However, our preliminary experiments show that applying a global relabeling at the beginning of the main while loop of G-PR leads significant performance improvements. Hence, we set iterGR to 0 (line ?? of Algorithm ??) which stores the iteration number in which the next global relabeling will be performed. Experiments in the literature show that the performance of the algorithm changes drastically with the strategy to set iterGR [?], [?]. In our preliminary experiments, we also observed this sensitivity with different implementations of GETITERGR (line ?? of Algorithm ??). In practice, a global relabeling is performed after every k × (n + m) pushes where the suggested k values are between 1 and 2 [?], [?]. Unlike the sequential pushrelabel algorithms, it is very expensive to count the number of pushes performed during kernel executions. Hence, for the GPU-based implementation, we need a different strategy. We initially experimented with fixed intervals and performed a global relabeling after every k iterations of the main while loop. We then developed an adaptive strategy where the next global relabeling is done after k × maxLevel loop iterations and experimented with several k values. The design rationale behind this strategy is based on the following theorem:

Theorem 2 ([?]): If M is a matching with a deficiency d there exists a set of d vertex disjoint augmenting paths.

Let S be the set of paths in Theorem ??. Since the paths are vertex disjoint, their total length can be used as a lower bound for m+n. Since the BFS in a global relabeling uses alternating paths and these paths stop at unmatched columns, a fraction of the maxLevel is a reasonable estimate on the average path length in S. Thus, following the sequential practice, using k × maxLevel kernel executions with d active columns can be promising. Our experiments show that the adaptive strategy is superior to global relabeling with fixed intervals.

B. Push relabel operations on GPUs

After a global relabeling in Algorithm ??, the kernel G-PR-KRNL processes the active columns and performs concurrent push-relabel operations. The first implementation of this kernel is given in Algorithm ??. As described in Section ??, there are three cases while processing an active column v ∈ V C : Let u ∈ V R be the row in Γ(u) with a minimum ψ value. If ψ(u) < m + n a single push matches v with u and consumes it (lines ?? and ??). If µ(u) was previously matched with a v ∈ V C a double push produces a new active column v . Note that we do not actually implement the double-push operation and fix µ(v ) since it has never been solely used, i.e., we always check µ(µ(v )) while we want to know if the column is matched or unmatched. If ψ(u) = m + n then v is set to inactive and consumed (line ??).

Concurrent push-relabel operations in a single G-PR-KRNL execution do not create a problem since when two active columns v and v perform pushes simultaneously with a row u, we have µ(v) = µ(v ) = u and µ(u) ∈ {v, v }. Hence, although one of the µ(v) and µ(v ) is wrong, the inconsistency in the matching can be captured by the algorithm. A problem might arise if there is an inconsistency in the neighborhood and matching invariants of the algorithm. We claim that the invariants are not violated neither. If both v and v select u at the same time for the push operation, they will both store the same ψ(u) as ψ min . If v completes the push operation earlier than v , first ψ(v) and ψ(u) will be updated for v. Although the matching of v will be invalidated by v , the update of ψ(v) does not violate the neighborhood invariant as ψ min = ψ(u) ≥ ψ(v) -1. Therefore, although v has another neighbor u with ψ(u ) = ψ min , both invariants will still be valid. The double update of ψ(u) does not create a problem either, as both of the threads set it to the same value. A ψ(u) entry can be set more than once with different values within a single execution of G-PR-KRNL. If this is the case, these push operations can be considered as two consecutive pushes in a sequential execution as they have different ψ min values stored at the beginning of push operation.

The while loop of the main algorithm G-PR iterates and continues to call the kernel G-PR-KRNL until all the active columns are consumed. Once the maximum cardinality matching is found, the algorithm terminates and the inconsistencies in µ is fixed FIXMATCHING kernel which implements: µ(v) ← -1 for any v with µ(µ(v)) = v. These inconsistencies do not prevent the algorithm obtaining a maximum cardinality matching. When all GPU kernels are done, such inconsistencies can only exist at the column entries of µ and the row matching will be correct. In fact, if there exist a perfect matching in G all the entries in µ are correct and no fix operation is required.

Algorithm 6 G-PR-KRNL

Input: G = (VR ∪ VC , E), µ, ψ 1: actExists ← false 2: for all column vertices v ∈ VC in parallel do

3: if µ(u) = -1 or µ(µ(v)) = v then 4: actExists ← true 5: ψmin ← m + n 6: u ← -1 7:
for all u ∈ Γ(v) do if ψmin < m + n then 14:

µ(u) ← v 15: µ(v) ← u 16: ψ(v) ← ψmin + 1 17: ψ(u) ← ψmin + 2 18: else 19: µ(v) ← -2

C. Reducing the number of GPU threads

In the G-PR-KRNL kernel, the proposed parallel pushrelabel-based algorithm uses n GPU threads to traverse the active column vertices. However, thanks to the initial matching heuristic, the number of unmatched columns is much smaller than n. Furthermore, we observed that most of the G-PR-KRNL executions are performed only with a few number of active columns, since it is harder to find augmenting paths when they are longer. This is usually the case when the deficiency gets smaller, which happens towards the end of the execution. In order to reduce the number of threads, we propose two main modifications on the algorithm.

1) Keeping the list of active columns: The first modification is to keep a list of active columns and have the threads work on those active columns (instead of all columns). Assume that we have an array A. Let v = A(i) be an active column and v ∈ V C be the column which becomes active once v is consumed. Since there can be at most one such column, in a sequential execution, setting A(i) ← v is sufficient to keep the id of the new active columns produced during the pushrelabel operations. When no new active column is produced after a consumption, one can set the corresponding value to -1. However, in a concurrent execution, some of the pushrelabel operations may need to be rolled back due to the conflicts.

As explained above, we need to keep the id of each consumed column to roll it back in case of conflicts. Hence, a single array is not sufficient to support correct and concurrent execution of push operations. In our implementation, we used two arrays A c and A p to maintain such a list and support concurrent operations on it at the same time. Initially, both arrays contain the unmatched column indices. In the course of the algorithm, the entries will contain either the ids of the active columns or -1 which appear when an augmentation (only a single push) takes place without producing new active columns.

Algorithm 7 G-PR: GPU-based Push-Relabel Matching

Input: G = (VR ∪ VC , E), µ 1: loop ← 0 2: actExists ← true 3: shrink ← false 4: iterGR ← 0 5: while actExists do 6:

if loop = iterGR then 7: maxLevel, ψ ← G-GR(G, µ) 8:
iterGR ← GETITERGR(maxLevel, loop, iterGR)

loop ← loop + 1 20: µ ← FIXMATCHING(µ) 21: return µ
In a single loop of the new G-PR, given in Algorithm ??, push-relabel operations are performed in two steps. In the first step, the algorithm fixes A c by using A p from the previous iteration with the kernel G-PR-INITKRNL given in Algorithm ??. The fix operation is done in two stages: first the conflicting operations in the previous execution are detected (line ??) and then they are rolled back (line ??). In addition to A c and A p , we used an additional array iA of size n where iA(v) = loop if and only if v is an active column. While fixing A c , G-PR-INITKRNL also sets the corresponding entries to loop (line ??). This array will be used to avoid duplicate processing of an active column by two different threads while performing push-relabel operations. Note that the number of GPU threads for this kernel execution is |A p | which is expected to be much smaller than n after the greedy matching heuristic.

Algorithm 8 G-PR-INITKRNL

Input: G = (VC ∪ VR, E), µ, ψ, loop, Ac, Ap, iA 1: actExists ← false 2: for all i ∈ {1, . . . , |Ap|} in parallel do

3: v ← Ap(i) 4: if v = -1 then 5: if µ(v) = -1 or v = µ(µ(v)) then 6: Ac(i) ← Ap(i) 7: else 8: v ← Ac(i) 9: if v = -1 then 10: iA(v) ← loop 11:

actExists ← true

Once A c is fixed, the push-relabel operations are performed by G-PR-PUSHKRNL given in Algorithm ??. The kernel works along the same lines with the previous GPU-based concurrent push-relabel kernel G-PR-KRNL. However, it uses less threads. The main difference is keeping the ids of the active columns produced throughout the execution (line 18)

and the extra check iA(µ(u))

? = loop at line 13. This check is necessary due to concurrency; an already active column v which is being processed can be inserted to A p in the same kernel execution. Although, this does not create a problem for the current execution, if G-PR-INITKRNL roll backs the operation on v, in the next execution v will be processed by two different threads. To avoid this problem, a column, which is already active at the beginning of the iteration is forbidden in A p . After completing G-PR-PUSHKRNL, at the end of each iteration, we swap the current and previous active column arrays for the next iteration (line 18 of Algorithm ??).

Algorithm 9 G-PR-PUSHKRNL

Input: G = (VR ∪ VC , E), µ, ψ, loop, Ac, Ap, iA 1: for all i ∈ {1, . . . , |Ap|} in parallel do 2: v ← Ac(i) 3: if v = -1 then 4: ψmin ← m + n 5: u ← -1 6: for all u ∈ Γ(v) do 7: if ψ(u ) < ψmin then 8: ψmin ← ψ(u ) 9: u ← u 10: if ψmin = ψ(v) -1 then 11: break 12:
if ψmin < m + n then Ap(i) ← -1

2) Dynamic compression of active-column lists: The size of the active column arrays used in Algorithms ?? and ?? is equal to the number of unmatched columns at the beginning. When the deficiency is smaller, we observed that the main loop of the algorithm iterates more in order to increase the matching size. Hence, the impact of a reduction on the size of the active column arrays will be significant. Thus, as a second improvement, we periodically shrink these arrays after each global-relabel operation by using G-PR-SHRKRNL (line 12). The kernel works along the same lines with G-PR-INITKRNL. But instead of directly writing the new active columns to A c , an initial pass on the previous and the current active columns is performed in order to obtain the number of local active columns within each thread. Then the threads perform a prefix sum operation on these counts. After the prefix sum, each thread has a private write region on A c . They traverse the previous active column array again, and instead of counting, this time they write the ids to their private regions. Hence, after G-PR-SHRKRNL the size of A c and A p is reduced to the exact number of active columns. We observed that although the shrinking can yield significant improvements, after some number of active columns the improvement does not compensate its overhead. Hence, we decide to use the shrink operation only when |A c | ≥ 512 (line 11 of Algorithm ??).

IV. EXPERIMENTS

The running time of the proposed implementations are compared against the sequential PR implementation (PR) [?], the multicore parallel implementations P-DBFS [?], and a GPU implementation of HKDW [?]. For the sequential PR implementation, we tried a set of k values to set the frequency of global-relabeling to k × (m + n). For our data set k = 0.5 was slightly better. Hence, we used it in our experiments. The CPU implementations are tested on a computer with 2.27GHz dual quad-core Intel Xeon CPUs with 2-way hyper-threading and 48GB main memory. The algorithms are implemented in C++ and OpenMP. The GPU implementations are tested on NVIDIA Tesla C2050 with usable 2.6GB of global memory. Tesla C2050 is equipped with 14 multiprocessors each containing 32 CUDA cores, totaling 448 CUDA cores. The implementations are compiled with gcc-4.4.4, cuda-4.2.9 and -O2 optimization flag. For the multicore algorithms, 8 threads are used. A standard heuristic called the cheap matching (see [?],

[?]) is used to initialize all tested algorithms. We compare the running time of the matching algorithms after this common initialization.

The algorithms are run on bipartite graphs corresponding to 28 matrices from variety of classes at UFL matrix collection [?]. The graphs are chosen from a large dataset given in [?]. The graphs for which all of the sequential algorithms, Pothen-Fan-Plus [?], Hopcroft-Karp [?], and the push-relabel algorithm (PR [?]), take less than one second are filtered out.

We first compare the performance of the proposed GPU algorithm with different parameters. Figure ?? shows the geometric means of the runtimes on different sets. In the figure, G-PR-First corresponds to our first implementation explained in Algorithm ??. G-PR-Shr corresponds to the implementation presented in Algorithm ??, while G-PR-NoShr is the one that omits the shrinking mechanism (G-PR-SHRKRNL functions omitted). Different strategies are used to specify the frequency of the global relabeling. Here, (adaptive, k) refers to the the adaptive global-relabeling strategy explained in Section ??. That is, the next global relabeling is performed after k × maxLevel push-relabel kernel calls. For the (fix, k) strategy, the next global-relabeling operation is performed after k push-relabel kernel executions. As the figure shows, we obtain the best results with (adaptive, 0.3) and (adaptive, 0.7). Having a fix frequency for all graphs throughout the execution is outperformed by the adaptive heuristic for almost all configurations. The results verify the necessity of adjusting the frequency of the global relabeling with respect to the graph structure and the current matching. The figure also shows that the proposed G-PR-active algorithm improves the performance of each configuration by 14% to 84%, as it decreased the divergence of the GPU threads. Shrinking the list of active columns throughout the execution improves the performance of the algorithms by another 2-8%. Since G-PR-shrink obtains the best performance with (adaptive, 0.7) configuration, we only compare the performance of this configuration with other implementations in the literature.

Figure ?? shows the speedup profiles of G-PR, G-HKDW, and P-DBFS on the experiment set. The speedups are calculated with respect to the sequential PR algorithm. A point (x,y) in the plots corresponds to the probability y of obtaining at least x speedup. As the plot shows, G-PR has a better speedup profile: with 39% probability it obtains a speedup at least 5. This probability is 21% and 14% for G-HKDW and P-DBFS, respectively. Furthermore, the proposed algorithm is faster than the sequential PR algorithm for 82% of the graphs. Although the performance of G-HKDW is close to that of G-PR, both G-HKDW and P-DBFS are outperformed by the proposed G-PR algorithm. On the average, the proposed algorithm is 1.30 and 2.82 times faster than G-HKDW and P-DBFS, respectively, when geometric means of the runtimes are compared.

Figure ?? shows the performance profiles of multicore and manycore 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 plot clearly shows the separation among G-PR and the other algorithms, thus marking G-PR as the fastest in most cases. For 75% of the cases, G-PR is at most 1.5 times worse than the best algorithm. These percentages are 46% and 14% for G-HKDW and P-DBFS, respectively. In particular, G-PR obtains the best performance in 61% of the selected graphs.

Figure ?? shows the indivual speedups of G-PR on the each graph. The x axis gives the graph ids given in Table ?? and the graphs are ordered with the increasing number of rows. The proposed algorithm obtains a speedup on 23 (out of 28) of the graphs. The maximum speedup achieved is on delaunay n24 as 12.60, while the minimum speedup is obtained as 0.31 on hugetrace-00000 graph. The overall average speedup is 3.05. Table ?? gives the actual running times of the best GPU and multicore algorithms, together with the sequential PR algorithm. As seen from this table, on 15 of the graphs, the GPU algorithm is faster than the other algorithms.

V. CONCLUDING REMARKS

We proposed a parallel push-relabel-based maximum cardinality matching algorithm for bipartite graphs on GPUs. We presented a lock-and atomic-free implementation and propose a set of modifications to obtain a better performance. We investigated the effect of the global-relabeling frequency and developed a strategy to amortize the cost of global and push relabeling. We presented experiments on various graphs and compared the performance of the proposed implementation against sequential, multicore, and manycore algorithms from the literature. The experiments showed that the proposed GPU implementation is faster than the existing implementations. We obtained speedups varied from 0.31 to 12.60, averaging 3.05, on a set of 28 graphs with respect to a recent sequential pushrelabel-based implementation.

As a future work, we intend to investigate more architectural and algorithmic advancements to improve the performance. One idea we want to study is the concurrent execution of global-relabeling and push-relabel kernels [?]. As far as we know, there is no GPU implementations of these techniques. After the Fermi architecture, NVIDIA GPUs support concurrent kernel executions with streams. Hence, it may be promising to occupy the device with two kernels when there is not enough parallelism with a single kernel.

ACKNOWLEDGMENT

This work was supported in parts by the NSF grants OCI-0904809 and OCI-0904802.

  Algorithm 4 G-GR: GPU-based Global RelabelingInput: G = (VR ∪ VC , E), µ 1: ψ ← INITRELABEL(µ) 2: uAdded ← true 3: cLevel ← 0 4: while uAdded do

Fig. 1 .

 1 Fig. 1. The comparison of the proposed G-PR algorithms with different global-relabeling strategies.

Fig. 2 .Fig. 3 .

 23 Fig.2. Speedup profiles of the parallel algorithms. A point (x,y) in the plots corresponds to the probability of obtaining at least x speedup is y.

Fig. 4 .

 4 Fig. 4. Individual speedups of G-PR algorithm on each graph instance.

  Algorithm 1 PR: Push-Relabel based Bipartite Matching

	Input: A bipartite graph G = (VR ∪ VC , E) and a (possibly empty)
		matching M
	Output: A maximum cardinality matching M
	1: Set ψ(u) = 0 for all u ∈ VR
	2: Set ψ(v) = 1 for all v ∈ VC
	3: Set all v ∈ VC unmatched by M to active
	4: while an active column v exists do
	5:	Find a row u ∈ Γ(v) of minimum ψ(u)
	6:	if ψ(u) < m + n then
	7:	if {u, w} ∈ M then
	8:	M ← M \ {u, w}/*Single push*/
	9:	Set w active
	10:	M ← M ∪ {u, v} /*Double push*/
	11:	ψ(v) ← ψ(u) + 1 /*Relabels v*/
	12:	ψ(u) ← ψ(u) + 2 /*Relabels u*/
	13:	Set v inactive
	14: return M = M

TABLE I ACTUAL

 I RUNNING TIME OF EACH ALGORITHM FOR ALL GRAPHS TOGETHER WITH THE NUMBER OF ROWS, COLUMNS, AND EDGES IN EACH GRAPH. IM AND MM GIVES THE CARDINALITY OF THE INITIAL AND MAXIMUM MATCHING FOR EACH GRAPH, RESPECTIVELY. THE GEOMETRIC MEANS OF THE RUN TIMES OF G-PR, G-HKDW, P-DBFS, AND PR ARE GIVEN IN THE BOTTOM ROW.

	ID Graph	#Row	#Cols	#Edges	IM	MM	G-PR	G-HKDW	P-DBFS	PR
	1	amazon0505	410,236	410,236	3,356,824	332,972	395,397	0.09	0.18	22.70	0.52
	2	coPapersDBLP	540,486	540,486 15,245,729	510,992	540,226	0.62	0.42	6.27	0.59
	3	amazon-2008	735,323	735,323	5,158,388	587,877	641,379	0.12	0.11	0.18	0.93
	4	flickr	820,878	820,878	9,837,214	285,241	367,147	0.13	0.22	0.35	0.99
	5	eu-2005	862,664	862,664 19,235,140	642,027	652,328	0.40	1.54	0.94	0.80
	6	delaunay n20	1,048,576	1,048,576	3,145,686	993,174	1,048,576	0.06	0.04	0.09	0.32
	7	kron g500-logn20	1,048,576	1,048,576 44,620,272	431,854	513,334	0.38	0.60	8.19	1.24
	8	roadNet-PA	1,090,920	1,090,920	1,541,898	916,444	1,059,398	0.33	0.14	0.29	0.59
	9	in-2004	1,382,908	1,382,908 16,917,053	781,063	804,245	0.58	1.44	2.16	0.56
	10	roadNet-TX	1,393,383	1,393,383	1,921,660	1,158,420	1,342,440	0.45	0.14	0.33	0.69
	11	Hamrle3	1,447,360	1,447,360	5,514,242	1,211,049	1,447,360	0.94	1.36	2.70	0.56
	12	as-Skitter	1,696,415	1,696,415 11,095,298	891,280	1,035,521	0.34	0.49	1.89	1.13
	13	GL7d19	1,911,130	1,955,309 37,322,725	1,904,144	1,911,130	0.24	0.58	0.38	1.38
	14	roadNet-CA	1,971,281	1,971,281	2,766,607	1,668,268	1,913,589	0.68	0.34	0.53	1.55
	15	delaunay n21	2,097,152	2,097,152	6,291,408	1,987,326	2,097,152	0.18	0.13	0.21	1.06
	16	kron g500-logn21	2,097,152	2,097,152 91,042,010	812,883	964,679	0.68	0.99	1.50	2.77
	17	wikipedia-20070206	3,566,907	3,566,907 45,030,389	1,623,931	1,992,408	0.62	1.09	5.24	3.11
	18	patents	3,774,768	3,774,768 14,970,767	1,892,820	2,011,083	0.54	0.88	0.84	3.65
	19	com-livejournal	3,997,962	3,997,962 34,681,189	2,577,642	3,608,272	2.08	4.58	22.46	9.67
	20	hugetrace-00000	4,588,484	4,588,484	6,879,133	4,581,148	4,588,484	2.71	1.96	0.83	0.84
	21	soc-LiveJournal1	4,847,571	4,847,571 68,993,773	2,831,783	3,835,002	1.35	3.32	14.35 12.66
	22	ljournal-2008	5,363,260	5,363,260 79,023,142	3,941,073	4,355,699	1.54	2.37	10.30 10.01
	23	italy osm	6,686,493	6,686,493	7,013,978	6,438,492	6,644,390	5.46	5.86	1.20	6.84
	24	delaunay n23	8,388,608	8,388,608 25,165,784	7,950,070	8,388,608	0.81	0.96	1.26	8.86
	25	wb-edu	9,845,725	9,845,725 57,156,537	4,810,825	5,000,334	2.00	33.82	8.61	3.94
	26	hugetrace-00020	16,002,413 16,002,413 23,998,813	15,535,760 16,002,413	14.19	7.90	393.13 28.69
	27	delaunay n24	16,777,216 16,777,216 50,331,601	15,892,194 16,777,216	1.83	1.98	2.41 23.01
	28	hugebubbles-00000	18,318,143 18,318,143 27,470,081	18,303,614 18,318,143	13.65	13.16	3.55	13.51
							GEOMEAN	0.70	0.92	1.99	2.15