Mehmet Deveci
email: mdeveci@bmi.osu.edu

Kamer Kaya

Bora Uçar
email: bora.ucar@ens-lyon.fr

Ümit V C ¸atalyürek

GPU accelerated maximum cardinality matching algorithms for bipartite graphs

Keywords: GPU, maximum cardinality matchings, bipartite graphs, breadth-first search

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 the 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 [START_REF] Kim | 3-D object recognition using bipartite matching embedded in discrete relaxation[END_REF], chemical structure analysis [START_REF] John | Kekulé patterns and Clar patterns in bipartite plane graphs[END_REF], and bioinformatics [START_REF] Azad | Multithreaded algorithms for matching in graphs with application to data analysis in flow cytometry[END_REF] (see also another two discussed by Burkard et al. [START_REF] Burkard | Assignment Problems[END_REF]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 [START_REF] Duff | Direct Methods for Sparse Matrices[END_REF]Chapter 6].

Achieving good parallel performance on graph algorithms is challenging: they are memory bounded; there are poor localities of the memory accesses; and the dependencies among the computations are irregular. 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 [START_REF] Vasconcelos | Bipartite graph matching computation on GPU[END_REF] propose a GPU implementation of an algorithm for the maximum weighted matching problem on bipartite graphs. Fagginger Auer and Bisseling [START_REF] Fagginger Auer | A GPU algorithm for greedy graph matching[END_REF] study an implementation of a greedy graph matching on GPU. Halappanavar et al. [START_REF] Halappanavar | Approximate weighted matching on emerging manycore and multithreaded architectures[END_REF] also study approximate matching on GPU. C ¸atalyürek et al. [START_REF] Deveci | Multithreaded clustering for multi-level hypergraph partitioning[END_REF] propose different greedy graph matching algorithms for multicore architectures. Azad et al. [START_REF] Azad | Multithreaded algorithms for maximum matching in bipartite graphs[END_REF] 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

A bipartite graph G = (V 1 ∪ V 2 , E) consists of a set of vertices V 1 ∪ V 2 where V 1 ∩ V 2 = ∅,
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. [START_REF] Duff | Design, implementation, and analysis of maximum transversal algorithms[END_REF]). Push-relabelbased algorithms form a second class [START_REF] Goldberg | A new approach to the maximum flow problem[END_REF]. A third class, pseudoflow algorithms, is based on a more recent work [START_REF] Hochbaum | The pseudoflow algorithm and the pseudoflow-based simplex for the maximum flow problem[END_REF]. There are O( √ nτ ) algorithms in the first two classes (e.g., Hopcroft-Karp algorithm [START_REF] Hopcroft | An n 5/2 algorithm for maximum matchings in bipartite graphs[END_REF] and a variant of the push relabel algorithm [START_REF] Goldberg | Global price updates help[END_REF]), 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 the ones in the third class have the runtime complexity of O(nτ ). These three classes of algorithms are described and compared in a recent study [START_REF] Kaya | Push-relabel based algorithms for the maximum transversal problem[END_REF]. 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 [START_REF] Duff | Design, implementation, and analysis of maximum transversal algorithms[END_REF][START_REF] Kaya | Push-relabel based algorithms for the maximum transversal problem[END_REF] and the original resources cited in those papers for other algorithms).

Augmenting-path-based algorithms follow the same common pattern: given an initial matching M (possibly empty), they search for an M-augmenting path P. If none exists, then M is maximum by a theorem of Berge [START_REF] Berge | Two theorems in graph theory[END_REF]. Otherwise, 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 edges of P were unmatched in M, we have |M ⊕ E(P)| = |M| + 1. The augmenting-path-based algorithms differ in the way the augmenting paths are found and the associated augmentations are realized. They mainly use either breadth-first-search (BFS), or depth-first-search (DFS), or a combination of them to locate and perform the augmenting paths.

Multicore counterparts of a number of augmenting-path based algorithms are proposed in a recent work [START_REF] Azad | Multithreaded algorithms for maximum matching in bipartite graphs[END_REF]. The parallelization of these algorithms is achieved by using atomic operations at BFS and/or DFS steps of the algorithm. Although atomic operations might not harm the performance on a multicore machine, they might cause a significant performance degradation on a GPU due to the possible serialization 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 [START_REF] Hopcroft | An n 5/2 algorithm for maximum matchings in bipartite graphs[END_REF] and HKDW [START_REF] Duff | Remarks on implementation of O(n 1/2 τ ) assignment algorithms[END_REF]. HK has the best known worst-case runtime 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 runtime while having the same 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. 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 alternates 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 runtime with a GPU.

Methods

We propose two GPU-based algorithms which find the augmenting paths via BFS, speculatively perform some of them, and fix any inconsistencies that can be resulting from speculative augmentations. The proposed algorithms exploit the GPU's vast number of thread parallelisms by assigning each vertex to a thread. Then, the threads concurrently perform independent operations for vertices in each kernel call, even though actual work is done for a portion of the vertices. Therefore, the GPU algorithm differs from a multi-core algorithm in which a shared data structure is used with atomic operations.

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 FixMatching functions. Algorithm 1 is given the adjacency list of the bipartite graph with its number of rows and columns. Any prior matching is given in rmatch and cmatch arrays as follows: 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.

rmatch[r] = c and cmatch[c] = r, if the row r is matched to the column c; rmatch[r] = -1, if r is unmatched; cmatch[c] = -1, if c is unmatched.
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.

We propose two implementations of the BFS kernel. Algorithm 2 is the first one. The BFS kernel is responsible from a single level BFS expansion. That is, it 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 column vertices in a single dimension. Each thread with id tid is assigned a number of columns which is obtained via the following function:

getP rocessCount(nc) = nc tot thread num if tid < nc mod tot thread num, nc tot thread num otherwise.
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 gives the description of the Alternate function. This kernel alternates the matched edges with the unmatched 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 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 set to -2 in the BFS kernel, only the threads whose row vertices' matches are -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 FixMatching kernel which implements:

rmatch[r] ← -1 for any r satisfying cmatch[rmatch[r]] = r.
Algorithm 4 gives the description of a slightly different BFS kernel function. This function takes a root array as an extra argument. Initially, the root array is filled with root

[c] = 0 if cmatch[c] > -1, and root[c] = c if cmatch[c] = -1.
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 the 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 proposed implementations are compared against the sequential HK and PFP implementations [START_REF] Duff | Design, implementation, and analysis of maximum transversal algorithms[END_REF], and against the multicore implementations P-PFP, P-DBFS, and P-HK obtained from [START_REF] Azad | Multithreaded algorithms for maximum matching in bipartite graphs[END_REF]. The CPU implementations are tested on a computer with 2.27GHz dual quad-core Intel Xeon CPUs with 2-way hyperthreading 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. 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 [START_REF] Duff | Design, implementation, and analysis of maximum transversal algorithms[END_REF]) is used to initialize all algorithms. We compare the runtime of the matching algorithms after this common initialization. The execution times of the GPU algorithms exclude memory copy time. But including memory copy time decreases the reported mean speedups across all data set by at most 6%.

The two main algorithms APFB and APsB can use different BFS kernel functions (GPUBFS and GPUBFS-WR). Moreover, 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.

We used 70 different matrices from variety of classes at UFL matrix collection [START_REF] Davis | The University of Florida sparse matrix collection[END_REF]. 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 [START_REF] Duff | Design, implementation, and analysis of maximum transversal algorithms[END_REF]. 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 runtime. 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 due to the unnecessary BFSs 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 BFSs not only increase the 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. As clearly seen from the figures, APFB variants converge 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 set of graphs, although the augmenting path exploration of APsB converges in larger number of iterations, the numbers of the BFS levels in the 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 is almost always the best algorithm, 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. )"

x" GPU" P1PF" P1DBFS" P1HK" HK"

(b) Permuted graphs Fig. 3. Log-scaled speedup profiles. show the performance profiles of the GPU and multicore algorithms. A point (x, y) in the plots means that with y probability, the corresponding algorithm obtains a performance that is at most x times worse than the best runtime. The plots clearly mark the GPU algorithm as the fastest in most cases, especially for the original graphs and for x ≤ 7 for the permuted ones. In particular, the GPU algorithm obtains the best performance in 61% of the original graphs, while this ratio increases to 74% for the permuted ones.

