Patrick Amestoy
email: patrick.amestoy@enseeiht.fr

Cleve Ashcraft

Olivier Boiteau
email: olivier.boiteau@edf.fr

Alfredo Buttari
email: alfredo.buttari@enseeiht.fr

Excellent Jean-Yves L'

Clément Weisbecker
email: clement.weisbecker@enseeiht.fr

Improving multifrontal methods by means of block low-rank representations *

Keywords: sparse direct methods, multifrontal method, low-rank approximations, elliptic PDEs AMS subject classifications. 05C50, 65F05, 65F50

Matrices coming from elliptic Partial Differential Equations (PDEs) have been shown to have a low-rank property: well defined off-diagonal blocks of their Schur complements can be approximated by low-rank products. Given a suitable ordering of the matrix which gives to the blocks a geometrical meaning, such approximations can be computed using an SVD or a rank-revealing QR factorization. The resulting representation offers a substantial reduction of the memory requirement and gives efficient ways to perform many of the basic dense algebra operations. Several strategies have been proposed to exploit this property. We propose a low-rank format called Block Low-Rank (BLR), and explain how it can be used to reduce the memory footprint and the complexity of direct solvers for sparse matrices based on the multifrontal method. We present experimental results that show how the BLR format delivers gains that are comparable to those obtained with hierarchical formats such as Hierarchical matrices (H matrices) and Hierarchically Semi-Separable (HSS matrices) but provides much greater flexibility and ease of use which are essential in the context of a general purpose, algebraic solver.

Utilisation de représentations de rang faible par blocs dans les méthodes multifrontales

Résumé : Il a été démontré que les matrices provenant d'équations aux dérivées partielles elliptiques ont une propriété de rang faible: certains blocs hors diagonaux de leurs compléments de Schur peuvent être approchés par des produits de rang faible. Etant donnée une permutation de la matrice qui donne un sens géométrique aux blocs, de telles approximations peuvent être calculées avec une décomposition en valeurs singulières ou une factorisation QR. La représentation correspondante offre une réduction significative en termes de besoins mémoire tout en fournissant des moyens plus efficaces d'effectuer les opérations d'algèbre linéaire dense. Plusieurs stratégies ont été proposées dans la littérature pour exploiter cette propriété. Nous proposons un format Block Low Rank (BLR) et montrons comment il peut être utilisé pour réduire les besoins mémoire et la complexité des solveurs directs pour matrices creuses basées sur la méthode multifrontale. Des résultats expérimentaux montrent que le format BLR donne des gains comparables à ceux obtenus avec des formats hiérarchiques, tout en fournissant une plus grande flexibilité et facilité d'utilisation, qui sont essentielles dans un solveur algébrique général.

Mots-clés : matrices creuses, méthodes directes, méthode multifrontale, approximation de rang faible, EDP elliptique 1. Introduction. We are interested in efficiently computing the solution of large sparse linear systems which arise from various applications such as mechanics and fluid dynamics. Modern applications commonly require the solution of linear systems stemming from the discretization of partial differential equations with several millions of unknowns. A sparse linear system is usually referred to as:

Ax = b , (1.1)
where A is a sparse matrix of order n, x is the unknown vector of size n and b is the right-hand side vector of size n.

Two types of methods [START_REF] Golub | Matrix computations[END_REF] are commonly used to solve (1.1). Iterative methods are cheap in memory consumption and provide better scalability in parallel environments but their effectiveness strongly depends on the numerical properties of the problem. Direct methods are more robust and reliable but costly (in terms of flops and memory). As the problems become larger and larger, time and memory complexity become critical in order to be able to tackle challenging applicative problems.

This work focuses on a well known direct approach called the multifrontal method [START_REF] Amestoy | The multifrontal method[END_REF]. The objective is to compute a Cholesky LL T factorization of A:

A = LL T , (1.2)
where L is a lower triangular matrix. Without loss of generality, here and in the rest of the paper we only discuss symmetric, positive definite systems for the sake of simplicity. In practice and in the experimental section, the presented methods and algorithms have been applied to the cases of symmetric, indefinite (LDL T factorization) and unsymmetric (LU factorization) matrices. We aim to improve the multifrontal method by means of low-rank approximation techniques, which give the opportunity to reduce the memory requirements and the complexity of dense factorizations.

After an introduction to the multifrontal method in the context of nested dissection, we show how it is possible to exploit low-rank properties of matrices to decrease their storage requirements as well as the complexity of some of the basic algebraic operations they are involved in. Then, we present the relevant details of some low-rank formats proposed in the literature and propose a new one called Block Low-Rank (BLR). Through a brief experimental comparison, we show that all these formats have comparable efficiency on a set of test problems. We then show how the BLR format can be integrated within a general purpose linear multifrontal solver. Lastly, we present and discuss experimental results in order to show the efficiency of our implementation of a low-rank sequential multifrontal solver. The study is carried out on a set of large-scale artificial and real life applicative problems.

2. The multifrontal method. The multifrontal method was first introduced by Duff and Reid [START_REF] Duff | The multifrontal solution of indefinite sparse symmetric linear systems[END_REF][START_REF]The multifrontal solution of indefinite sparse symmetric linear systems[END_REF] in 1983 and, since then, has been the object of numerous studies and the method of choice for several, high-performance, software packages such as MUMPS [START_REF] Amestoy | MUMPS: a general purpose distributed memory sparse solver[END_REF], UMFPACK [START_REF] Davis | Algorithm 832: UMFPACK V4.3 -an unsymmetric-pattern multifrontal method[END_REF], WSMP [START_REF] Gupta | WSMP: A high-performance shared-and distributed-memory parallel sparse linear equation solver[END_REF], HSL [START_REF] Gould | A numerical evaluation of HSL packages for the direct solution of large sparse, symmetric linear systems of equations[END_REF], PSPASES [START_REF] Gupta | PSPASES: An efficient and scalable parallel sparse direct solver[END_REF] and DSCPACK [START_REF] Raghavan | DSCPACK: Domain-separator codes for the parallel solution of sparse linear systems[END_REF].

This section provides a brief presentation of the multifrontal algorithm based on the nested dissection method which better suits the context of our work and allows for an easier understanding of how low-rank approximation techniques can be used within sparse, multifrontal solvers; more general and formal presentations are available in the literature [START_REF] Amestoy | The multifrontal method[END_REF][START_REF] Duff | Direct Methods for Sparse Matrices[END_REF]. We assume that the reader is familiar with the basics of graph theory.

It is well known that, despite their good reliability and numerical robustness, sparse direct solvers pose heavy requirements in terms of computational and memory resources. This is mostly due to the fact that the factors of a sparse matrix are, in general, much denser and this, in turn, is due to the appearance of fill-in, i.e., new nonzero coefficients introduced by the factorization process. Many methods and algorithms have been proposed in order to reduce the amount of fill-in and, among them, one of the most commonly used is the nested dissection method [START_REF] George | Nested dissection of a regular finite element mesh[END_REF].

Given a sparse, symmetric matrix A of size n, its structure can be represented with an adjacency graph, i.e., a graph G(V, E) containing n vertices (one for each unknown in A) and edges (i, j) for all a ij = 0. The fill-in can then be easily modeled using this graph. Specifically, eliminating a variable of A amounts to removing the associated vertex from G along with all incident edges and adding new edges that connect the neighbors of the eliminated vertex that were not already connected to each other; these newly introduced edges represent fill-in coefficients. Now assume that a vertex separator S of G is computed, i.e., a subset of vertices which, if removed, splits the graph into two subgraphs D 1 and D 2 and assume that A is permuted in such a way that all the variables in D 1 are eliminated first, all those in D 2 are eliminated second and all those in S are eliminated last. Because all the neighbors of vertices in D 1 are either in D 1 or in S, no fill-in will be generated inside the submatrix that connects the variables in D 1 to those in D 2 . This is illustrated in Figure 2.1 for the simple case of a matrix from a 5-point stencil operator. Because fill-in occurring inside the diagonal blocks A 1,1 and A 2,2 may still be excessive, this procedure can be recursively and independently applied to the two subgraphs corresponding to D 1 and D 2 until the fill-in inside diagonal blocks can be considered negligible or not worth being reduced. This procedure essentially describes the nested dissection method although modern nested dissection ordering tools [START_REF] Karypis | MeTiS: A Software Package for Partitioning Unstructured Graphs, Partitioning Meshes, and Computing Fill-Reducing Orderings of Sparse Matrices Version 4[END_REF][START_REF] Pellegrini | Scotch and libscotch 5.0 User's guide[END_REF] employ much more sophisticated algorithms.

• • • • • • • • • • • • • • • • • • • • D 1 D 2 S (a)
Applying the nested dissection algorithm to a graph generates a separators tree which implicitly defines a fill-in reducing permutation of the input matrix. This is shown in Figure 2.2, where, for the leaves of the tree, we assumed S i = D i which is equivalent to saying that S i is a separator that splits D i into two empty subdomains.

The separators tree can also be regarded as an elimination or assembly tree -a data structure introduced by Schreiber [START_REF] Schreiber | A new implementation of sparse Gaussian elimination[END_REF] which is at the base of all modern sparse factorization algorithms. This tree defines the dependencies between the variables of a matrix and thus, implicitly, the order in which they have to be eliminated. Specifically, the elimination tree states that the elimination of the variables associated with a node (or with the corresponding separator, in our case) only affects variables associated with ancestor nodes, that is, nodes along the path that connects the eliminated node to the root of the tree. Therefore, all the pivotal orders defined by topological traversals of the tree are equivalent in the sense that they produce the same amount of fill-in.

The multifrontal method [START_REF] Amestoy | The multifrontal method[END_REF][START_REF] Duff | The multifrontal solution of indefinite sparse symmetric linear systems[END_REF][START_REF]The multifrontal solution of indefinite sparse symmetric linear systems[END_REF] achieves the factorization of a sparse matrix through a sequence of operations on relatively small dense matrices called frontal matrices or, simply, fronts. One of these matrices is associated with each node of the tree and is formed by coefficients related to the variables associated with that node as well as their neighbors. Note that the neighbors include those variables that have become such due to fill-in introduced by the elimination of previous variables.

The multifrontal factorization consists in a topological order (i.e., bottom-up) traversal of the tree where, each time a node is visited, two operations are performed:

• assembly: the frontal matrix is formed by summing the coefficients in the rows and columns of the variables associated with the tree node with coefficients produced by the factorization of its children. • factorization: once the frontal matrix is formed, a partial Cholesky factorization is performed on it in order to eliminate the variables associated with the tree node. The result of this operation is a set of rows of the global L factor and a Schur complement, also commonly referred to as contribution block (CB), containing coefficients that will be assembled into the parent node. As shown in Figure 2.3, the set of variables in a frontal matrix can be split into two subsets:

• fully-summed (FS) variables: these are the variables associated with the tree node or, equivalently, with the corresponding separator. They are called fully-summed because the rows and columns in the front are up to date with respect to previously eliminated variables, that is, those that have been eliminated at descendant nodes. • non fully-summed (NFS) variables: this subset contains the neighbors of the fully-summed variables and is formed by pieces of separators belonging to ancestor nodes which form a border around the fully-summed variables. For the frontal matrix at node S 7 of the tree in Figure 2.2, this border is formed by half of the S 1 separator and the entire S 3 separator; variables in this border have become neighbors of those in S 7 due to fill-in coefficients introduced by the eliminations of variables in S 14 and S 15 . Once a frontal matrix is factorized, the resulting factors are stored apart while the contribution block is held in a temporary memory area which is then freed when the parent front is processed. For this reason the global amount of memory needed to run the multifrontal factorization is, in general, bigger than the size of the resulting factors. The temporary memory area is managed as a stack which is very convenient if the tree is traversed in postorder. In the rest of this paper we will refer to this area simply as CB stack.

This paper describes a technique based on the usage of low-rank approximations to represent frontal matrices and to perform all the operations they are involved in. This allows to reduce both the memory consumption as well as the total volume of computations of the multifrontal method.

Low-rank approximations.

In this section, we present how low-rank properties of matrices can be revealed and exploited. Several matrix representations are described and discussed.

Low-rank property.

A low-rank matrix can be represented in a form which decreases its memory requirements and the complexity of basic linear algebra operations it is involved in, such as matrix-matrix products and triangular solves. This is formalized by Definition 3.1.

Definition 3.1 ([7]

Low-rank matrix). Let A be a matrix of size m × n. Let k ε be the approximated numerical rank of A at accuracy ε. A is said to be a low-rank matrix if there exist three matrices U of size m × k ε , V of size n × k ε and E of size m × n such that :

A = U • V T + E ,
where E 2 ≤ ε and k ε (m + n) < mn.

The numerical rank at precision ε, k ε , can be computed together with U and V with a Singular Value Decomposition (SVD), in which case U and V are unitary, or, less precisely but much faster, with a Rank-Revealing QR (RRQR) factorization, in which case only U is unitary. Low-rank approximation techniques are based upon the idea to ignore E and simply represent A as the product of U and V T . Thus, A = U • V T is said to be the low-rank form of A and it is an approximation whose precision can be controlled through the parameter ε called low-rank threshold. In the rest of this paper, the subscript ε in k ε will be omitted for the sake of clarity. Definition 3.1 simply states that k has to be small enough to guarantee that the low-rank form of A takes less memory that the standard form. Analogous criteria can be defined in the case where the objective is to reduce the complexity of operations involving A. Table 3.1 shows the requirements on the rank k in order to reduce the operation count of the matrix-matrix product and of the triangular system solve, when the resulting matrix is stored in dense form.

= U 1 V T 1 and A 2 = U 2 V T 2 are n × n, L is lower triangular, U 1 , V 1 , U 2 and V 2 are n × k.
The bracketing is critical to exploit the smaller dimension of each matrix.

operation type dense low-rank rank requirement

Cholesky factorization

LL T n 3 /3 - - RRQR compression U 1 V T 1 6kn 2 -6k 2 n + 10 3 k 3 - - triangular solve U 1 (V T 1 L -T ) n 3 3kn 2 k < n/3 matrix-matrix product U 1 (V T 1 U 2 )V T 2 2n 3 2kn 2 + 4k 2 n k < n/2
In practice, matrices coming from application problems are not low-rank, which means that they cannot be directly written in low-rank form. However, Börm [START_REF] Börm | Efficient Numerical Methods for Non-local Operators[END_REF] and Bebendorf [START_REF] Bebendorf | Hierarchical Matrices: A Means to Efficiently Solve Elliptic Boundary Value Problems[END_REF] show that low-rank approximations can be performed on submatrices defined by an appropriately chosen partitioning of matrix indices. Assuming I = {1, . . . , n} is the set of row (and column) indices of A, a set of indices σ ⊂ I is called a cluster. Then, a clustering of I is a disjoint union of clusters which equals I. b = σ × τ ⊂ I × I is called a block cluster based on clusters σ and τ . A block clustering of I × I is then defined as a disjoint union of block clusters which equals I × I. Let b = σ × τ be a block cluster, A b = A στ is the subblock of A with row indices σ and column indices τ .

A considerable reduction of both the memory footprint and of the operations complexity can be achieved if the block clustering defines blocks whose singular values have a rapid decay (for instance, exponential) in which case each block can be accurately represented by a lowrank product U • V T . Clearly, this condition cannot be directly used in practice to define the matrix block clustering. In many practical cases it is, however, possible to exploit the knowledge of the mathematical problem or the geometrical properties of the domain where the problem is defined in order to define an admissibility condition, i.e., a heuristic rule that can be cheaply checked to establish whether or not a block is (likely to be) low-rank or that can be used to guide the block clustering computation. For instance, in the case of matrices deriving from discretized elliptic PDEs one such admissibility condition is presented in Definition 3.

diam(σ) + diam(τ ) ≤ 2η dist(σ, τ ) ,
where diam() and dist() are classical geometric diameter and distance, respectively, and η is a problem dependent parameter.

This admissibility condition follows the intuition that variable sets that are far away in the domain are likely to have weak interactions which translates into the fact that the corresponding block has a low rank; this idea is depicted in Figure 3.1(a). Figure 3.1(b) shows that the rank of a block A στ is a decreasing function of the geometric distance between clusters σ and τ . This experiment has been done on a top-level separator of a 3D 128 3 wave propagation problem called Geoazur128 (and further described in Table 6.1), with square clusters of dimension 16 × 16, so that each subblock has size 256. It shows that depending on the distance between clusters, there is potential for compression which can be exploited. The dashed line at y = 128 shows the cutoff point where it pays to store the subblock using a low-rank representation. An admissibility condition being given, it is possible to define an admissible block clustering, i.e., a block clustering with admissible blocks, which can be used to efficiently represent a dense matrix by means of low-rank approximations with an accuracy defined by the low-rank threshold ε in Definition 3.1.

It has to be noted that, in practice, admissibility conditions may be even simpler than the one presented above and may also take into account aspects that, for example, pertain to the efficiency of basic linear algebra kernels on the U and V matrices of the low-rank form; further details about admissibility conditions are given in Section 4.

Matrix representations.

In order to exploit low-rank blocks arising from admissible block clusterings, several matrix representations have been proposed in the literature. In this section, we very briefly describe Hierarchical matrices [START_REF] Bebendorf | Hierarchical Matrices: A Means to Efficiently Solve Elliptic Boundary Value Problems[END_REF] (H-matrices) and Hierarchically Semi-Separable matrices [START_REF]Fast algorithms for hierarchically semiseparable matrices[END_REF] (HSS matrices) because they have been used in the context of a multifrontal solver. The Hierarchically Block-Separable (HBS) format introduced by Gillman et al. [START_REF] Gillman | A direct solver with O(N) complexity for integral equations on one-dimensional domains[END_REF] is highly relevant too but will not be discussed in details as it is essentially equivalent to HSS. All these structures can be used to approximate fronts and thus globally reduce the memory consumption as well as the flop count of the multifrontal process.

Hierarchical matrices (H-matrices).

The H-matrix [START_REF] Bebendorf | Hierarchical Matrices: A Means to Efficiently Solve Elliptic Boundary Value Problems[END_REF][START_REF] Börm | Efficient Numerical Methods for Non-local Operators[END_REF] format, where H stands for Hierarchical, is historically the first low-rank format for dense matrices. This format is based on an admissible block clustering which is obtained by a recursive subdivision of I × I. Algorithm 1 shows how the H-matrix format of a matrix A is built.

The result of this recursive procedure is a block cluster tree built in a top-down fashion. Note that a set of siblings of the tree defines a block clustering of the block cluster associated with their parent node and that the final admissible clustering is defined by the leaves of the tree.

Figure 3.2 shows an example of a H-matrix with the corresponding block cluster tree. In this case the block clustering has been defined by splitting the set of matrix indices as

I = I 7 = I 3 ∪ I 6 = {I 1 ∪ I 2 } ∪ {I 4 ∪ I 5 }.
In a more general case, the structure of a H-matrix is not necessarily as regular as in the provided example which means that the diagonal blocks may not have the same size. For more details about H-matrices, please refer to Bebendorf [START_REF] Bebendorf | Hierarchical Matrices: A Means to Efficiently Solve Elliptic Boundary Value Problems[END_REF], Börm [START_REF] Börm | Efficient Numerical Methods for Non-local Operators[END_REF] or Hackbusch [START_REF] Hackbusch | A sparse matrix arithmetic based on H-matrices. Part I: Introduction to H-matrices[END_REF].

Algorithm 1 H-matrix construction

Input A, a matrix defined on row and column indices I. Output A, the H-matrix form of A.

1: initialize list with [I × I] 2: while list is not empty do if b is admissible then 5:

A b ← low-rank form of A b 6:
else if b is large enough to split then end if 12: end while

A =      D 11 U 12 V T 12 U 36 V T 36 U 21 V T 21 D 22 U 63 V T 63 D 44 U 45 V T 45 U 54 V T 54 D 55      I 7 × I 7 I 3 × I 3 I 3 × I 6 I 6 × I 3 I 6 × I 6 I 1 × I 1 I 1 × I 2 I 2 × I 1 I 2 × I 2 I 4 × I 4 I 4 × I 5 I 5 × I 4 I 5 × I 5 Figure 3.
2. H-matrix structure and associated block cluster tree 3.2.2. Hierarchically Semi-Separable matrices (HSS matrices). HSS matrices are also based on a hierarchical blocking defined by a cluster tree which is built by recursively splitting the set of matrix indices.

Definition 3.3 ([33] HSS matrix).

Assume A is an n × n dense matrix, and I = {1, 2, . . . , n}. I is recursively split into smaller pieces following a binary tree T with k nodes, denoted by j = 1, 2, . . . , k ≡ root( T ). Let t j ⊂ I be a cluster associated with each node j of T . (T is used to manage the recursive partition of A.) We say A is an HSS form with the corresponding postordered HSS tree T if :

