Grégoire Pichon

Mathieu Faverge

Pierre Ramet

Jean Roman

Reordering strategy for blocking optimization in sparse linear solvers

Keywords: Sparse linear solver, nested dissection, sparse matrix ordering, heterogeneous architectures algèbre linéaire creuse, dissection emboitée, renumérotation, architectures hétérogènes Chevron

Solving sparse linear systems is a problem that arises in many scientic applications, and sparse direct solvers are a time consuming and key kernel for those applications and for more advanced solvers such as hybrid direct-iterative solvers. For this reason, optimizing their performance on modern architectures is critical. The preprocessing steps of sparse direct solvers, ordering and block-symbolic factorization, are two major steps that lead to a reduced amount of computation and memory and to a better task granularity to reach a good level of performance when using BLAS kernels. With the advent of GPUs, the granularity of the block computation became more important than ever. In this paper, we present a reordering strategy that increases this block granularity. This strategy relies on the block-symbolic factorization to rene the ordering produced by tools such as Metis or Scotch, but it does not impact the number of operations required to solve the problem. We integrate this algorithm in the PaStiX solver and show an important reduction of the number of o-diagonal blocks on a large spectrum of matrices. This improvement leads to an increase in eciency of up to 20% on GPUs.

Introduction

Many scientic applications, such as electromagnetism, astrophysics, and computational uid dynamics, use numerical models that require solving linear systems of the form Ax = b. In those problems, the matrix A can either be considered as dense (almost no zero entries) or sparse (mostly zero entries). Due to multiple structural and numerical dierences that appear in those problems, many dierent solutions exist to solve them. In this paper, we focus on problems leading to sparse systems with a symmetric pattern and more specically on direct methods which factorize the matrix A in LL t , LDL t or LU , with L, D and U respectively unit lower triangular, diagonal, and upper triangular according to the problem numerical properties. Those sparse matrices appear mostly when discretizing Partial Dierential Equations (PDEs) on 2D and 3D nite element or nite volume meshes. The main issue with such factorizations is the ll-in zero entries becoming non-zero that appears in the factorized form of A during the execution of the algorithm. If not correctly considered, the ll-in can transform the sparse matrix into a dense one which might not t in memory. In this context, sparse direct solvers rely on two important preprocessing steps to reduce this ll-in and control where it appears.

The rst one nds a suitable unknown ordering that aims at minimizing the ll-in to limit the memory overhead and the oating point operations (Flop) required to complete the factorization.

The problem is then transformed into (P AP t )(P x) = P b where P is an orthogonal permutation matrix. A wide array of literature exists on solutions to graph reordering problems; the most commonly used for sparse direct factorization being the nested dissection recursive algorithm introduced by George [START_REF] George | Nested dissection of a regular nite element mesh[END_REF].

The second preprocessing step of sparse direct solvers is the block-symbolic factorization [START_REF] Charrier | Algorithmic study and complexity bounds for a nested dissection solver[END_REF]. This step analytically computes the block-structure of the factorized matrix from the reordering step and from a supernode partition of the unknowns. It allows the solver to create the data structure that will hold the nal matrix instead of allocating it at runtime. The goal of this step is also to block the data in order to eciently apply matrix-matrix operations, also known as BLAS Level 3 [3], on those blocks instead of scalar operations. For this purpose, extra ll-in, and by extent extra computations, might be added in order to reduce the time to solution. However, the size of those blocks might not reach the sucient size to extract all the performance from the BLAS kernels.

Modern architectures, whether based on CPUs, GPUs, or Intel Xeon Phi may be ecient with a performance close to the theoretical peak. This can be achieved only if the data size is large enough to take advantage of caches, vector units, and provides a larger ratio of computation per byte. Accelerators such as GPUs or Intel Xeon Phi require even larger blocking sizes than the ones for CPUs due to their particular architectural features.

In order to provide more suited block sizes to kernel operations, we propose in this paper an algorithm that reorders the unknowns of the problem to increase the average size of the odiagonal blocks in block-symbolic factorization structures. The major feature of this solution is that, based on an existing nested dissection ordering for a given problem, our solution will keep constant the amount of ll-in generated during the factorization. So, the amount of memory and computation to store and compute the factorized matrix is invariant. The consequence of this increased average size is that the number of o-diagonal blocks is largely reduced, diminishing the memory overhead of the data structures used by the solver and the number of tasks required to compute the solution in task-based implementations [START_REF] Lacoste | Taking Advantage of Hybrid Systems for Sparse Direct Solvers via Task-Based Runtimes[END_REF][START_REF] Agullo | Multifrontal QR factorization for multicore architectures over runtime systems[END_REF], increasing the performance of BLAS kernels.

Section 2 gives a brief background on the block-symbolic factorization for sparse direct solvers and introduces the problem when classical reordering techniques are used. Section 3 states the problem, describes our reordering strategy, and gives an upper bound of its theoretical cost RR n°8860

for the class of graphs with bounded degree, that is, graphs from real-life 2D or 3D numerical simulations. The quality of the symbolic structure and its impact on the performance of a sparse direct solver, here PaStiX [START_REF] Hénon | PaStiX: A High-Performance Parallel Direct Solver for Sparse Symmetric Denite Systems[END_REF], is studied in section 4. Finally, we conclude and present some future opportunities for this work.

Background

This section provides some background on sparse direct solvers and the associated preprocessing steps. A detailed example showing the impact of the ordering on the block-structure for the factorized matrix is also presented.

Sparse direct solvers

The common approach to sparse direct solvers is composed of four main steps: 1) ordering and computation of a supernodal partition for the unknowns; 2) block-symbolic factorization; 3) numerical block-factorization; and 4) triangular system solve. The rst step exploits the property of the Fill-in Characterization theorem (see theorem 2.1) to minimize the ll-in, zeros becoming non-zeros during factorization, and the second one predicts the block-structure that will facilitate ecient numerical factorization and solve. Steps 3) and 4) perform the actual computation.

Theorem 2.1. (Fill-in Characterization Theorem from [START_REF] Rose | Algorithmic aspects of vertex elimination in directed graphs[END_REF]) Given an n × n sparse matrix A, and its adjacency graph G = (V, E), any entry a i,j = 0 from A will become a non-zero entry in the factorized matrix if and only if there is a path in G from vertex i to vertex j that only goes through vertices with a lower index than i and j.

Among the ordering techniques known to eciently reduce the matrix ll-in are the Approximate Minimum Degree (AMD) algorithm [START_REF] Amestoy | An approximate minimum degree ordering algorithm[END_REF] and the Minimum Local Fill (MF) algorithm [START_REF] Tinney | Direct solutions of sparse network equations by optimally ordered triangular factorization[END_REF][START_REF] Luce | On the minimum ops problem in the sparse cholesky factorization[END_REF].

However, these orderings fail to expose a lot of parallelism in the block computation during the factorization. In order to both reduce ll-in and exhibit parallelism, an ordering algorithm based on nested dissection [START_REF] George | Nested dissection of a regular nite element mesh[END_REF] has been introduced and is now the most widely used in sparse direct solvers. This class of algorithms works on the symmetric undirected graph associated with the matrix and recursively partitions the graph to expose independent subproblems that can be solved in parallel while reducing the ll-in of the matrix.

