Ümit V C ¸atalyürek

Kamer Kaya

On shared-memory parallelization of a sparse matrix scaling algorithm

Keywords: Shared-memory parallelization, sparse matrices, hypergraphs, matrix scaling

We discuss efficient shared memory parallelization of sparse matrix computations whose main traits resemble to those of the sparse matrix-vector multiply operation. Such computations are difficult to parallelize because of the relatively small computational granularity characterized by small number of operations per each data access. Our main application is a sparse matrix scaling algorithm which is more memory bound than the sparse matrix vector multiplication operation. We take the application and parallelize it using the standard OpenMP programming principles. Apart from the common race condition avoiding constructs, we do not reorganize the algorithm. Rather, we identify associated performance metrics and describe models to optimize them. By using these models, we implement parallel matrix scaling algorithms for two wellknown sparse matrix storage formats. Experimental results show that simple parallelization attempts which leave data/work partitioning to the runtime scheduler can suffer from the overhead of avoiding race conditions especially when the number of threads increases. The proposed algorithms perform better than these algorithms by optimizing the identified performance metrics and reducing the overhead.

I. INTRODUCTION

Sparse matrix computations working exclusively on nonzero entries are difficult to parallelize. The difficulty arises from the fact that the amount of computation per nonzero access is very small. Among the most well-known computations are sparse matrix-dense vector multiplies (Sp-MxV) of the form y Ax. In SpMxV, each nonzero is accessed, multiplied with an x-vector entry, and the result is added to a y-vector entry only once. That is, there are two reads and one read-and-write per each nonzero accessed. In this work, we parallelize an iterative matrix scaling algorithm where the computational core of each iteration is an operation that performs three reads and two read-and-writes per each nonzero accessed. In this light, the computations that we are targeting in this paper are harder to efficiently parallelize than the SpMxV operation. Our aim is to identify metrics to optimize parallel execution of the scaling algorithm and describe models for this purpose.

We assume that the algorithms will be parallelized using the standard OpenMP parallelization techniques, without a huge effort in reorganizing the algorithms. In particular, we take the point of view of a programmer that adopts loop-level parallelism and the so-called low-level, singleprogram multiple-data parallelism [1, p.32]. In both cases, to avoid race conditions, the programmer uses the nowstandard technique of OpenMP directives, locks, atomic instructions, and/or large amount of private memory for each thread. These techniques incur overhead as additional work, a burden for CPU, or private read/write memory, a burden for cache efficiency. We propose models and methods to reduce the private memory usage, number of locks, atomic operations, and extra parallel work that is incurred. In all of these cases, we use hypergraphs and relate the objective of reducing the overhead to a suitable cost function in hypergraph partitioning problem.

Hypergraph partitioning-based models have been widely used to efficiently parallelize SpMxV and similar operations for distributed memory processor systems [START_REF] Bisseling | Mondriaan sparse matrix partitioning for attacking cryptosystems by a parallel block Lanczos algorithm-A case study[END_REF], [START_REF] Aykanat | Hypergraph-partitioning based decomposition for parallel sparse-matrix vector multiplication[END_REF], [START_REF] Aykanat | On twodimensional sparse matrix partitioning: Models, methods, and a recipe[END_REF], [START_REF] Uc | Hypergraph partitioning[END_REF], [START_REF] Uc | Revisiting hypergraph models for sparse matrix partitioning[END_REF], [START_REF] Vastenhouw | A two-dimensional data distribution method for parallel sparse matrix-vector multiplication[END_REF]. These models identify the load balance and the total volume of communication as two relevant metrics for efficient parallelization. Although there are other communication cost metrics (such as the total number of messages or the maximum number or volume of messages per processors), the total volume of messages is deemed to be important to minimize first. Other communication cost metrics can be addressed during or after optimizing the partitioning for the total communication volume.

As should be evident, the sparse matrix storage format and the performance metrics are related. We focus on sparse matrix algorithms that use two of the most common formats. The coordinate format (COO) keeps an array of Z triplets of the form ha ij , i, ji for a sparse matrix A with Z entries. Each triplet contains a nonzero entry a ij and its row and column indices (i, j). In computations that use COO, it is therefore natural to assign works to threads based on triplets. The compressed row storage format (CRS) uses three arrays to store a m⇥n sparse matrix A with Z nonzeros. One array of size Z keeps the values of nonzeros where the nonzeros in a row are stored consecutively. Another array parallel to the first one keeps the column index of each nonzero. The third array keeps the starting index of the nonzeros at a given row where the ending index of the nonzeros at a row is one less than the starting index of the next row. In computations that use CRS, it is natural to assign works to threads by rows. We identify metrics and models related to the optimization of algorithms that uses those formats.

We view the scaling algorithm that we parallelize as a model parallelization problem posing somewhat more challenges than the SpMxV operations. Nonetheless, most of the methods proposed for efficient parallelization of the SpMxV operations should be applicable to our problem. Some of the existing studies [START_REF] Buluc | Parallel sparse matrix-vector and matrixtranspose-vector multiplication using compressed sparse blocks[END_REF], [START_REF] Shantharam | Exploiting dense substructures for fast sparse matrix vector multiplication[END_REF], [START_REF] Williams | Optimization of sparse matrix-vector multiplication on emerging multicore platforms[END_REF] reorganize the data, manipulate the data structures, and update the algorithms for parallel settings. It should be possible to apply our methods to at least some of these improved algorithms. This will require careful analyze of many different algorithms and surpass the limits and aims of this paper. Some other studies [START_REF] Kotakemori | Performance evaluation of parallel sparse matrix-vector products on SGI Altix3700[END_REF], [START_REF] Williams | Optimization of sparse matrix-vector multiplication on emerging multicore platforms[END_REF], [START_REF] Wu | Efficient PageRank and SpMV computation on AMD GPUs[END_REF] apply architecture-aware techniques to greatly improve the performance of parallelized algorithms. It seems harder to add our techniques on top of these techniques. Yet some other studies [START_REF] Yzelman | Cache-oblivious sparse matrix-vector multiplication by using sparse matrix partitioning methods[END_REF], [START_REF]Two-dimensional cache-oblivious sparse matrix-vector multiplication[END_REF] improve the performance of the serial SpMxV computations by reordering the matrix and optimizing the cache performance. Our methods can be combined with these methods by reordering the submatrices assigned to each thread.

When compared with hypergraph partitioning, the execution time of matrix scaling is small. However, if one has several matrices with the same nonzero pattern and different nonzeros, which is the case for several applications, the reduction on the scaling time may compensate the cost for preprocessing. Furthermore, we chose scaling only as a model application. As the improvements proposed for SpMxV in the literature can be used for scaling, the models and techniques proposed in this paper can also be useful for SpMxV and related sparse-matrix problems. We refer the reader to [START_REF] Uc | A library for parallel sparse matrixvector multiplies[END_REF], [START_REF] Akbudak | Hypergraphpartitioning-based models and methods for exploiting cache locality in sparse-matrix vector multiplication[END_REF] for a comparison of the time spent for hypergraph partitioning with respect to the SpMxV time.

The organization of the paper is as follows. In Section II, we present the scaling algorithm and background material on hypergraph partitioning. In Section III, we present the shared-memory implementation of the scaling algorithm that uses the CRS and COO storage formats. We discuss two different standard OpenMP parallelization approaches and relate their performance metrics to well-known partitioning objectives in the literature. We evaluate the parallel algorithms experimentally in Section IV in order to show that the theoretical improvements due to the discussed models lead to improved running times in practice. We conclude the paper with a discussion in Section V.

II. BACKGROUND

We present the scaling algorithm below. We then describe the hypergraph partitioning problem and three common objective functions. Hypergraph models to represent sparse matrices used in this paper are also presented.

A. The scaling algorithm

Matrix scaling algorithms have been used to improve the numerical stability in solving linear systems [START_REF] Duff | Direct Methods for Sparse Matrices[END_REF]. Scaling a matrix A consists in pre-and postmultiplying A by two diagonal matrices, say D r and D c , yielding a scaled matrix b A = D r AD c . We describe a scaling algorithm (originally proposed by Ruiz [START_REF] Ruiz | A scaling algorithm to equilibrate both rows and columns norms in matrices[END_REF]; see also [START_REF] Ruiz | A symmetry preserving algorithm for matrix scaling[END_REF]) which finds scaling matrices D r and D c such that the rows and columns of b A = D r AD c has the same length in some norm. Common choices of the norm include 1-norm and 1-norm scaling. We use kxk 1 to denote the 1-norm of the vector x where kxk

1 = P |x i | (similarly kxk 1 is the 1-norm of x where kxk 1 = max |x i |).
Experiments with a sparse solver suggest the application of this algorithm both with 1-and 1-norms [START_REF] Ruiz | A symmetry preserving algorithm for matrix scaling[END_REF], [START_REF] Amestoy | MUMPS[END_REF]. Here we focus on the 1-norm scaling, as the algorithm in 1norm converges very fast [START_REF] Ruiz | A symmetry preserving algorithm for matrix scaling[END_REF]. Since the 1-norm scaling is possible only for square matrices, we describe the algorithm with 1-norm on an n ⇥ n, real, sparse matrix A.

For i = 1, . . . , n, let r i = a T i⇤ 2 R n⇥1 denote the ith row of A as a vector, and for j = 1, . . . , n, let c j = a ⇤j 2 R n⇥1 denote the jth column of A. Let diag(•) be an operator that takes a vector of size n and returns an n⇥n diagonal matrix having the entries of the vector in the main diagonal, i.e., for D = diag(x) it holds that d ii = x i and d ij = 0 for 1  i 6 = j  n. With these definitions, the scaling algorithm is formulated as shown in Algorithm 1.

✓ q kr

(k)

i k1 ◆ i=1,...,n 2 D2 diag ✓ q kc (k) j k1 ◆ j=1,...,n b A (k+1) D 1 1 b A (k) D 1 2 D (k+1) r D (k) r D 1 1 D (k+1) c D (k) c D 1 2
In this algorithm, b A (k) denotes the scaled matrix at the kth iteration whose ith row and jth column are represented, respectively, by r

d (k) r (i) ⇥ a ij ⇥ d (k) c (j) while computing r (k) i
and c

(k) j . The diagonal matrices D 1 and D 2 are created anew at each iteration (using the square root of the norms of the rows and columns, respectively). The output, scaling matrices D r and D c , are updated at the end of each iteration, where those resulting from the last iteration are returned. The convergence is obtained if

max i n |1 kr (k) i k 1 | o  " and max j n |1 kc (k) j k 1 |
o  " for a given value of " > 0. The computational core of the algorithm is the lines 1 and 2. In a proper sequential implementation of this algorithm, these two steps should be carried out as follows:

v d (k) r (i) ⇥ a ij ⇥ d (k) c (j) and d 1 (i) d 1 (i) + v and d 2 (j) d 2 (j) + v.
In other words, one accesses a nonzero entry, retrieves the components of the D r and D c (three reads), performs the scalar multiplications, and then adds the result to the components of D 1 and D 2 (two read-and-writes). This high number of memory accesses per nonzero renders the overall algorithm hard to efficiently parallelize. Furthermore, our preliminary experiments verify that handling read/write conflicts during parallelization is crucial since the scaling algorithm does not self-stabilize in case of a conflict.

A distributed memory parallelization of this algorithm using MPI is reported by Amestoy et al. [START_REF] Amestoy | A parallel matrix scaling algorithm[END_REF]. In that study, it has been shown that the communication volume requirements of such a parallel implementation is related to the communication volume requirements of a sparse matrixvector multiply operation.

B. Hypergraphs and hypergraph models for sparse matrices

A hypergraph H = (V, N ) is defined as a set of vertices V and a set of nets (hyperedges) N among those vertices. A net n 2 N is a subset of vertices, and the vertices in n are called its pins. The size of a net is the number of its pins, and the degree of a vertex is equal to the number of nets that contain it. Graph is a special instance of hypergraph such that each net has size two. Vertices can be associated with weights, denoted with w[•], and nets can be associated with costs, denoted with c[•].

A K-way partition of a hypergraph H is denoted as

⇧ = {V 1 , V 2 , . . . , V K } where • parts are pairwise disjoint, i.e., V k \ V `= ; for all 1  k < ` K, • each part V k is a nonempty subset of V, i.e., V k ✓ V and V k 6 = ; for 1  k  K, • union of K parts is equal to V, i.e., S K k=1 V k = V.
In a partition ⇧, a net that has at least one pin (vertex) in a part is said to connect that part. The number of parts connected by a net n, i.e., connectivity, is denoted as n . A net n is said to be uncut (internal) if it connects exactly one part, and cut (external), otherwise (i.e., n > 1).