1. T is a full binary tree in its postordering, or, each node j is either a leaf or a non-leaf node with two children j 1 and j 2 which satisfy j 1 < j 2 < j ; 2. The index sets satisfy t j1 ∪ t j2 = t j and t j1 ∩ t j2 = ∅ for each non-leaf node j, with t k ≡ I ; 3. For each node j, there exist matrices

D i , U i , V i , R i , W i , B i (called HSS generators),
which satisfy the following recursions for each non-leaf node j:

D j ≡ A| tj ×tj = D j1 U j1 B j1 V T j2 U j2 B j2 V T j1 D j2 , U j = U j1 R j1 U j2 R j2 , V j = V j1 W j1 V j2 W j2 , where U k , V k , R k , W k and B k are not needed (since D k ≡ A
is the entire diagonal block without a corresponding off-diagonal block). Due to the recursion, only the D j , U j , V j generators associated with a leaf node j of T are stored. An example of a HSS matrix is given in Figure 3.3. The construction of an HSS matrix is achieved through a topological order traversal of the T tree, also referred to as HSS tree. Each time a node j is visited, the block-row and block-column corresponding to the related cluster (i.e., A tj ,: = A tj ,(I\tj ) and A :,tj = A (I\tj ),tj , respectively) are compressed into a low-rank form. Note that, apart from the leaves, the compressed block-row or block-column is formed by combining the result of previous

A =      D 1 U 1 B 1 V T 2 U 1 R 1 B 3 W T 4 V T 4 U 1 R 1 B 3 W T 5 V T 5 U 2 B 2 V T 1 D 2 U 2 R 2 B 3 W T 4 V T 4 U 2 R 2 B 3 W T 5 V T 5 U 4 R 4 B 6 W T 1 V T 1 U 4 R 4 B 6 W T 2 V T 2 D 4 U 4 B 4 V T 5 U 5 R 5 B 6 W T 1 V T 1 U 5 R 5 B 6 W T 2 V T 2 U 5 B 5 V T 4 D 5      D 1 D 2 D 4 D 5 D 3 D 6 Figure 3.3.
Structure of a HSS matrix and its basic construction scheme compressions. The bottom of Figure 3.3 shows this process when T has four leaves (seven nodes, total). Although it should be possible to have a different blocking for the upper and lower triangular parts of the matrix as for the H-matrices, we have not found any example of this case in the literature.

For a more detailed presentation and study of HSS matrices, please refer to Xia et al. [START_REF] Xia | Efficient structured multifrontal factorization for general large sparse matrices[END_REF][START_REF]Fast algorithms for hierarchically semiseparable matrices[END_REF][START_REF] Xia | Superfast multifrontal method for large structured linear systems of equations[END_REF] and Wang et al. [START_REF] Wang | On 3D modeling of seismic wave propagation via a structured parallel multifrontal direct Helmholtz solver[END_REF][START_REF] Wang | Efficient scalable algorithms for hierarchically semiseparable matrices[END_REF].

Block Low-Rank matrices (BLR matrices).

We propose a new structure based on a non-hierarchical blocking of the matrix.

Definition 3.4 (Block Low-Rank matrix (BLR)). Given an admissibility condition, let P be an admissible block clustering with p 2 clusters. Let A be an n × n matrix.

A =      B 11 B 12 • • • B 1p B 21 B 22 • • • B 2p . . . . . . . . . . . . B p1 B p2 • • • B pp      is called a "Block Low-Rank matrix" if ∃(σ, τ ) ∈ {1, 2, . . . , p} 2 , there exists k στ such that B στ is a low-rank block with rank k στ .
The structure of a BLR matrix is not hierarchical: a flat block matrix structure is used, as illustrated in Figure 3

A =      D 1 U 12 V T 12 U 13 V T 13 U 14 V T 14 U 21 V T 21 D 2 U 23 V T 23 U 24 V T 24 U 31 V T 31 U 32 V T 32 D 3 U 34 V T 34 U 41 V T 41 U 42 V T 42 U 43 V T 43 D 4      . (3.1)
Figure 3.4 shows the global structure of the BLR representation of a dense Schur complement of order 128×128 corresponding to the top level separator of a 128×128×128 Laplacian problem, with a low-rank threshold set to 10 -14 . The numbering scheme illustrated in 3.4(a) is recursive although this is only required for H and HSS matrices.

BLR matrices can be viewed as a particular case of H-matrices where all the subblocks have been subdivided identically, i.e., where all the branches of the block cluster tree have the same depth. (b) Structure of a BLR matrix. The darkness of a block is proportional to its storage requirement (the lighter a block is, the smaller is the memory needed to store it). Each block in the matrix is of size 512 × 512.

Comparative study.

A preliminary comparative study of the three formats described above (both in terms of memory and computational costs) is presented in this section in order to validate the potential of our BLR format in the context of the development of algebraic methods for exploiting low-rank approximations within a general purpose, multifrontal solver. The H, HSS and BLR formats are used for compressing the Schur complement associated with a separator of a cubic mesh (or graph) of size 128; the separator is a 128 × 128 surface lying in the middle of the cubic domain and the corresponding Schur complement is a dense matrix of order 16384. This is assessed for two different problems: the Geoazur problem (see Section 6.1) and a Laplacian operator discretized with a 11-point stencil. For all three formats, the clustering was defined by a 8 × 8 recursive checkerboard partitioning of the separator into blocks of size 256, using the same approach as in the 4 × 8 case from Figure 3.4(a); this clustering is essentially what is referred to as the weak admissibility condition in Börm [START_REF] Börm | Efficient Numerical Methods for Non-local Operators[END_REF]Remark 3.17]. Other block sizes have been experimented and give comparable results for both problems. These results are in compliance with what was observed in Wang et al. [START_REF] Wang | Massively parallel structured multifrontal solver for time-harmonic elastic waves in 3D anisotropic media[END_REF].

In terms of memory compression, Figures 3.5(a) and 3.6(a) show that the three formats appear to be roughly equivalent at low precision: the truncation is so agressive that almost no information is stored anymore. In the case of the Geoazur problem, the three formats provide comparable gains with BLR being slightly worse than the other two at high accuracy and better when the approximation threshold is bigger than 10 -10 . For the Laplacian problem, BLR is consistently better than the hierarchical formats although all three provide, in general, considerable gains. These preliminary results suggest that the BLR format delivers gains that are comparable to those obtained with the hierarchical ones on the test problems. show the cost of computing the H, HSS and BLR formats starting from a standard dense matrix. This is computed as the cost of the partial RRQR factorizations; because for H and HSS matrices these operations are performed on much larger blocks, the conversion to BLR format is much cheaper with respect to the cases of hierarchical formats and also with respect to the cost of an LU factorization of the same front. The cost of a full rank factorization is indeed 10 (Geoazur128 problem) or 100 (Laplacian problem) times larger than the BLR compression cost at full accuracy (ε = 10 -14 ). Moreover, because a truncated rank-revealing QR factorization is used to compress the blocks, the compression cost decreases when the accuracy decreases. This property is extremely useful as it allows to switch from full-rank to low-rank and vice versa at an acceptable cost compared to the dense LU factorization cost.

A deep and detailed comparison of low-rank formats, both theoretical and experimental, is out of the scope of this paper; nonetheless, it is the subject of ongoing research work and collaborations. 3.3. Motivation to use BLR format in a multifrontal solver. H and HSS matrices have been successfully used to accelerate multifrontal solvers and to reduce their memory consumption [START_REF] Xia | Efficient structured multifrontal factorization for general large sparse matrices[END_REF]. These approaches, however, either rely on the knowledge of mathematical properties of the problem or on the geometric properties of the domain on which it is defined (i.e., the discretization mesh) or lack some of the features that are essential for a general purpose, sparse, direct solver like, for instance, robust threshold pivoting or dynamic load balancing in a parallel environment. Our objective is, instead, to exploit low-rank approximations within a general purpose, multifrontal solver which aims at providing robust numerical features in a parallel framework and where no knowledge of the problem can be assumed except the matrix itself. In such a context, hierarchical structures are likely to be too hard to handle and may severely limit the necessary flexibility:

• In a parallel multifrontal solver, frontal matrices may be statically or dynamically partitioned, in an irregular way in order to achieve a good load and memory balance; hierarchical formats may pose heavy constraints or may be complex to handle in this case.

• The HSS format achieves good compression rates only if each node of the HSS tree defines an admissible cluster. This means that the diagonal blocks cannot be permuted in an arbitrary way but have to appear along the diagonal in a specific order which depends on the geometry of the separator. This information may not be available or too complex to extrapolate in a general purpose, algebraic solver. • In HSS matrices the compressions are achieved through SVD or RRQR decompositions on block-rows and block-columns, i.e., strongly over or under-determined submatrices. This has a twofold disadvantage. First it is based on the assumption that all the columns in a block-row (or, equivalently, rows in a block-column) lie in a small space which has been proven only for some classes of problems and may not be true in general [START_REF] Chandrasekaran | On the numerical rank of the off-diagonal blocks of Schur complements of discretized elliptic PDEs[END_REF]. Second, due to the shape of the data, these operations are often inefficient and difficult to parallelize; this has led researchers to consider the use of communication-avoiding RRQR factorizations or the use of randomization techniques. • It is not natural, and does not seem to have been done in practical implementations, to use the H or HSS formats on the L 21 part of a frontal matrix (see Figure 2.

3(b)).

The approach presented by Xia [START_REF] Xia | Efficient structured multifrontal factorization for general large sparse matrices[END_REF], for instance, handles this submatrix as a whole low-rank block which may result in sub-optimal compression and which again requires inefficient operations to compute the low-rank form. • The assembly of a frontal matrix is very difficult to achieve if the contribution blocks are stored in a hierarchical, low-rank format. For this reason, intermediate full-rank representations are used in practice [START_REF] Xia | Efficient structured multifrontal factorization for general large sparse matrices[END_REF][START_REF] Wang | On 3D modeling of seismic wave propagation via a structured parallel multifrontal direct Helmholtz solver[END_REF][START_REF] Xia | Superfast multifrontal method for large structured linear systems of equations[END_REF]. This, however, comes at a very significant cost (see Figures 3.5(b) and 3.6(b)). The analysis presented in Section 3.2.4, although restricted to a limited number of problems, suggests that a simpler format such as BLR delivers benefits comparable to the H and HSS ones; BLR, though, represents a more suitable candidate for exploiting low-rank techniques within a general-purpose, multifrontal solver as it presents several advantages over hierarchical formats:

• The matrix blocking is flat, i.e., not hierarchical, and no relative order is imposed between the blocks of the clustering. This is a very valuable feature because it allows to compute the index clustering more easily and because it delivers much greater flexibility for distributing the data in a parallel environment. • The size of blocks is homogeneous and is such that the SVD or RRQR operations done for computing the low-rank forms can be efficiently executed concurrently with sequential code (for example, the LAPACK GESVD or GEQP3 routines). This property is extremely valuable in a distributed memory environment where reducing the volume of communications is important. • The assembly of frontal matrices is relatively easy if the contribution blocks are in BLR format. We also have the option of switching to an intermediate full-rank format because compression is relatively cheap (see Figures 3.5(b) and 3.6(b)). Motivated by the previous observations, we focus in this paper on the BLR format for exploiting, in an algebraic setting, low-rank techniques. The remainder of this paper describes how block clustering is performed and the main algorithmic issues of the low-rank factorization phase.

4. Block clustering. We describe in this section how block clusterings can be computed in order to obtain efficient BLR representations of fronts.

Admissibility condition.

The admissibility condition presented in Definition 3.2 requires geometric information in order to properly compute diameters and distances. In the context of an algebraic solver, this is not conceivable because the only available information is the matrix.

Thus, another condition must be used and a natural idea is to use the graph of the matrix. In Börm [START_REF] Börm | Efficient Numerical Methods for Non-local Operators[END_REF] and Grasedyck [START_REF] Grasedyck | Parallel black box H -LU preconditioning for elliptic boundary value problems[END_REF]

diam G (σ) + diam G (τ ) ≤ 2η dist G (σ, τ ) ,
where diam G and dist G are classical graph diameter and distance in G, respectively, and η a problem dependent parameter.

This admissibility condition may still be unpractical because, first, it is not clear how to choose the η parameter in an algebraic context and, second, because computing diameters and distances of clusters of a graph may be costly. This condition, however, can be further simplified and complemented with other practical considerations in order to define a clustering strategy suited for the BLR format:

1. For efficiency reasons, the size of the clusters should be chosen between a minimum value that ensures a good efficiency of the BLAS operations performed later with U and V matrices and a maximum value that allows to perform the SVDs or RRQRs on the blocks with sequential routines as well as an easy and flexible distribution of blocks in a parallel environment. 2. The distance between any two clusters σ and τ has to be greater than zero (dist G (σ, τ ) > 0) which amounts to saying that all the clusters are disjoint. Note that this point only is equivalent to the weak admissibility condition proposed by Börm [8, Remark 3.17]. 3. For a given size of clusters (and consequently a given number of clusters in G) the diameter of each cluster should be as small as possible in order to group within a cluster only those variables that are likely to strongly interact with each other. For example, in the case where G is a flat surface, it is better off to define clusters by partitioning G in a checkerboard fashion rather than cutting it into slices. Once the size of the clusters is chosen, the objectives 2 and 3 can be easily achieved by feeding the graph G to any modern graph partitioning tool such as METIS [START_REF] Karypis | MeTiS: A Software Package for Partitioning Unstructured Graphs, Partitioning Meshes, and Computing Fill-Reducing Orderings of Sparse Matrices Version 4[END_REF] or SCOTCH [START_REF] Pellegrini | Scotch and libscotch 5.0 User's guide[END_REF], as discussed in the next two sections which show how to compute the clustering of the fully assembled and non-fully assembled variables in a frontal matrix. Note that, in practice, partitioning tools take the number of partitions as input and not their size; however the size of the resulting clusters will differ only slightly because partitioning methods commonly try to balance the weight of the partitions. Because of the simple nature of the BLR format which does not require a relative order between clusters, simpler and cheaper partitioning or clustering techniques may be employed instead of complex tools such as METIS or SCOTCH. Experimental results in Section 6 show that the cost induced by the separators clustering using METIS is however acceptable.

Separator clustering.

Remember that two sets of variables are present in a front. An admissible partition of both sets must be computed. We consider G S , the graph induced by the nodes of the separator S corresponding to front F . A graph partitioning tool such as METIS or SCOTCH with a given partition size will compute BLR admissible blocks.

However, the connectedness of G S cannot be guaranteed. Figure 4.1(a) indeed shows the worst case of disconnected separator.

In this case, any partitioner cannot do anything else than computing ten groups of size one because only singletons would be given in input. This leads to a blocking made of 1 × 1 blocks, which are not compressible. In this extreme case, the clustering computed by the partitioning tool is completely useless.

Although in reality this would rarely occur, it is quite commonly the case where a separator is formed by multiple connected components that are close to each other in the global graph G. This may lead to a sub-optimal clustering because variables that strongly interact due to their adjacency will end up in different clusters.

This problem may be overcome by reconnecting G S in a way that takes into account the geometry and shape of the separator: variables close to each other in the original graph G have to be still close to each other in the reconnected G S . We describe an approach that achieves this objective by extending the subgraph induced by each separator with a halo formed by a relatively small number of level sets. The graph partitioning tool is therefore run on this extended graph and the resulting partitioning projected back on the original subgraph G S . Figure 4.1 shows how this is done on the example above using just one level set. For a limited number of level sets, the extended graph preserves the shape of the separator, keeps the cost of computing the clustering limited and allows to compute clusterings that better comply with the strategy presented in the previous section. In practice with complex physical domains, more layers may be needed to effectively reconnect the separators. On regular grids, we observed that two layers are enough to reconnect the separator and to obtain good performance.

(a) separator sub- graph G S (b) halo subgraph G H • • • • • • • • • (c) 3 partitions of G H • • • (d) 3 partitions of G S

Explicit border clustering.

The previous section shows how the halo method can be employed to compute clusterings of the separators. This method can obviously be used for clustering any set of variables and, therefore, also those in the non fully-assembled part of the front that form a border around the corresponding separator (see Section 2). Although this approach will lead to an efficient clustering of the non fully-summed variables, it may suffer performance issues due to its high cost.

Inherited border clustering.

As one can observe on Figure 4.2(b), one particular variable belongs to one separator which, in turn, may be (partially) included in multiple borders. Inversely, as Figure 2.2 illustrates, a border can be viewed as a union of parts of separators lying on the path that connects the related node to the root of the elimination tree. As a consequence, clustering the variables of all the separators is sufficient to obtain, by induction, a clustering of the borders.

A top-down traversal of the separator tree is performed and the variables of each separator are clustered with the halo method proposed in Section 4.2. Then, the clustering of a given border is inherited from the clustering of all the separators it is made of, as shows Figure 4.2(b). Depending on where the separators intersect, small clusters, which are not admissible, may indeed be formed; as a result, the CB may include blocks which are too small to be effectively compressed and to achieve a good BLAS efficiency for the related operations. Note that this problem also affects the blocking of the L 21 submatrix but to a lower extent because the effect is damped by the good clustering computed for the separator (see Figure 4.3).

(a) L 11 blocking : optimal blocks because any block interconnects two optimal groups from the corresponding original separator clustering (b) L 21 blocking : close to optimal block because any block interconnects at least one optimal cluster from the corresponding original separator clustering (c) CB blocking : not always optimal because a block possibly interconnects two non-optimal small clusters from the inherited clustering of the border. Braces show a reclustering possibility. To avoid having small blocks, reclustering strategies are being investigated. Figure 4.3(c) shows a basic strategy which merges neighbor clusters together in order to improve the BLAS efficiency as well as the compression gains.

It has to be noted that the inherited clustering also provides another convenient property: the blocking of a frontal matrix is compatible with the blocking of its parent front. This translates into the fact that one block of its Schur complement will be assembled into exactly one block of the parent front. This may considerably ease the assembly of frontal matrices especially in the parallel case where frontal matrices are distributed.

We will report in Section 6.3 on the effects of explicit and inherited clustering strategies and show related experimental results.

5. Block Low-Rank multifrontal method. This section describes how the BLR format can be integrated within a multifrontal solver in order to reduce its complexity (and, thus, its time to completion) as well as its memory footprint for solving problems with low-rank properties.

The previous sections describe how the BLR format possesses properties and features that are extremely favorable in a parallel fully functional general multifrontal solver. However, presenting and analysing such a solver is out of the scope of this work. The main objective of this article is, instead, to evaluate the potential of the BLR format and to show that, when used within a multifrontal solver, it is a valid alternative to other approaches such as those based on hierarchical formats. BLR matrices can be efficiently employed to improve the multifrontal method. In particular, we will show how the fronts factorization (Section 5.2), where most of the computational time of the global execution is spent, can benefit from the BLR representation of the fronts. Clearly, using the BLR format also affects the way the fronts are assembled. This is described in Section 5.1 5.1. Assembly. The assembly phase of the BLR multifrontal method depends on how the BLR factorization (see Section 5.2) is done.

Because of their nature, assembly operations are relatively inefficient [START_REF] Ashcraft | The influence of relaxed supernode partitions on the multifrontal method[END_REF] and their cost can significantly affect the speed of the overall multifrontal factorization. In the experimental code we developed, we decided not to perform the assembly operations in BLR format in order to avoid an excessive amount of indirect addressing which may seriously impact performance. An analogous choice was also made in related work on the usage of HSS matrices within multifrontal solvers [START_REF] Xia | Efficient structured multifrontal factorization for general large sparse matrices[END_REF].

Frontal matrices are, thus, assembled in full-rank form and compressed progressively, by panels, as the partial front factorization proceeds, as shown in the next section. This, however, does not mean that contribution blocks are not compressed; despite the fact that compressing contribution blocks does not help reducing the overall number of floating point operations, storing them in low-rank form has a twofold advantage:

1. It allows to reduce the peak of active memory. In a sequential, multifrontal method, at any moment, the active memory is defined as the sum of the size of the CB stack(described in Section 2) plus the size of the front currently being factorized.

Stacking the contribution blocks in low-rank form reduces the active memory. Although this makes the relative weight of the current front in the active memory higher, it has to be noted that if the partial front factorization is done in a left-looking fashion, frontal matrices can be assembled panelwise which means that in the active memory only one panel at a time is stored in full-rank. We have not implemented this feature yet, and, as discussed in the next section, the front factorization is still performed in a right-looking fashion. 2. In a parallel environment, it reduces the communications volume. In a parallel solver, frontal matrices are generally mapped onto different nodes which means that assembly operations involve transferring contribution blocks from one node to another. By storing contribution blocks in low-rank form, the total volume of communications can be reduced.