The top-left part of gure 1 shows the adjacency graph of a 2D symmetric 5 × 5 grid with a possible 2-level partitioning of the graph. The goal of the nested dissection method is to recursively partition the graph G = (V, E) into A ∪ B ∪ C such that no edge directly connects a vertex from A to a vertex of B, and such that C is as small as possible. C is called the separator and corresponds to a supernode. This separation of A and B combined with theorem 2.1 guarantees that if all vertices of C are numbered with larger numbers than those of A and B, no ll-in appears between a vertex from A and a vertex from B. This partitioning and ordering process is then recursively applied on A and B, until a small enough size is reached for the subgraphs. Then, local ordering heuristics like AMD are used on these remaining subgraphs. A global supernode partition of the unknowns is obtained by merging the set of supernodes from the nested dissection process (all the separators) and the set of supernodes achieved from the reordered non-separated subgraphs (by using the algorithm introduced in [START_REF] Liu | On nding supernodes for sparse matrix computations[END_REF][START_REF] Liu | The role of elimination trees in sparse factorization[END_REF]). This partitioning and ordering operation is usually performed through an external tool such as Metis [START_REF] Karypis | A fast and high quality multilevel scheme for partitioning irregular graphs[END_REF] or Scotch [START_REF] Pellegrini | Scotch and libScotch 5.1 User's Guide[END_REF]. Adjacency graph (G). Elimination tree (T).

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
Factorized matrix (L). Given this supernodal partition, one can compute the block-symbolic data structure of the factorized matrix, as presented on the right part of gure 1. The goal is to predict the blockdata structure of the nal L matrix for the numerical factorization and to gather information in blocks that will enable the use of ecient kernels as BLAS Level 3 operations [3]. This blockdata structure is composed of N column blocks, one for each supernode of the partition, with a dense diagonal block (in gray in the gure) and with several dense o-diagonal blocks (in green or in red in the gure) corresponding to interactions between supernodes; some additional ll-in is accepted to form dense blocks in order to be more CPU-ecient. The block-symbolic factorization computes this block-data structure with Θ(N ) space and time complexities [START_REF] Charrier | Algorithmic study and complexity bounds for a nested dissection solver[END_REF]. From this structure, one can deduce the quotient graph which describes all the interactions between supernodes during the factorization (for example, supernode 1 will contribute to supernodes 3 and 7), and the elimination tree which describes the amount of parallelism in the computations as a supernode will contribute only to supernodes belonging to its ancestors. Finally, before distributing the column blocks on the processors, the biggest column blocks corresponding to the top most supernodes in the tree are split in order to exploit the parallelism inside the dense computations [START_REF] Hénon | PaStiX: A High-Performance Parallel Direct Solver for Sparse Symmetric Denite Systems[END_REF].

The rst two steps of a direct solver are preprocessing stages, independent from numerical values. Note that those steps can be computed once to solve the same problem several times with dierent numerical values. Steps 3) and 4) are numerical. Figure 2 presents how the elimination of a column block is divided into three stages:

1. Factorization of the dense diagonal block;

2. Application of an in place Solve on the o-diagonal blocks;

3. Update of the underlying matrix.

Usually, the solve stage (stage 2) is done through one or multiple calls to BLAS kernels according to the data distribution used by the solver. In the PaStiX solver, a 1D distribution is RR n°8860 used and all blocks are stored contiguously in dense storage in order to perform the solve stage in only one BLAS call. The update stage (stage 3) can be done in multiple ways. The rst option is to do it similarly to dense factorization with one matrix-matrix multiply per couple of o-diagonal blocks (red updates in gure 2). However, the granularity of the tasks in sparse solvers is often so small before reaching the top levels of the elimination tree that it is inecient.

The most adopted solution for supernodal methods is to compute a matrix-matrix multiply for each column block that requires updates, meaning one per blue o-diagonal block (only the rst two are represented in gure 2). The temporary result is then scattered and added to the target column block. The last solution, similar to what is done in a multifrontal solver, consists in a single matrix-matrix multiplication that is followed by a 2D scatter of the updates. In the last two options, if the updates are too discontinuous and spread all over the updated submatrix, this can lead to memory bound updates while the operation is originally compute bound. It is then interesting to consider an ordering solution, compatible with the nested dissection method, that will limit the number of o-diagonal blocks to have more compact updates. It will also reduce the memory bound aspect of the update operation, which is the most time consuming for the factorization and solve steps.

Intra-node reordering

Let us now illustrate the problem of current ordering solutions and how to overcome this problem.

For this purpose, we consider a regular 3D cube of n 3 vertices presented in gure 3. We apply the nested dissection process to this cube. Naturally, the rst separator, in gray, is a plane of n 2 vertices cutting the cube into two halves of balanced parts. Then, by recursively applying the nested dissection process, we partition the two-halves' subparts with the two red separators, and again dissect the resulting partitions by the four third-level green separators giving us eight nal partitions. We know from this process that each separator will be ordered with higher indices than those in lower levels.

Inside each separator, vertices have to be ordered as well, and it is common to use techniques such as the Reverse Cuthill-McKee [START_REF] George | Computer Solution of Large Sparse Positive Denite[END_REF] (RCM) algorithm in order to have an internal separator ordering as continuous as possible to limit the number of o-diagonal blocks in the associated column block. This strategy works with only the local graph induced by the separator. It starts Inria from a peripheral vertex and orders, consecutively, vertices at distance 1, then at distance 2, and so on, giving indices in reverse order. It is close to a Breadth-First Search (BFS) algorithm.

However, such an algorithm uses only interactions within a supernode, without taking into account contributing supernodes. On the quotient graph of gure 1, it means that this will reorder unknowns inside a node of this graph without considering interactions with other nodes of this graph. However, these interactions are the ones related to o-diagonal blocks in the factorized matrix. Therefore, it is important to note that the ordering inside a supernode can be rearranged to take into account interactions with vertices outside its local graph without changing the nal ll-in of the L block-structure used by the solver. Then, we can expect that a complete knowledge of the local graph and of its outer interactions will lead to a better quality in terms of the number of o-diagonal blocks.

Figure 3: Three-levels of nested dissection on a regular cube. Based on this, we propose an optimal ordering, in gure 4b, computed by hand, as opposed to an RCM algorithm, in gure 4a. This ordering is considered optimal as it minimizes the number of o-diagonal blocks to one per column. One can note that RCM will not order consecutively vertices that will receive contributions from the same supernodes, leading to a substantially larger number of o-diagonal blocks than the optimal solution. For instance, the four blue vertices in the top right of RCM ordering will RR n°8860 create four dierent o-diagonal blocks. The general idea is that some projections will be cut by RCM following the neighborhood, while those vertices could have been ordered together to reduce the number of o-diagonal blocks. On the right, the optimal ordering tries to consider this rule by ordering vertices with similar connections in a contiguous manner. This leads to a smaller number of o-diagonal blocks as shown in the block-data structure computed by the block-symbolic factorization for these two orderings in gures 5a and 5b. The ordering proposed in gure 4b is optimal in terms of number of o-diagonal blocks as long as is it impossible to exhibit a block-symbolic structure will less than 14 o-diagonal blocks: there is no more than one o-diagonal block per column block. We have demonstrated with this simple example that RCM does not fulll the correct objective in a more global view of the problem. This is especially true in the context of 3D graphs, where the separator is a 2D structure, receiving contributions from 3D structures on both sides.

With 2D graphs, the separator is a 1D structure and in such a case, RCM will provide generally a good solution by following the neighborhood in the BFS algorithm. However, it often happens that the separators found by generic tools such as Metis or Scotch are disconnected graphs making this previous statement incorrect as long as it is impossible to recover and follow the spatial neighborhood of the vertices (issued from the associated mesh).

Note that if it is quite easy to manually compute the optimal ordering on our example, it is harder in practice. Indeed, given an initial partition V = A ∪ B ∪ C, nothing guarantees that subparts A and B will be partitioned in a similar fashion, and that the resulting projection will match. For instance, gure 6 presents the projection of level-1 (in red) and level-2 (in green)

supernodes on the rst separator of a 40 × 40 × 40 laplacian partitioned with Scotch. One can note that there are crossed contributions, meaning that subparts A and B are partitioned dierently.

In the next section, we propose a new reordering strategy that permutes the rows to compact the o-diagonal information. Note that such a reordering strategy will not impact the global ll-in as long as the diagonal blocks are considered as dense blocks. The rst solution that appears on gure 6 would be to cluster vertices by common connections to nodes of the quotient graph. However, in most cases, that would result in clusters of O(1) size that would still need to be ordered correctly, taking into account their level in the elimination tree of the connected supernodes. The solution we propose to remedy this problem relies on the computed blockdata structure. Our objective is to express an algorithm providing the optimal solution before proposing a heuristic with a reasonable complexity.

Inria

Related Works

Studying the structure of o-diagonal blocks was used in dierent contexts. In [START_REF] George | On the application of the minimum degree algorithm to nite element systems[END_REF], the purpose was to reduce the overhead associated with each single o-diagonal block. The authors proposed a reordering strategy that renes the ordering provided by the minimum degree algorithm. Their experiments have been applied to 2D graphs and successfully reduce the number of o-diagonal blocks. However, the authors did not provide a theoretical study of their reordering algorithm and their solution did not apply in the context of 3D graphs.

In HSL [?, 17], the authors have introduced reordering techniques in the context of both supernodal and multifrontal solvers. The objective is to create larger o-diagonal blocks to enhance data locality and reduce factorization time in the MA87 supernodal solver. Each diagonal block is reordered according to its set of children. Given a child and one of its ancestors, section 2.2

showed that the set of rows from the ancestor connected to the child should be ordered contiguously to create a single o-diagonal block and avoid scattering operations. HSL reordering strategy starts by sorting children accordingly to their contribution size to the studied supernode, in other words with the number of rows that is connects the two nodes. Then, the rows connected to the larger child are numbered continuously to create a single o-diagonal block coming from this child, and the process is repeated on the remaining rows until all of them are reordered.