0.0# 0.2# 0.4# 0.6# 0.8# 1.0# 1# 2# 3# 4# 5# 6# 7# 8# 9# 10# Frac1on#of#Test#Cases# GPU$ P%PF$ P%DBFS$ P%HK$ (a) Original graphs 0.0# 0.2# 0.4# 0.6# 0.8# 1.0# 1# 2# 3# 4# 5# 6# 7# 8# 9# 10# Frac1on#of#test#cases# GPU$ P%PF$ P%DBFS$ P%HK$ (b) Permuted graphs
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. Moreover, the runtimes obtained by the GPU algorithm are robust among the repeated executions on a graph instance. We calculated the standard deviation of the execution times for different repetitions. For the graphs in O S1 set, the ratios of the standard deviations to the average time are observed to be less than 10%, 18%, and 47% for 20, 5, and 3 graphs, respectively. 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 GPU implementation of a BFS-based maximum cardinality matching algorithm for bipartite graphs. We compared the performance of the proposed implementation against sequential and multicore algorithms on various datasets. The experiments showed that the GPU implementation is faster than the existing multicore implementations. The speedups achieved with respect to well-known sequential implementations varied from 0.03 to 629.19, averaging 9.29 w.r.t. the fastest sequential algorithm on a set of 20 hardest problems. A GPU is a restricted memory device. Although, an out-of-core or distributedmemory 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 matching algorithms for extreme-scale bipartite graphs on GPUs.

Algorithm 1 : 5 vertex inserted ← true; 6 while vertex inserted do 7 predecessor

 1567 Shortest augmenting paths (APsB) Data: cxadj, cadj, nc, nr, rmatch, cmatch 1 augmenting path f ound ← true; 2 while augmenting path f ound do 3 bf s level ← L0; 4 InitBfsArray(bf s array, cmatch, L0); ←Bfs(bf s level, bf s array, cxadj, cadj, nc, rmatch, 8 vertex inserted, augmenting path f ound); 9 if augmenting path f ound then 10 break; 11 bf s level ← bf s level + 1; 12 cmatch, rmatch ← Alternate (cmatch, rmatch, nc, predecessor); 13 cmatch, rmatch ← FixMatching (cmatch, rmatch);

3

  col vertex ← i × tot thread num + tid; 4 if bf s array[col vertex] = bf s level then 5 for j from cxadj[col vertex] to cxadj[col vertex + 1] do 6 neighbor row ← cadj[j]; 7 col match ← rmatch[neighbor row]; 8 if col match > -1 then 9 if bf s array[col match] = L0 -1 then 10 vertex inserted ← true; 11 bf s array[col match] ← bf s level + 1; 12 predeccesor[neighbor row] ← col vertex; 13 else 14 if col match=-1 then 15 rmatch[neighbor row] ← -2; 16 predeccesor[neighbor row] ← col vertex; 17 augmenting path f ound ← true; takes the set of vertices at a BFS level and adds the union of the unvisited neighbors of those vertices as the next level of vertices. Initially, the input bf s array is filled with bf s array[c] = L 0 -1 if cmatch[c] > -1 and bf s array[c] = L 0 if cmatch[c] = -1 by a simple InitBfsArray kernel (L 0 denotes BFS start level).

3Fig. 1 .

 1 Fig. 1. Vertices r1 and c2 are matched; others are not. Two augmenting paths starting from c1 are possible.

Algorithm 4 : 18 if

 418 BFS Kernel Function-2 (GPUBFS-WR)Data: bf s level, bf s array, cxadj, cadj, nc, rmatch, root vertex inserted, augmenting path f ound 1 process cnt ← getP rocessCount(nc); 2 for i from 0 to process cnt -1 do3 col vertex ← i × tot thread num + tid; 4 if bf s array[col vertex] = bf s level then 5 myRoot ← root[col vertex]; 6 if bf s array[myRoot] < L0 -1 then7 continue; 8 for j from cxadj[col vertex] to cxadj[col vertex + 1] do 9 neighbor row ← cadj[j]; 10 col match ← rmatch[neighbor row]; 11 if col match > -1 then 12 if bf s array[col match] = L0 -1 then 13 vertex inserted ← true; 14 bf s array[col match] ← bf s level + 1; 15 root[col match] ← myRoot; 16 predeccesor[neighbor row] ← col vertex; 17 else col match=-1 then 19 bf s array[myRoot] ← L0 -2; 20 rmatch[neighbor row] ← -2; 21 predeccesor[neighbor row] ← col vertex; 22 augmenting path f ound ← true;

Fig. 2 .

 2 Fig.2. The BFS ids and the number of kernel executions for each BFS in APsB and APFB variants for two graphs. The x axis shows the id of the while iteration at line 2 of APsB. The y axis shows the number of the while iterations at line 6 of APsB.

Fig. 4 .

 4 Fig. 4. Performance profiles.