Let W k denote the total weight in V k (i.e.,

W k = P v2V k w[v]
) and W avg denote the weight of each part when the total vertex weight is equally distributed (i.e.,

W avg = ( P v2V w[v])/K). If each part V k 2 ⇧ satisfies the balance criterion W k  W avg (1 + "), for k = 1, 2, . . . , K (1)
we say that ⇧ is balanced where " represents the maximum allowed imbalance ratio.

The set of external nets of a partition ⇧ is denoted as N E . Let (⇧) denote the cost, i.e., cutsize, of a partition ⇧. There are various cutsize definitions [START_REF] Lengauer | Combinatorial Algorithms for Integrated Circuit Layout[END_REF] such as:

cut (⇧) = X n2N E c[n] , (2)
conn (⇧) = X n2N c[n]( n 1) , (3)
SOED (⇧) = X n2N E c[n] n . (4)
The cutsize metric given in (2) will be referred to here as cut-net metric, the one in (3) will be referred as connectivity-1 metric, and the one in (4) will be referred to as the SOED metric (widely used in the VLSI domain [23, p.10], [START_REF] Karypis | Multilevel k-way hypergraph partitioning[END_REF], and recently found applications in the scientific computing domain [START_REF] Yamazaki | Combinatorial problems in a parallel hybrid linear solver[END_REF]). Given " and an integer K > 1, the hypergraph partitioning problem can be defined as the task of finding a balanced partition ⇧ with K parts such that (⇧) is minimized. The hypergraph partitioning problem is NPhard [START_REF] Lengauer | Combinatorial Algorithms for Integrated Circuit Layout[END_REF] with any of the above objective functions. We used a state-of-the-art partitioning tool PaToH [START_REF] Aykanat | PaToH: A Multilevel Hypergraph Partitioning Tool[END_REF] which already has options to perform hypergraph partitioning in order to optimize cut-net and connectivity-1 metrics (while of course achieving balance). We have instrumented PaToH with the technique proposed by Yamazaki et al. [START_REF] Yamazaki | Combinatorial problems in a parallel hybrid linear solver[END_REF] to address the SOED metric (4).

There are three well-known hypergraph models for sparse matrices. These are the column-net [START_REF] Aykanat | Hypergraph-partitioning based decomposition for parallel sparse-matrix vector multiplication[END_REF], row-net [START_REF] Aykanat | Hypergraph-partitioning based decomposition for parallel sparse-matrix vector multiplication[END_REF], and finegrain models [START_REF]A fine-grain hypergraph model for 2D decomposition of sparse matrices[END_REF]. We describe these models below for a sparse matrix A of size m ⇥ n with Z nonzeros.

In the column-net model, A is represented as a unit-cost hypergraph H R = (V R , N C ) with |V R | = m vertices, |N C | =
n nets, and Z pins. In H R , there exists one vertex v i 2 V R for each row i. Weight w[v i ] of a vertex v i is equal to the number of nonzeros in row i. There exists one unit-cost net n j 2 N C for each column j. Net n j connects the vertices corresponding to the rows that have a nonzero in column j. That is, v i 2 n j if and only if a ij 6 = 0. The row-net model is the column-net model of the transpose of A.

In the fine-grain model, A is represented as a unit-weight and unit-cost hypergraph

H Z = (V Z , N RC ) with |V Z | = Z vertices, |N RC | =
m+n nets and 2Z pins. In V Z , there exists one unit-weight vertex v ij for each nonzero a ij . In N RC , there exist one unit-cost row-net r i for each row i and one unit-cost column-net c j for each column j. The row-net r i connects the vertices corresponding to the nonzeros in row i, and the column-net c j connects the vertices corresponding to the nonzeros in column j. That is, v ij 2 r i and v ij 2 c j if and only if a ij 6 = 0.

III. MULTITHREADED ALGORITHMS FOR MATRIX

SCALING

The multithreaded implementations of Algorithm 1 differ in the way they store the matrix A and hence, in the way they compute the row and column norms. We use d 1 [i] and d 2 [j] to store the norm of row i and column j in the algorithms (instead of the notation D 1 and D 2 ) and call d 1 and d 2 as vectors; similarly d r and d c are used (instead of D r , D c ).

As mentioned above, we do not update the entries of A. Instead, for each nonzero a ij , we compute the corresponding scaled nonzero In all the algorithms that we present, we distinguish three phases: init, put, and get. We will focus on these phases in which the algorithms compute the norms of the rows and columns of the iteration matrix b A (k) . The implementation of other parts of the algorithms is almost the same-they are parallelized with OpenMP parallel-for directives. In the init phase, each thread t initializes its own private memory. In the put phase, the nonzeros of the matrix are processed. After reading a ij from the memory and computing the corresponding scaled nonzero val, a thread t adds val to either shared

d r [i] ⇥ a ij ⇥ d c [j]. At each iteration, each scaled nonzero d r [i] ⇥ a ij ⇥ d c [j]
d 1 [i] or private d t 1 [i] (and similarly to d 2 [j] or d t 2 [j]
). In the get phase, the values in private arrays are transferred from the private arrays to the shared array.

A naive approach which is suitable for the CRS format is given in Algorithm 2. While processing the nonzeros in Algorithm 2, a parallelization in the row-level is promising since the CRS format is used to store A. That is, all the nonzeros in the ith row of A is processed by the same thread t which is the only thread writing to d 1 [i]. Thus there is no conflict and hence no need for a private memory for these writes. However, the accesses to d 2 need to be synchronized carefully. In Algorithm 2, each thread t uses its private d t 2 array to store the partial norms of the columns and in get, these partial norms are added to the shared d 2 array.

Algorithm 3 parallelizes the scaling process in nonzerolevel when the matrix is stored in COO format. With this approach, in addition to the column norms, a row norm d 1 [i] can also be updated by different threads at the same time. Hence, Algorithm 3 uses a private d t 1 array for each thread t to store the partial row norms. Besides, a thread cannot pack a row norm in a private variable sum t while processing the nonzeros of a row. Similar to Algorithm 2, the partial row and column rows are computed in put phase and added back to the shared norm arrays in the get phase.

There are two overheads of conflict resolution by using Algorithm 2: Simple parallel scaling with CRS Input: A: n ⇥ n input matrix in CRS format Output: dr, dc: row and column scaling vectors

for i = 1 to n in parallel do dr[i] 1 dc[i] 1 while not converged do for i = 1 to n in parallel do d1[i] 0 d2[i] 0 init for t = 1 to ⌧ in parallel do for i = 1 to n do d t 2 [i] 0 put for i = 1 to n in parallel do It is the current thread id sum t 0 for each nonzero aij in row i do val dr[i] ⇥ aij ⇥ dc[j] add val to sum t and d t 2 [j] d1[i] sum t get for t = 1 to ⌧ do for i = 1 to n in parallel do d2[i] d2[i] + d t 2 [i] error max (maxi(|1 d1[i]|), maxi(|1 d2[i]|)) if error < " then converged true else for i = 1 to n in parallel do dr[i] dr[i]/ p d1[i] dc[i] dc[i]/ p d2[i]
Algorithm 3: Simple parallel scaling with COO Input: A: n ⇥ n input matrix in COO format Output: dr, dc: row and column scaling vectors

• • • while not converged do • • • init for t = 1 to ⌧ in parallel do for i = 1 to n do d t 1 [i] 0 d t 2 [i] 0 put for each nonzero aij in parallel do It is the current thread id val dr[i] ⇥ aij ⇥ dc[j] add val to d t 1 [i] and d t 2 [j] get for t = 1 to ⌧ do for i = 1 to n in parallel do d1[i] d1[i] + d t 1 [i] for i = 1 to n in parallel do d2[i] d2[i] + d t 2 [i] • • •
private arrays: computation and memory. For the CRS based parallel algorithm, the total computation overhead due to init and get phases are 2⌧ n at each iteration. For the COO based algorithm, this overhead is 4⌧ n. The total private memory used by these algorithms are ⌧ n and 2⌧ n, respectively.

A. Conflict resolution with partitioning for parallelization in the row-level

In Algorithm 2, each thread is given a set of rows during the execution of the scaling process. Hence, the rows are implicitly and dynamically partitioned among the threads by the scheduler. The decisions are given with respect to a scheduling policy without considering any metric or restriction other than load balancing. Here, we show that the memory and computation overhead can be reduced by modeling the problem as a hypergraph partitioning problem and using a static partition.

Let ⇧ = {R 1 , R 2 , . . . , R ⌧ } be a partition of the row set {1, 2, . . . , n} in the column-net hypergraph model. Assume that thread t will process the nonzeros of the rows in R t . Hence, a conflict will occur only if there exists a column j such that a ij and a i 0 j are nonzero and rows i and i 0 are in different parts. The number of such columns is equal to the number of external nets with respect to ⇧, i.e., cut (⇧). By using this fact, one can greatly reduce the size of private d 2 arrays by partitioning the hypergraph with the metric cut (⇧) as in [START_REF] Bisseling | Mondriaan sparse matrix partitioning for attacking cryptosystems by a parallel block Lanczos algorithm-A case study[END_REF]. Given a partition ⇧, we first permute the rows and columns of the matrix with respect to ⇧. The rows are permuted in the partial order of (R 1 R 2 • • • R ⌧ ). The external columns are put before the other columns, for the rest of the columns, column j is put before column j 0 if j is internal to t, j 0 is internal to t 0 and t < t 0 .

Let us call a nonzero a ij external if column j is in the cut. Otherwise, we say that a ij is internal. We also call a row i internal if all a i⇤ s are internal. Otherwise, we call that row external. Let r ext be the number of these external rows. Algorithm 4, CRS-Cut, parallelizes the scaling process in the row-level with static partitioning. To partition the rows among threads, we use ⌧ pointers to keep the id of the first row of each part. The init and get phases are similar to the simple approach described in Algorithm 4. Note that the private array sizes are cut (⇧) instead of n. Hence, the total private memory used is ⌧ cut (⇧) and the total computation overhead is 2⌧ cut (⇧). To access to the internal/external rows in an efficient way, we separate them in the row lists of the parts. We then use an additional pointer set of size ⌧ to keep the location of the first internal row for each part. Again for further efficiency, we partition the nonzero list of each external row into two by separating the internal and external nonzeros. We keep a set of pointers of size r ext to store the location of the first internal nonzero in each list.

In the init and get phases of Algorithm

d1[i] sum t get for t = 1 to ⌧ do for i = 1 to cut in parallel do d2[i] d2[i] + d t 2 [i] • • • that i 2 R t .
In hypergraph partitioning, this corresponds to the connectivity of net j with part t. For each column j, the number of nonzero d t 2 [j]s is equal to its connectivity j . Hence, at each iteration, the total number of nonzeros in private arrays is P n j=1 j = SOED (⇧) as defined in 4. By using the part-to-net connectivity information while traversing the private arrays in init and put phases, we propose a modified version of CRS-Cut as shown in Algorithm 5.

Algorithm 5, CRS-SOED, uses the same amount of private memory as CRS-Cut to store the partial column norms. To reduce the computation overhead, CRS-SOED also keeps the connectivity information for each part. The size of this additional information is ⌧ + SOED (⇧). Note that when ⌧ = 2, whole private memory will be traversed by CRS-SOED since each net in the cut is connected to both parts. Hence, considering its memory overhead over CRS-Cut, using CRS-SOED can be advantageous when ⌧ > 2.

B. Conflict resolution with partitioning for parallelization in the nonzero-level

When A is stored in COO format, a parallelization in the nonzero-level is more suitable as shown in Algorithm 3. In the simple approach, the scheduler assigns a set of nonzeros

[i] d2[i] + d t 2 [i] • • •
to each thread during the execution of the scaling process. This is an implicit partitioning of the nonzeros among the threads where the assignments are done with respect to a scheduling policy by considering the load balance. Similar to the CRS-based versions, we show that COO-based storage can benefit from a static hypergraph partitioning but this time with the fine-grain model (see Section II-B).

Let ⇧ = {Z 1 , Z 2 , . . . , Z ⌧ } be a partition of the nonzeros of the matrix A. Let us call a nonzero a ij row/col-external if both of its nets, which correspond to ith row and jth column, are external. If only row i (or column j) is external we call a ij row-external (or col-external). Otherwise, we call a ij internal. Algorithm 6, COO-SOED, uses a partition to distribute the nonzeros among threads when the matrix is in the COO format. Similar to the previous algorithms, for each external nonzero, the partial row and column norms are written to the private arrays. Otherwise, if the nonzero is internal, shared norm arrays d 1 and d 2 are used. Hence, the number of private memory accesses in the init and get phases is equal to SOED (⇧). To distribute the work, given ⇧, we permute the nonzeros as in

(Z 1 Z 2 • • • Z ⌧ )
and use ⌧ pointers to keep the location of the first nonzero of each part. For efficiency, we partition each Z i into four to separate row/col-external, row-external, col-external, and internal nonzeros. For efficient access to these nonzero classes, we keep 3 additional pointers for each part to store the location of the first nonzero of each class.

C. Other approaches for shared-memory parallelization of scaling

Using private arrays is not the only way for conflict-free parallelization of the scaling process. A simple alternative is using locks and/or atomic operations for concurrency. The overheads of these approaches such as the number of locks can also be related with a partitioning metric. Consider a matrix stored in CRS format. Given a partition ⇧, the number of external nonzeros is equal to SOED (⇧) given

[i] d1[i] + d t 1 [i] for each external col. i conn. to Zt in parallel do d2[i] d2[i] + d t 2 [i] • • •
in 2 where the cost of each net is equal to the number of nonzeros in the corresponding column. Note that the same is also true for the COO format and parallelization in the nonzero-level. Hence, if one uses an atomic operation for each external nonzero the number of atomic operations is equal to SOED (⇧).

For a further reduction on the number of nonzeros in private arrays, the algorithms can be modified as follows. We assumed that an external net is external to all parts. However, such nets can be assumed to be internal to only one part which permits only one thread to update the corresponding shared norm value in the put phase. Hence, there will be no conflicts in this phase. When this implementation is used the number of private memory accesses init and get phases will be equal to conn (⇧) given in [START_REF] Aykanat | Hypergraph-partitioning based decomposition for parallel sparse-matrix vector multiplication[END_REF]. Although this yields to a reduction on the computational overhead for parallelization, it may have an undesirable effect. For example, in the current implementation of CRS-SOED and COO-SOED, the threads do not need to synchronize between put and get phases since the sets of shared memory locations accessed in these phases are disjoint. Hence, one can use the nowait directive of OpenMP in the implementation to avoid a thread synchronization after the put phase. However, if the conn (⇧) metric and the related implementation is used, a synchronization is necessary since the sets of memory locations will not be disjoint anymore.

IV. EXPERIMENTAL RESULTS

We used two different architectures to test the algorithms. The first one has a 2.27GHz dual quad-core Intel Xeon (Bloomfield) CPU and 48GB main memory. Each core in a socket has 32KB L1 and 256KB L2 caches. Each socket has an 8MB L3 cache shared by 4 cores. The second architecture has a dual quad-core AMD Opteron (Shanghai) processor and 32GB main memory. Each core in a socket has 64KB L1 and 512KB L2 caches. Each socket has a 6MB L3 cache shared by 4 cores. All of the algorithms are implemented in C and OpenMP. The compiler is icc version 12.0 and 11.1 for the first and the second architecture, respectively, and -O3 optimization flag is used. None of our optimization methods is architecture specific. We therefore expect similar behavior from the algorithms on the two architectures.

We used a set of real life matrices from different application domains that are available at the University of Florida (UFL) Sparse Matrix Collection (http://www.cise.ufl.edu/research/sparse/matrices). The properties of the matrices used in the experiments are given in Table I. In the experiments, we do not compute the error on the norm values of the scaled matrix after each iteration. Instead, we iterate the process for 100 times to measure the execution time. To generate the partitions used by the algorithms, we used PaToH [START_REF] Aykanat | PaToH: A Multilevel Hypergraph Partitioning Tool[END_REF] with the appropriate objective for each algorithm and imbalance ratio 0.05.

We tested the validity of our motivation for minimizing the objectives cut and SOED as follows. For each ⌧ , we partitioned the rows/nonzeros of the matrices with a loadbalancing heuristic. To do this, by mimicking the longest processing time first rule of Graham [START_REF] Graham | Bounds on multiprocessing timing anomalies[END_REF], we first sort the rows in decreasing order by number of their nonzeros. We then visit each row in this order and assign it to the part with the least nonzeros assigned so far. For CRS based algorithms, we measured that the imbalances are always smaller than 10 5 for all of the matrices and ⌧ values except when trans5 is partitioned to 8 threads. For COO based algorithms, it is clear that the heuristic partitions the nonzeros in the best way possible. Hence, the partitions we obtain provide almost perfect load balancing among the threads. After running a scaling algorithm with this partition and measuring the execution time, we do a normalization by dividing it to the time spent for scaling the same matrix with the partition generated by PaToH with the appropriate metric. The comparisons among the running time would therefore show that the difference in the practical execution times is related to the optimized cutsize metrics. Figure 1 shows the averages of these normalized values.

Algorithm 7: COO-SOED-Atom: get for t = 1 to ⌧ in parallel do for each external row i conn. to Zt do (atomic)

d1[i] d1[i] + d t 1 [i] for each external column i conn. to Zt do (atomic) d2[i] d2[i] + d t 2 [i]
As Figure 1 shows, the algorithms perform bad when they use the balanced-only partitions instead of the ones generated by PaToH. For CRS-Cut, CRS-SOED and COO-SOED, the ratio of the performances are 1.64, 2.13, and 1.31 for 8 threads. These ratios show that minimizing the cut is an effective approach for sparse matrix computations in sharedmemory architectures. For this experiment, we also implemented an experimental variant of COO-SOED by using atomic operations. This variant, COO-SOED-Atom, differs from COO-SOED only in the get phase. The implementation of this phase is given in Algorithm 7. As expected, COO-SOED-Atom is slower than COO-SOED in practice since atomic operations are costlier. In Figure 1, it can be seen that the normalized performance of COO-SOED-Atom with balanced partitions is worse (9%, 9%, and 6% more for 2, 4, and 8 threads, respectively) than those of COO-SOED with balanced partition. Since, the rest of the algorithms is exactly the same, we can argue that when PaToH partition is used instead of balanced partition, the reduction on COO-SOED-Atom's atomic operations improves the runtime more than COO-SOED. Note that for a partition ⇧, the number of atomic operations is equal to SOED (⇧). Thus, Figure 1 shows that minimizing the cut-size is helpful in practice.

Figures 2 and3 show the average speedup values for the Intel-based and AMD-based architectures, respectively. In the figures, CRS and COO correspond to OpenMP versions of Algorithms 2 and 3 with guided scheduling policy and chunksize 50. These two algorithms do not actually scale after 4 threads. When ⌧ is increased to 8, the execution time of these algorithms also increase. This can be expected especially for matrices with a few number of nonzeros per row/column. Note that in simple CRS-and COO-based algorithms, each additional thread decreases the number of nonzeros per thread but also requires an additional private memory of size ⇥(n). In addition to the computation overhead due to init and get phases, the pressure on the cache will increase in the put phase. As Table II shows, the execution time of COO-Simple is decreased only for three matrices when ⌧ is increased from 4 to 8. These matrices, cage14, Chebyshev4, and torso1 have 18.02, 78.94, and 73.32 nonzeros per row/column on the average (they are the top three matrices with respect to this criteria). For the remaining 9 matrices, we compute the increase in the execution time and pair these values with the average numbers of nonzeros. The scatter plot for these pairs is given in Figure 4. As the figure shows, there is an inverse correlation between the average number of nonzeros of each matrix and performance decrease when 8 threads are used.

Our parallel scaling algorithms based on staticpartitioning obtain much better speedups compared with the simple approaches. In the Intel-based architecture, Figure 2, for 8 threads, the average speedups obtained by CRS-Cut, Scatter plot of the matrices for which an increase on the execution time of COO-Simple is observed when ⌧ is increased from 4 to 8. Each matrix is represented by an (x, y) in the plot where x is the average number of nonzeros per row/column in the matrix and y is the percentage of the increase on the execution time of COO-Simple.

CRS-SOED, and COO-SOED are 3.69, 3.95, and 4.02, respectively. For the AMD-based architecture, Figure 3, these values are 3.49, 3.52, and 3.67. Although we are not aware of any multithreaded scaling algorithm, there exist several SpMxV algorithms designed for shared-memory architectures. In a recent one, the architectures used are similar to the ones in this paper [START_REF] Buluc | Parallel sparse matrix-vector and matrixtranspose-vector multiplication using compressed sparse blocks[END_REF]. Buluc ¸et al. obtained average speedups between 3-3.5 for 8 threads. They also observed that on an Intel Harpertown architecture, which is a dual-socket quad-core system with 12MB of L2 cache, the average Mflops/sec on 8 threads is only 3.5% higher than on 4 processors. Similar to our experiments, the Mflops/sec even decreased for some matrices when ⌧ is increased from 4 to 8. Note that in SpMxV, there is one input and one output vector. On the other hand, for scaling we have two of each. Hence the scaling problem can be considered as hard as SpMxV, or even more harder.

We used different sequential algorithms while computing the speedups for CRS-based algorithms and COO-SOED. Hence the speedup results do not imply that COO-SOED is the fastest parallel scaling algorithm in this paper. As the relative performances in Figures 5 and6, and the execution times in Table II show, this is not true. The figures show the average relative performance of each algorithm for 1, 2, 4, and 8 threads. The relative performance is computed as follows: for each matrix, the execution time of the algorithm is divided to that of the best one. That is, for each matrix the best algorithm is considered as 1 and the other algorithms are evaluated with that value. The averages of these values are computed per each ⌧ value. As Figures 5 and6 show, COO-based algorithms are slower than the CRSbased ones. This is expected since the CRS format implicitly solves the conflict problem while computing the row norms. Furthermore, in a CRS-based algorithm, the partial results can be hold in a private variable (sum t in the algorithms) until the final result is obtained. With this approach each location in d 1 is accessed only once at each iteration.

For both of the architectures CRS-Cut is better than CRS-SOED on the average when ⌧ is 1, 2, or 4. As mentioned before, the when ⌧ is 2, CRS-SOED cannot do any reductions. However, its additional memory overhead will still be there. When ⌧ increases to 8, the reduction on init and get phases starts to show itself and it compensates the overhead of CRS-SOED due to the storage of connectivity information. This experiment shows that minimizing SOED (⇧) can be promising when ⌧ is large.

V. CONCLUSION AND DISCUSSION

We discussed the relevance of hypergraph partitioningbased methods to minimize the parallelization overhead associated with the standard OpenMP parallelization techniques. In particular, we have argued that three common objective functions in the hypergraph partitioning problem correspond exactly to minimizing the parallelization overhead. We note that the use of objective functions of cut-net and connectivity-1 metrics correspond closely to their use in distributed memory parallelization problem. On the other hand, the use of the SOED metric in the shared memory parallelization, to the best of knowledge, does not have such a correspondence.

We have achieved significant speedups by making use of the discussed parallelization techniques. It is assumed that the parallelization would be performed using only OpenMP directives. Therefore, we did not consider hardware aware optimizations techniques (these can help for further performance optimization). Optimizing the codes that use those techniques would require revisiting the models and encapsulate the techniques. Although we believe this is a viable task, it can also be rather intricate, especially for architecture aware optimizations and partitioning which require mapping and partitioning to be performed together.

We did not put a lot of effort to optimize the the cache efficiency of a single thread. The methods proposed for optimizing cache performance in SpMxV operations can also be used to reorganize computations of a single thread (see [START_REF] Yzelman | Fast sparse matrix-vector multiplication by partitioning and ordering[END_REF]Chapter 5]). Here, the issues of partitioning for threads and optimizing a single thread's performance are separate issues and therefore they can be used together.

Algorithm 1 :

 1 1-norm scaling algorithmInput: A: input matrix, square of size n ⇥ n Output: Dr, Dc: row and column scaling matrices b A In for k = 0, 1, 2, . . . until convergence do 1 D1 diag

  , one does not need to store the iterated b A (k) , rather one can access it through left and right multiplications, respectively, with the current scaling matrices D other words, one uses

Algorithm 5 :

 5 Part. based scaling with CRS-SOEDInput: A: n ⇥ n input matrix in CRS format and a partition ⇧ = {R1, R2, . . . , R⌧ } of rows Output: dr, dc: row and column scaling vectors • • • while not converged do • • • init for t = 1 to ⌧ in parallel do for each external column i connected to Rt do d t 2 [i] 0 put • • • IAs same as CRS-Cut get for t = 1 to ⌧ do for each external column i of Rt in parallel do d2

Algorithm 6 :

 6 Part. based scaling with COO-SOEDInput: A: n ⇥ n input matrix in COO format and a partition ⇧ = {Z1, Z2, . . . , Z⌧ } of nonzeros Output: dr, dc: row and column scaling vectors

Figure 1 .

 1 Figure 1. Average normalized execution times of the parallel algorithms using a perfectly balanced partition. For each matrix-⌧ pair, the execution time is normalized by dividing it to the execution time of the algorithm with the partition generated by PaToH. The figure shows that the cut-size minimized partitions improve the performance.

Figure 2 .

 2 Figure 2. Average speedups of scaling algorithms on the Intel-based architecture. For CRS-based algorithms, the speedups are computed by using the execution time of the CRS-based sequential algorithm. For the COO-based ones, we used the COO-based sequential algorithm.

Figure 3 .

 3 Figure 3. Average speedups of scaling algorithms on the AMD-based architecture. For CRS-based algorithms, the speedups are computed by using the execution time of the CRS-based sequential algorithm. For the COO-based ones, we used the COO-based sequential algorithm.

Figure 4 .

 4 Figure 4.Scatter plot of the matrices for which an increase on the execution time of COO-Simple is observed when ⌧ is increased from 4 to 8. Each matrix is represented by an (x, y) in the plot where x is the average number of nonzeros per row/column in the matrix and y is the percentage of the increase on the execution time of COO-Simple.

Figure 5 .

 5 Figure 5. Average relative performances of parallel scaling algorithms on the Intel-based architecture. The relative performance of each algorithm on a particular problem is computed by dividing the average execution time of the algorithm by the best average time for the same problem.

Figure 6 .

 6 Figure 6. Average relative performances of parallel scaling algorithms on the AMD-based architecture. The relative performance of each algorithm on a particular problem is computed by dividing the average execution time of the algorithm by the best average time for the same problem.

  needs to be added to d 1 [i] and d 2 [j]. In a multithreaded setting, these write operations can result in race conditions on d 1 or d 2 . With a naive approach, the threads can evade from write conflicts by using private d 1 s and/or d 2 s. Let ⌧ be the number of the threads and d t 1 and d t 2 denote the private row and column norm arrays for a thread t 2 {1, . . . , ⌧ }, respectively.

Table I PROPERTIES

 I OF THE MATRICES USED IN THE EXPERIMENTS.

	Matrix	n	Z	Avg. deg
	atmosmodd	1,270,432	8,814,880	6.94
	atmosmodl	1,489,752 10,319,760	6.93
	cage14	1,505,785 27,130,349	18.02
	Chebyshev4	68,121	5,377,761	78.94
	Hamrle3	1,447,360	5,514,242	3.81
	NotreDame	325,729	929,849	2.85
	pre2	659,033	5,834,044	8.85
	rajat21	411,676	1,876,011	4.56
	rajat30	643,994	6,175,244	9.59
	Stanford Berk.	683,446	7,583,376	11.10
	torso1	116,158	8,516,500	73.32
	trans5	116,835	749,800	6.42

Table II AVERAGE

 II RUNNING TIMES OF 100 EXECUTIONS OF THE ALGORITHMS IN SECONDS FOR 100 ITERATIONS.

				CRS		COO
	Matrix	T	Simp.	Cut SOED Simp. SOED
	atmosmodd	1	5.54	5.37	5.55	6.72	6.71
		2	3.20	2.90	2.85	4.25	3.38
		4	2.45	1.82	1.84	3.53	2.14
		8	2.76	1.42	1.55	4.62	2.04
	atmosmodl	1	6.95	6.71	6.93	8.33	8.27
		2	3.83	3.65	3.56	5.21	4.30
		4	3.03	2.13	2.28	4.27	2.70
		8	3.27	1.72	1.93	5.52	2.13
	cage14	1	17.80 17.38	18.41	21.04	22.94
		2	9.10	9.77	9.70	11.97	10.87
		4	5.65	5.88	5.96	7.70	6.44
		8	4.93	5.70	4.59	7.69	5.20
	Chebyshev4	1	1.93	1.91	1.92	2.86	3.16
		2	1.22	0.95	1.04	1.55	1.60
		4	1.10	0.67	0.65	0.88	0.79
		8	1.23	0.56	0.55	0.79	0.65
	Hamrle3	1	6.13	5.96	6.19	6.76	6.62
		2	3.50	3.30	3.33	4.49	3.48
		4	2.67	1.86	2.01	3.53	2.50
		8	2.93	1.69	1.42	4.92	1.80
	NotreDame	1	1.18	1.09	1.06	1.17	1.03
		2	0.68	0.60	0.60	0.69	0.64
		4	0.47	0.34	0.31	0.65	0.39
		8	0.57	0.27	0.20	0.88	0.25
	pre2	1	4.02	3.99	4.03	4.67	4.84
		2	2.23	2.05	2.05	2.75	2.40
		4	1.39	1.23	1.17	1.96	1.48
		8	1.45	0.88	0.92	2.48	1.24
	rajat21	1	1.90	1.83	1.88	2.04	2.03
		2	1.08	1.02	1.01	1.25	1.04
		4	0.74	0.58	0.60	1.01	0.70
		8	0.79	0.59	0.52	1.30	0.58
	rajat30	1	4.14	4.20	4.33	4.75	5.04
		2	2.33	2.29	2.48	2.86	2.54
		4	1.61	1.66	1.61	2.08	1.51
		8	1.80	1.41	1.36	2.41	1.32
	Stanford Berk.	1	4.04	3.97	4.01	4.49	4.53
		2	2.24	2.22	2.08	2.69	2.39
		4	1.40	1.15	1.17	1.93	1.40
		8	1.46	0.93	0.83	2.16	1.13
	trans5	1	0.57	0.57	0.58	0.65	0.66
		2	0.31	0.34	0.33	0.36	0.32
		4	0.20	0.19	0.21	0.28	0.17
		8	0.25	0.19	0.20	0.34	0.11
	torso1	1	2.59	2.60	2.60	4.38	4.64
		2	1.32	1.38	1.43	2.25	2.29
		4	0.85	0.94	0.94	1.40	1.21
		8	0.78	0.78	0.78	1.22	1.01

ACKNOWLEDGMENT

This work was supported in parts by the DOE grant DE-FC02-06ER2775 and by the NSF grants CNS-0643969, OCI-0904809, and OCI-0904802.


Ümit V C ¸atalyürek

Mehmet Deveci
email: mdeveci@bmi.osu.edu

Kamer Kaya

Multithreaded Clustering for Multi-level Hypergraph Partitioning

Keywords: Multi-level hypergraph partitioning, coarsening, multithreaded clustering algorithms, multicore programming

Requirements for efficient parallelization of many complex and irregular applications can be cast as a hypergraph partitioning problem. The current-state-of-the art software libraries that provide tool support for the hypergraph partitioning problem are designed and implemented before the game-changing advancements in multi-core computing. Hence, analyzing the structure of those tools for designing multithreaded versions of the algorithms is a crucial tasks. The most successful partitioning tools are based on the multi-level approach. In this approach, a given hypergraph is coarsened to a much smaller one, a partition is obtained on the the smallest hypergraph, and that partition is projected to the original hypergraph while refining it on the intermediate hypergraphs.

The coarsening operation corresponds to clustering the vertices of a hypergraph and is the most time consuming task in a multilevel partitioning tool. We present three efficient multithreaded clustering algorithms which are very suited for multi-level partitioners. We compare their performance with that of the ones currently used in today's hypergraph partitioners. We show on a large number of real life hypergraphs that our implementations, integrated into a commonly used partitioning library PaToH, achieve good speedups without reducing the clustering quality.

I. INTRODUCTION

Hypergraph partitioning is an important problem widely encountered in parallelization of complex and irregular applications from various domains including VLSI design [START_REF] Lengauer | Combinatorial Algorithms for Integrated Circuit Layout[END_REF], parallel scientific computing [START_REF] Aykanat | Hypergraph-partitioning based decomposition for parallel sparse-matrix vector multiplication[END_REF], [START_REF] Aykanat | On twodimensional sparse matrix partitioning: Models, methods, and a recipe[END_REF], sparse matrix reordering [START_REF] Aykanat | Hypergraph partitioning-based fill-reducing ordering for symmetric matrices[END_REF], static and dynamic load balancing [START_REF] Boman | A repartitioning hypergraph model for dynamic load balancing[END_REF], software engineering [START_REF] Bisseling | Partitioning a call graph[END_REF], cryptosystem analysis [START_REF] Bisseling | Mondriaan sparse matrix partitioning for attacking cryptosystems by a parallel block lanczos algorithm-A case study[END_REF], and database design [START_REF] Shekhar | Efficient join-index-based spatial-join processing: A clustering approach[END_REF], [START_REF] Liu | A hypergraph based approach to declustering problems[END_REF], [START_REF] Ozdal | Hypergraph models and algorithms for data-pattern-based clustering[END_REF]. Being such an important problem, considerable effort has been put in providing tool support, see hMeTiS [START_REF] Karypis | hMeTiS: A hypergraph partitioning package[END_REF], MLpart [START_REF] Caldwell | Improved algorithms for hypergraph bipartitioning[END_REF], Mondriaan [START_REF] Vastenhouw | A two-dimensional data distribution method for parallel sparse matrix-vector multiplication[END_REF], Parkway [START_REF] Trifunovic | Parkway 2.0: A parallel multilevel hypergraph partitioning tool[END_REF], PaToH [START_REF] Korpeoglu | PaToH: A Multilevel Hypergraph Partitioning Tool[END_REF], and Zoltan [START_REF] Boman | Zoltan 3.0: Parallel Partitioning, Load Balancing, and Data-Management Services; User's Guide[END_REF].

All the tools above follow the multi-level approach. This approach consists of three phases: coarsening, initial partitioning, and uncoarsening. In the coarsening phase, the original hypergraph is reduced to a much smaller hypergraph after a series of coarsening levels. At each level, vertices that are deemed to be similar are grouped to form vertex clusters, and a new hypergraph is formed by unifying a cluster as a single vertex. That is, the clusters become the vertices for the next level. In the initial partitioning phase, the coarsest hypergraph is partitioned. In the uncoarsening phase, the partition found in the second phase is projected back to the original hypergraph where the partition is locally refined on the hypergraphs associated with each coarsening level.

The coarsening phase is the most important phase of the multi-level approach. This is for the following three reasons. First, the worst-case running time complexity of this phase is higher than the other two phases (initial partitioning and uncoarsening phases have, in most common implementations, linear worst-case running time complexity). Second, as the uncoarsening level performs local improvements, the quality of a partition is highly affected by the quality of the coarsening phase. For example, given a hypergraph, a coarsening algorithm, a conventional initial partitioning algorithm and a refinement algorithm based on the most common ones, very slight variations on vertex similarity metrics can effect the performance quite significantly (see for example the start of Section 5.1 of [START_REF] Uc ¸ar | A matrix partitioning interface to PaToH in Matlab[END_REF]). Third, it is usually the case that the better the coarsening, the faster the uncoarsening phase is. Therefore, the coarsening phase also affects the running time of the other phases.

Our aim in this paper is to efficiently parallelize the coarsening phase of PaToH, a well-known and commonly used hypergraph partitioning tool. The algorithmic kernel of this phase is a clustering algorithm that marks similar vertices to be coalesced. There are two classes of clustering algorithms in PaToH. The algorithms in the first class allow at most two vertices in a cluster. These algorithms are called matching-based or matching algorithms in short. The algorithms in the second class, called agglomerative algorithms, allow any number of vertices to become together to form a cluster. The most effective clustering algorithms in PaToH are agglomerative ones whereas the fastest ones are matching based. We propose efficient parallelization of these two classes of algorithms (Section III). We report practical experiments with PaToH (and its coarsening phase alone) on a recent multicore architecture (Section IV). Our techniques are easily applicable to some other sequential hypergraph partitioners, since they use the same multilevel approach and have similar data structures.

II. BACKGROUND, PROBLEM FORMULATION AND RELATED WORK

A. Hypergraph Partitioning

A hypergraph H = (V, N ) is defined as a set of vertices V and a set of nets (hyperedges) N among those vertices. A net n ∈ N is a subset of vertices, and the vertices in n are called its pins. The size of a net is the number of its pins, and the degree of a vertex is equal to the number of nets that contain it. Graph is a special instance of hypergraph such that each net has size two. We use pins[n] and nets [v] to represent the pins of a net n, and the set of nets that contain vertex v, respectively. Vertices can be associated with weights, denoted with w[•], and nets can be associated with costs, denoted with c[•].

A K-way partition of a hypergraph H is denoted as

Π = {V 1 , V 2 , . . . , V K } where • parts are pairwise disjoint, i.e., V k ∩ V = ∅ for all 1 ≤ k < ≤ K, • each part V k is a nonempty subset of V, i.e., V k ⊆ V and V k = ∅ for 1 ≤ k ≤ K, • union of K parts is equal to V, i.e., K k=1 V k = V.
In a partition Π, a net that has at least one pin (vertex) in a part is said to connect that part. The number of parts connected by a net n, i.e., connectivity, is denoted as λ n . A net n is said to be uncut (internal) if it connects exactly one part (i.e., λ n = 1), and cut (external), otherwise (i.e., λ n > 1).

Let W k denote the total vertex weight in V k (i.e., W k =

v∈V k w[v]
) and W avg denote the weight of each part when the total vertex weight is equally distributed (i.e., W avg = ( v∈V w[v])/K). If each part V k ∈ Π satisfies the balance criterion

W k ≤ W avg (1 + ε), for k = 1, 2, . . . , K (1)
we say that Π is balanced where ε represents the maximum allowed imbalance ratio. The set of external nets of a partition Π is denoted as N E . Let χ(Π) denote the cost, i.e., cutsize, of a partition Π. There are various cutsize definitions [START_REF] Lengauer | Combinatorial Algorithms for Integrated Circuit Layout[END_REF] such as:

χ(Π) = n∈N E c[n] (2)
χ(Π) = n∈N c[n](λ n -1) . (3)
In ( 2) and (3), each cut net n contributes c[n] and c[n](λ n -1) to the cutsize, respectively. The cutsize metric given in (2) will be referred to here as cut-net metric and the one in (3) will be referred as connectivity metric. Given ε and an integer K > 1, the hypergraph partitioning problem can be defined as the task of finding a balanced partition Π with K parts such that χ(Π) is minimized. The hypergraph partitioning problem is NP-hard [START_REF] Lengauer | Combinatorial Algorithms for Integrated Circuit Layout[END_REF].

B. Clustering algorithms for hypergraph partitioning

As said before, there are two classes of clustering algorithms: matching-based ones and agglomerative ones. The matching-based ones put at most two similar vertices in a cluster whereas the agglomerative ones allow any number of similar vertices. There are various similarity metrics-see for example [START_REF] Aykanat | Hypergraph-partitioning based decomposition for parallel sparse-matrix vector multiplication[END_REF], [START_REF] Alpert | Recent directions in netlist partitioning: A survey[END_REF], [START_REF] Karypis | Multilevel hypergraph partitioning[END_REF]. All these metrics are defined on two adjacent vertices (one of them can be a vertex cluster). Two vertices are adjacent if they share a net, i.e., the vertices u and v are matchable if

N uv = nets[u] ∩ nets[v] = ∅.
In order to find a given vertex u's adjacent vertices, one needs to visit each net n ∈ nets[u] and then visit each vertex v ∈ pins [n]. Therefore, the computational complexity of the clustering algorithms is at least in the order of

n∈N |pins[n]| 2 .
As mentioned in the introduction, the other two phases of the multi-level approach have a linear time worst case time complexity. As n∈N |pins[n]| 2 is most likely to be the bottleneck, an effective clustering algorithm's worst case running time should not pass this limit for the algorithm to be efficient as well.

The sequential implementations of the clustering algorithms in PaToH proceed in the following way to have a running time proportional to n∈N |pins[n]| 2 . The vertices are visited in a given (possibly random) order and each vertex u, if not clustered yet, is tried to be clustered with the most similar vertex or cluster. In the matching-based ones, the current vertex u if not matched yet, chooses one of its unmatched adjacent vertices according to a criterion. If such a vertex v exists, the matched pair u and v are marked as a cluster of size two. If there is no unmatched adjacent vertex of u, then vertex u remains as a singleton cluster. In the agglomerative ones, the current vertex u, if not marked to be in a cluster yet, can choose a cluster to join (thus forming a cluster of size at least three), or can create another cluster with one of its unmatched adjacent vertices (thus forming a cluster of size two). Hence in agglomerative clustering, vertex u never remains as a singleton, as long as it is not isolated (i.e., not connected to any net).

For the clustering algorithms in this paper, there exists a representative vertex for each cluster. When a vertex u ∈ V is put into a cluster, we set rep[u] to the representative of this group. When a singleton vertex u chooses another one v, we choose one as the representative and set rep[u] and rep[v] accordingly. For all the algorithms, we assume that rep[u] is initially null for all u ∈ V. This will also be true if u remains singleton at the end of the algorithm.

Algorithm 1 presents one of the matching-based clustering algorithms that are available in PaToH. In this algorithm, the vertex u (if not matched yet) is matched with currently unmatched neighbor v with the maximum connectivity, where the connectivity refers to the sum of the costs of the common nets. This matching algorithm is called as Heavy Connectivity Matching (HCM) in PaToH [START_REF] Aykanat | Hypergraph-partitioning based decomposition for parallel sparse-matrix vector multiplication[END_REF], and Inner Algorithm 1: Sequential greedy matching (HCM)

Data: H = (V, N ), rep 1 for each vertex u ∈ V in a given order do if rep[u] = null then 2 adju ← {} for each net n ∈ nets[u] do for each vertex v ∈ pins[n] do if rep[v] = null then 3 pins[n] ← pins[n] \{v} else 4 if v / ∈ adju then 5 adju ← adju ∪ {v} 6 conn[v] ← conn[v] + c[n] v * ← u conn * ← 0 for each vertex v ∈ adju do if conn[v] > conn * and v = u then conn * ← conn[v] v * ← v conn[v] ← 0 if u = v * then rep[v * ] ← u rep[u] ← u
Product Matching (IPM) in Zoltan [START_REF] Devine | Parallel hypergraph partitioning for scientific computing[END_REF] and Mondriaan [START_REF] Vastenhouw | A two-dimensional data distribution method for parallel sparse matrix-vector multiplication[END_REF].

One can have different variations of this algorithm by changing the vertex visit order (line 1) and/or using different scaling schemes while computing the contribution of each net to its pins (line 6). The array conn[•] of size |V| is necessary to compute the connectivity of the vertex u and all its adjacent vertices in time linearly proportional to the number of adjacent vertices. The operation at line 3 is again for efficiency. It removes the matched vertices from pins[n], hence the next searches on pins[n] will take less time. Algorithm 2 presents one of the agglomerative clustering algorithms that are available in PaToH. Similar to the sequential HCM algorithm, vertices are again visited in a given order. If a vertex u has already been clustered, it is skipped. However, an unclustered vertex u can choose to join an existing cluster, can start a new cluster with a vertex, or stay as a singleton cluster. Therefore, compared to the previous algorithm, all adjacent vertices of the current vertex u are considered for selection. In order to avoid building an extremely large cluster (which would cause load balance problem in initial partitioning and refinement phases), we also enforce that weight of a cluster must be smaller than a given value maxW . Our experience shows that such restriction is not needed in matching based clustering, since at each level only at most two vertices can be clustered together.

In Algorithm 2, we use the total shared net cost (heavy connectivity clustering) as the similarity metric. In practice Algorithm 2: Sequential agglomer. clustering (HCC)

Data: H = (V, N ), maxW , rep for each vertex u ∈ V in a given order do if rep[u] = null then adju ← {} for each net n ∈ nets[u] do for each vertex v ∈ pins[n] do if v / ∈ adju then adju ← adju ∪ {v} conn[v] ← conn[v] + c[n] v * ← u conn * ← 0 for each vertex v ∈ adju do if u = v then continue v r ←rep[v] if v r = null then v r ← v if v r = v then conn[v r ] ← conn[v r ] + conn[v] conn[v] ← 0 adju ← adju ∪ {v r } \ {v} 1 totW ← w[u] + w[v r ] if conn[v r ] > conn * then if totW < maxW then conn * ← conn[v r ] v * ← v r for each vertex v ∈ adju do conn[v] ← 0 if u = v * then rep[v * ] ← v * rep[u] ← v * w[v * ] ← w[v * ] + w[u]
(and in our experiments), we use the absorption clustering metric (implemented in PaToH) which divides the contribution of each net to the number of its pins. That is, a net

n contributes c[n]/|pins[n]| to the similarity value instead of c[n]
. This metric favors clustering vertices connected via nets of small sizes. The sequential code in PaToH also divides the overall similarity score between two vertices by the weight of the cluster which will contain u (the value totW at line 1). Hence, to compare the performance of our multithreaded clustering algorithms with PaToH, we also use this modified similarity metric in our implementations. However, for simplicity, we will continue to use the heavy connectivity clustering metric in the paper.

C. Metrics

We define the metrics of cardinality and quality to compare different clustering methods. The cardinality is defined as the number of clustering decisions taken by an algorithm, i.e., cardinality

= u∈V,rep[u] =null (|{v ∈ V : rep[v] = u}| -1) .
In the multi-level framework, this represents the reduction on the number of vertices between two consecutive coarsening levels. The quality of a clustering is defined as the sum of the similarities between each vertex pair which resides in the same cluster, i.e., quality

= u∈V v∈Cu n∈Nuv c[n] 2 ,
where

C u = {v ∈ V \ {u} : rep[v] = rep[u] and rep[v] =
null} is the set of vertices which are in the same cluster with u, and N uv is the set of nets shared by u and v.

Although the definitions are generic to be used for both matching-based and agglomerative clustering, we do not use these criteria to compare a matching-based clustering heuristic with an agglomerative one, since the latter has an obvious advantage.

D. Related work

For a given hypergraph H = (V, N ), let A be the vertexnet incidence matrix, i.e., the rows of A correspond to the vertices of H, and the columns of A correspond to the nets of H such that a vn = 1 iff v ∈ pins [n]. Consider now the symmetric matrix M = AA T -diag(AA T ). The matrix M can be effectively represented by an undirected graph G(M ) with |V| vertices and having an edge of weight m uv between two vertices u and v if m uv = 0. That is, there is a one-toone correspondence between the vertices of H and G(M ). As m uv = 0 iff the vertices u and v of H are adjacent, the correspondence implies that a matching among the vertices of H corresponds to a matching on the vertices of G(M ). Therefore, various matching algorithms and heuristics for graphs can be used on G(M ) to find a matching among the vertices of H.

Bisseling and Manne [START_REF] Manne | A parallel approximation algorithm for the weighted maximum matching problem[END_REF] propose a distributed memory, 1/2-approximate algorithm to find weighted matchings in graphs. Building on this work, C ¸atalyürek et al. [START_REF] Wasniewski | Distributed-memory parallel algorithms for matching and coloring[END_REF] present efficient distributed-memory parallel algorithms and scalable implementations. Halappanavar et al. [START_REF] Halappanavar | Approximate weighted matching on emerging manycore and multithreaded architectures[END_REF] present an efficient shared-memory implementations for computing 1/2approximate weighted matchings. For maximum cardinality matching problem, in a recent work, Patwary et al. [START_REF] Patwary | Parallel greedy graph matching using an edge partitioning approach[END_REF] propose a distributed memory, sub-optimal algorithm.

There are a number of reasons why we cannot use aforementioned algorithms. First and foremost, storing the graph G(M ) requires a large memory. The time required to compute this graph is about as costly as computing a matching in H in a sequential execution. Second, it is our experience (with the coarsening algorithms within the multilevel partitioner PaToH) that one does not need a maximum weighted matching, nor a maximum cardinality one, nor an approximation guarantee to find helpful coarsening schemes. Third, while matching the vertices of a hypergraph, we sometimes need to avoid matched vertices become too big, or favor vertex clusters with smaller weights (due to multilevel nature of the partitioning algorithm), and vertices that are mostly related via nets of smaller size. These modifications can be incorporated into the graph G(M ) by adjusting the edge weights (or by leaving some edges out). This will help reduce the memory requirements of the graph matching based algorithms. However, the computational cost of constructing the graph remains the same. Almost all of the most effective sequential clustering algorithms implemented in PaToH for coarsening purposes has the same worst case time complexity but are much faster in practice. We, therefore, cannot afford building the graph G(M ) or its modified versions and call existing graph matching algorithms. Furthermore, agglomerative clustering algorithms cannot be accomplished by using the aforementioned algorithms.

We highlight that the matching-based clustering algorithms considered in this paper can be perceived as a graph matching algorithm adjusted to work on an implicit representation of the graph G(M ) or its modified versions. However, to the best of our knowledge, there is no immediate parallel graph-matching based algorithm that is analogous to the agglomerative clustering algorithms considered in this work (although variants of agglomerative coarsening for graphs exist, see for example [START_REF] Abou-Rjeili | Multilevel algorithms for partitioning power-law graphs[END_REF] and [START_REF] Chevalier | Comparison of Coarsening Schemes for Multilevel Graph Partitioning[END_REF]).

As a sanity check, we implemented a modified version of the sequential 1/2-approximation algorithm of Halappanavar et al. [START_REF] Halappanavar | Approximate weighted matching on emerging manycore and multithreaded architectures[END_REF] which works directly on hypergraphs. In other words, instead of explicitly constructing the graph G(M ), adjacencies of vertices are constructed on the fly, like Algorithm 1. We compared the quality and cardinality of this algorithm with that of the greedy sequential matching HCM. The approximation algorithm obtained matchings with better quality by 14% while the cardinalities were the same. However, this good performance comes with a significant execution time overhead, yielding 6.5 times slower execution. When we integrated the 1/2-approximation algorithm into the coarsening phase of PaToH, we observed that better matching quality helps the partitioner to obtain better cutsize. For example, when the partitioner is executed 10 times with random seeds, the average cutsize of 1/2approximation algorithm is 8% better than the one obtained by using HCM. However, when we compare the minimum of these cutsizes, HCM outperforms the approximation algorithm by 2%. Moreover, the difference between the minimum cut obtained by using HCM with the average cut obtained by using the approximation algorithm is 14% in favor of HCM. After these preliminary experiments, we decided to parallelize HCM and HCC, since they are much faster, and one can obtain better cutsize by using them.

III. MULTITHREADED CLUSTERING FOR COARSENING

In this section, we will present three novel parallel greedy clustering algorithms. The first two are matching-based and the third one is a greedy agglomerative clustering method.

A. Multithreaded matching algorithms

To adapt the greedy sequential matching algorithm for multithreaded architectures, we use two different approaches. In the first one, we employ a locking mechanism which prevents inconsistent matching decisions between the threads. In the second approach, we let the algorithm run almost without any modifications and then use a conflict resolution mechanism to create a consistent matching.

The lock-based algorithm is given in Algorithm 3. The structure of the algorithm is similar to the sequential one except the lines 2 and 5, where we use the atomic CHECKAND-LOCK operation. To lock a vertex u, this operation first checks if u is already locked or not. If not, it locks it and returns true. Otherwise, it returns false. Its atomicity guarantees that a locked vertex is never considered for a matching. That is, both the visited vertex u (at line 1) and the adjacent vertex v must not be locked to consider the matching of u and v. If they are, and if the similarity of v is bigger than the current best (at line 3), the algorithm keeps v as the best candidate v * . When a better candidate is found, the old one is UNLOCKed to make it available for other threads (line 6), and to construct better matchings in terms of cardinality and quality.

Algorithm 3: Parallel lock-based matching

Data: H = (V, N ), rep 1 for each vertex u ∈ V in parallel do 2 if CHECKANDLOCK(u) then adju ← {} for each net n ∈ nets[u] do for each vertex v ∈ pins[n] do if v / ∈ adju then adju ← adju ∪ {v} conn[v] ← conn[v] + c[n] v * ← u conn * ← 0 3 for each vertex v ∈ adju do 4 if conn[v] > conn * then 5 if CHECKANDLOCK(v) then if u = v * then 6 UNLOCK(v * ) conn * ← conn[v] v * ← v conn[v] ← 0 if u = v * then rep[u] ← u rep[v * ] ← u
As a different approach without a lock mechanism, we modify the sequential code slightly and execute it in a multithreaded environment. If the for loop at line 1 of Algorithm 1 is executed in parallel, different threads may set rep[u] to different values for a vertex u. Hence, the rep array will contain inconsistent decisions. To solve this issue, one can make each thread use a private rep array and store all of its matching decisions locally. Then, a consistent matching can be devised from this information in another phase. Another idea is keeping the sequential code (almost) as is, letting the threads create conflicts, and resolving the conflicts later. Our preliminary experiments show that there is not much difference between the performances of these two approaches in terms of the cardinality and the quality. However, the first one requires more memory: one rep array per thread compared to a shared one. Hence, we followed the second idea and use a conflict resolution scheme with O(|V|) complexity. Algorithm 4 shows the pseudocode of our parallel resolution-based algorithm.

Algorithm 4: Parallel resolution-based matching

Data: H = (V, N ), rep for each vertex u ∈ V in parallel do if rep[u] = null then adju ← {} for each net n ∈ nets[u] do for each vertex v ∈ pins[n] do if v / ∈ adju then adju ← adju ∪ {v} conn[v] ← conn[v] + c[n] v * ← u conn * ← 0 for each vertex v ∈ adju do 1 if rep[v] = null then if conn[v] > conn * then if u = v then conn * ← conn[v] v * ← v conn[v] ← 0 2 if rep[u] = rep[v] = null then 3 rep[u] ← v * 4 rep[v * ] ← u 5 for each vertex u ∈ V in parallel do v ← rep[u] if v = null and u = rep[v] then rep[u] ← null for each vertex u ∈ V in parallel do v ← rep[u] if v = null and u < v then rep[u] ← u
Our conflict resolution scheme starts at line 5 of Algorithm 4. Note that instead of setting a fixed representative for the matched vertices u and v * , we set their rep values to each other. This bidirectional information is then used in our resolution scheme to check if the rep array contains inconsistent information. That is if u = rep[rep [u]] for a vertex u, we know that at least two threads matched either u or v with different vertices. If this is the case, the resolution scheme acts greedily and aggressively sets rep[u] to null indicating u will be unmatched. After the first loop at line 5, which is executed in parallel, the rep array will contain the matching decisions consistent with each other. Then, with the last parallel loop, we set the representatives for each matched pair.

The proposed resolution scheme is sufficient to obtain a valid matching in the multithreaded setting. However, we slightly modify Algorithm 1 to obtain better matchings. Since each conflict will probably cost a pair, and losing a pair reduces matching cardinality and hence, quality, we desire lesser conflicts. To avoid them, at line 1 of Algorithm 4, we check if a vertex v adjacent to u is already matched. If we detect an already matched candidate, we do not consider it as a possible mate for u. Furthermore, at line 2, we again verify if u and v are already matched right before matching them. This check is necessary since either of them could have been matched by another thread after the current one starts considering them.

B. Multithreaded agglomerative clustering

To adapt the sequential agglomerative clustering algorithm to multithreaded setting, we use the same lock-based approach integrated to Algorithm 3. The pseudocode of the parallel agglomerative clustering algorithm is given in Algorithm 5. The algorithm visits the vertices in parallel, and when a thread visits a vertex u, it tries to lock u. If u is already locked the thread skips u and visits the next vertex. If u is not locked but is already a member of a cluster, the thread unlocks u. Since a cluster cannot be the source of a new cluster, this is necessary. On the other hand, if u is a singleton vertex, the thread continues by computing the similarity values for each adjacent vertex and then traverses the adjacency list adj u along the same lines as the sequential algorithm. The main difference here is the locking request for v r (line 1) which is either set to v if v is singleton, or to the representative of the cluster that v resides in. Before considering v r as a matching candidate, this lock is necessary. However, if v r is already the best candidate, it is not so (since the thread has already grabbed v r ). When the lock is granted, the thread checks if the adjacent vertex v, which was singleton before, has been assigned to a cluster by another thread. If this is the case, the thread unlocks the representative and continues with the next adjacent vertex. Otherwise, it recomputes the total weight of u and v r (line 3) since new vertices might have been inserted to v r 's cluster by other threads. Since insertions can only increase w[v r ] and conn[v r ], we do not need to compare conn[v r ] with conn * again. On the other hand, since we cannot construct clusters with large weights, we need to check if totW is still smaller than maxW (line 4). When the best candidate v * is found, we put u in the cluster v * represents and update the rep and w arrays accordingly. Unlike the matching based algorithms, a cluster is allowed to be a candidate more than once throughout the execution. Hence, at the end of the iteration (lines 5 and 6) we unlock all the vertices that are locked during this iteration.

C. Implementation Details

To obtain lock functionality for the multithreaded clustering algorithms described in the previous section, we use the compare and exchange CPU instruction which exists in x86 and Itanium architectures. We first allocate a lock array of length |V| and initialize all entries to 0. For each call of the corresponding function in C, __sync_bool_compare_and_swap, the entry related with the lock request is compared with 0. In case of equality, it is set to 1, and the function returns true. On the other hand, if the entry is not 0 then it returns false. To unlock a vertex, we simply set the related entry in the lock array to 0.

Although this function provides great support and flexibility for concurrency, our preliminary experiments show that it can also reduce the efficiency of a multithreaded algorithm. To alleviate this, we try reduce the number of calls on this function by adding an if statement before each lock request which helps us to see if the lock is really necessary. We observe significant improvements on the execution times due to these additional if statements. For example, the parallel lock-based matching algorithm described in the previous section should be implemented as in Algorithm 6 to make it much faster. We stress that the if statements at lines 1 and 3 do not change anything in the execution flow. That is if a vertex is not locked, it cannot be also matched since a matched vertex always stays locked. Hence everything that can pass the lock requests (lines 2 and 4) can also pass the previous if statements. However, the opposite is not true.

We used the same hypergraph data structure with PaToH. We store the ids of the vertices of each net n, that is its pins, consecutively in an array ids of size In the actual implementation of Algorithm 1, the set adj u corresponds to an array of maximum size |V|, and an integer which keeps the number of adjacent vertices in the array. With this pair, the vertex addition (line 5) and reset (line 2) operations take constant time. Furthermore, to find if a vertex v is a member of adj u (line 4), we use conn[v] since the edge costs are positive, and conn[v] > 0 if and only if v ∈ adj u . The implementation of these lines are the same for other algorithms.

IV. EXPERIMENTAL RESULTS

The algorithms are tested on a computer with 2.27GHz dual quad-core Intel Xeon CPUs with 2-way hyper-threading enabled, and 48GB main memory. All of the algorithms are implemented in C and OpenMP. The compiler is icc version 11.1 and -O3 optimization flag is used.

To generate our hypergraphs, we used real life matrices from the University of Florida (UFL) Sparse Matrix Collection (http://www.cise.ufl.edu/research/sparse/matrices). We randomly choose 70 large, square matrices from the library and create corresponding hypergraphs using the column-net hypergraph model [START_REF] Aykanat | Hypergraph-partitioning based decomposition for parallel sparse-matrix vector multiplication[END_REF]. An overall summary of the properties of these hypergraphs is given in Table I. The complete list of matrices is at http://bmi.osu.edu/ ∼ kamer/multi coarse matrices.txt.

A. Individual performance of the algorithms

We first compare the performance of the multithreaded clustering algorithms with respect to the cardinality, quality and speedup as standalone clustering algorithms.

1) Performance on cardinality and quality: For matching based clustering, a matching with the best quality can be found by first constructing the matrix M = AA Tdiag(AA T ) where A is the vertex-net incidence matrix. Then a maximum weighted matching in G(M ), the associated weighted graph of M , is also the maximum quality matching.

We use Gabow's maximum weighted matching algorithm [START_REF] Gabow | An efficient implementation of Edmonds' algorithm for maximum matching on graphs[END_REF]-implemented by Rothberg and available as a part of The First DIMACS Implementation Challenge (available at ftp://dimacs.rutgers.edu/pub/netflow). This algorithm has a time complexity of O(|V| 3 ) for a graph on |V| vertices and finds a maximum weighted matching in the graph, not necessarily among the maximum cardinality ones.

Due to the running time complexity of the maximum quality matching algorithm, it is impractical to obtain the relative performance of the clustering algorithms on our original dataset. We therefore use an additional data set containing considerably small matrices. The new dataset contains all 289 square matrices in the UFL sparse matrix collection with at least 3, 000 and at most 10, 000 rows. We construct the hypergraphs for each of these matrices and find the maximum quality matching on them. The relative performance of an algorithm is computed by dividing its cardinality and quality scores to those of Gabow's quality matching algorithm. Table II shows the minimum, the maximum, and the geometric mean of all 289 relative performance for each algorithm. As the table shows, the sequential algorithm and its parallel lock-based variant are only 17-19% far from the optimal in terms of quality and almost equal in terms of cardinality. Considering the difference in computational complexities, we can argue that their relative performance is reasonably good. The lockbased algorithm performs slightly better than the sequential one. This demonstrates that while reducing the execution time, the proposed lock-based parallelization does not hamper the performance of the sequential matching algorithm in terms of both cardinality and quality. For this experiment, the parallel algorithms were executed with 8 threads.

Figure 1 shows the performance profiles generated to Figure 1. Performance profiles for sequential and multithreaded matching algorithms with respect to the maximum quality matchings. A point (x, y) in the profile graph means that with y probability, the quality of the matching found by an algorithm is larger than max/x where max is the maximum quality for that instance.

analyze the results in more detail. A point (x, y) in the profile graph means that with y probability, the quality of the matching found by an algorithm is more than max/x where max is the maximum quality for that instance. The figure shows that for this data set, the resolution-based algorithm performs worse than the other algorithms. The lock-based algorithm obtains matchings which are at most 1.25 times less than max with 74% probability. However, the resolution-based matching has this performance only for 45% probability. While obtaining matchings having at most 15% worse quality than the optimum, the probabilities are 56% and 28% for the lock-and resolution-based algorithms, respectively. Hence, the former performs two times better than the latter. For the original data set with large hypergraphs, the relative performance of the multithreaded algorithms are given with respect to that of their sequential versions. Table III shows the statistics for this experiment. Table III shows that the multithreaded versions of lockbased and agglomerative algorithms perform as good as their sequential versions in terms of cardinality and quality. Hence, once again, we can conclude that parallelization via locks does not diminish the performance of the multithreaded algorithms. On the other hand, the resolution-based algorithm is outperformed by other algorithms in terms of quality. Speedups for matching and clustering: For the lock-and resolution-based algorithms, speedup is computed by using the execution times of sequential greedy matching algorithm. For the parallel agglomerative one, its sequential version is used.

Table IV shows the average numbers of matched vertices and conflicts for the proposed resolution-based algorithm with respect to the number of threads. To compute the averages, we execute the algorithm on each hypergraph ten times and report the geometric mean of these results. As expected, the number of the conflicts increases with the number of threads. However, when compared with the cardinality of the matching, the conflicts are at most 0.7% of the total match count for a single graph instance. This shows that the probability of a conflict is very low even with 8 threads. 2) Speedup: Figure 2 shows the speedups achieved by the multithreaded algorithms. On the average, the algorithms obtain decent speedups compared to their sequential versions up to 8 threads. In 8-thread experiments, the resolution-based algorithm has the highest speedup of 5.87, which is followed by the parallel agglomerative algorithm with a speedup of 5.82. The lock-based algorithm has the least speedup of 5.23 in this category. However, this is still decent especially when we consider the 5-7% overhead due to OpenMP and atomic operations.

To analyze the scalability of the algorithms from a closer point of view, we draw the speedup profiles in Figure 3. The resolution-based algorithm has better scalability in general. For example, with 4 threads, the probability that the resolution-based algorithm obtains a speedup of at least 3.2 is 83%. The same probabilities for the lock-based and parallel agglomerative algorithms are 54% and 65%, respectively. With 8 threads, the resolution-based version achieves at least 6.6 speedup for 33% the hypergraphs. However, the lock-based and parallel agglomerative algorithms achieve the same speedup only in for 16% and 20% of them. Hence, the resolution-based algorithm is the best among the ones proposed in this paper in terms of scalability.

B. Multi-level performance of the algorithms

As mentioned in the introduction, we integrated our algorithms into the coarsening phase of PaToH [START_REF] Korpeoglu | PaToH: A Multilevel Hypergraph Partitioning Tool[END_REF]. In this section, we first investigate how our algorithms scale for the clustering operations inside PaToH. The overall performance of a clustering algorithm in such a setting can be different from the standalone performance of the same algorithm, since in the multi-level framework, the hypergraphs are coarsened until the coarsest hypergraph is considerably small (for example, until the number vertices reduces below 100).

Figure 4 shows that the speedups on the multi-level clustering part are slightly worse than that of the standalone clustering. For example, the average speedups for the 8( 4 Since we only parallelize the clustering operations inside the partitioner, the speedups we obtain on the overall execution time cannot be equal to the number of threads, even in the ideal case. To find the ideal speedups, we use Amdahl's law [START_REF] Amdahl | Validity of the single processor approach to achieving large-scale computing capabilities[END_REF]:

speedup ideal = 1 (1 -r) + r #threads (4)
where r is the ratio of the total clustering time to the total time of a sequential execution. To find the ideal speedup on the average, we compute (4) for each hypergraph and then take the geometric mean since we do the same for actual speedups.

Figure 5 shows the ideal and actual speedups of the multithreaded algorithms. Since the ideal speedup lines (in dashed style) are drawn by using different sequential algorithms for the matching-based and agglomerative clustering, we separated these two cases and draw two different charts. On the average, all the algorithms obtain speedups close to the ideal. Among the matching-based algorithms, the lockbased one is more efficient since its speedup is closer to the ideal. This is interesting since according to Figure 4, it has less speedup on the total clustering time. At first sight, this looks like an anomaly because this is the only part that has been parallelized. However, the since lock-based algorithm's quality is better (Table III), there is probably less work remaining for the refinement heuristics in the uncoarsening phase. Hence, in total, one can achieve better speedup by using the lock-based algorithm rather than the resolutionbased one.

We also obtain good speedups with the parallel agglomerative algorithm. For 2, 4, and 8 threads, the algorithm makes PaToH only 6%, 6% and 10% slower, respectively, than the best possible parallel execution time.

Parallelization of both matching-based and agglomerative clustering algorithm reduces the total execution time significantly. As mentioned before, matching-based algorithms are faster than the agglomerative ones. However, Figure 6 shows that PaToH is 20-30% faster when an agglomerative algorithm is used in the coarsening phase rather than a matching based one. According to our experiments, the coarsening phase is indeed 25% slower with the agglomerative clustering algorithm. However, the total execution time is 13% less. The difference comes from the reduction on the time of initial partitioning and uncoarsening/refinement phases. The initial partitioning takes less time because the coarsest hypergraph has fewer number of vertices with an agglomerative algorithm. In addition, the agglomerative clustering results in 25% less cutsize compared to the matching-based clustering. Hence, we can claim that it is more suitable for the cutsize definition given in [START_REF] Aykanat | On twodimensional sparse matrix partitioning: Models, methods, and a recipe[END_REF].

When equipped with the multithreaded clustering algorithms, the cutsize of the partition found by PaToH is almost equal to the original cutsize obtained by using the sequential versions. For the agglomerative case, the cutsize changes only up to 1%. This is also true with the lock-based matching algorithm when compared with sequential greedy matching. For the resolution-based matching algorithm, there is at most 3% percent increase in the cutsize on average. On the other hand, as shown above the algorithms scale reasonably well.

V. CONCLUSION

Clustering algorithms are the most time consuming part of the current-state-of-the-art hypergraph partitioning tools that follow the multi-level framework. We have investigated the matching-based and agglomerative clustering algorithms. We have argued that the matching-based clustering algorithms can be perceived as a matching algorithm on an implicitly represented undirected, edge weighted graph, whereas there is no immediate equivalent algorithm for the agglomerative ones.

We have proposed two different multithreaded implementations of the matching-based clustering algorithms. The first one uses atomic lock operations to prevent inconsistent matching decisions made by two different threads. The second one lets the threads perform matchings as they would do in a sequential setting and then later on resolves the conflicts that would arise. We have also proposed a multithreaded agglomerative clustering algorithm. This algorithm also uses locks to prevent conflicts.

We have presented different sets of experiments on a large number of hypergraphs. Our experiments have demonstrated that the multithreaded clustering algorithms perform almost as good as their sequential counter parts, sometimes even better in terms of clustering quality and cardinality. The experiments have also shown that our algorithms achieve decent speedups (the best was 5.87 with 8 threads).

We integrated our algorithms to a well-known hypergraph partitioner PaToH. This integration makes PaToH 1.85 times faster where the ideal speedup is 2.07. In addition, it does not worsen the cutsizes obtained. We observed that clusterings with better quality helps the partitioner to obtain better cuts. Fortunately, the multi-level framework may tolerate slower algorithms which generate better clusterings in terms of cardinality and quality. This is because of the fact that such clusterings will reduce the time required for the initial partitioning and uncoarsening phases. However, there is a limit with this tolerance. If the algorithm is too slow, one can execute the partitioner several times with a faster, parallelizable algorithm that generates clusterings with acceptable quality and achieve even better cutsizes.

  n∈N |pins[n]|. We also keep another array xids of size |N | + 1, which stores the start index of each net's pins. Hence, in our implementation, the pins of a net n, denoted by pins[n] in the pseudocodes, are stored in ids[xids[n]] through ids[xids[n + 1] -1]. With the data structures above, the computational complexity of the clustering algorithms in this paper are in the order of n∈N |pins[n]| 2 , since all non-loop lines in their pseudocodes have O(1) complexity. For example, in Algorithm 1, to remove matched vertices from pins[n] (line 3), we keep a pointer array netend of size N where netend[n] initially points to the last vertex in pins[n] for all n ∈ N . Then, to execute pins[n] ← pins[n] \{v}, we only decrease netend[n] and swap v with the vertex in the new location. In this way, we also keep the list of vertices connected to each net unchanged since we only reorder them.

Figure 2 .

 2 Figure 2.

  ) threads case, are 5.25(3.22), 5.56(3.40), and 5.47(3.23) for the lock-based, resolution-based, and parallel agglomerative algorithms, respectively.

8 Figure 3 .

 83 Figure 3. Speedup profiles for the multithreaded algorithms: A point (x, y) in the profile graph means that with y probability, the speedup obtained by the parallel algorithm will be more than x.

Figure 4 .

 4 Figure 4. Speedups on the time spent by the clustering algorithms in the multi-level approach.

Figure 5 .Figure 6 .

 56 Figure 5. Overall speedup on the total execution time of PaToH. The ideal speedup line is drawn by using Amdahl's law.

if rep[v] = null then

if CHECKANDLOCK(v) then • • • • • • • • •

ACKNOWLEDGMENT

This work was supported in parts by the DOE grant DE-FC02-06ER2775 and by the NSF grants CNS-0643969, OCI-0904809, and OCI-0904802.

Algorithm 5: Parallel agglomerative clustering

Data: H = (V, N ), maxW , rep

Algorithm 6: Parallel lock-based matching: modified