This strategy has led to performance gains in a multi-threaded context [?]. However, from construction, this algorithm gives priority to the largest branch of the elimination tree neglecting the other ones. In practice, we observe that leaves from both sides of the elimination tree might be connected to the same unknowns of their ancestors. Thus, we can expect that an algebraic view of the problem will provide better results by using a global view of the connections of one supernode to all its descendants in the elimination tree.

In [START_REF] Wissam | Scaling the solution of large sparse linear systems using multifrontal methods on hybrid shared-distributed memory architectures[END_REF], a reordering strategy is proposed for the multifrontal solver Mumps [START_REF] Amestoy | MUMPS[END_REF]. The objective is to provide a row ordering and the associated mapping on a set of processors to minimize the total volume of communications. The strategy studied minimizes communication between a parent and its set of children. As opposed to HSL and our solutions, this algorithm, designed for multifrontal solver only, dynamically reorders the rows at each level of the elimination tree.

It considers only interactions from children to direct parent to minimize scattering operations at each level when updating the frontal matrix. This way, rows from a parent node can have dierent orderings at each level of the elimination tree, or even between dierent child branches.

Thus, it is out of scope of the global ordering proposed in this paper for a supernodal solver.

RR n°8860

3 Improving the blocking size

As presented in section 2.2, the RCM algorithm widely used to order supernodes generates many extra o-diagonal blocks by not considering supernode interactions, which leads to an increased number of less ecient block operations. In this section, we present an algorithm that intends to reorder supernodes thanks to a global view of the nested dissection partition.

We expect that considering contributing supernodes will lead to a better quality a smaller number of larger blocks. Our main idea is to consider the set of contributions for each row of a supernode, before using a distance metric to minimize the creation of o-diagonal blocks when permuting rows.

Problem modeling

The strategy is to rely on the block-symbolic factorization of L instead of the original graph of A. Indeed, it allows us to take into account ll-in elements that were computed thanks to the block-symbolic factorization process instead of re-computing those elements with the matrix graph. Let us consider the th diagonal block C of the factorized matrix that corresponds to a supernode, and the set of supernodes C k with k < corresponding to the supernodes in lower levels of the elimination tree than C . Note that we refer to N as the total number of diagonal blocks appearing in the structure of the factorized matrix, as opposed to n for the total number of unknowns.

We dene for each supernode C :

row ik = 1 if vertex i from C is connected to C k 0 otherwise , k ∈ 1, -1 , i ∈ 1, |C | . ( 1
)
row ik is then equal to 1 when the vertex i, or row i, of the supernode C is connected to any vertex of the supernode k belonging to a lower level in the elimination tree. It is equal to 0, if not, meaning that no non-zero element connects the two in the initial matrix, or no ll-in will create that connection. Let's now dene for each vertex the binary vector B i = (row ik ) k∈ 1, -1 . We can then dene w i , the weight of a row i, as in equation [START_REF] Charrier | Algorithmic study and complexity bounds for a nested dissection solver[END_REF], that represents the number of supernodes contributing to that row i, and the distance between two rows i and j, d i,j , as in equation [START_REF] Dongarra | A set of level 3 basic linear algebra subprograms[END_REF]. It is known as the Hamming distance [START_REF] Hamming | Error Detecting and Error Correcting Codes[END_REF] between two binary vectors, and allows for measuring the number of o-diagonal blocks induced by the succession of two rows i and j. Indeed, d i,j represents the number of o-diagonal blocks that belongs to only one of the two rows, which can be seen as the number of blocks that end at row i or start at row j.

w i = -1 k=1 row ik , (2)
d i,j = d(B i , B j ) = -1 k=1 row ik ⊕ row jk , ( 3
)
where ⊕ is the exclusive or operation. Thus, the total number of o-diagonal blocks, odb , contributing to the diagonal block C can be dened as:

odb = 1 2 (w 1 + |C |-1 i=1 d i,i+1 + w |C | ), (4)

Inria

where the Hamming weights of the rst and last row of the supernode C correspond respectively to the number of blocks in the rst row and in the last one, and the distances between two consecutive rows gives the evolution in the number of blocks when travelling through them.

Figure 7 illustrates the computation of the number of o-diagonal blocks with the previous equation ( 4) on an example of four rows. The computation of the distance from the second row to the third one is illustrated on the left: there are two dierences making it a distance of 2. Table 7b summarizes the distances between each couple of rows in this example. With those information and the weight of the rst and last rows, respectively 3 and 2, one can compute the number of o-diagonal blocks, odb, from the formula: Thus, to reduce the total number of o-diagonal blocks in the nal structure, the goal is to minimize this metric odb for each supernode, by computing a minimal path visiting each node, with a constraint on the rst and the last node. This problem is known as the Shortest

1 2 (w 1 +d 1,2 +d 2,3 +d 3,4 +w 4 ) = 1 2 (3+3+2+2+2) = 6.
Hamiltonian Path Problem, and is an NP-hard problem.

Proposed heuristic

We rst propose to introduce an extra virtual vertex, S 0 , for which B 0 is the null set. Thus, we have:

∀i ∈ 1, |C | , d 0,i = d i,0 = w i , (5)
The problem can now be transformed in a Travelling Salesman Problem [START_REF] Applegate | The Traveling Salesman Problem: A Computational Study[END_REF] (TSP):

|C | i=0 d i,(i+1) , (6)
which is also an NP-Hard problem, but for which multiple heuristics have been proposed in the literature [START_REF] Johnson | The Traveling Salesman Problem: A Case Study in Local Optimization[END_REF], as opposed to Shortest Hamiltonian Path Problem. Furthermore, our problem presents properties that make it suitable for better heuristics and theoretical models that guarantee the maximum distance to the optimal solution. Firstly, our problem is symmetric as:

d ij = d ji , ∀(i, j) ∈ 1, |C | 2 , (7)
and secondly, respects the triangular inequality:

d ij ≤ d ik + d kj , ∀(i, j, k) ∈ 1, |C | 3 . (8)
This sets our problem as an Euclidean TSP, and so heuristics for these specic cases can be used.

Dierent TSP heuristics that can be used to solve this problem, with their respective cost and quality with respect to the optimal, are presented in table 1.

RR n°8860

To keep a global complexity below that of the numerical factorization, we explain in the following section 3.3 that the complexity of the TSP algorithm has to remain equal or lower to Θ(p 2 ), where p is the number of vertices in the cycle. Thus, it prevents advanced algorithms such as the Christodes algorithm [START_REF] Christodes | Worst-case analysis of a new heuristic for the travelling salesman problem[END_REF] from being used. Furthermore, as p might reach several hundred or more, the use of Nearest neighbor or Clarke and Wright heuristics might provide low quality results. From the remaining options, we decided to use the nearest insertion method which is a quadratic algorithm and guaranties a maximal distance to the optimal of 2 [START_REF] Rosenkrantz | An analysis of several heuristics for the traveling salesman problem[END_REF]. A quality comparison of our algorithm over 3D Laplacian matrices and real matrices against the Concorde [25] TSP solver that returns optimal solutions have shown that our nearest insertion algorithm provides results in less than 10% from the optimal.

Algorithm Complexity for p nodes Quality (wrt optimal) Nearest neighbor

Θ(p 2 ) 1 2 (1 + log(p)) Nearest insertion Θ(p 2 ) 2
Clarke and Wright

Θ(p 2 log(p)) Θ(log(p))
Cheapest insertion

Θ(p 2 log(p)) 2
Minimum spanning tree

Θ(p 2 ) 2
Christodes

Θ(p 3 ) 1.5
Table 1: Complexity and quality of dierent TSP algorithms.

Our nal algorithm is then decomposed in three stages presented in algorithm 1 that are applied to each separator of the nested dissection. Note that it is not applied on the leaves of the elimination tree since they will not receive contributions from other supernodes. The rst step is to compute the B i vectors for each row i of the current separator. Then, it computes the distances matrix of the separator: D = (d i,j ) (i,j)∈ 0,|C | 2 . Finally, the TSP algorithm is executed using this matrix to produce the local ordering of the supernode that minimizes the equation ( 6).

Inria

The rst stage of this algorithm builds the vector B i for each row i. In fact, to minimize the storage, only contributing supernodes (row ik = 1) are stored for B i . In order to do so, we rely on the structure of the block-symbolic factorization that provides a compressed storage of the information similar to compressed sparse row (CSR) format. Given a supernode, one can easily access the o-diagonal blocks contributing to this supernode, and due to the sparse property, the number of these blocks is much smaller than ( -1). The accumulated operations for all the supernodes in the matrix is in Θ(n). Note that we store the contributing supernode numbers in an ordered fashion for faster computation of the distances. Furthermore, the memory overhead of this operation is limited by the fact that each supernode is treated independently.