Figures 4 (

 4 Figures 4(a) and 4(b) show the performance profiles of the GPU and multicore algorithms. A point (x, y) in the plots means that with y probability, the corresponding algorithm obtains a performance that is at most x times worse than the best runtime. The plots clearly mark the GPU algorithm as the fastest in most cases, especially for the original graphs and for x ≤ 7 for the permuted ones. In particular, the GPU algorithm obtains the best performance in 61% of the original graphs, while this ratio increases to 74% for the permuted ones.Figure5gives 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. Moreover, the runtimes obtained by the GPU algorithm are robust among the repeated executions on a graph instance. We calculated the standard deviation of the execution times for different repetitions. For the graphs in O S1 set, the ratios of the standard deviations to the average time are observed to be less than 10%, 18%, and 47% for 20, 5, and 3 graphs, respectively.Table2gives the actual runtime for O Hardest20 set for the best GPU and multicore algorithms, together with the sequential algorithms. The compete set of runtimes can be found at http://bmi.osu.edu/hpc/software/ matchmaker2/maxCardMatch.html. As table shows, except six instances among

Fig. 5 .

 5 Fig. 5. Overall speedup of the proposed GPU algorithm w.r.t. PFP (left bars) and HK (right bars) algorithms.

Table 1 .

 1 Geometric mean of the runtime (in seconds) of the GPU algorithms on different sets of instances. Table 1 compares the proposed GPU implementations on different sets. As we see from the table, using a constant number of threads (CT) always increases

		APFB	APsB
		GPUBFS GPUBFS-WR GPUBFS GPUBFS-WR
		MT CT MT	CT MT CT MT	CT
	O S1	2.96 1.89 2.12	1.34 3.68 2.88 2.98	2.27
	O Hardest20	4.28 2.70 3.21	1.93 5.23 4.14 4.20	3.13
	RCP S1	3.66 3.24 1.13	1.05 3.52 3.33 2.22	2.14
	RCP Hardest20 7.27 5.79 3.37	2.85 12.06 10.75 8.17	7.41

Table 2

 2 gives the actual runtime for O Hardest20 set for the best GPU and multicore algorithms, together with the sequential algorithms. The compete set of runtimes can be found at http://bmi.osu.edu/hpc/software/ matchmaker2/maxCardMatch.html. As table shows, except six instances among

		12%
		10%
	Speedup&	6% 8%
		4%
		2%
		0%	O_S1%	O_Hardest20%	RCP_S1%	RCP_Hardest20%
		PFP%	4.26%	5.58%	3.54%	9.29%
		HK%	3.61%	3.96%	7.40%	9.62%

Table 2 .

 2 The actual runtime of each algorithm for the O Hardest20 set.

			Original graphs		Permuted graphs
	Matrix name	GPU P-DBFS	PFP	HK GPU P-DBFS PFP	HK
	roadNet-CA	0.34	0.53	0.95	2.48 0.39	1.88	3.05	4.89
	delaunay n23	0.96	1.26	2.68	1.11 0.90	5.56	3.27 14.34
	coPapersDBLP	0.42	6.27	3.11	1.62 0.38	1.25 0.29	1.26
	kron g500-logn21	0.99	1.50	5.37	4.73 3.71	4.01 64.29 16.08
	amazon-2008	0.11	0.18	6.11	1.85 0.41	1.37 61.32	4.69
	delaunay n24	1.98	2.41	6.43	2.22 1.86	12.84	6.92 35.24
	as-Skitter	0.49	1.89	7.79	3.56 3.27	5.74 472.63 29.63
	amazon0505	0.18	22.70	9.05	1.87 0.24	15.23 17.59	2.23
	wikipedia-20070206 1.09	5.24	11.98	6.52 1.05	5.99	9.74	5.73
	Hamrle3	1.36	2.70	0.04	12.61 3.85	7.39 37.71 57.00
	hugetrace-00020	7.90 393.13	15.95	15.02 1.52	9.97	8.68 38.27
	hugebubbles-00000 13.16	3.55	19.81	5.56 1.80	10.91 10.03 38.97
	wb-edu	33.82	8.61	3.38	20.35 17.43	20.10 9.49 51.14
	rgg n 2 24 s0	3.68	2.25	25.40	0.12 2.20	12.50	5.72 31.78
	patents	0.88	0.84	92.03	16.18 0.91	0.97 101.76 18.30
	italy osm	5.86	1.20	1.02 122.00 0.70	3.97	6.24 18.34
	soc-LiveJournal1	3.32	14.35 243.91	21.16 3.73	7.14 343.94 20.71
	ljournal-2008	2.37	10.30 360.31	17.66 6.90	7.58 176.69 23.45
	europe osm	57.53	11.21	14.15 1911.56 7.21	37.93 68.18 197.03
	com-livejournal	4.58	22.46 2879.36	34.28 5.88	17.19 165.32 29.40


Kamer Kaya

Bora Uçar
email: bora.ucar@ens-lyon.fr

Ümit V C ¸atalyürek

Analysis of Partitioning Models and Metrics in Parallel Sparse Matrix-Vector Multiplication

Keywords:

  L'archive ouverte pluridisciplinaire

Introduction

Repeated sparse matrix-vector (SpMxV) and sparse matrix-transpose-vector multiplies that involve the same large, sparse matrix are the kernel operations in various iterative algorithms involving sparse linear systems. Such iterative algorithms include solvers for linear systems, eigenvalues, and linear programs. Efficient parallelization of SpMxV operations is therefore very important in virtually all large scale scientific computing applications. A number of partitioning methods and models based on hypergraphs have been used to enable efficient parallelization of SpMxV. These partitioning methods address different communication cost metrics for some variants of parallel SpMxV operations. In general, the importance of the communication cost metrics, such as the total volume of communication, the total number of messages and these two quantities on per processor basis, depends on the machine architecture, problem size, and the underlying parallel algorithm. In this study, we investigate the effects of the partitioning methods in order to identify the most relevant metrics and quantify their effects in various configurations. Our aims are to help the partitioner developers identify the important metrics, and to help the users of those partitioners to identify the most suitable partitioning method for their use case.

The standard hypergraph based models minimize the total volume of communication explicitly [START_REF] Aykanat | Hypergraph-partitioning-based decomposition for parallel sparse-matrix vector multiplication[END_REF]. Some more recent variants do that while imposing a limit on the total number of communication by a 2D partitioning approach [START_REF] Aykanat | A hypergraph-partitioning approach for coarsegrain decomposition[END_REF][START_REF] Aykanat | On two-dimensional sparse matrix partitioning: Models, methods, and a recipe[END_REF]. More sophisticated approaches [START_REF] Bisseling | Communication balancing in parallel sparse matrixvector multiplication[END_REF][START_REF] Uçar | Minimizing communication cost in fine-grain partitioning of sparse matrices[END_REF][START_REF] Uçar | Encapsulating multiple communication-cost metrics in partitioning sparse rectangular matrices for parallel matrix-vector multiplies[END_REF] minimize different communication cost metrics on top of the total volume. Experimental investigations in these studies demonstrate that different communication cost metrics and their interplay can be important to achieve scalable parallel algorithms. It is therefore important to understand the effects of different metrics (optimized by different partitioning models) on the running time of applications under different configurations.

The contribution of this paper is two-fold. We designed and conducted several experiments in a system with 512 processors to show the effects of partitioning models and metrics on SpMxV performance. As far as we know, this is the first work which compares the existing partitioning models and metrics in modern architectures with modern software following message-passing paradigm. Our experiments confirm that it is difficult, if not impossible, to define the correct partitioning model and metric without analyzing the characteristics of the input matrices and the SpMxV library being used. We experimented with three existing libraries, PETSc [START_REF] Balay | PETSc users manual[END_REF][START_REF] Balay | Efficient management of parallelism in object oriented numerical software libraries[END_REF], Trilinos [START_REF] Heroux | An overview of the trilinos project[END_REF], and an in-house library SpMV [START_REF] Uçar | A library for parallel sparse matrix-vector multiplies[END_REF]. In order to overcome the mentioned difficulty, we carefully analyze the results using regression analysis techniques and relate the execution time of SpMxV implementations to different partitioning metrics. We portray this analysis, which forms out the second contribution, in detail so as to suggest improved objective functions for partitioning software and a guideline to choose partitioning methods for practitioners. Although, we only had an access to a 512-processor machine, the experiments and their analysis show that to scale larger systems, one needs to be more careful while partitioning the matrix-in our experiments the fact that the communication metrics greatly related to the execution time is observable starting from 64 processors.

Parallel SpMxV operation and software

Consider the sparse matrix-vector multiply operation of the form y ← Ax, where the nonzeros of the m × n matrix A are partitioned among K processors such that each processor P k owns a mutually disjoint subset of nonzeros, A (k) where A = k A (k) . The vectors y and x are also partitioned among processors, where the processor P k holds x (k) , a dense vector of size n k , and it is responsible for computing y (k) , a dense vector of size m k .

The standard parallel SpMxV algorithm [START_REF] Aykanat | On two-dimensional sparse matrix partitioning: Models, methods, and a recipe[END_REF][START_REF] Uçar | A library for parallel sparse matrix-vector multiplies[END_REF][START_REF] Vastenhouw | A two-dimensional data distribution method for parallel sparse matrix-vector multiplication[END_REF] based on the described nonzero and vector entry partitioning is called the row-column-parallel algorithm. In this algorithm, each processor P k executes the following steps:

1. Expand: Send entries of x (k) that are needed by others. Receive entries of

x that are needed but owned by others.

2. Scalar multiply-adds: Perform ȳ ← A (k) x, where x contains x (k) and the received entries of x. 3. Fold: Send partial results from ȳ to the responsible processors. Receive contributions to the y (k) vector.

If A is distributed columnwise, and the x-vector entries are partitioned conformably with the column partition of A, then the expand operation is not needed. Similarly, if A is distributed rowwise, and the y-vector entries are partitioned conformably with the rows of A, then the fold operation is not needed.

Libraries

There are different implementations of the above algorithm. We summarize three implementations with which we have experimented. Two of the implementations are in the well-known general libraries Trilinos [START_REF] Heroux | An overview of the trilinos project[END_REF] and PETSc [START_REF] Balay | PETSc users manual[END_REF][START_REF] Balay | Efficient management of parallelism in object oriented numerical software libraries[END_REF]; the third one, SpMV, is an in-house library [START_REF] Uçar | A library for parallel sparse matrix-vector multiplies[END_REF]. Trilinos provides an implementation which can be described as in Algorithm 1 from the point of view of the processor P k . In this implementation, the expand operations are finished before doing any computation. Then, all the scalar multiply-add operations are performed. Later on, the fold operations are completed. Trilinos uses Irecv/Isend and waitall communication primitives to handle the communications at steps 1 and 2 of Algorithm 1. It issues Irecvs, performs Isends and then before commencing the computations ensures that all in the incoming data is received by using the waitall operation. PETSc provides an implementation of the above algorithm only for the rowparallel case. Algorithm 2 summarizes that implementation from the point of view of P k . First, the expand operation is initiated using Irecv and Isend primitives. Then, instead of waiting the reception of all necessary x-vector entries, it performs some local computations so as to overlap communication and computations. In particular, the processor P k performs scalar multiply-add operations using local a ij 's for which µ(x j ) = P k and there is no a i with µ(x ) = P k . Then, upon verifying the reception of all needed x-vector entries using waitall, P k continues with scalar multiply-add operations with the nonzeros on the rows that has at least one nonzero in a column j for which µ(x ) = P k . The implementation can also be seen in an earlier technical report [START_REF] Saad | P-SPARSLIB: A portable library of distributed memory sparse iterative solvers[END_REF]. Figure 1 describes the algorithm pictorially. After issuing Isends and Irecvs (for x(k) ), processor P k performs the computations associated with the horizontally shaded matrix zone. Then, waitall is executed to have all x (k) before continuing with the rows that are below the horizontal ones. Note that the local matrices are actually permuted into the displayed form (local rows and the interface rows). The advantage of this implementation with respect to the Algorithm 1 is that it allows overlap between the reception of messages for the expand operation and scalar multiply-add operations with the nonzeros in local rows. Consider again the matrix A (k) of processor P k as shown in Fig. 1. Before executing the waitall operation, are some more scalar multiply-add operations that P k can perform before the reception of any x(k) . These operations are related to the nonzeros that are in the hatched zone in the figure. In order to exploit the hatched zone for communication computation overlap, one can store that zone in the compressed column storage (CCS) format. This way, one can delay the invocation of the waitall operation for some more time. In fact, we can get rid of the waitall operation and maximize the communication computation overlap by performing all scalar multiply-operations that involve a received x-vector entry before waiting the reception of any other message. This requires storing the vertically shaded zones of the matrix in Fig. 1 in CCS format: with this, when P k receives x(k) , it can visit the respective column and perform all operations. This way of storing the vertically shaded and hatched zones in CCS maximizes the amount of overlap in the strict sense (optimal amount of overlap) when a processor receives a single message from each sender (as should be the case in a proper SpMxV code). The third library that we investigate in this work, SpMV [START_REF] Uçar | A library for parallel sparse matrix-vector multiplies[END_REF], implements this approach for row-parallel and row-column parallel algorithms (see descriptions in the accompanying technical report [START_REF] Kaya | On analysis of partitioning models and metrics in parallel sparse matrix-vector multiplication[END_REF]).

x (k) ∧ 1 ... x (k) ∧ k ... x (k) ∧ K y (k) ∧

Investigated partitioning metrics and methods

We study the practical effects of five different metrics: the maximum number of nonzeros assigned to a processor (MaxNnz) which defines the load balance; the total communication volume (TotVol); the maximum send volume of a processor (MaxSV); the total number of messages (TotMsg); and the maximum number of messages sent by a processor (MaxMS). Our investigations are necessarily experimental, yet some a priori results can be told about these metrics, see the accompanying technical report [START_REF] Kaya | On analysis of partitioning models and metrics in parallel sparse matrix-vector multiplication[END_REF] and the references therein.

We investigated the rowwise and columnwise partitioning methods CN and RN (the naming convention corresponds to the hypergraph models in the original paper [START_REF] Aykanat | Hypergraph-partitioning-based decomposition for parallel sparse-matrix vector multiplication[END_REF]) among the one-dimensional partitioning approaches. Among the twodimensional ones we investigated the fine grain (FG), and checkerboard (CB) partitioning models (see [START_REF] Aykanat | On two-dimensional sparse matrix partitioning: Models, methods, and a recipe[END_REF] and the references therein). These four methods try to reduce the total volume of communication and obtain load balance; the 2D methods implicitly reduce the total and maximum number of messages per processor. We also used a block partitioning (BL) model whose aim is just to have load balance in rowwise partitioning of the matrices. In this model, we traverse the rows from 1 to m, generate a part with approximately τ /K nonzeros, and continue with the next part when this number is exceeded. For matrices based on 2D meshes, we used another rowwise partitioning model called MP. This model tiles the 2D plane with a diamond-like shape [4, Section 4.8] and associates each shape (corresponding to a set of rows) with a processor. This approach balances the volume and the number of messages the processors send and receive.

Experimental investigations

We carried our experiments on a 64-node cluster where each node has a 2.27GHz dual quad-core Intel Xeon (Bloomfield) CPU and 48GB main memory. Each core in a socket has 64KB L1 and 256KB L2 caches, and each socket has an 8MB L3 cache shared by 4 cores. The interconnection network is 20Gbps DDR InfiniBand. For parallelism, mvapich2 version 1.6 is used. We built SpMV, PETSc, and Trilinos with gcc 4.4.4 and used optimization flag -O3. For PETSc experiments, we used the matrix type MPIAIJ and the multiplication routine MatMult. We used PaToH [START_REF] Aykanat | PaToH: A Multilevel Hypergraph Partitioning Tool[END_REF] with default setting quality for partitioning the matrices (this allows 0.03 imbalance). Since each node has 8 cores, we have 512 processors in total. In the experiments, we use K ∈ {1, 8, 16, 32, 64, 128, 256, 512}. For an experiment with K = 1 processors, we fully utilize K/8 nodes of the cluster. To measure the time of one SpMxV operation (in secs), we do 500 multiplications for each execution. The tables and figures show the averages of these 500 runs.

We used seven large real-life square matrices from different application domains that are available at the University of Florida (UFL) Sparse Matrix Collec-tion (http://www.cise.ufl.edu/research/sparse/matrices) and three synthetically generated matrices corresponding to [START_REF] Aykanat | PaToH: A Multilevel Hypergraph Partitioning Tool[END_REF] In the experiments on real-life matrices, we use PETSc with rowwise models CN and BL, and SpMV and Trilinos with all models except MP. For meshes, we added MP to each library's model set. The reason is technological: PETSc provides SpMxV routines only for rowwise partitioning.

We designed a set of experiments to show the effect of different partitioning metrics in the actual running time of SpMxV computations. Our first two sets of experiments (Fig. 2 and Table . 2 of the accompanying technical report [START_REF] Kaya | On analysis of partitioning models and metrics in parallel sparse matrix-vector multiplication[END_REF]) showed clearly that the total volume of communication, the total number of messages, the maximum number and volume of messages sent by a processor affect significantly the running time of SpMxV. Which one of these four communication metric is the most important in different libraries? To what extent? What about their combinations? In order to answer these questions we carried out the following regression analysis.

Regression analysis

To evaluate the performance of the libraries with respect to the partitioning metrics, we use linear regression analysis techniques and solve the nonnegative least squares problem (NNLS). In NNLS, given a variable matrix V and a vector t, we want to find a dependency vector d which minimizes Vdt s.t. d ≥ 0. In our case, V has five columns which correspond to the partitioning metrics MaxNnz, TotVol, MaxSV, TotMsg, and MaxSM. Each row of V corresponds to an SpMxV execution where the execution time is put to the corresponding entry of t. Hence, we have the same V but a different t for each library. We apply a well-known technique in regression analysis and standardize each entry of V by subtracting its column's mean and dividing it to its column's standard deviation so that the mean and the standard deviation of each column become 0 and 1, respectively. This way, the units are removed and each column becomes equally important throughout the analysis. We then used MATLAB's lsqnonneg to solve NNLS. Each entry of the output d i shows the dependency of the execution time to the partitioning metric corresponding to the ith column of V. Tables 2,3, and 4 show the dependency values found in various settings.

We first apply regression analysis to each library with all matrices and rowwise partitioning models CN (column-net) and BL (block). The analysis shows that when K ≤ 64, SpMxV performance depends rigorously on the maximum number of nonzeros assigned to a processor. In this case, the dependency values for MaxNnz are 8.02, 7.81, and 6.80 for SpMV, PETSc, and Trilinos, respectively. As Table 2 shows, the next important metric is MaxSV with values 1.66, 1.53, and 2.20. The latency-based (TotMsg, MaxSM) partitioning metrics do not effect the performance for K ≤ 64. However, when K gets larger, these metrics are of utmost importance. Furthermore, the importance of MaxNnz decreases drastically for all the libraries. For SpMV and PETSc, MaxNnz becomes the 3rd important variable, whereas for Trilinos, it is the 4th. This shows that SpMV and PETSc handle the increase in the communication metrics better than Trilinos. When K ≥ 128, the dependency of Trilinos to TotMsg is much less than that of SpMV and PETSc. On the contrary, Trilinos' MaxSM dependency is almost 1.75 times more than SpMV and PETSc. This is expected since Trilinos uses Algorithm 1 which has no communication-computation overlap due to the use of waitall primitive. Such primitives can cause close coupling among the processors. When MaxNnz and the variance on the number of messages per processor are large, the overhead due to the bottleneck processor can result in poor SpMxV performance. Note that the dependency profiles of SpMV and PETSc, which are similar due to the communication-computation overlap, do not point out a similar bottleneck.

We extend the regression analysis to all matrices and all partitioning models and show the results in Table 3. The performance of SpMV and Trilinos rigorously depend on MaxNnz if K ≤ 32, and on TotMsg and MaxSM when K ≥ 256. Once again, Trilinos' MaxSM dependency is higher than that of SpMV due to the waitall primitive. To see the effect of matrix structure on regression analysis, we use only mesh-based matrices in the next experiment. As Table 4 shows, we observe that for these matrices, the performance of SpMV and Trilinos mostly depend on MaxNnz even when K ≥ 64. Note that these matrices are banded and the communication metrics have relatively lower values compared to those of real-life matrices. Hence, the most dominant factor is MaxNnz.

In the light of the regression analysis experiments, we note that the partitioning metrics effect the performance of parallel SpMxV libraries. The best metric (or function of metrics) that needs to be minimized depends on the number of processors, the size and structure of the matrix, which we are planning to investigate in the future, and even the library itself. Although some of these variables are known while generating the partitions, predicting the others may need a preprocessing phase. For example, we already know that the libraries in this paper are employing point-to-point communication primitives which makes the connectivity metric suitable. However, if collective communication primitives, e.g., MPI ALLGATHER, had been used, it would be better to minimize the cut-net metric as the main partitioning objective (however, we should note that such collective operations introduce unnecessary synchronization and messages especially for large K values). On the other hand, the matrix structure can be different for each input and a partitioner needs either a manual direction or a preprocessing to predict the best metric for each matrix.

Summary of further results

We provide summary of some further results and refer the interested reader to the accompanying technical report [START_REF] Kaya | On analysis of partitioning models and metrics in parallel sparse matrix-vector multiplication[END_REF]. We observed that for all libraries and partitioning methods, minimizing TotMsg is more important than minimizing MaxSV for reducing the execution time, especially when K is large. Starting from K = 64, the difference becomes obvious in favor of TotMsg which is concordant with the regression analyses. For K ∈ {8, 16}, minimizing MaxSV or TotMsg are equally important.

The checkerboard method (CB) which is demonstrated (see Fig. 4 in [START_REF] Kaya | On analysis of partitioning models and metrics in parallel sparse matrix-vector multiplication[END_REF]) to reduce most of the communication cost metrics seems to be the method of choice when K is not small (when K is small, there is no much difference between the models as also revealed by the regression analysis before). The relative performances of the partitioning methods do not change with the increasing K. However, their difference tend to increase and hence, the model used for partitioning becomes more important as the parallel matrix-vector multiplication times of the libraries show in Figure 2. When K is 256, the only significant reduction on the execution time is obtained by SpMV with the CB model.

Conclusion

We have carried out a detailed study to understand the importance of partitioning models and their effects in parallel SpMxV operations. As mentioned in the experiments, minimizing the right metric with the right partitioning model is crucial to increase throughput. For example, for the real-life matrices in our test set, CB model is the only one which can obtain a significant reduction on the SpMxV time when K is increased from 128 to 256 (after that we did not see any speed up). It is obvious that the other models fail to obtain such a reduction since the gain by dividing MaxNnz by two does not compensate the communication overhead induced by multiplying K by two. Hence, assuming the communication overhead is doubled on the average, doubling K increases the relative importance of communication on SpMxV four times.

Matrices from today's scientific and industrial applications can be huge. If one has only a few processors, partitioning may not matter, since the contribution of communication to the execution time will be low and the overall improvement on SpMxV via a good partitioning will be insignificant. However, as the regression analyses of Section 3.1 show, after a number of processors, the communication overhead will start to dominate the SpMxV time. For our experiments, this number is somewhere between 32 and 64, and it depends on the characteristics of the matrix, the library, and the architecture used for SpMxV operations. Although it may be more than 64, considering the advancements on CPU hardware, we can easily argue that this number will remain highly practical and partitioning will matter more for systems that are larger than those considered here.

Algorithm 1 :3

 1 ParSpMxV-Trilinos variantInput: A, x, µ Output: y 1 Send and Receive x vector entries so that each processor has the requiredx-vector entries2 Compute y (k) i← aij xj for the local nonzeros, i.e., the nonzeros for which µ(aij) = P k Send and Receive local nonzero partial results y

4

  µ(yi) = P k , for all nonzero y (k) i Compute yi ← y i for each yi with µ(yi) = P k

Algorithm 2 :

 2 ParSpMxV-Overlap-PETSc variant Input: A, x, µ Output: y 1 Send local xj (i.e., µ(xj) = P k ) to those processors that have at least one nonzero in column j 2 Compute y k i ← aij xj for the local nonzeros and local xj, i.e., the nonzeros for which µ(aij) = P k and µ(xj) = P k 3 Receive all non-local xj (i.e., µ(xj) = P k ) 4 Compute y k i ← y k i +aij xj for the local nonzeros and non-local xj, i.e., the nonzeros for which µ(aij) = P k and µ(xj) = P k

Fig. 1 :

 1 Fig. 1: The zones of the matrix A (k) of processor P k with respect to the vector x assuming a row-parallel algorithm.

Fig. 2 :

 2 Fig. 2: Mean SpMxV times on real-life matrices in log scale for each library with respect to partitioning model.

Table 1 :

 1 -point stencil meshes in 2D with sizes 1024 × 1024, 2048 × 2048, and 4096 × 4096. The properties of the matrices are given in Table 1. Properties of the experiment matrices.

	Matrix	Description	n	τ	Matrix	Description	n	τ
	atmosmodl Atmosp. model. 1,489,752 10,319,760 cage15	DNA electrop. 5,154,859 99,199,551
	TSOPF RS Opt. pow. flow	38,120 16,171,169 HV15R	3D engine fan 2,017,169 283,073,458
	Freescale1 Semicon. sim.	3,428,755 17,052,626 mesh-1024 5-point stencil 1,048,576 5,238,784
	rajat31	Circuit sim.	4,690,002 20,316,253 mesh-2048 5-point stencil 4,194,304 20,963,328
	RM07R	Comp. fluid dyn. 381,689 37,464,962 mesh-4096 5-point stencil 16,777,216 83,869,696

Table 2 :

 2 Regression analysis of SpMV, PETSc and Trilinos with all matrices and models CN and BL.

			8 ≤ K ≤ 64		128 ≤ K ≤ 512
	Metric	SpMV PETSc Trilinos SpMV PETSc Trilinos
	MaxNnz	8.02	7.81	6.80	0.49	0.44	0.83
	TotVol	0.18	0.38	1.00	0.39	0.36	1.06
	MaxSV	1.66	1.53	2.20	0.00	0.00	0.11
	TotMsg	0.15	0.28	0.00	7.90	8.03	4.51
	MaxSM	0.00	0.00	0.00	1.22	1.18	3.49
		8 ≤ K ≤ 32	64 ≤ K ≤ 128	256 ≤ K ≤ 512
	Metric	SpMV Trilinos SpMV Trilinos SpMV Trilinos
	MaxNnz	8.43	7.54	2.75	2.52	0.00	0.02
	TotVol	0.23	0.89	0.52	1.94	0.38	0.98
	MaxSV	1.35	1.57	1.57	1.69	0.04	0.50
	TotMsg	0.00	0.00	4.66	2.38	6.24	3.06
	MaxSM	0.00	0.00	0.49	1.47	3.34	5.44

Table 3 :

 3 Regression analysis of SpMV and Trilinos with all matrices and partitioning models. PETSc is not shown in this table because it cannot handle all the schemes.

		8 ≤ K ≤ 32	64 ≤ K ≤ 128	256 ≤ K ≤ 512
	Metric	SpMV Trilinos SpMV Trilinos SpMV Trilinos
	MaxNnz	8.97	9.38	8.83	9.05	5.10	5.47
	TotVol	0.00	0.00	0.00	0.24	0.00	0.00
	MaxSV	0.72	0.48	0.43	0.09	0.92	0.52
	TotMsg	0.00	0.00	0.42	0.07	0.42	0.99
	MaxSM	0.31	0.14	0.33	0.55	3.55	3.02

Table 4 :

 4 Regression analysis of SpMV and Trilinos with mesh-based matrices and all partitioning models.


Emil Björnson
email: emil.bjornson@supelec.fr

Agisilaos Papadogiannis
email: agisilaos.papadogiannis@chalmers.se

Michail Matthaiou
email: michail.matthaiou@chalmers.se

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

ON THE IMPACT OF TRANSCEIVER IMPAIRMENTS ON AF RELAYING

Keywords:

Recently, it was shown that transceiver hardware impairments have a detrimental impact on the performance of communication systems, especially for high-rate systems. The vast majority of technical contributions in the area of relaying assume ideal transceiver hardware. This paper quantifies the impact of transceiver hardware impairments in dual-hop Amplify-and-Forward (AF) relaying, both for fixed and variable gain relays. The outage probability (OP) in this practical scenario is a function of the instantaneous end-toend signal-to-noise-and-distortion ratio (SNDR). This paper derives closed-form expressions for the exact and asymptotic OPs under Rayleigh fading, accounting for hardware impairments at both the transmitter and the relay. The performance loss is small at low spectral efficiency, but can otherwise be very substantial. In particular, it turns out that for high signal-to-noise ratio (SNR), the instantaneous end-to-end SNDR converges to a deterministic constant, called the SNDR ceiling, which is inversely proportional to the level of impairments. This stands in stark contrast to the ideal hardware case for which the end-to-end SNDR grows without bound in the high SNR regime.

INTRODUCTION

The use of relays for improving coverage, throughput, and qualityof-service in wireless systems has been a hot research topic over the past decade, both in academia [START_REF] Laneman | Cooperative diversity in wireless networks: Efficient protocols and outage behavior[END_REF][START_REF] Hasna | A performance study of dualhop transmissions with fixed gain relays[END_REF] and in industry [START_REF] Hua | Theories and methods for advanced wireless relays-Issue I[END_REF]. This is due to the fact that, unlike base stations, relays are low-cost nodes that can be easily deployed and, hence, enhance the network agility. The vast majority of works in the context of relaying networks make the assumption of ideal transceiver hardware. However, in practice, hardware suffers from impairments; for example, phase noise, IQ imbalance, and amplifier non-linearities [START_REF] Costa | m-QAM-OFDM system performance in the presence of a nonlinear amplifier and phase noise[END_REF][START_REF] Schenk | RF Imperfections in High-Rate Wireless Systems: Impact and Digital Compensation[END_REF][START_REF] Studer | MIMO transmission with residual transmit-RF impairments[END_REF]. These have a deleterious impact on the achievable performance [START_REF] Zetterberg | Experimental investigation of TDD reciprocitybased zero-forcing transmit precoding[END_REF][START_REF] Björnson | Capacity limits and multiplexing gains of MIMO channels with transceiver impairments[END_REF][START_REF] Björnson | Optimal coordinated beamforming in the multicell downlink with transceiver impairments[END_REF][START_REF] Samuel | Impact of RF impairments in cellular wireless metropolitan area networks[END_REF][START_REF] Riihonen | BEP analysis of OFDM relay links with nonlinear power amplifiers[END_REF][START_REF] Qi | Analysis and compensation of I/Q imbalance in amplify-and-forward cooperative systems[END_REF][START_REF] Björnson | Optimal resource allocation in coordinated multi-cell systems[END_REF]. This effect is more pronounced in high-rate systems, especially those employing inexpensive hardware [START_REF] Schenk | RF Imperfections in High-Rate Wireless Systems: Impact and Digital Compensation[END_REF]. For instance, some recent works have demonstrated that non-ideal hardware severely affects single-hop multiantenna systems; [START_REF] Studer | MIMO transmission with residual transmit-RF impairments[END_REF][START_REF] Björnson | Capacity limits and multiplexing gains of MIMO channels with transceiver impairments[END_REF] proved that there is a finite capacity limit at high SNR, while [START_REF] Björnson | Optimal coordinated beamforming in the multicell downlink with transceiver impairments[END_REF][START_REF] Björnson | Optimal resource allocation in coordinated multi-cell systems[END_REF] showed that existing signal processing algorithms need to be re-designed to account for these impairments.

Despite the importance of transceiver impairments, their impact on relaying has only been partially investigated; in particular, bit error rate simulations were conducted in [START_REF] Samuel | Impact of RF impairments in cellular wireless metropolitan area networks[END_REF], while [START_REF] Riihonen | BEP analysis of OFDM relay links with nonlinear power amplifiers[END_REF][START_REF] Qi | Analysis and compensation of I/Q imbalance in amplify-and-forward cooperative systems[END_REF] derived expressions for the bit/symbol error rates, considering only non-linearities or IQ imbalance, respectively. In this paper, we try to bridge this gap by analytically evaluating the performance of dualhop AF relaying networks in the presence of aggregate transceiver impairments, both for fixed and variable gain relaying. In particular, we obtain the SNDRs and tractable closed-form expressions for the exact and asymptotic OPs of the system. This enables us to characterize the impact of impairments. To the best of our knowledge, this is the first paper presenting an analytical study of relaying with transceiver impairments under the generalized system model of [START_REF] Schenk | RF Imperfections in High-Rate Wireless Systems: Impact and Digital Compensation[END_REF]. One important observation is that for both types of relaying, an increase in the SNRs lead to an "SNDR ceiling" effect, which is explicitly quantified. The value of the SNDR ceiling is inversely proportional to the level of impairments. This manifests that transceiver impairments have a significant impact on AF relaying-especially at high SNRs and when high spectral efficiencies are desirable.

SIGNAL AND SYSTEM MODEL

This paper revisits the classical dual-hop relaying in Fig. 1a, where a source communicates with a destination through a relay using AF.

Ideal Hardware

In the classical complex baseband system model (see Fig. 1a), the received signals at the relay and destination are yi

= hisi + νi, i = 1, 2, (1)
where s1, s2 ∈ C are the transmitted signals from the source and relay, respectively, with average signal power Pi = Es i {|si|2 }, for i = 1, 2. Observe that E{•} denotes the expectation operator. In addition, νi ∼ CN (0, Ni) represents the additive circularly-symmetric complex Gaussian noise terms for i = 1, 2. The random channel fading coefficients are modeled as Rayleigh fading1 such that hi ∼ CN (0, Ωi), where Ωi E h i {|hi| 2 } is the average power of the fading channel. Recall that, in this case, the cumulative distribution functions of the channel gains, ρi = |hi| 2 , are given by

Fρ i (x) = 1 -e -x Ω i , x ≥ 0, i = 1, 2. (2)
The quantity P i Ω i N i will be referred to as the average SNR, for i = 1, 2. A high SNR is achieved by high signal and/or fading power.

The transmitted signal s2 at the relay is simply an amplified version of the signal y1 received at the relay: s2 = G y1. If the relay has instantaneous knowledge of the fading channel, h1, it can apply variable gain relaying with G v P2/Es 1 ,ν 1 {|y1| 2 } [START_REF] Emamian | Multi-user spatial diversity in a shadow-fading environment[END_REF]. Otherwise, fixed gain relaying with G f P2/E s 1 ,ν 1 ,h 1 {|y1| 2 } can be applied using only statistical information [START_REF] Hasna | A performance study of dualhop transmissions with fixed gain relays[END_REF]. 2 With ideal (id) hardware, the amplification factors can be expressed as

G f id = P2 P1 Ω1 + N1 , G v id = P2 P1 ρ1 + N1 . (3)
The received signal at the destination is

y2 = h2 Gid(h1s1 + ν1) + ν2. (4)
After some simple algebraic manipulations, the end-to-end SNDRs of (1) for variable and fixed gain relays are obtained as

γ f id = P1 ρ1 ρ2 ρ2 N1 +N2/(G f id ) 2 , γ v id = P1 P2 ρ1 ρ2 P2 ρ2 N1 + (P1ρ1 +N1)N2 , (5) respectively
. Similar expressions were previously derived in [START_REF] Hasna | A performance study of dualhop transmissions with fixed gain relays[END_REF][START_REF] Emamian | Multi-user spatial diversity in a shadow-fading environment[END_REF].

Non-Ideal Hardware: Transceiver Impairments

The model in (1) implicitly assumes ideal transceiver hardware. However, physical radio-frequency (RF) transceivers suffer from a variety of impairments that create a mismatch between the intended signal for transmission, si, and the signal that is actually generated and emitted. The impact of each impairment (e.g., phase-noise, IQ imbalance, and amplifier non-linearities) can be modeled separately for a given hardware setup [START_REF] Schenk | RF Imperfections in High-Rate Wireless Systems: Impact and Digital Compensation[END_REF][START_REF] Samuel | Impact of RF impairments in cellular wireless metropolitan area networks[END_REF]. Interestingly, the combined influence is well-modeled by the generalized system model

yi = hi(si + η tx i ) + η rx i + νi, i = 1, 2, (6)
where η tx i , η rx i are distortion noises in the transmitter and the receiver, respectively [START_REF] Schenk | RF Imperfections in High-Rate Wireless Systems: Impact and Digital Compensation[END_REF]. Experiments and theoretical investigations (e.g., [START_REF] Costa | m-QAM-OFDM system performance in the presence of a nonlinear amplifier and phase noise[END_REF][START_REF] Studer | MIMO transmission with residual transmit-RF impairments[END_REF][START_REF] Zetterberg | Experimental investigation of TDD reciprocitybased zero-forcing transmit precoding[END_REF]) have verified that these are well-modeled as

η tx i ∼ CN 0, (κ tx i ) 2 P , η rx i ∼ CN 0, (κ rx i ) 2 |hi| 2 P , (7)
where the Gaussianity can be explained by the aggregate effect of many impairments. 3 Note that the distortion depends on the intended signal for transmission/reception, in the sense that the variance of ηi is proportional to its average power, here denoted as P.

For convenience, we define ηi η tx i + η rx i /hi ∼ CN (0, κ 2 i P) as the aggregate distortion noise, where

κ 2 i = (κ tx i ) 2 + (κ rx i ) 2 .
The parameter κi > 0 describes the level of impairments on the ith channel and is an important parameter in this paper. Note that κi can be interpreted as the total error vector magnitude (EVM)-a common measure for quantifying the quality of RF transceivers. 4 This parameter is essentially constant while operating in the dynamic range of the power amplifier, but can otherwise increase very rapidly [START_REF] Zetterberg | Experimental investigation of TDD reciprocitybased zero-forcing transmit precoding[END_REF][START_REF] Björnson | Optimal coordinated beamforming in the multicell downlink with transceiver impairments[END_REF]. The generalized system model of ( 6) is illustrated in Fig. 1b.

For the non-ideal (ni) hardware counterpart of (4), the received signal at the destination is now obtained as y2 = h2 Gni h1 (s1 + η1)

+ ν1 + h2 η2 + ν2 (8)
= Gni h1 h2 s1 + Gni h1 h2 η1 + Gni h2 ν1 + h2 η2 + ν2.

For fixed and variable gain relaying, Gni reads respectively as

G f ni P2 P1 Ω1(1 + κ 2 1 ) + N1 , G v ni P2 P1 ρ1(1 + κ 2 1 ) + N1 . ( 9
)
The aggregate distortion noise is equal to η1 ∼ CN 0, κ 2 1 P1 for the source and η2 ∼ CN 0,

κ 2 2 G 2 ni Es 1 ,ν 1 {|y1| 2 } for the relay. The latter becomes η f 2 ∼ CN 0, G f ni 2 κ 2 2 (P1 ρ1(1 + κ 2 1 ) + N1) and η v 2 ∼ CN 0, κ 2 2 P2
for fixed and variable gain relays, respectively. After some algebraic manipulations (using the expressions for G v ni ), the end-to-end SNDRs for fixed and variable gain relaying are

γ f ni = P1 ρ1 ρ2 P1 ρ1 ρ2 d + ρ2 N1(1 + κ 2 2 ) + N 2 (G f ni ) 2 , (10)
γ v ni = P1 P2 ρ1 ρ2 P1 P2 ρ1 ρ2 d+ρ1P1N2(1+κ 2 1 )+ρ2P2N1(1+κ 2 2 )+N1N2 , (11) respectively
, where d κ 2 1 + κ 2 2 + κ 2 1 κ 2 2 will appear repeatedly in the paper. These SNDR expressions reduce to (5) for κ1 = κ2 = 0.

OUTAGE PROBABILITY ANALYSIS

This section derives new analytical expressions for the exact and asymptotic OPs with AF under the presence of transceiver impairments, for fixed and variable gain relaying. These results are compared with the standard results in the literature based on the assumption of ideal hardware. The OP, Pout(x), is defined as the probability that the channel fading makes the SNDR fall below a certain threshold x; that is,

Pout(x) = Pr{γ ≤ x}, (12)
where γ is the end-to-end SNDR of the system.

Exact SNR Analysis: Arbitrary Fading Distributions

This subsection considers the OP with arbitrary channel fading. The following lemma plays an important role when deriving Pout(x).

Lemma 1. Let c1, c2, c3 be strictly positive constants and let ρ be a non-negative random variable with cdf Fρ(•). Then,

Pr c1ρ c2ρ + c3 ≤ x = Fρ c 3 x c 1 -c 2 x , 0 ≤ x < c 1 c 2 , 1, x ≥ c 1 c 2 . ( 13
)
Suppose c2 = 0 instead, then (13) simplifies to

Pr c1ρ c3 ≤ x = Fρ c3x c1 . ( 14
)
Proof. The left-hand side of ( 13) is equal to

Pr c1ρ ≤ (c2ρ + c3)x = Pr ρ ≤ c3x (c1 -c2x) (15)
after some algebra. The last expression is exactly Fρ( c 3 x c 1 -c 2 x ). Note that if (c1 -c2x) ≤ 0, then the outage event always occurs.

The effective SNDRs under AF relaying are given in ( 5) and ( 10)- [START_REF] Riihonen | BEP analysis of OFDM relay links with nonlinear power amplifiers[END_REF] for ideal and non-ideal transceiver hardware, respectively. The OPs can be computed using Lemma 1.

Theorem 1. Suppose ρi is an independent non-negative random variable with cdf Fρ i (•) and pdf fρ i (•) for i = 1, 2. For fixed and variable gain relaying, the OPs with ideal hardware are

P id,f out (x) = ∞ 0 Fρ 1 x N1ρ2 + N 2 (G f id ) 2 P1ρ2 fρ 2 (ρ2)dρ2, (16)
P id,v out (x) = 1- ∞ 0 1-Fρ 1 (x+x 2 )+ P 2 x N 2 z P 1 P 2 N 1 N 2 z fρ 2 z + N2x P2 dz, (17)
respectively. The corresponding OPs with transceiver hardware impairments satisfy P ni,f out (x) = P ni,v out (x) = 1 for x ≥ 1 d and

P ni,f out (x) = ∞ 0 Fρ 1 (b1ρ2 + a1)x ρ2(1 -xd) fρ 2 (ρ2)dρ2, (18)
P ni,v out (x) = 1- ∞ 0 1-Fρ 1 x+x 2 1-xd + b 1 xz a 2 (1-xd) a 2 z fρ 2 z + xb2 1-xd dz, (19)
for x < 1 d . In these expressions, we have a1

N 2 P 1 (G f ni ) 2 , a2 N 1 N 2 P 1 P 2 , b1 N 1 (1+κ 2
2 ) P 1 , and b2

N 2 (1+κ 2 1 ) P 2 .
Proof. The OP in [START_REF] Qi | Analysis and compensation of I/Q imbalance in amplify-and-forward cooperative systems[END_REF] depends on the end-to-end SNDR γ, which is a function of both ρ1 and ρ2. Using the law of total probability to condition on ρ2, we can write Pr{γ ≤ x} = ∞ 0 Pr{γ ≤ x|ρ2}fρ 2 (ρ2)dρ2. The conditional probability Pr{γ ≤ x|ρ2} can then be evaluated using Lemma 1, which immediately gives ( 16) and [START_REF] Lau | Downlink scheduling and rate adaptation design of multi-user, multiple-antenna base station with imperfect CSIT[END_REF]. Finally, ( 17) is obtained by the additional change of variable z = ρ2 -N 2 x P 2 , while (19) is achieved by setting z = ρ2 -xb 2 1-xd .

This theorem provides integral expressions for the OPs. In the following, we particularize these expressions to the case of Rayleigh fading. Interestingly, the OP with transceiver impairments always equals 1 when x ≥ 1 d . This has an intuitive explanation since γ ≤ 1 d for any transmit power and fading distribution; see Section 3.3.

Exact SNR Analysis: Rayleigh Fading Channels

Under Rayleigh fading and ideal hardware, the OPs with fixed and variable gain relaying were obtained in [2, Eq. ( 9)] and [14, Eq. ( 14)], respectively. Using our notation, we get

P id,f out (x) = 1 -2e - N 1 x P 1 Ω 1 c f id x K1 2 c f id x , (20)
P id,v out (x) = 1 -2e -x N 1 P 1 Ω 1 + N 2 P 2 Ω 2 x + x 2 c v id K1 2 x + x 2 c v id , (21)
where

c f id N 2 P 1 (G f id ) 2 Ω 1 Ω 2 , c v id P 1 P 2 Ω 1 Ω 2 N 1 N 2
, and K1(•) denotes the first-order modified Bessel function of the second kind.

The non-ideal case with transceiver hardware impairments is more involved, as can be seen from the following corollary.

Corollary 1. Suppose h1, h2 are independent and hi ∼ CN (0, Ωi) for i = 1, 2. The OPs under transceiver hardware impairments are

P ni,f out (x) =    1-2e -B f ni x 1-xd A f ni x 1-xd K1 2 A f ni x 1-xd , x < 1 d , 1, x ≥ 1 d , (22)
P ni,v out (x) =    1-2e -B v ni x 1-xd √ A v ni (x+x 2 ) 1-xd K1 2 √ A v ni (x+x 2 ) 1-xd , x < 1 d , 1,
x ≥ 1 d , (23) for fixed and variable gain relaying, respectively. In these expressions, we have

A f ni N 2 P 1 (G f ni ) 2 Ω 1 Ω 2 , B f ni N 1 (1+κ 2 2 ) Ω 1 P 1 , A v ni N 1 N 2 P 1 P 2 Ω 1 Ω 2 , and B v ni N 1 (1+κ 2 2 ) Ω 1 P 1 + N 2 (1+κ 2 1 ) Ω 2 P 2 .
Proof. The proof consists of evaluating ( 18) and ( 19) using the fading distribution in [START_REF] Hasna | A performance study of dualhop transmissions with fixed gain relays[END_REF]. For fixed gain relaying,

P ni,f out (x) becomes ∞ 0 1 -e - (b 1 ρ 2 +a 1 )x Ω 1 ρ 2 (1-xd) e - ρ 2
Ω 2 dρ2 = 1 -e -B f ni x 1-xd ∞ 0 e - A f ni x ρ 2 (1-xd) e - ρ 2 Ω 2 dρ2 (24)
for x < 1 d and is equal to 1 otherwise. The remaining integral is evaluated using [START_REF] Gradshteyn | Table of Integrals, Series, and Products[END_REF]Eq. (3.324.1)], leading to (22). The expression for P ni,v out (x) is obtained in the same way. This corollary provides new and tractable closed-form expressions for the OPs in case of transceiver impairments.

Asymptotic SNR Analysis

To obtain some insights on the fundamental impact of impairments, we now elaborate on the high-SNR regime. For the ease of presentation, we assume that P1, P2 grow large with P1 = βP2 for some β > 0, such that the relay gain remains finite and strictly positive.

Corollary 2. Suppose P1, P2 grow large with a finite non-zero ratio and consider any random distributions on ρ1, ρ2 that are strictly positive (with probability one). The asymptotic OP with transceiver impairments and either fixed or variable gain relaying satisfies lim

P 1 ,P 2 →∞ Pout(x) = 0, x ≤ 1 κ 2 1 +κ 2 2 +κ 2 1 κ 2 2 , 1, x > 1 κ 2 1 +κ 2 2 +κ 2 1 κ 2 2 . ( 25
)
Proof. Referring back to [START_REF] Samuel | Impact of RF impairments in cellular wireless metropolitan area networks[END_REF] and by taking the limit P1, P2 → ∞ (with P1 = βP2), we can easily see that the instantaneous end-toend SNDR, for fixed gain relaying, converges to lim

P 1 ,P 2 →∞ γ f ni = γ * = 1 d = 1 κ 2 1 + κ 2 2 + κ 2 1 κ 2 2 . (26)
Therefore, (25) is trivially obtained. The proof for the case of variable gain relaying follows a similar line of reasoning.

A number of conclusions can be drawn from Corollary 2. First, an SNDR ceiling effect appears in the high-SNR regime, which significantly limits the performance of AF relaying systems. This phenomenon is fundamentally different from the ideal hardware case, in which an increasing SNR makes the end-to-end SNDR grow without bound and Pout(x) → 0 for any x. Note that this ceiling effect is independent of the fading distributions of the two hops.

Moreover, the SNDR ceiling γ * in ( 26) is inversely proportional to the squares of κ1, κ2. We also note that the impact of κ1 and κ2 on the value of γ * is symmetric. Recall that these parameters describe the level of transceiver impairments and equal the EVMs of the source and relay hardware, respectively. This validates that transceiver impairments dramatically affect performance of AF relaying, and, hence, they should be taken into account when designing and evaluating relaying networks. To this end, we now assume that we have a fixed total EVM constraint κtot, such that κ1 +κ2 = κtot, and present the following insightful corollary. Corollary 3. Suppose κ1+κ2 = κtot for some constant 0 < κtot < 2. The SNDR ceiling in (26) is maximized if κ1 = κ2 = κtot 2 . Proof. The maximization of (26) is equivalent to a minimization of f (κ1) = κ 2 1 +(κtot -κ1) 2 +κ 2 1 (κtot -κ1) 2 . Simple differentiation reveals that the function is convex for any 0 < κtot < 2 and that κ1 = κtot 2 is the only real-valued solution to ∂ ∂κ 1 f (κ1) = 0. Corollary 3 shows that it is better to have the same hardware quality on the two hops of the relay channel, than having one highquality node and one low-quality node. We finally point out that the SNDR ceiling can be trivially bounded as

1 (κ 1 +κ 2 ) 2 ≤ γ * ≤ 1 κ 2 1 +κ 2 2 .

NUMERICAL ILLUSTRATIONS

This section provides numerical illustrations of our analytical results. We consider a symmetric Rayleigh fading scenario where the average SNR, P i Ω i N i , is the same for i = 1, 2. This enables us to focus on the impact of transceiver hardware impairments.

The OP, Pout(x), with fixed and variable gain relaying is shown in Fig. 2 as a function of the average SNR. We compare the system performance with ideal transceivers and non-ideal transceivers with κ1 = κ2 = 0.1. Two thresholds are considered: x = 2 2 -1 = 3 and x = 2 5 -1 = 31. These correspond to spectral efficiencies of 2 and 5 bits/s/Hz, respectively. The curves were generated by the analytical expressions in (20)-(23), while the marker symbols show the results of Monte-Carlo simulations. As seen in Fig. 2, there is a minor performance loss caused by transceiver impairments in the low threshold case of x = 3. However, there is a substantial performance loss of around 5 dB in SNR when the threshold is increased to x = 31. This effect applies to both variable and fixed gain relaying. Note that variable gain relaying outperforms the fixed one in most scenarios of interest, which is in line with [START_REF] Hasna | A performance study of dualhop transmissions with fixed gain relays[END_REF]. Next, we consider the average rate of data correctly received under fading (in bits/s/Hz) [17, Chapter 4.2.3], or the so-called goodput [START_REF] Lau | Downlink scheduling and rate adaptation design of multi-user, multiple-antenna base station with imperfect CSIT[END_REF]. This is defined as the product between the gross spectral efficiency and the corresponding probability of successful communication. For any given setup, the maximum goodput is defined as is the SNDR ceiling. This is consistent with Corollary 2, because γ * is the highest SNDR for which the OP approaches zero as P1, P2 → ∞. The saturation level increases when κ1, κ2 are decreased. Finally, Fig. 4 shows the value of the SNDR ceiling for different sets of κ1, κ2 for which κ1 + κ2 = 0.15. As stated in Corollary 3, the maximum is achieved by κi = 0.15 2 for i = 1, 2.

Rgoodput = max

CONCLUSIONS

Physical transceiver hardware introduces impairments that distort the emitted signals. Despite the practical importance of these impairments, little was known about their impact on the achievable performance of relaying systems. This paper has analytically shown that the performance of dual-hop AF relaying is highly affected by these impairments, particularly when high spectral efficiencies are required. Closed-form expressions for the exact and asymptotic OPs were derived for fixed and variable gain relaying. These expressions effectively characterize the impact of impairments and manifest the existence of an SNDR ceiling that cannot be crossed by increasing transmit powers or improving fading conditions. The engineering insight is that both hops should have the same hardware quality.

Fig. 1 :

 1 Fig. 1: Block diagram of AF relaying with (a) ideal hardware or (b) transceiver impairments with aggregate distortion noise ηi, i = 1, 2.

Fig. 2 :

 2 Fig. 2: Outage probability Pout(x) with AF relaying with ideal and non-ideal transceiver hardware.

2 Fig. 3 : 2 Fig. 4 :

 2324 Fig.3: Maximum goodput Rgoodput(x) with AF relaying with ideal and non-ideal transceiver hardware.

x≥0 log 2

 2 (1 + x) (1 -Pout(x)) . (27)The analytical expressions for Pout(x) in (20)-(23) allow simple numerical computation of the maximum goodput in (27). This is illustrated in Fig.3for different levels of impairments: κ1 = κ2 ∈ {0.05, 0.1, 0.15}. This figure demonstrates the deleterious impact of non-ideal hardware; the goodput with transceiver impairments saturates at high SNR and approaches log 2 (1 + γ * ), where γ * =

Note that the following analysis is generic and applies for any fading distribution. We focus on the case of Rayleigh fading for the sake of tractability.

[START_REF] Hasna | A performance study of dualhop transmissions with fixed gain relays[END_REF] The relay then has a long-term power constraint P 2 = E{|G f id y 1 | 2 } where expectation is taken over signal, noise and channel fading

realizations.[START_REF] Hua | Theories and methods for advanced wireless relays-Issue I[END_REF] The Gaussianity holds in particular for the residual distortion when compensation algorithms are applied to handle multiplicative signal errors[START_REF] Studer | MIMO transmission with residual transmit-RF impairments[END_REF].

The total EVM is defined as Eη i {|η i | 2 }/Es i {|s i | 2 }. 3GPP LTE has EVM requirements for transmitters in the range κ tx i ∈ [0.08, 0.175], where smaller values are needed to support high spectral efficiency [15, Sec. 14.3.4].