Because fronts are assembled in full-rank form, if a contribution block is in low-rank form, it has to be decompressed before the associated assembly takes place. This only increases the overall number of floating point operations by a small amount because the cost of converting to and from the BLR format is contained, as shown in Figures 3.5 and 3.6. Because it is not necessary to systematically compress the contribution blocks, in the experimental results of Section 6 the gains provided by this operation and the associated cost are always reported separately.

Factorization.

Once a front has been assembled, a partial factorization is performed in order to compute the corresponding part of the global factors. These computations have a large influence on the global performance of the software. Moreover, the storage of these factor parts becomes larger and larger as matrix size grows. It is thus critical to improve this phase. We assume that clusterings C F S for the FS variables and C N F S for the NFS variables have been computed, which give a blocking of the front. In the full-rank standard case, four fundamental tasks must be performed in order to compute the blocked incomplete factorization of the front F :

1. factor (F) L ρ,ρ L T ρ,ρ = F ρ,ρ for ρ ∈ C F S . 2. solve (S) L τ,ρ = F τ,ρ L -T ρ,ρ for ρ ∈ C F S , τ ∈ C F S C N F S , where τ > ρ (τ > ρ if F τ,τ appears after F ρ,ρ on the diagonal). 3. internal update (U) F τ,σ = F τ,σ -L τ,ρ L T σ,ρ for ρ, σ ∈ C F S , τ ∈ C F S C N F S , where τ ≥ σ > ρ. 4. external update (U) F τ,σ = F τ,σ -L τ,ρ L T σ,ρ for ρ ∈ C F S , τ, σ ∈ C N F S , where τ ≥ σ.
We add a new task which corresponds to the compression of a block defined by the clusterings:

5. compress (C) L τ,ρ U τ,α1 V T ρ,α2 for ρ ∈ C F S , τ ∈ C F S C N F S ,
where

|α 1 | = |α 2 | = numerical rank of L τ,ρ .
Due to the flexibility of the BLR format, several versions of the BLR factorization of a front can be implemented based on these five tasks, depending on the position of the compression. Figure 5.1 presents five different algorithms for the factorization of a front. Note that for the sake of simplicity, the order relations between σ, ρ and τ are now ignored. Similarly, the subscripts 1 and 2 in α 1 and α 2 will be omitted. Remember that U and V are both unitary only if a SVD is used for the compression. The first algorithm is the conventional full-rank partial factorization [START_REF] Golub | Matrix computations[END_REF] and is called FSUU, which stands for Factor, Solve, internal Update, external Update. The other four are based on the BLR representation of the front and differ on when the compression is performed. As we move down the figure, from FSUUC to FCSUU, the computational cost decreases as we are performing the compression earlier and thus making more and more of an approximation to the factorization.

We decided to implement and study Algorithm FSCUU as it gives the best compromise between savings (for both memory and flops) and robustness of the solver:

• The number of operations needed for the solve task is much lower than for the internal and external updates so that we can focus on these two latter tasks. • The approximation done in the factors is only due to approximations done at lower levels of the elimination tree. The solve task is done accurately. The algorithms in Figure 5.1 target different objectives in the context of applicative solvers. For instance, FSUUC can be used to speed-up the solve phase based on a full off-line compression of the accurate factors, which, for example, may be beneficial in Newton type solvers where multiple solves have to be done based on the same factorization.

Note that a sixth task can be added which corresponds to the compression of the CB:

6. external compress (C e ) F τ,σ U τ,α V T σ,α for τ, σ ∈ C N F S .
This task can be performed at the end of any of the algorithms in Figure 5.1 in order to decrease the active memory size as well as the amount of communication necessary to assemble the parent front in a parallel context, as explained in the previous section.

Front selection.

It is not necessary, nor efficient, to use low-rank approximations on all the frontal matrices. Indeed, to obtain good compression rates which overcome the compression costs, it is better to consider only larger fronts. How much work and memory is done and consumed in fronts larger than a given size can be easily assessed on a regular 9-point stencil by reusing some work and notations by Alan George [START_REF] George | Nested dissection of a regular finite element mesh[END_REF] on nested dissection; the results of this analysis are shown in Figure 5.2. The elimination tree directly based on a complete nested dissection is used (i.e. a nested dissection where the smallest separators have size 1). However, in practice, the nested dissection is stopped before the separators reach size 1 and the corresponding elimination tree is post-processed with different techniques, such as amalgamation, which aim at merging fronts in order to increase efficiency. The resulting tree is called the assembly tree and contains basically less fronts but larger ones. For this reason, in a practical context, the graphs shown in Figure 5.2(a) would be translated upwards. This remark should be taken into account when reading the analysis below. FSUU : no compression for ρ ∈ C F S in ascending order

F factor F ρ,ρ = L ρ,ρ L T ρ,ρ S solve L τ,ρ = F τ,ρ L -T ρ,ρ for τ ∈ C F S C N F S U internal update F τ,σ = F τ,σ -L τ,ρ L T σ,ρ for σ ∈ C F S , τ ∈ C F S C N F S U external update F τ,σ = F τ,σ -L τ,ρ L T σ,ρ
for τ, σ ∈ C N F S FSUUC : compress after updates for ρ ∈ C F S in ascending order

F factor F ρ,ρ = L ρ,ρ L T ρ,ρ S solve L τ,ρ = F τ,ρ L -T ρ,ρ for τ ∈ C F S C N F S U internal update F τ,σ = F τ,σ -L τ,ρ L T σ,ρ for σ ∈ C F S , τ ∈ C F S C N F S U external update F τ,σ = F τ,σ -L τ,ρ L T σ,ρ for τ, σ ∈ C N F S C compress L τ,ρ U τ,α V T ρ,α for τ ∈ C F S C N F S
FSUCU : compress after internal and before external updates for ρ ∈ C F S in ascending order

F factor F ρ,ρ = L ρ,ρ L T ρ,ρ S solve L τ,ρ = F τ,ρ L -T ρ,ρ for τ ∈ C F S C N F S U internal update F τ,σ = F τ,σ -L τ,ρ L T σ,ρ for σ ∈ C F S , τ ∈ C F S C N F S C compress L τ,ρ U τ,α V T ρ,α for τ ∈ C F S C N F S U external update F τ,σ = F τ,σ -U τ,α V T ρ,α V ρ,β U T σ,β
for τ, σ ∈ C N F S FSCUU : compress before updates for ρ ∈ C F S in ascending order

F factor F ρ,ρ = L ρ,ρ L T ρ,ρ S solve L τ,ρ = F τ,ρ L -T ρ,ρ for τ ∈ C F S C N F S C compress L τ,ρ U τ,α V T ρ,α for τ ∈ C F S C N F S U internal update F τ,σ = F τ,σ -U τ,α V T ρ,α V ρ,β U T σ,β for σ ∈ C F S , τ ∈ C F S C N F S U external update F τ,σ = F τ,σ -U τ,α V T ρ,α V ρ,β U T σ,β for τ, σ ∈ C N F S
FCSUU : compress before solve for ρ ∈ C F S in ascending order show that most of the factor entries are computed and almost all of the floating point operations are done within fronts of relatively large size. This is particularly interesting considering that these larger fronts only account for a very small fraction of the total number of fronts in the elimination tree, as shown in Table 5. 1. This shows that there is no need to compress small fronts for two reasons: Proportions of entries of the factors and of number of operations in the top levels of the multifrontal tree. The problem studied is a 2D Laplacian on a 9-point stencil. The corresponding matrix size is the square of the mesh size. Note that the largest fronts we look at in these plots represent very few fronts compared to the total number of fronts, as shows Table 5.1. It would be too much work for a small gain (in small fronts we can only target a few entries and a few computations, and moreover the SV D compression has a cost) ; 2. It is more critical to focus on the top of the tree because the corresponding large fronts represent most of the work to be done in the multifrontal process, in terms of computations and memory.

F factor F ρ,ρ = L ρ,ρ L T ρ,ρ C compress F τ,ρ U τ,α Z T ρ,α for τ ∈ C F S C N F S S solve L τ,ρ = U τ,α (Z T ρ,α L -T ρ,ρ ) for τ ∈ C F S C N F S i.e, L τ,ρ = U τ,α V T ρ,α with V T ρ,α = Z T ρ,α L -T ρ,ρ U internal update F τ,σ = F τ,σ -U τ,α V T ρ,α V ρ,β U T σ,β for σ ∈ C F S , τ ∈ C F S C N F S U external update F τ,σ = F τ,σ -U τ,α V T ρ,α V ρ,β U T σ,β for τ, σ ∈ C N F S
6. Experiments. This section presents experimental results on the set of problems described in Section 6.1. Our methods have been developed and incorporated in the general purpose symmetric and unsymmetric sparse linear solver MUMPS [START_REF] Mumps | MUltifrontal Massively Parallel Solver[END_REF] (in complex and real arithmetic) but can be applied to any multifrontal solver.