The second stage computes the distance matrix. When computing the distance d ij between rows i and j from C , we take advantage of the sorted sets B i and B j to realize this computation in Θ(

|B i | + |B j |) operations.
The third stage executes the nearest insertion heuristics to solve the TSP problems on the vertices of the supernode based on the previously computed distance matrix. As stated previously, this step is computed in Θ(|C | 2 ) operations. It is known that the solution given is not optimal but will be at a distance 2 of the optimal in the worst case. Figure 8 presents the block-symbolic factorization of a 3D Laplacian of size 8×8×8 reordered with the Scotch nested dissection algorithm. In gure 8a, our reordering algorithm has not been applied, and supernode ordering results only from the local RCM applied by Scotch. One can notice that some rows can be easily aggregated to reduce the number of o-diagonal blocks.

In gure 8b, our algorithm has to reorder unknowns within each supernode. The nal structure exhibits more compact blocks that are larger. Note that the ll-in of the matrix has not changed due to the dense storage of the diagonal blocks. Our algorithm does not impact the ll-in outside those diagonal blocks.

Complexity study

For this study we consider graphs issued from nite element mesh coming from real-life simulations of 2D or 3D physical problems. From a theoretical point of view, the majority of those graphs have a bounded degree and are specic cases of bounded-density graphs [START_REF] Miller | Density graphs and separators[END_REF]. In this RR n°8860 section, we provide a complexity study of our reordering algorithm in the context of a nested dissection partitioning strategy for this class of graphs.

Good separators can be built for bounded-density graphs or more generally for overlap graphs [START_REF] Miller | A unied geometric approach to graph separators[END_REF]. In d-dimension, such n-node graphs have separators whose size grows as Θ(n (d-1)/d ).

In this study, we consider the general framework of separator theorems introduced by Lipton

and Tarjan [START_REF] Lipton | A separator theorem for planar graphs[END_REF] for which we will have σ = d-1 d .

Denition 3.1. A class ϕ of graphs satises an n σ -separator theorem, 1 2 ≤ σ < 1, if there are constants 1 2 ≤ α < 1, β > 0 for which any n-vertex graph in ϕ has the following property: the vertices of G can be partitioned into three sets A, B, and C such that: no vertex in A is adjacent to any vertex in B, |A| ≤ αn, |B| ≤ αn, and C ≤ βn σ where C is the separator of G.

Theorem 3.2. (From [START_REF] Charrier | Algorithmic study and complexity bounds for a nested dissection solver[END_REF]) The number of o-diagonal rows in the block data structure for the factorized matrix L is at most Θ(n).

This result comes from [START_REF] Charrier | Algorithmic study and complexity bounds for a nested dissection solver[END_REF]. In this paper, the authors demonstrated that the number of o-diagonal blocks is at most Θ(n) and this was achieved by proving that this upper bound is in fact true for the total number of rows inside the o-diagonal blocks, leading to theorem 3.2.

Using this theorem, we demonstrate theorem 3.3. Theorem 3.3. (Reordering Complexity) For a graph of bounded degree satisfying a n σ -separation theorem, the reordering algorithm complexity is bounded by Θ(n σ+1 ).

Proof

The main cost of the reordering algorithm is issued from the distance matrix computation. As presented in section 3.2, we compute a distance matrix for each supernode. This matrix is of size |C |, and each element of the matrix, D , is the distance between two rows of the supernode.

The overall complexity is then given by:

C = N =1 |C | i=1 (row ik ) k∈ 1, -1 × (|C | -1). (9)
More precisely, for a supernode C , the complexity is given by the number of o-diagonal rows that contribute to it multiplied by the number of comparisons: (|C | -1). For instance, given gure 5a, one can note that the complexity will be proportional to the colored surface (blue, green, and red blocks), where row ik = 1, as well as in the number of rows. Using the compressed sparse information (colored blocks) only instead of the dense matrix is important for reaching a reasonable theoretical complexity, as long as this number of o-diagonal blocks is bounded in the context of nite element graphs. Given theorem 3.2, we know that the number of o-diagonal contributing rows in the complete matrix L is in Θ(n). In addition, the largest separator is asymptotically smaller than the Inria maximum size of the rst separator that is Θ(n σ ). The complexity is then bounded by:

C ≤ max 1≤ ≤N (|C | -1) Θ(n σ ) × N =1 ( |C | i=1 (row ik ) k∈ 1, -1 ) Θ(n) = Θ(n σ+1 ).
For graphs of bounded degree, this result leads to: for the graph family admitting an n 1 2 -separation theorem (2D meshes), the reordering cost is bounded by Θ(n √ n), and is at worst as costly as the numerical factorization;

for the graph family admitting an n 2 3 -separation theorem (3D meshes), the reordering cost is bounded by Θ(n

3 ), and is cheaper than the numerical factorization, which grows as Θ(n 2 ).

Analysis

Note that this complexity is as said before larger than the complexity of the TSP nearest insertion heuristic. For a subgraph of size p respecting the p σ -separation theorem, this heuristic complexity is in Θ(p 2σ ). Using [START_REF] Charrier | Algorithmic study and complexity bounds for a nested dissection solver[END_REF], we can compute the overall complexity as a recursive function depending on the complexity on one supernode. It leads to an overall complexity in Θ(n log(n)) for 2D graphs and Θ(n 3 ) for 3D graphs, and is then less expensive than the complexity of computing the distance matrix.

The reordering is as costly as the numerical factorization for 2D meshes, but RCM is usually giving a good ordering on 2D graphs, as long as the separators are contiguous lines. For the 3D cases, the reordering strategy is cheaper than the numerical factorization. Thus, this reordering strategy is interesting for any graph with 1 2 < σ < 1, including graphs with a structure between 2D and 3D meshes. This algorithm can easily be parallelized since each supernode is an independent subproblem, and the distance matrix computation can also be computed in parallel. Thus, the sequential cost of this reordering step can be lowered and should be negligible compared to the numerical factorization.

Figure 9 presents the complexity study on 3D Laplacian matrices. We computed the practical complexity of our reordering algorithm with respect to the upper bound we demonstrated. The red curve presents the sequential time taken by our reordering algorithm. It is compared to the theoretical complexity demonstrated previously, but scaled to match on the middle point (size 150 3 ) to ease the read and check, so we can conrm that the trends of both curves are identical to a constant factor. Finally, in green we also plotted the practical complexity: total number of comparisons performed during our reordering algorithm, to see if the theoretical complexity was of the same order. This curve is also scaled to match on the middle point. One can note that the three curves are quite close, which conrms that we found a good upper bound complexity for a large set of sizes.

Note that this complexity seems to be signicant with respect to the factorization complexity.

Nevertheless, the nature of operations (simple comparisons between integers) is much cheaper than the numerical factorization operations. In addition, if we use the partitioner to obtain large enough supernodes, it will reduce by a notable factor the complexity of our algorithm, as long as we operate on a column block and not on each element contributing to each row. This parameter can be set in ordering tools as Metis and Scotch, and has an impact on the global ll-in of the matrix. As presented before, the reordering stage takes part of preprocessing steps, and can be used for many numerical steps and it enhances both factorization and solve steps.

Strategies to reduce computational cost

As we have seen, the total complexity of the reordering step can still reach the complexity of the numerical factorization. We now introduce heuristics to reduce the computational cost of our reordering algorithm.

Multi-Level: partial computation based on the elimination tree

As presented in section 2, the obtained partition allows us to decompose contributing supernodes according to the elimination tree. With the characterization theorem, we know that when we consider one row, if this row receives a contribution from a supernode in the lowest levels, then it will receive contributions from all its descendants to this node. This helps us divide our distance computation into a two dimension distance d i,j = d high i,j +d low i,j to reduce its cost. Given a split level parameter, we rst compute the high-levels distance, d high i,j , by considering only the contributions from the supernodes in the split level levels directly below the studied supernode. This distance gives us a minimum of the distance between two rows. Indeed, if considering all supernodes, the distance will be necessarily equal or larger to the high-levels distance by construction of the elimination tree. Then, we only compute the low-levels distance, d low i,j , only if the rst one is equal to 0.

In practice, we observed that for a graph of bounded degree, not especially regular, a ratio of 3 to 5 between the number of lower and upper supernodes largely reduces the number of complete distances computed while conserving a good quality in the results. The split level parameter is then adjusted to match this ratio according to the part of the elimination tree considered. It is important to notice that it is impossible to consider the distances level by level, since the goal here is to group together the rows which are connected to the same set of leaves in the elimination tree. It means that they will receive contributions from nodes on identical paths in this tree. The partial distances consider only the beginning of those paths and not their potential reconnection Inria further down the tree. That is why it is important to take multiple levels at once to keep a good quality.

Stopping criteria: partial computation based on distances

The second idea we used to reduce the cost of our reordering techniques is to stop the computation of a distance if it exceeds a threshold parameter. This solution helps to quickly disregard the rows that are far-away from each other. This limits the complexity of the distance computation, reducing the overall practical complexity. In most cases, a small value such as 10 can already provide good quality improvement. However, it depends on the graph properties and the average number of dierences between rows. Unfortunately, if this heuristic is used alone, this improvement is not always guaranteed and it might lead to a quality worsening. In association with the previous multi-level heuristic, the results are always improved, as we will see in the following section.

Experimental study

In this section, we present experiments with our reordering strategy, both in terms of quality (number of o-diagonal blocks) and impact on the performance of numerical factorization. We compare here three dierent strategies to the original ordering provided by the Scotch library.

Two are based on our strategy, namely TSP, with the full distance computation or with the multilevel approximation. The last one, namely HSL, is the one implemented in the HSL library for the MA87 supernodal solver (optimize_locality routine from MA87). Note that those three reordering strategies are applied to the partition found by Scotch, and they do not modify the ll-in.

Context

We used a set of large matrices arising from real-life applications originating from the University of Florida's [START_REF] Davis | The university of orida sparse matrix collection[END_REF] Sparse Matrix Collection. For that experiment, we took all matrices from this collection with a size between 500.000 and 10.000.000. From this large set, we extracted matrices that are applicants for solving linear systems. Thus, we remove matrices originating from the Web and from DNA problems. This nal set is composed of 104 matrices, sorted by families. We also conduct some experiments with a matrix of 10 7 unknowns, taken from a CEA simulation, an industrial partner in the context of the PaStiX project.

We utilized the Plafrim 1 supercomputer for our experiments. For the performance experiments on an heterogeneous environment, we used the mirage cluster, where nodes are composed of two Intel Westmere Xeon X5650 hexa-core CPUs running at 2.67 GHz with 36 GB of memory, and enhanced by three Nvidia GPUs, M2070. We used Intel MKL 2016.0.0 for the BLAS kernels on the CPUs, and we used the Nvidia Cuda 7.5 development kit to compile the GPU kernels. For the scalability experiments in a multi-threaded context, we used the miriel cluster.

Each node is equipped with two Intel Xeon E5-2680 v3 runtime system and presented in [START_REF] Lacoste | Taking Advantage of Hybrid Systems for Sparse Direct Solvers via Task-Based Runtimes[END_REF].

For the initial ordering step, we used Scotch 5.1.11 with the congurable strategy string from PaStiX to set the minimal size of non separated sub-graphs, cmin, to be 20 as in [START_REF] Lacoste | Scheduling and memory optimizations for sparse direct solver on multicore/multi-gpu cluster systems[END_REF]. We 1 https://plafrim.bordeaux.inria.fr RR n°8860 also set the f rat parameter to 0.08, meaning that columns aggregation is allowed by Scotch as long as the ll-in introduced does not exceed 8% of the original matrix. It is important to use such parameters to increase the width of the column blocks and reach a good level of performance using accelerators. Even if it increases the ll-in, the nal performance gain is usually more important and makes the memory overhead induced by the extra ll-in acceptable.

Reordering quality and time

First, we study the quality and the computational time of the three reordering algorithms, the two versions of our TSP and HSL, compared to the original ordering computed by Scotch that is known to be in Θ(n × log(n)). Note that sequential implementation is used for all algorithms, except in subsection Parallelism.

For the quality criteria, the metric we use is the number of o-diagonal blocks in the matrix.

We always use Scotch to provide the initial partition and ordering of the matrix, thus the number of o-diagonal blocks only reects the impact of the reordering strategy. Another related metric we could use is the number of o-diagonal blocks per column block. In ideal cases, it would be, respectively, 4 and 6 for 2D and 3D meshes. However, since the partition computed by scotch is not based on the geometry, this optimum is never reached and varies a lot from one matrix to another, so we stayed with the global number of o-diagonal blocks and its evolution compared to the original solution given by Scotch.

Quality

Figure 10 presents the quality of reordering strategies in terms of the number of o-diagonal blocks with respect to the Scotch ordering. We recall that Scotch uses RCM to order unknowns within each supernode. Three metrics are represented: one with HSL reordering, one for our multi-level heuristic, and nally one for the full distance computation heuristic. We can see that our algorithm reduces the number of o-diagonal blocks on all test cases. On the 3D Inria problems, our reordering strategy improves the metric by 50 to 60%, while on the 2D problems, the improvement is of 20 to 30%. Furthermore, we can observe that the multi-level heuristic does not signicantly impact the quality of the ordering. It only reduces it by a few percent in 11 cases over the 104 matrices tested, while giving the same quality in all other cases. The HSL heuristic improves the initial ordering in average by approximately 20%, and up to 40%, but is outperformed by our TSP heuristic regardless of the multi-level distances approximation on all cases.

In addition, we observed that for matrices not issued from meshes with an unbalanced elimination tree, and not presented here, using the multi-level heuristic can deteriorate the solution.

Indeed, in this case, the multi-level heuristic is unable to distinguish close interactions from far interactions with only the rst levels, leading to incorrect choices.

Time

Figure 11 presents the cost of the three reordering strategies in sequential with respect to the cost of the initial ordering performed by Scotch. The reordering is in fact an extra step in the preprocessing stage of sparse direct solvers. One can note that despite the higher theoretical complexity, the reordering step of our TSP heuristic is 2 to 10 times faster than Scotch. Thus, adding the reordering step creates an overhead of no more than 10% to 50% in most cases when called in sequential. However, on specic matrix structures, with a lot of connections to the last supernode, the reordering operation can be twice as expensive as Scotch. In those cases, the overhead is largely diminished by the multi-level heuristic, which reduces the time of the reordering step to the same order as Scotch. We observe that the multi-level heuristic is always benecial to the computational time. For the second matrix an optimization problem with a huge density on the left of the gures , we can observe a quality gain of more than 95%, while the cost is more than 600 times larger than the ordering time. This problem illustrates the limitation of our heuristic using a global view compare to the local heuristic of the HSL algorithm that is still faster than the Scotch ordering but gives only 20% improvement. This problem is typically not suited for sparse linear solver, due to its large number of non-zeros as RR n°8860

well as its consequent ll-in.

Strategy

Stopping criteria

Table 2 shows the impact of the stopping criteria on a large test case issued from a 10 million unknowns matrix from the CEA. The rst line presents the results without reordering and the time of the Scotch step. We compared this to the number of o-diagonal blocks and the time obtained with our reordering algorithm when using dierent heuristics. The STOP parameter refers to the criteria introduced in section 3.4 and denes after how many dierences a distance computation must be stopped. One can notice that with all congurations the quality is within 39% to 42% of the original, which means that those heuristics have a low impact on the quality of the result. However, this can have a large impact on the time to solution, since a small STOP criterion combined with the multi-level heuristic can divide the computational time by more than

2.

In conclusion, we can say that for a large set of sparse matrices, we obtain a resulting number of o-diagonal blocks between two and three times smaller than the original Scotch RCM ordering, while the HSL heuristic reduces them in average only by a fth. It is interesting as it should reduce by the same factor the overhead associated to the tasks management in runtimes, and should improve the kernel eciency of the solver. Up to our experiments, we reach a practical complexity close to Scotch ordering process, leading to a preprocessing stage that is not too costly compared to the numerical factorization. Furthermore, it should accelerate the numerical factorization and solve steps to hide this extra cost when only one numerical step is made, and give some global improvement when multiple factorizations or solves are performed. While HSL reordering overhead might be much smaller than our heuristic, we hope that the dierence in the quality gain, as well as the fact that our strategy improves all children instead of giving advantage to the largest one, will benet the factorization step by a larger factor.

Parallelism