• For each test matrix, we study a linear system whose right-hand side is a vector of ones. To measure the effectiveness of our approach, several metrics are presented in Section 6.2. • The global ordering used for MUMPS is METIS [START_REF] Karypis | MeTiS: A Software Package for Partitioning Unstructured Graphs, Partitioning Meshes, and Computing Fill-Reducing Orderings of Sparse Matrices Version 4[END_REF], although experiments have shown comparable results with other global orderings such as AMD [START_REF] Amestoy | An approximate minimum degree ordering algorithm[END_REF]. • Only fronts larger than 100 are candidate for compression. The halo depth used during the clustering phase is 2. Inherited clustering with cluster size set up to 200 is used, except in Section 6.3 where the clustering strategy is studied. • The compression is achieved through a truncated rank-revealing QR factorization (modified from Lapack's CGEQP3) with an absolute low-rank threshold. We are currently investigating the efficiency of a relative dropping in terms of compression, stability and accuracy. • Results are obtained with a sequential code and only from the factorization phase.

For the solve phase, the flops reduction is equal to the factor memory reduction. The supercomputer used for these experiments is a SMP Altiv UV (ccNUMA architecture) equipped with 48 WESTMERE EX octo-core processors @2.67GHz and 3TB of RAM.

6.1. Set of problems. The experiments presented throughout this paper have been run on a set of problems coming from different physics applications. These matrices will be used all along this section to illustrate each aspect of our work and are described in Table 6 Curl5000 is an academic problem and corresponds to the curl-curl operator, which is widely used in computational electromagnetism. Geoazur128 matrix is a complex-valued impedance matrix resulting from the finite-difference discretization of the heterogeneous Helmholtz equation which is the second-order visco-acoustic time-harmonic wave equation for pressure p. The aim is the modeling of visco-acoustic wave propagation in a 3D viscoacoustic medium parameterized by wavespeed, density and quality factor. For more details about these matrices, see Operto [START_REF] Operto | 3D finite-difference frequency-domain modeling of visco-acoustic wave propagation using a massively parallel direct solver: A feasibility study[END_REF]. The last two matrices are thermal and mechanical simulations from the French Electricity Company (EDF). For matrix ME RAFF12, we will also study different mesh refinements of the model to analyse the impact on the LR properties of the matrix. Finally, various sizes of a 3D Laplacian problem defined on an 11-point stencil are also used to demonstrate the potential of a BLR based factorization. 6.2. Metrics. The effectiveness of the proposed techniques will be measured with different metrics mostly related to the reduction of memory and flops:

• The factor compression |L| is defined as the ratio of the number of entries needed to store the factor computed with the BLR approach over the number of entries needed to store the regular factor. • The maximum size of CB stack compression |CB| is the ratio of the maximum size of the CB stack with the BLR approach over the maximum size of the CB stack with a regular full rank approach. For a definition of the CB stack, please refer to Section 2. • For each task (or tasks combination) defined in Section 5.2, we will indicate either the corresponding absolute low-rank flop count (execution time), either the corresponding flops (execution time) compression as a percentage of the full-rank factorization (FR facto) flop count (execution time). For instance, a column called "F+S" shows absolute data or percentages related to flops or execution time (it is always indicated) related to the "Factorization" and "Solve" tasks of the algorithm.

6.3. Clustering strategies. Table 6.2 reports experimental results related to the clustering strategies. As the problem size grows, the explicit border clustering time becomes excessively time consuming. Experiments confirm that the inherited clustering is 2 to 8 times faster than the explicit border clustering, with a comparable efficiency. The overhead due to this specific phase is very low compared to the full rank analysis time. In pratice, a basic reclustering strategy (as explained in Section 4) can be performed in a negligible time. It improves slightly the compression rates and substantially the BLAS operations, which was our target. 6.4. Influence of the cluster size. By definition of the BLR admissibility condition described in Section 4.1, a target block size has to be chosen in order to define suitable admissible blocks. For a general algebraic multifrontal solver, having such a parameter can be a constraint if it has to be tuned finely to ensure a good efficiency of the method for a given problem. Although block sizes of 100 or 150 are too small to achieve good compression rates, results presented in Table 6.3 show that the dependence between compression and block size is usually quite small, so that a variation in the block size does not have a critical impact on the global efficiency of the method. This property is very desirable because it will ease the unsymmetric (and symmetric indefinite) pivoting, which can in some cases modify dynamically the number of fully summed variables within a front. From a more general point of view, block sizes can also be modified by reclustering (see Section 4), or because of out-of-core and parallelism, which are already strong algorithmic constraints in a general purpose sparse solver. Thus, the flexibility in the block sizes is very desirable to suppress a new constraint. 6.5. Memory, flops and accuracy. In this section we present global results related to memory, flops and accuracy with different low-rank thresholds. Figure 6.1 shows the cost and memory occupancy of the BLR based solver relative to the full-rank one; for instance, for problem TH RAFF7 with ε = 10 -14 , the number of operations needed for the factorization is divided by 5, the maximum size of CB stack by 4 and the factor memory by 2. In the figure, also the accuracy of the computed solution is reported, measured in terms of Componentwise Scaled Residual [START_REF] Skeel | Scaling for numerical stability in Gaussian elimination[END_REF][START_REF] Higham | Componentwise perturbation theory for linear systems with multiple right-hand sides[END_REF][START_REF] Arioli | Solving sparse linear systems with sparse backward error[END_REF], defined as CSR = max

i |b -Ax| i (|b| + |A| |x|) i .
The accuracy of the solution follows the low-rank threshold used for the compression. No propagation is observed and ε gives to the user a good control over the numerical behaviour of the solver. As far as memory compression is concerned, the method shows good efficiency at high to middle accuracy. At low accuracy, the results are even better. With the ME RAFF12 and TH RAFF7 matrices, good compression rates are also obtained at full accuracy. At full accuracy, the compression of the maximum size of CB stack (|CB|) can be quite interesting even when we do not have a global gain in terms of factor compression (|L|) and flops reduction (F+S+C+U+U flops), see Figure 6.1(a). This is probably due to the fact that the full-rank maximum size of CB stack is reached on large fronts on which the compression of the CB block is quite large. It is an interesting feature which could be exploited to limit the memory requirement for the factorization or to compress the multifrontal CB stack on demand.

In terms of computational cost compression, gains are almost always higher than those due to memory compression for two reasons:

• as shows Figure 5.2, for a given front size, the top level fronts concentrate a higher flops ratio than entries in factor ratio (with respect to the global flops and entries in factor). • operations on blocks have a cubic complexity whereas memory has a square complexity. For matrix Curl5000, with ε = 10 -14 and ε = 10 -12 , we do slightly more operations than in a usual multifrontal solver because of the additional cost for QR factorizations performed to reveal uncompressable subblocks (|L| = 100%).

6.6. Scalability with respect to the problem size. The efficiency of low-rank approximation techniques on multifrontal factorizations can strongly depend on the size of the problem. Table 6.4 illustrates this idea on different mesh refinements of the applicative prob-lem class ME RAFF*. The efficiency of the low-rank methods increases when refining the mesh. In order to further investigate this behavior, a full set of tests were conducted on matrices of increasing size from the Laplacian operator discretized with a 3D 11-point stencil. Table 6.5 shows how the effectiveness of low-rank techniques obtained with the BLR format increases with the mesh size, in terms of memory, flops and time reductions. Table 6.5 Scalability of the BLR factorization with respect to the mesh size M . Results are given for a matrix coming from the Laplacian operator discretized with a 3D 11-point stencil. The low-rank threshold is set up to 10 -14 so that a full accuracy is kept. FR facto is a full rank factorization obtained with MUMPS. The percentages are given with respect to FR facto. Ce can be viewed as the price to pay to obtain the memory reduction |CB|. See Section 6.2 for more details about the metrics. For a mesh size of 160, a factor of almost 20 is obtained in terms of flops reduction. The flop counts reveal an O(N 4/3 ) complexity, which is comparable to what is obtained in 3D with an HSS solver [START_REF] Wang | Massively parallel structured multifrontal solver for time-harmonic elastic waves in 3D anisotropic media[END_REF][START_REF] Xia | Efficient structured multifrontal factorization for general large sparse matrices[END_REF].

The overhead due to the compression of CBs appears to be globally very small, which is critical to decrease the active memory consumption. The same behavior is observed in terms of memory: for mesh size of 160, the size of the factors has been reduced by a factor of 3 and the maximum size of CB stack by a factor of almost 10.

The second part of Table 6.5 shows how the flops reduction is translated into time reduction. Results show that a large part of the flops reduction is converted into time reduction. When the flops compression rate is important, the benefits of BLAS 3 operations are partly lost because of the very low rank of the blocks, which explains why the efficiency decreases slightly. However, for large enough mesh sizes, very interesting time reductions are obtained (a factor of 7 is already obtained for mesh size 160), that are increasing with the size of the problem.

Conclusion.

A new format for low-rank multifrontal solvers has been presented. More flexible than other low-rank matrix structures, the BLR format shows a high efficiency on applicative problems, both in terms of memory and flops reduction. The low-rank threshold ε gives a full control on the global accuracy of the solver and no error propagation has been observed. We showed that this format is very suitable for the algorithmic requirements of a robust, general multifrontal solver, and presents some distinctive features over hierarchical formats such as HSS and H that make its usage more immediate and, likely, more effective. By means of an experimental comparison, we showed that the gains achieved with the BLR format are comparable to those obtained with hierarchical formats. Through this comparison we also showed that the cost of constructing the low-rank representation is much lower in the case of the BLR format with respect to the HSS and H approaches; this provides an additional property which is very favorable in the context of a sparse, multifrontal solver involving operations with complex data access patterns. We proposed a method for computing the blocking of frontal matrices for the BLR format; thanks to the properties of BLR, this technique does not require any knowledge of the geometry of the problem and can be run at a cost which is marginal relative to the cost of the analysis phase of a multifrontal solver.

In order to evaluate the potential of the proposed techniques, we integrated the BLR format and the blocking technique in the sequential version of the MUMPS solver and we ran experiments on a number of problems from real world applications as well classical, textbook problems. The experimental results show that considerable gains can be achieved already at very accurate approximation levels, in which case the final solution backward error is comparable to the one obtained with a standard, full-rank solver. Experimental results also show that the BLR format is tolerant with respect to variations in the size of the blocksa property that may be used to accommodate the BLR format of frontal matrices to data distribution in a parallel environment as well as to classical pivoting techniques.

The work presented in this document lays the foundations of an ongoing research effort which aims at exploiting the BLR format within a parallel, fully featured multifrontal solver.

  One level of nested dissection based on a 5×5 square grid mesh. A is reordered with respect to the nested dissection.

Figure 2 . 1 .

 21 Figure 2.1. One level of nested dissection.

Figure 2 . 2 .

 22 Figure 2.2. A general nested dissection and its corresponding separator tree.

  Before elimination: the front is assembled. FS variables are ready to be eliminated. NFS variables will be updated. After elimination: FS variables are eliminated, CB is the Schur complement that will assembled into the parent front.

Figure 2 . 3 .

 23 Figure 2.3. Structure of a front F .

2 .

 2 Definition 3.2 ([7,[START_REF] Börm | Efficient Numerical Methods for Non-local Operators[END_REF] Admissibility condition for elliptic PDEs). Let b = σ × τ be a block cluster. b is admissible if

  Correlation between graph distance and full accuracy block rank.

Figure 3 . 1 .

 31 Figure 3.1. Two illustrations of the admissibility condition for elliptic PDEs.

add b 1 , b 2 , b 3 and b 4 to list 9

 129

. 4 .

 4 Based on the same block clusters as in example Figures(3.

  2) and (3.3), we show in Equation (3.1) an example of a 4 × 4 BLR matrix. Remember that U and V are not necessarily unitary.

Figure 3 . 4 .

 34 Figure 3.4. Illustration of a BLR matrix of a dense Schur complement of a 128 × 128 × 128 Laplacian problem with a low-rank threshold set up to 10 -14 . The corresponding clustering of its 128 × 128 plan graph into 4 × 8 = 32 blocks is also given.

  Figures 3.5(b) and 3.6(b)

10

 10

Figure 3 . 5 .

 35 Figure 3.5. Comparison between H, HSS and BLR formats on a 128 3 Geozur problem. The matrix considered is the factorized Schur complement associated with the top level separator of a perfect hand generated nested dissection tree. The clustering is hand-computed and geometrical.

Figure 3 . 6 .

 36 Figure 3.6. Comparison between H, HSS and BLR formats on a 128 3 Laplacian problem. The matrix considered is the factorized Schur complement associated with the top level separator of a perfect handcomputed nested dissection tree. The clustering is hand-computed and geometrical.

Figure 4 . 1 .

 41 Figure 4.1. Halo-based partitioning of G S with depth 1. G H is partitioned using METIS.

Figure 4 .

 4 2(a) illustrates the way the clusters are computed using the halo method for both the separator and the border.

Figure 4 . 2 .

 42 Figure 4.2. Two different ways to perform the clustering of a border

Figure 4 . 3 .

 43 Figure 4.3. Relation between inherited clustering and front blocking. The small clusters correspond to clusters located in corners in Figure 4.2(b). The large clusters correspond to optimal clusters which are integrally kept in the current front. Note that not all the large clusters of Figure 4.2(b) are represented here.

Figure 5 . 1 .

 51 Figure 5.1. Standard and BLR factorizations of a front. The flexibility of the BLR format allows to define 4 different BLR factorization algorithms. When two different compressed blocks have to be used, variable β is used in addition to α.

Figures 5 .

 5 Figures 5.2(a) and5.2(b) show that most of the factor entries are computed and almost all of the floating point operations are done within fronts of relatively large size. This is particularly interesting considering that these larger fronts only account for a very small fraction of the total number of fronts in the elimination tree, as shown in Table5.1. This shows that there is no need to compress small fronts for two reasons:

  of the e ntrie s of the factor compute d in fronts large r than N Proportion of entries of the factors computed in the top levels of the multifrontal tree. of ope rations pe rforme d in fronts large r than Proportion of the total number of operations performed for the partial factorization of the top levels of the multifrontal tree.

Figure 5 .

 5 Figure 5.2.Proportions of entries of the factors and of number of operations in the top levels of the multifrontal tree. The problem studied is a 2D Laplacian on a 9-point stencil. The corresponding matrix size is the square of the mesh size. Note that the largest fronts we look at in these plots represent very few fronts compared to the total number of fronts, as shows Table5.1.

Figure 6 . 1 .

 61 Figure 6.1. Global results of the BLR methods implemented within MUMPS. The solver is fully algebraic.Compression results are shown for the memory needed to store the factors and the contribution blocks, and for the number of operations needed to perform the factorization (left-hand side vertical axis). Results related to the accuracy of the computed solution (right-hand side vertical axis) are also given. For the exact definition of the metrics used in these graphs, please refer to Section 6.2.

Table 3 . 1

 31 Cost of dense and low-rank basic linear algebra operations. A, L, A 1

  • The L 21 submatrix can be easily represented in BLR format with a block clustering induced by those of the L 11 and L 22 submatrices.

• Pivoting techniques seem hard and inefficient to use when factorizing a matrix stored in a hierarchical format; this has led researchers to employ different types of factorizations, e.g. a U LV factorization

[START_REF] Chandrasekaran | A fast ULV decomposition solver for hierarchically semiseparable representations[END_REF][START_REF]Fast algorithms for hierarchically semiseparable matrices[END_REF]

. BLR format is more naturally suited for applying partial threshold pivoting techniques within a standard LU factorization.

  , other admissibility conditions have been studied. They only need the graph of the matrix G and are thus called black box methods. Definition 4.1 details this graph version of Definition 3.2. Definition 4.1 ([8, 19] Black box admissibility condition). Let b = σ × τ a block cluster where σ and τ are sets of graph nodes. b is admissible if

Table 5 . 1

 51 Proportion of fronts larger than N with respect to the total number of fronts in the elimination tree, for different mesh sizes.

				mesh size
		2 10	2 11	2 12	2 13	2 14	2 15
	N = 200	1.61	1.68	1.71	1.73	1.74	1.74
	N = 400	0.37	0.40	0.42	0.42	0.43	0.43
	N = 600	0.17	0.20	0.21	0.22	0.23	0.23
	N = 800	0.07	0.09	0.10	0.11	0.11	0.11
	N = 1000 0.05	0.06	0.06	0.06	0.06	0.06
	1.

Table 6 . 1

 61 .1. They correspond to four important classes of applications within the field of elliptic PDEs. Set of problems used for the experimentations. They are finite-difference or finite elements methods simulations. CSR = Componentwise Scaled Residual = max

	Name	Prop.	Arith.	N	NZ	factors	flops CSR	application
				(×10 6 ) (×10 9 ) storage (×10 12 ) (full rank)	field
	Curl5000	2D/sym.	D	50	0.2	27 GB	5 2 × 10 -15	electromag.
	Geoazur128 3D/unsym.	Z	2		55	46 GB	62 2 × 10 -12	wave prop.
	TH RAFF7	3D/sym	D	8	118 138 GB	100 8 × 10 -15	thermal
	ME RAFF12	2D/sym.	D	134		1 215 GB	200 4 × 10 -15	mechanical
	precision complex.				i	|b -Ax| i (|b| + |A| |x|) i	. D=double precision real, Z=double

Table 6 . 2

 62 Comparison of the efficiency of the explicit and inherited clustering strategies. Each inh column gives results with the clustering done during the analyis. Each exp column gives results with the clustering done during the factorization. The clustering time is the total time spent for the clustering of all the fronts. The low-rank threshold is 10 -8 . The full rank analysis time is the time spent in the analysis phase for a regular full rank multifrontal factorization (excluding clustering time).

			memory		flops			time
		|L|		|CB|		F+S+C+U+U	clustering	full rank
	clustering	inh	exp	inh	exp	inh	exp	inh	exp	analysis
	Curl5000	63.7% 62.4%	7.0%	5.5%	10.9% 11.1%	13 s	27 s	897 s
	Geoazur128	79.0% 77.0%	47.0% 45.0%	60.8% 59.1%	5 s	42 s	62 s
	TH RAFF7	34.1% 30.7%	17.5% 16.2%	7.2%	6.6%	34 s 206 s	387 s
	ME RAFF12	52.9% 51.1%	4.8%	4.1%	6.1%	6.1%	121 s 239 s	1971 s

Table 6 . 3

 63 Influence of the block size used for the BLR compression on the efficiency of the BLR multifrontal method. Study for Geoazur128 with basic reclustering and low-rank threshold set to 10 -8 .

		memory		flops
	block size	|L| |CB| F+S+C+U+U F+S+U+U	C
	100 88%	55%	75%	4.6E+13	9.8E+11
	150 83%	50%	68%	4.1E+13	1.5E+12
	200 79%	46%	63%	3.7E+13	1.9E+12
	250 76%	45%	60%	3.4E+13	2.4E+12
	300 75%	44%	58%	3.3E+13	2.7E+12
	350 74%	44%	56%	3.2E+13	3.1E+12
	400 73%	45%	56%	3.1E+13	3.4E+12
	450 73%	44%	56%	3.1E+13	3.7E+12
	500 73%	44%	57%	3.1E+13	4.4E+12

Table 6 . 4

 64 Overview of the influence of the mesh refinement degree on the memory and flops reductions. The low-rank threshold is set up to 10 -8 .

			memory	flops
	Refinement	N	|L| |CB| F+S+C+U+U
	ME RAFF8	2, 101, 258 71%	27%	32%
	ME RAFF11	33, 570, 826 59%	11%	12%
	ME RAFF12	134, 250, 506 52%	7%	6%

RR n°8199

Inria

Acknowledgment. We wish to thank Xiaoye S. Li (LBNL) and Artem Napov (LBNL) for the helpful discussions on low-rank techniques and for the experimental HSS code used in Section 3. This work was granted access to the HPC resources of CALMIP under the allocation 2012-p0989 and is funded by EDF under collaboration contract 8610-AAP-5910070284.

RESEARCH CENTRE GRENOBLE -RHÔNE-ALPES


Sébastien Zamith

Grégoire De Tournadre

Pierre Labastie

Jean-Marc L'hermite

Attachment cross-sections of protonated and deprotonated water clusters

Attachment cross-sections of water molecules onto size selected protonated and deprotonated water clusters have been measured in the size range n=30-140 for 10 eV kinetic energy of the clusters in the laboratory frame. Within our experimental accuracy, the attachment cross-sections are found to have the same magnitude and size dependence for both species. It is shown that electrostatic interactions are likely to play a role even for the largest sizes investigated.

I.

Introduction

Charged particles play a crucial role in the nucleation of droplets in the atmosphere [1].

Due to electrostatic interactions, charged particles promote the formation of droplets. Since the cloud chamber experiment of Wilson back in 1897 [2], ions of different sign have been shown to induce significantly different nucleation rates [3,4,5,6,7,8,9]. Negatively charged water clusters seem to nucleate more efficiently than positively charged species; the origin of this sign preference for nucleation is still controversial [10,11,12,13,14,15]. Are higher attachment cross sections for anions than for cations responsible for this sign preference?

The present paper brings a contribution to this issue. We compare the attachment crosssection of protonated and deprotonated water clusters, for sizes ranging from n=30 to 140. Attachment rate and cross-sections are measured for clusters colliding with water molecules at 10 eV collision energy in the laboratory frame. We find the cross-sections to be essentially identical for anions and cations.

The paper is organized as follow: Section II gives a rapid overview of the experimental setup. Section III details how the cross-sections and rates are deduced from our measures and Section IV presents and discusses the experimental results. Conclusions are drawn in Section V.

II. Experimental setup

The experimental setup and procedure used here to measure attachment cross-sections of deprotonated water clusters are essentially the ones used in ref. 16. Briefly, both protonated and deprotonated water clusters are produced in a gas aggregation source.

Depending on the polarity of the voltages used in the setup one can study either positively or negatively charged clusters. After their production, the clusters seeded in helium carrier gas are thermalized at 25 K in a heat bath. The clusters are then mass selected and slowed down to 10 eV kinetic energy in the laboratory frame [17]. They enter a collision cell with a controlled pressure of water vapor. At the output of the collision cell, collision products are mass analyzed with a time of flight mass spectrometer. Our raw experimental data consist of time of flight mass spectra with one peak corresponding to the parent cluster and several additional peaks corresponding to the attachment of water molecules onto the parent cluster. We have shown in ref. 16 that no evaporation occurs after attachment and until the clusters reach the detector, in the present experimental conditions. As will be discussed below, absolute attachment rates or cross-sections can be directly deduced from such mass spectra.

In order to be able to compare as accurately as possible the cross-sections for protonated and deprotonated water clusters, measurements were performed for some sizes alternately for both species in the same runs, without changing anything in the setup but the polarity of the high voltages. The measurements were reproduced on different days under different source conditions and pressures in the cell. We thus achieve a much better precision for the measurement of the ratios of the cross sections of anions and cations than for the cross sections themselves.

III. Attachment rate and cross-section

From our measured mass spectra one can easily deduce the attachment rates or collision cross-sections. The attachment rate is simply given by:

Eq. 1 where is the cluster velocity in the laboratory frame, is the density in the cell and is the cell length. is the intensity of the peak of intact parent cluster and is the total intensity integrated over all the peaks.

The attachment rate is related to the cross section by: Eq. 2

where [START_REF] Chapman | The mathematical theory of non-uniform gases[END_REF]:

Eq. 3
is the distribution of relative velocities between the cluster and the colliding molecules. is the temperature of the cell, the Boltzmann constant and is the mass of the water molecule. Furthermore, if one assumes that the cross-section is a slowly varying function of , the cross-section deduced from Eq. 1 and Eq. 2 is:

Eq. 4

where is the average relative velocity given by [START_REF] Chapman | The mathematical theory of non-uniform gases[END_REF]:

Eq. 5

with .

In the following both the attachment rates and cross-sections will be presented. It is advantageous to compare directly the experimental attachment rates to theoretical ones since no assumption has to be made about the cross-section.

IV. Experimental Results and Discussion

In Figure 1 are presented the measured attachment cross-sections for protonated and deprotonated water clusters, and Figure 2 shows the corresponding attachment rates.

At first we should mention that compared to our previous results [16,[START_REF] Zamith | [END_REF] we find crosssections roughly 30% higher than previously measured. A slightly different, more accurate data processing has been performed and it turned out that results of ref. Fig. 1 and Fig. 2 show that the attachment cross sections and the attachment rates, respectively, of protonated water clusters are close to those measured for deprotonated clusters. As expected, cross-sections and attachment rates increase with the cluster size. The attachment probabilities of protonated clusters seem to be slightly lower than those of deprotonated clusters at small sizes, but this trend is not marked enough, within our experimental accuracy, to draw a definite conclusion about a possible charge effect. In any case, the tiny difference that might exist between the attachment probabilities of cations and anions discard higher attachment cross sections of anions as a possible cause of the sign preference for nucleation (in the frame of classical nucleation theory indeed, the nucleation rate depends only weakly on the attachment cross section [20,21]).

Before analyzing the effect of electrostatic interactions between the impinging molecule and the cluster, we will first compare our experimental results with the simplest model, considering the clusters and the water molecules as hard-spheres. The geometric crosssection is given by [START_REF] Levine | Molecular Reaction Dynamics and Chemical Reactivity[END_REF]:

Eq. 6
where is the cluster radius and is the molecular radius. The rate is then:

Eq. 7
It is difficult to have a reliable estimate of both the molecular and cluster radius. One can assume that the clusters have the same density as the bulk and deduce the radius of a cluster of size as . The density can be taken to be the density of hexagonal ice, liquid water or amorphous ice for instance. These three densities yield values for in the range 1.93-1.98 Å. However, clusters might have densities differing from the bulk values. In the size range considered here, clusters are expected to be amorphous [START_REF] Torchet | [END_REF]24]. In ref. [24], the authors deduce from theoretical calculations a diameter of 2.7 nm for a cluster containing 293 molecules. This yields a value =2 Å. Furthermore from the distribution of molecular distances from the cluster center of mass (figure 24 of ref. 24), one can deduce values for the molecular radius around 1.9-2.0 Å in the size range n=22-931.

From our measurements, if we assume that for the bigger sizes the attachment cross-section is given by the hard-sphere model, we deduce a molecular radius of 2.28 Å (red line in figure 1).

Recently Lengyel et al [25] have measured attachment cross-sections for neutral water clusters. For an average size of n=260, they have found a value (955 Å 2 ) very similar to the one measured in our group for protonated water clusters (980 Å 2 for n=250 [16]) for similar collision energies. However in [25] the authors claim that they observe a cross-section which is about 2.5 times bigger than the geometric cross-section, taking a value of 1.77 Å for and calculating the geometric cross-section as

. The corresponding geometric cross-section is plotted as a dotted line in Fig. 1. The correct definition for the geometric cross-section should be the one given by Eq. 6

and we believe that the value taken by Lengyel et al is probably a low estimate.

In reference [25] the collisions occur between neutrals whereas we studied ion-neutral collisions: One should therefore expect higher cross-sections in our experiment, due to iondipole and ion-induced dipole attractive interactions. The fact that similar cross-sections are found in both experiments, for neutral and charges species, would tend to indicate that, in this size range and this collision energy range, electrostatic interactions do not play a major role. In this case, provided that clusters are approximately spherical, the attachment crosssections can be approximated to a good accuracy by hard sphere cross-sections.

However, Lengyel et al observe a variation of the cross-section when using different collision partners. These observations are rather in favor of a cross-section still dominated by long range interactions. On the other hand in [16,[START_REF] Zamith | [END_REF] we observed that as the size increases, the cross-section no longer depend on the collision energy. This observation is in favor of our argument that the cross-section becomes dominated by the geometric size of the clusters for the higher sizes.

In order to evaluate the enhancement of the cross-sections (or the rates) due to electrostatic interactions, we calculated these cross sections in the frame of a Langevin model, using the average dipole orientation (ADO) theory [26]. In this approach the interaction between the charged cluster and the neutral molecule is given by: Eq. 8

The first term represents the charge-induced dipole interaction, the second term is the charge-permanent dipole interaction, and the third term is the centrifugal barrier. is the distance between the center of the cluster and the molecule, E c is the collision energy in the center of mass frame and is the impact parameter. The polarizability of water is = 9.93 a.u. [27] and the dipole moment of water equals 0.73 a.u. [28]. is the dipole locking constant introduced by Su et al [26] in ADO theory. In order to introduce size effects, we used the models developed by Kummerlöwe et al [29]. In this paper two models are considered: the hard sphere average dipole orientation (HSA) and the surface charge capture (SCC) model. In the SCC model the charge is mobile at the surface of the cluster whereas in the HSA model it is localized at the center of the cluster. Those two models give lower and upper bounds for the effect of the electrostatic interactions, which are likely to be underestimated in the HSA model and overestimated in the SCC model. The attachment rate in the frame of these models is calculated as:

Eq. 9

where the expressions for are given in ref. 29. The rates obtained from a numerical integration of Eq. 9 are plotted together with the experimental rates in figure 2.

The theoretical rates have been calculated assuming a molecular radius of 1.98 Å and a dipole locking constant 215 . 0  l C [26]. For large sizes, the overall result of these calculations are not in contradiction with our experimental observations, since experimental attachment rates fall in between the predictions of the HSA and the SCC models, which are expected to underestimate and overestimate, respectively, the electrostatic interactions thus the attachment rates. But none of those two models can explain our observed variation of the attachment rates with size. The only way we have found so far to reproduce our results is to invoke dynamical effects [16]: when the collision duration becomes shorter than the time required for redistributing the collision energy among the internal degrees of freedom of the cluster, the attachment probability is reduced.

Our experimental cross-sections are about 1.3-1.4 higher than the geometric crosssections assuming Å. So electrostatic interactions probably still play a role even at the largest cluster sizes investigated here. However, assuming a slightly bigger molecular radius would reduce the deduced role of electrostatic interactions. Finally, despite our effort to measure reliable cross-sections, it is hard to completely dismiss the possibility that the pressure reading is off by an amount that would suffice to explain the deviation from the hard sphere cross-section at large sizes.

The modeling of particles growth requires the knowledge of attachment cross-sections [20,21]. Concerning water nucleation in the atmosphere, the cross sections of the particles is often estimated in the frame of the hard sphere model using the ice density [30]. Our measures indicate that more realistic cross-sections should be used. However, these cross sections cannot be directly deduced from the present work since our measurements are not made at temperatures usually present in the atmosphere. Indeed, our lowest collision energy in the present study corresponds to mean kinetic energy of water molecules at a temperature of about 900 K. Nevertheless, we can try to extrapolate our results to small cluster sizes at room temperature. To do so, we assume that the cross-section is given by either the HSA or SCC model multiplied by a factor due to dynamical effects, namely (see ref. 16). Fig. 3 presents the ratio of cross-sections calculated in this way to geometric cross-section. This crude model predicts an enhancement of the cross-section compared to the geometric cross-section, with a maximum for sizes around 10. As the size increases, the cross-section is a factor in between 1.2 and 2 larger than the geometric crosssection. Quite surprisingly, the size dependence of the calculated cross sections does not strongly depart, except for very small clusters, from geometric cross sections. This behavior is due to two conjugate effects that almost counterbalance: on the one hand, the electrostatic enhancement of the cross sections at low energy, on the other hand their reduction by the dynamical effect presented in reference [START_REF] Zamith | [END_REF].

V. Conclusion

We have made a direct comparison of attachment cross-sections for protonated and deprotonated water clusters. Within our experimental conditions, we observe no significant effect of the sign of the charge of the clusters. The predominant effect is due to previously invoked dynamical effects [16,[START_REF] Zamith | [END_REF]. At large sizes, the attachment cross sections are higher than geometric cross sections, which seems to indicate that electrostatic interactions enhance the attachment probability. On the other hand, our cross sections being close to the ones obtained for neutral species by J. Lengyel and coworkers [25] suggests that electrostatic interactions do not play a major role in the attachment process. This contradiction could not be removed by analyzing our results in the frame of a Langevin model; clearly, the problem of the attachment of water molecules onto water clusters still needs further investigations. Ratios of calculated attachment cross-section to geometric cross-section. The attachment crosssection is calculated using the HSA (dashed line) and SCC (full line) cross-sections and by taking into account dynamical effects.

Figure captions

16 should be multiplied by a factor 1 . 3 .

 13 The conclusions of ref. 16 remain the same however: the model presented in ref. 16 still reproduces satisfactorily our results.

Figure 1 :

 1 Figure 1: (color online)

Figure 2 :

 2 Figure 2: (color online)

Figure 3 :

 3 Figure 3: (color online)

Figure 3

 3 Figure 1