As previously stated, the reordering algorithm is largely parallel as each supernode can be reordered independently of the others. The rst level of parallelism is a dynamic bin-packing that distributes supernodes in reverse order of their sizes. However, some supernodes are too large and take too long to be reordered compared to all others. They represent almost all the computational requirements. We then divided the set of supernodes into two parts. For the smaller set, we just reorder dierent supernodes in parallel, and for the larger set, we parallelize the distance matrix computation. Figure 12 shows the speedup obtained with 24 threads over the best sequential version on a miriel node. This simple parallelization accelerates the algorithm by 10 on average and helps to totally hide the cost of the reordering step in a multi-threaded context, where ordering tools are hard to parallelize. Note that for many matrices, the parallel implementation of our reordering strategy has an execution time smaller than 1s. In a few cases, the speed-up is still limited to 5 because the TSP problem on the largest supernode remains sequential and may represent a large part of the sequential execution.

4.3

Impact on supernodal method: PaStiX

In this section, we measure the performance gain brought by the reordering strategies. For these experiments, we extracted 6 matrices from the previous collection, and we use the number of operations (Flops) that a scalar algorithm would require to factorize the matrix with the ordering returned by Scotch to compute the performance of our solver. We recall that this number is stable with all reordering heuristics.

Figure 13 presents the performance on a single mirage node, for three algorithms based on the original ordering from Scotch. The rst one leaves the Scotch ordering untouched. The HSL heuristic is applied on the second one, and nally the third one includes the TSP heuristic with full distances. For each matrix and each ordering, scalability of the numerical factorization is presented with all the 12 cores of the architecture enhanced by 0 to 3 GPUs. All results are an average performance on ve runs.

To explain the dierent performance gains, we rely on table 3, which presents the average number of rows in the o-diagonal blocks with and without reordering, and not the total number of blocks to give an insight on the size of the updates. The block width is dened by the Scotch partition and is the same for all experiments on each matrix. This number is important, as it is especially benecial to enlarge blocks when the original solution provides small data blocks.

For multi-threaded runs, both reordering strategies give a slight benet up to 7% on the performance. Indeed, on the selected matrices, the original o-diagonal block height is already large enough to get a good CPU eciency since the original solver already runs at up to 67% of the theoretical peak of the node (128.16 GFlop/s). It is also true for HSL reordering. In general, when the solver exploits GPUs, the benet is more important and can reach up to 20%.

In gure 13, we can see that with the afshell10 matrix, extracted from a 2D application, reordering strategies have a low impact on the performance, and the accelerators are also not helpful for this lower computation case. For the Flan1565 matrix, the gain is not important for If we compare with HSL reordering strategy, we can see that our reordering is helpful to improve slightly the performance of the solver. Hence, choosing between both strategies depends on the number of factorizations that are performed.

Table 3 also presents our TSP reordering strategy overhead with the parallel implementation and the resulting gain on the numerical factorization time using the mirage architecture with the three GPUs. Those numbers reect that it is interesting to use our reordering strategy in many cases with a small overhead that is immediately recovered by the performance improvement of the numerical factorization. As long as several problems presenting the same structure are solved, this small overhead is again diminished. Similarly, if GPUs are used, the gain during the factorization is higher, completely hiding the overhead of the reordering. The cost of the HSL strategy being really small with respect to Scotch, it is always recommended to apply it for a single factorization or for homogeneous computations. However, if GPUs are involved, HSL reordering impact on the performance of the numerical factorization is really slight and is goes from slight slow down to slight speed up (around 3%). This validates the use of more complex heuristics as the proposed TSP. Figure 14 presents a scalability study on one miriel node with 24 threads, with and without our reordering stage on the 10 million unknowns matrix from the CEA. This matrix, despite being a large 3D problem, presents a really small average block size of less than 5 when no reordering is applied. The reordering algorithm raises up to 12.5, explaining the larger average gain of 8 -10% that is observed. In both cases, we notice that the solver manages to scale correctly over the 24 threads, and even a little better when the reordering is applied. A slight drop in the performance on 14 threads is explained by the overow on the second socket.

Conclusion

We presented a new reordering strategy, that according to our experiments succeeds in reducing the number of o-diagonal blocks in the block-symbolic factorization. It allows one to signicantly improve the performance of GPU kernels, and the BLAS CPU kernels in smaller RR n°8860 ratios, as well as reduces the number of tasks when using a runtime system. The resulting gain can be up to 20% on heterogeneous architectures enhanced by Nvidia Fermi architectures. Such an improvement is signicant, as long as it is dicult to reach a good level of performance with sparse algebra on accelerators. This gain can be observed on both the factorization and the solve steps. It is particularly well working for graphs issued from nite element meshes of 3D problems.

In the context of 2D graphs, partitioner tools can be sucient, as long as separators are close to 1D structures and can easily be ordered by following the neighborhood. For other problems, the strategy enhances the number of o-diagonal blocks, but might be costly on graphs where vertices have large degrees.

Furthermore, we proposed a parallel implementation of our reordering strategy, leading to a computational cost that is really low with respect to the numerical factorization and that is counterbalanced by the gain on the factorization. In addition, if multiple factorizations are applied on the same structure, this benets the multiple factorization and solve steps at no extra cost. We proved that such a preprocessing stage is cheap in the context of 3D graphs of bounded degree, and showed that it works well for a large set of matrices. We compared with HSL reordering, which targets the same objective of reducing the overall number of o-diagonal blocks. While the TSP heuristic is often more expensive, the quality is always improved, leading to better performance. In the context of multiple factorizations, or when using GPUs, the TSP overhead is recovered by performance improvement, while it may be better to use HSL for the other cases.

For future work, we plan to study the impact of our reordering strategy in a multifrontal context with the Mumps [START_REF] Amestoy | MUMPS[END_REF] solver and compare it with the solution studied in [START_REF] Wissam | Scaling the solution of large sparse linear systems using multifrontal methods on hybrid shared-distributed memory architectures[END_REF] that performs the permutation during the factorization. The main dierence with the static ordering heuristics studied in this paper, is that the MUMPS heuristic is applied dynamically at each level of the elimination tree. Such a reordering technique is also important in the objective of integrating variable size batched operations currently under development for the modern GPU architectures.

Finally, one of the most important perspectives is to exploit this result to guide matrix compression methods in diagonal blocks for using hierarchical matrices in sparse direct solvers. Indeed, considering the diagonal block by itself for compression without external contributions leads to incorrect compression schemes. Using the reordering algorithms to guide the compression helps to gather contributions corresponding to similar far or close interactions.

Inria

Figure 1 :

 1 Figure 1: Nested dissection and block-data structure for L.

Figure 2 :

 2 Figure 2: Steps to factorize the matrix -symmetric positive denite case.

Figure 4 :

 4 Figure 4: Projection of contributing supernodes and ordering on the rst separator (gray in gure 3).

Figure 4

 4 Figure 4 presents the vertices of the gray separator from the 3D cube case with n = 5. The projection of contributing supernodes on this separator is shown. The blue parts are the vertices connected only to the leaves of the elimination tree. Thanks to the nested dissection process, the nodes of the gray separator have the largest numbers and their connections to other supernodes represent the o-diagonal contributions. Based on this, we propose an optimal ordering, in

Figure 5 :

 5 Figure 5: O-diagonal blocks contributing to the rst separator in gure 3.

Figure 6 :

 6 Figure 6: Projection of contributing supernodes on the rst separator of a 3D laplacian of size 40 × 40 × 40, using Scotch.

Figure 7 :

 7 Figure 7: Example of a symbolic structure and its associated distance matrix. Computation of the distance between rows 2 and 3 is illustrated.

Figure 8 :

 8 Figure 8: Block-symbolic factorization of 8 × 8 × 8 Laplacian initially reordered with Scotch.

4

 4

Figure 9 :

 9 Figure 9: Comparison of the time of reordering against theoretical and practical complexities on 3D Laplacians.

  after reordering / Number of blocks with Scotch HSL reordering TSP with multi-level distances TSP with full distances

Figure 10 :

 10 Figure 10: Impact of the heuristic used on the ratio of o-diagonal blocks over those produced by the initial Scotch ordering on the Florida set of matrices. The lower the better.

Figure 11 :

 11 Figure 11: Time of the sequential reordering step with respect to the initial Scotch ordering on the Florida set of matrices. The lower the better.

Figure 12 :

 12 Figure 12: Speedup of the reordering step (full distance heuristic) with 24 threads on the full set of matrices.

Figure 14 :

 14 Figure 14: Scalability on the CEA 10 million unknowns matrix with 24 threads.

  Reordering algorithm for each supernode C in the elimination tree do for each row i in the supernode C do for each contributing node k ∈ 1, -1 doSet row ik to 1 Build the structure B i end for end for for each row i in the supernode C do for each row j in the supernode C do Compute the distance between rows i and j

		Compute the distances
	end for end for
	Cycle = {S 0 , 1} for i ∈ 2, |C | do
	Insert row i in Cycle such that (6) is minimized end for	Order rows
	Split Cycle at S 0 end for

Algorithm 1

  12-cores running at 2.50 GHz and 128 GB of memory. The same version of the Intel MKL is used. The PaStiX version used for our experiments is the one implemented on top of the Parsec [?]

Table 2 :

 2 Number of o-diagonal blocks and reordering times on the CEA 10 million unknowns matrix.

		Number of blocks	Time (s)
		Full	Multi-level	Full	Multi-level
	No reordering	9760700	360
	Reordering / Stop= 10	4100616	4095986	33.2	31.1
	Reordering / Stop= 20	3896248	3897179	42.6	38.5
	Reordering / Stop= 30	3891210	3891262	50.7	43.3
	Reordering / Stop= 40	3891803	3891962	58.1	46.3
	Reordering / Stop= ∞	3891825	3892522	64.8	47.7

  Figure 13: Performance impact of the reordering algorithms on the PaStiX solver on top of the Parsec runtime with 1 node of the hybrid mirage architecture.

		Scotch 12 threads HSL 12 threads TSP 12 threads	+ 1 GPU + 1 GPU + 1 GPU	+ 2 GPUs + 2 GPUs + 2 GPUs	+ 3 GPUs + 3 GPUs + 3 GPUs
	Matrix afshell10	Avg. number of rows of o-diag. block Scotch TSP HSL 46.05 45.52 54.09	TSP times Overhead Gain 0.133 s 0 s
	FilterV2	8.794	11.33	19.91	0.164 s	1.23 s
	Flan1565	29.13	32.33	62.40	0.644 s	0.52 s
	audi	17.94	20.57	41.76	0.748 s	2.08 s
	MHD	16.86	17.04	27.64	1.16 s	4.42 s
	Geo1438	18.79	23.17	49.74	1.78 s	7.48 s

Table 3 :

 3 Impact of the reordering strategies on the number of rows per o-diagonal block and on the timings. The timings shows the overhead of the parallel TSP, and the gain it provides on the factorization step using the mirage architecture. Inria both reordering strategies because the original o-diagonal block height is already large enough for eciency. On other problems, issued from 3D applications, we observe signicant gains from 10 to 20% which reect the block height increase of 1.5 to 2.5 presented in table 3.

Pichon & Faverge & Ramet & Roman

RR n°8860

Acknowledgments

This material is based upon work supported by Bordeaux INP and the DGA under a DGA/Inria grant. Experiments presented in this paper were carried out using the PLAFRIM experimental platform.


Daniel Barlet

published or not. The documents may come

Note on the semi-continuity of the algebraic dimension


Egon W Stemle
email: egon.stemle@eurac.edu

2015 -A minimally-deep learning PoS-tagger (trained for German CMC and Web data)

This article describes the system that participated in the Part-of-speech tagging subtask of the EmpiriST 2015 shared task on automatic linguistic annotation of computer-mediated communication / social media.

The system combines a small assertion of trending techniques, which implement matured methods, from NLP and ML to achieve competitive results on PoS tagging of German CMC and Web corpus data; in particular, the system uses word embeddings and character-level representations of word beginnings and endings in a LSTM RNN architecture. Labelled data (Tiger v2.2 and EmpiriST) and unlabelled data (German Wikipedia) were used for training.

Introduction

Part-of-speech (PoS) tagging is an essential processing stage for virtually all NLP applications. Subsequent tasks, like parsing, named-entity recognition, event detection, and machine translation, often utilise PoS tags, and benefit (directly or indirectly) from accurate tag sequences. However, frequent phenomena in computer-mediated communication (CMC) and Web corpora such as emoticons, acronyms, interaction words, iteration of letters, graphostylistics, shortenings, addressing terms, spelling variations, and boilerplate [START_REF] Jannis | Neue Medien -neue Schriftlichkeit?[END_REF]Bernardini et al., 2008;[START_REF] Beißwenger | Das Dortmunder Chat-Korpus: ein annotiertes Korpus zur Sprachverwendung und sprachlichen Variation in der deutschsprachigen Chat-Kommunikation[END_REF] deteriorate the performance of PoS-taggers [START_REF] Giesbrecht | Is Part-of-Speech Tagging a Solved Task? An Evaluation of POS Taggers for the German Web as Corpus[END_REF]Baldwin et al., 2013).

To this end, the EmpiriST shared task (ST) invited developers of NLP applications to adapt their tokenisation and PoS tagging tools and resources for the processing of written German CMC and Web data (Beißwenger et al., to appear). The ST was divided into two subtasks, tokenisation and PoS tagging, and for each subtask two data sets were provided (see Subsection 4.1.3). The systems were evaluated by the organisers on raw data for the tokenisation subtask, and on unlabelled but pre-tokenised data for the PoS tagging subtask (both on the same approx. 14,000 tokens).

We participated in the PoS tagging subtask of the ST with our new minimally-deep learning PoS-tagger: We combine word2vec (w2v) word embeddings (WEs) with a single-layer Long Short Term Memory (LSTM) recurrent neural network (RNN) architecture; strictly speaking, w2v is shallow. Therefore we call the combination with a single hidden layer minimally-deep. The sequence of unlabelled w2v representations of words is accompanied by the sequence of n-grams of the word beginnings and endings, and is fed into the RNN which in turn predicts PoS labels.

The paper is organised as follows: We present our system design in Section 2, the implementation in Section 3, and its evaluation in Section 4. Section 5 concludes with an outlook on possible implementation improvements.

Design

Overall, our design takes inspiration from as far back as Benello et al. (1989) who used four preceding words and one following word in a feedforward neural network with backpropagation for PoS tagging, builds upon the strong foundation laid down by [START_REF] Collobert | [END_REF] for a NN architecture and learning algorithm that can be applied to various natural language processing tasks, and ultimately is a variation of Nogueira dos Santos and Zadrozny (2014) who trained a NN for PoS tagging, with character-level and WE representations of words.

Word Embeddings

Recently, state-of-the-art results on various linguistic tasks were accomplished by architectures using neural-network based WEs. [START_REF] Baroni | [END_REF] conducted a set of experiments comparing the popular w2v (Mikolov et al., 2013a;Mikolov et al., 2013b) implementation for creating WEs to other distributional methods with state-of-the-art results across various (semantic) tasks. These results suggest that the word embeddings substantially outperform the other architectures on semantic similarity and analogy detection tasks. Subsequently, [START_REF] Levy | [END_REF] conducted a comprehensive set of experiments and comparisons that suggest that much of the improved results are due to the system design and parameter optimizations, rather than the selected method. They conclude that "there does not seem to be a consistent significant advantage to one approach over the other".

Word embeddings provide high-quality low dimensional vector representations of words from large corpora of unlabelled data, and the representations, typically computed using NNs, encode many linguistic regularities and patterns (Mikolov et al., 2013b).

Character-Level Sub-Word Information

The morphology of a word is opaque to WEs, and the relatedness of the meaning of a lemma's different word forms, i.e. its different string representations, is not systematically encoded. This means that in morphologically rich languages with longtailed frequency distributions, even some WE representations for word forms of common lemmata may become very poor [START_REF] Kim | [END_REF].

We agree with Nogueira dos Santos and Zadrozny (2014) and [START_REF] Kim | [END_REF] that subword information is very important for PoS tagging, and therefore we augment the WE representations with character-level representations of the word beginnings and endings; thereby, we also stay language agnostic-at least, as much as possible-by avoiding the need for, often language specific, morphological pre-processing.

Recurrent Neural Network Layer

Language Models are a central part of NLP. They are used to place distributions over word sequences that encode systematic structural properties of the sample of linguistic content they are built from, and can then be used on novel content, e.g. to rank it or predict some feature on it. For a detailed overview on language modelling research see [START_REF] Mikolov | Statistical Language Models Based on Neural Networks[END_REF].

A straight-forward approach to incorporate WEs into feature-based language models is to use the embeddings' vector representations as features. Having said that, WEs are also used in neural network architectures, where they constitute (part of) the input to the network.

Neural networks (NNs) consist of a large number of simple, highly interconnected processing nodes in an architecture loosely inspired by the structure of the cerebral cortex of the brain (O'Reilly and [START_REF] Munakata | Computational Explorations in Cognitive Neuroscience Understanding the Mind by Simulating the Brain[END_REF]. The nodes receive weighted inputs through these connections and fire according to their individual thresholds of their shared activation function. A firing node passes on an activation to all successive connected nodes. During learning the input is propagated through the network and the output is compared to the desired output. Then, the weights of the connections (and the thresholds) are adjusted stepwise so as to more closely resemble a configuration that would produce the desired output. After all input cases have been presented, the process typically starts over again, and the output values will usually be closer to the correct values.

RNNs are NNs where the connections between the elements are directed cycles, i.e. the networks have loops, and this enables them to model sequential dependencies of the input. However, regular RNNs have fundamental difficulties learning long-term dependencies, and special kinds of RNNs need to be used [START_REF] Hochreiter | Untersuchungen zu dynamischen neuronalen Netzen[END_REF]; a very popular kind is the so called long short-term memory (LSTM) network proposed by [START_REF] Hochreiter | Long short-term memory[END_REF].

Implementation

We maintain the implementation in a source code repository at https://github.com/ bot-zen/. The version tagged as 0.9 comprises the version that was used to generate the results submitted to the ST. The version tagged as 1.0 is identical at its core but comes with ex-plicit documentation on how to download and install external software, and how to download and pre-process required corpora.

Our system feeds WEs and character-level subword information into a single-layer RNN with a LSTM architecture.

Word Embeddings

We incorporates w2v's original C implementation for learning WEs1 in an independent preprocessing step, i.e. we pre-compute the WEs. Then, we use gensim2 , a Python tool for unsupervised semantic modelling from plain text, to load the data, and to extract the vector representations of the embedded words as input to our NN.

Character-Level Sub-Word Information

Our implementation uses a one-hot encoding with a few additional features for representing subword information. The one-hot encoding transforms a categorical feature into a vector where the categories are represented by equally many dimensions with binary values. We convert a letter to lower-case and use the sets of ASCII characters, digits, and punctuation marks as categories for the encoding. Then, we add dimensions to represent more binary features like 'uppercase' (was uppercase prior to conversion), 'digit' (is digit), 'punctuation' (is punctuation mark), whitespace (is white space, except the new line character; note that this category is usually empty, because we expect our tokens to not include white space characters), and unknown (other characters, e.g. diacritics). This results in vectors with more than a single one-hot dimension.

Recurrent Neural Network Layer

Our implementation uses Keras, a minimalist, highly modular NNs library, written in Python and capable of running on top of either TensorFlow or Theano [START_REF] Chollet | Keras: Deep Learning library for Theano and TensorFlow[END_REF]. In our case it runs on top of Theano, a Python library that allows to define, optimize, and evaluate mathematical expressions involving multi-dimensional arrays efficiently (The Theano Development Team et al., 2016).

The input to our network are sequences of the same length as the sentences we process. During training we group sentences of the same length into batches. Each single word in the sequence is represented by its sub-word information and two WEs that come from two sources (see Section 4). Unknown words, i.e. words without a WE, are mapped to a randomly generated vector representation once, and this representation is reused later. In Total, each word is represented by 1, 800 features: two times 500 (WEs), and ten times 80 for two 5-grams (word beginning and ending). (If words are shorter than 5 characters their 5-grams are zero-padded.)

This sequential input is fed into a LSTM layer that, in turn, projects to a fully connected output layer with softmax activation function. We use categorical cross-entropy as loss function and backpropagation in conjunction with the RM-Sprop optimization for learning. At the time of writing, this was the Keras default-or the explicitly documented option to be used-for our type of architecture.

Case Study

We used our implementation to participate in the EmpiriST 2015 shared task. First, we describe the corpora used for training, and then the specific system configuration(s) for the ST. [START_REF] Brants | [END_REF] containing German newspaper texts. The corpus was semi-automatically PoS tagged, and is one of the standard corpora used for German PoS tagging. It contains 888,238 tokens in 50,472 sentences. For research and evaluation purposes, the TIGERCorpus can be downloaded for free.

German Wikipedia (w2v)

de.wiki'15 4 are user talk pages (messages from users to users, often questions and advice), article talk pages (questions, concerns or comments related to improving a Wikipedia article), and article pages of the German wikipedia from 2015, made available by the Institut für Deutsche Sprache 5 . The corpus contains 2 billion tokens (talk:379m, article talk:447m, article:1,1bn) in 79 million sentences (talk:15m, article talk:17m, article:47m), is well-sized for w2v, and also (partly) resembles or target data. It is available under the CC BY-SA 3.0 6 license.

EmpiriST 2015 shared task

For the ST we used one overall configuration for the system, but we used three different corpus configurations for training. Consequently, we participated in the ST with three runs: we used PoS tags from empirist (run 1), from Tiger v2.2 (run 2), and from both (run 3). For w2v we trained a 500-dimensional skip-gram model on empirist that ignored all words with less than 3 occurrences within a window size of 10; it was trained with negative sampling (value 5) and erroneously 9 also with hierarchical softmax. We also trained a 500dimensional continuous bag-of-words model on de.wiki'15 that ignored all words with less than 25 occurrences within a window size of 10; it was trained with negative sampling (value 3) and erroneously also with hierarchical softmax. 6 Creative Commons Attribution-ShareAlike 3.0 Unported, i.e. the data can be copied and redistributed, and adapted for any purpose, even commercially. See http:// creativecommons.org/licenses/by-sa/3.0/ for more details.

7 https://sites.google.com/site/ empirist2015/home/shared-task-data 8 For evaluation during the development phase we used empirist-trial. Unfortunately, we found out only later that the CMC part of the trial data is also part of the training data, i.e. for the CMC data our evaluation data was identical with the training data.

9 According to w2v's author, technically negative sampling and hierarchical softmax can be combined but one should avoid combining them (see https: //groups.google.com/forum/#!topic/ word2vec-toolkit/WUWad9fL0jU).

We had forgotten to deactivate an option in a data processing script.

The rational behind training the two models differently was that according to w2v author's experience10 a skip-gram model "works well with small amount[s] of the training data, [and] represents well even rare words or phrases", and a cbow model is "several times faster to train than the skip-gram, [and has] slightly better accuracy for the frequent words". The other w2v parameters were left at their default settings11 .

To optimize the system's output we ran a simple grid search for three parameters: the hidden LSTM layer's size, the dropout value for the projections from the LSTM to the output layer during training, and the number of epochs during training.

Conclusion & Outlook

We presented our submission to the EmpiriST 2015 shared task, where we participated in the PoS tagging sub-task with fair results on the CMC data and adequate results on the Web data. Still, our implementation, albeit following state-of-the art designs and methods, is quite unpolished, and can certainly gain performance with more detailed tuning. For example, adding special sequence start and sequence stop symbols to the input is typically done as a pre-processing step, which might improve the results at the beginning and the end of sentences; or we might gain some performance by adding additional hidden layers to enable the network to learn more intermediate abstractions.

A more profound design change could also help, e.g. Recurrent Memory Network are a novel recurrent architecture that have been shown to outperform LSTMs on some language modelling tasks. Finally, for learning the word embeddings we could use different corpora, or selectively extract parts from large web-corpora resemblingas much as possible-the type of data that is to be tagged.

4. 1

 1 Training Data for w2v and PoS Tagging 4.1.1 Tiger v2.2 (PoS) Tiger v2.2 3 is version 2.2 of the TIGERCorpus

  4.1.3 EmpiriST 2015 Data (PoS and w2v) empirist 7 is the CMC and Web data made available by the organizers of the ST. It contains data samples from different CMC genres and samples from text genres on the Web. The training corpus contains 10,053 tokens and was PoS tagged by two annotators (unclear cases were decided by a third person). The trial corpus contains around 3,600 tokens (2,100 CMC 8 , 1,500 Web) and was PoS tagged by one annotator (without systematic error checks). See Beißwenger et al. (to appear) for more details.

Table 1 :

 1 Official results of our PoS tagger for the three runs on the EmpiriST 2015 shared task data.

	The found values were size:1024, dropout:0.1,
	epochs:20.
		CMC Web
	(1) empirist	81.03 86.97
	(2) Tiger v2.2	73.56 89.73
	(3) empirist+Tiger v2.2 85.42 90.63
	Winning Team	87.33 93.55

https://code.google.com/archive/p/ word2vec/

https://radimrehurek.com/gensim/

http://www.ims.uni-stuttgart.de/ forschung/ressourcen/korpora/tiger.html

http://www1.ids-mannheim.de/kl/ projekte/korpora/verfuegbarkeit.html# Download

http://www.ids-mannheim.de

https://groups.google.com/d/ msg/word2vec-toolkit/NLvYXU99cAM/ E5ld8LcDxlAJ

-sample 1e-3 -iter 5 -alpha 0.025 for skip-gram and -alpha 0.05 for continuous bag-of-words