Antoine Gérard
email: antoine.gerard@inria.fr

Annabelle Collin

Yves Coudière

Validation of a sequential data assimilation method applied to cardiac electrophysiology

Improved electroanatomical recordings and imaging capabilities prompts for trying to personalize cardiac electrical models. Though it is challenging, data assimilation is a family of methods relevant to this problem. It consists in estimating the state of the system thanks to observations. Among the two main approaches, variationnal and sequential, we choose to investigate the possibilities of a sequential method based on an atrial model and atrial electroanatomical maps. The method combines a state observer and a Kalman filter to estimate some parameters in the model.

First, we simulate a sinus rythm coupled to different pacing scenarios using a validated realistic atria model, see [3], summarized as :

8 > > > > > > > < > > > > > > > : A m C m ∂ u ∂t (k) + f (u (k) ,t) ! div(s (k) u (k) ) = ( 1) k g(u (1) u (2) ) (1) ∂ w ∂t (k) = g CRN (u (k) , w (k) )
The anatomical geometry was segmented from clinical MRI, while the electrophysiology is based on the Courtemanche ionic model and monodomain equations on coupled surfaces (endo and epicardum, labeled k = 1, 2) so as to efficiently represent transmural effects. Activation times were then recorded in the left atria, and provide a number of in silico electroanatomical maps.

For the data assimilation work, these maps are projected on a coarser mesh of the left atria only. We applied to a simpler monodomain model the Luenberger observer devised by A. Collin & al. [2]. In addition the conductivity parameters were estimated using a Reduced-order-Unscented-Kalman-Filter (RoUKF) [4] in order for the simplified model to fit the recorder maps. Eq. 2 states both the observer and the filter :

8 > > > > > < > > > > > : A m ✓ C m ∂ ũ ∂t + f ( ũ,t) ◆ div(s -ũ) = G L (z, ũ) ∂ w ∂t + g MS ( ũ, w) = 0 (2) q (t) = G q (D(z, ũ)) FIGURE 1 -Reconstructed activation map
Firstly, we tried to estimate a unique parameter a such that the conductivity tensor s = a Id is isotropic. Secondly, we tried to estimate two conductivities parameters g l and g t such that the conductivity tensor is s (x) = g t Id +(g l g t ) f (x) ⌦ f (x), where the unit vector f (x) points towards the fiber direction.

In order to validate the data assimilation technique, we compared the activation maps coming from a monodomain simulation with estimated parameters, to the originally recorder activations maps.

All simulations were done using Verdandi library [1] developed by M3DISIM Research team and the software CEPS developed by Carmen Research team.

Aknowledgments :

The authors would like to thanks Inria research team M3DISIM for their help with Verdandi software. This study recieved financial support from the French program managed by the ANR, Grant reference ANR-10-IAHU-04.

1

Proceedings of the V Virtual Physiological Human conference: 5-7 September 2018, Zaragoza, Spain


H U Hsien-Tsung
email: hsien-tsung.hu@umontpellier.fr

Cabrolier Pierre

Alkadri Ahmad

Carlier Capucine

Gril Joseph

Langbour Patrick

Beyond « Clear Wood »: Exploring the Structure-Properties of Figured Woods and of Woods with Strong Gradients Selected in Craftsmanship

published or not. The documents may come

A paradox between the viewpoints of engineering sciences and of the artisanal usage for wood is posed:

• In natural, homogeneous wood is singularity

• From engineering viewpoint, a good quality wood is orthotropic(clear wood), with homogeneous character.

• To the artisanal usage, heterogeneity and/or singularity could be seemed as defect or quality.

The goal for this program is to fulfil the gap of knowledge of heterogeneous woods(e.g. figured wood and grain deviation).

•

  Beyond " clear wood " : exploring the structure-properties of figured woods and of woods with strong gradients selected in craftsmanship Hsien-Tsung Hu, Pierre Cabrolier, Ahmad Alkadri, Capucine Carlier, Joseph Gril, Patrick Langbour, Iris Brémaud To cite this version: Hsien-Tsung Hu, Pierre Cabrolier, Ahmad Alkadri, Capucine Carlier, Joseph Gril, et al.. Beyond " clear wood " : exploring the structure-properties of figured woods and of woods with strong gradients selected in craftsmanship. 7èmes journées du GDR Sciences du bois, Nov 2018, Cluny, France. hal-C.Exploration of botanical diversity and variability to determine the certain singularities and mechanical characteristics. D.Reducing the problem according to the typology of the observation scales relevant for the multiscale analysis. Experimental characterization (and modelling) of the structural mechanical repercussions in conditions with high gradients of properties.E. Study of the visio-tactile sensory perception by craftsmen, especially for figured woods.

  to present a beginning Ph.D. research for the study of wood with artisanal usage with non-straight or heterogeneous.

  Alkadri, A., C. Carlier, I. Wahyudi, J. Gril, P. Langbour and I. Brémaud (2018). "Relationships between anatomical and vibrational properties of wavy sycamore maple." IAWA Journal 39(1): 63-86 • Brémaud, I. (2016). Vers une vision « ethno-bio-mécanique » des bois -La rhéologie du matériau-bois à l'interface entre diversité botanique et usages en lutheries et artisanats. HDR Univ. Montpellier. 149pp • Cabrolier, P., Description et comportement mécanique des bois contrefilés, 2007, Université Montpellier 2 (Sciences et Techniques) EXPLORATION 1910 BORTHWICK AW: The production of adventitious roots and their relation to bird's-eye formation (Maserholz) in the wood of various trees. 1933 HOLMBERG LH: Is suppression a possible cause of bird's-eye in sugar maple? PILLOW MY: Detection of figured wood in standing trees. 1969 CHAFE SC: Dimpled grain in wood. CARLQUIST S: Wood anatomy of Lobelioideae(Campanulaceae)."Figured wood" appeared in journal of timber industries, veneer, furniture making, commercial brochures, etc., … 1952 BAUMER M: A propos des figures du bois.

  1977 BEALS HO & DAVIS TC: Figure in wood: An illustrated review. 1979 HEJNOWICZ Z: The common basis of wood grain figures is the systematically changing orientation of cambial fusiform cells. Anatomy observation, genetic relationship, appearance of figured for CAD are focused on… Measuring instrument system(Vybris) for the dynamic modulus properties, developed in LMGC for a decade's researches of the anisotropic characteristics for the wood with strong gradients(interlocked grain, wavy grain, curly grain...) : Description et comportement mécanique des bois contrefilés 2018 Alkadri A et al: Relationships between anatomical and vibrational properties of wavy sycamore maple Mechanical properties of microspecimens; comparing it with the multi-scale prediction and the further step: macro-scopic behavior Calculating the ratio of figured tissue and non-figured tissue Correlation for micro-, meso-and macro-scopic en.woodtec.co.jp/products/individuality-tree/ 7 e G D R d u B O I S | C l u n y , F r a n c e | 2 0 -2 2 , N o v , 2


Patrick R Amestoy
email: patrick.amestoy@enseeiht.fr

Jean-Yves L'excellent
email: yves.l.excellent@ens-lyon.fr

Gilles Moreau
email: gilles.moreau@ens-lyon.fr

ON EXPLOITING SPARSITY OF MULTIPLE RIGHT-HAND SIDES IN SPARSE DIRECT SOLVERS

Keywords: sparse linear algebra, sparse matrices, direct method, multiple sparse right-hand sides AMS subject classifications. 05C50, 65F05, 65F50

The cost of the solution phase in sparse direct methods is sometimes critical. It can be larger than that of the factorization in applications where systems of linear equations with thousands of right-hand sides (RHS) must be solved. In this paper, we focus on the case of multiple sparse RHS with different nonzero structures in each column. In this setting, vertical sparsity reduces the number of operations by avoiding computations on rows that are entirely zero, and horizontal sparsity goes further by performing each elementary solve operation only on a subset of the RHS columns. To maximize the exploitation of horizontal sparsity, we propose a new algorithm to build a permutation of the RHS columns. We then propose an original approach to split the RHS columns into a minimal number of blocks, while reducing the number of operations down to a given threshold. Both algorithms are motivated by geometric intuitions and designed using an algebraic approach, so that they can be applied to general systems. We demonstrate the effectiveness of our algorithms on systems coming from real applications and compare them to other standard approaches. Finally, we give some perspectives and possible applications for this work.

Introduction. We consider the direct solution of sparse systems of linear equations ( 1)

AX = B,
where A is an n × n sparse matrix with a symmetric structure and B is an n × m matrix of righthand sides (RHS). When A is decomposed in the form A = LU with a sparse direct method [START_REF] Duff | Direct Methods for Sparse Matrices[END_REF], the solution can be obtained by triangular solves involving L and U . In this study, we are interested in the situation where not only A is sparse, but also B, with the columns of B possibly having different structures, and we focus on the efficient solution of the forward system

(2)

LY = B,
where the unknown Y and the right-hand side B are n × m matrices. We will see in this study that the ideas developed for Equation [START_REF] Amestoy | On computing inverse entries of a sparse matrix in an out-of-core environment[END_REF] are indeed more general and can be applied in a broader context. In particular, they can be applied to the backward substitution phase, in situations where the system U X = Y must be solved for a subset of the entries in X [START_REF] Amestoy | On computing inverse entries of a sparse matrix in an out-of-core environment[END_REF][START_REF] Tz | Parallel triangular solution in the out-of-core multifrontal approach for solving large sparse linear systems[END_REF][START_REF] Yeung | Interactively cutting and constraining vertices in meshes using augmented matrices[END_REF][START_REF] Yeung | AMPS: an augmented matrix formulation for principal submatrix updates with application to power grids[END_REF]. In direct methods, the dependencies of the computations for factorization and solve operations can be represented by a tree [START_REF] Liu | The role of elimination trees in sparse factorization[END_REF], which plays a central role and expresses parallelism. The factorization phase is usually the most costly phase but, depending on the number of columns m in B or on the number of systems to solve with identical A and different B, the cost of the solve phase may also be significant.

As an example, electromagnetism, geophysics, or imaging applications can lead to systems with sparse multiple right-hand sides for which the solution phase is significantly more costly than the factorization [START_REF] Amestoy | Fast 3D frequency-domain full waveform inversion with a parallel Block Low-Rank multifrontal direct solver: application to OBC data from the North Sea[END_REF][START_REF] Shantsev | Large-scale 3D EM modeling with a Block Low-Rank multifrontal direct solver[END_REF]. Such applications motivate the algorithms presented in this study.

It is worth considering a RHS as sparse when doing so improves performance or storage compared to the dense case. The exploitation of RHS sparsity (later extended to reduce computations when only a subset of the solution is needed [START_REF] Tz | Parallel triangular solution in the out-of-core multifrontal approach for solving large sparse linear systems[END_REF][START_REF] Yeung | Interactively cutting and constraining vertices in meshes using augmented matrices[END_REF]) was formalised by Gilbert [START_REF] Gilbert | Predicting structure in sparse matrix computations[END_REF] and Gilbert and Liu [START_REF] Gilbert | Elimination structures for unsymmetric sparse LU factors[END_REF], who showed that the structure of the solution Y from Equation (2) can be predicted from the structures of L and B. From this structure prediction, one can design mechanisms to reduce computation. In particular, tree pruning suppresses nodes in the elimination tree involving only computations on zeros. When solving a problem with multiple RHS, the preferred technique is usually to process all the RHS in one shot. The subset of the elimination tree to be traversed is then the union of the different pruned trees (see Section 2). However, when the RHS have different structures, extra operations on zeros are performed. In order to limit these extra operations, several approaches may be applied. The one that minimizes the number of operations processes the RHS columns one by one, each time with a different pruned tree. However, such an approach is not practical and leads to a poor arithmetic intensity (e.g., it will not exploit level 3 BLAS [START_REF] Dongarra | A set of level 3 basic linear algebra subprograms[END_REF]). In the context of blocks of RHS with a predetermined number of columns per block, heuristics based on a specific order of the columns or on hypergraph partitioning have been proposed to determine which columns to include within which block. The function to minimize might be the volume of accesses to the factor matrices [START_REF] Amestoy | On computing inverse entries of a sparse matrix in an out-of-core environment[END_REF], or the number of operations [START_REF] Yamazaki | On partitioning and reordering problems in a hierarchically parallel hybrid linear solver[END_REF]. When possible, large sets of columns, possibly the whole set of m columns, may be processed in one shot. Thanks to the different sparsity structure of each column of B, it is then possible to work on less than m columns at most nodes in the tree. Such a mechanism has been introduced in the context of the parallel computation of entries of the inverse [START_REF] Amestoy | Parallel computation of entries of A -1[END_REF], where at each node, computations are performed on a contiguous interval of RHS columns.

In this work, we propose an approach to permute B that significantly reduces the amount of computation with respect to previous work. We then go further by dividing B into blocks of columns. However, instead of enforcing a constant block size, we aim at minimizing the number of blocks, which keeps flexibility for the implementation and may improve arithmetic intensity. Because of tree pruning, RHS sparsity limits the amount of available tree parallelism. Therefore, when possible, our heuristics aim at choosing an approach that maximizes tree parallelism.

This paper is organized as follows. Section 1 presents the general context of our study and Section 2 describes the tree pruning technique arising from Gilbert results, and node intervals, where different intervals of columns may be processed at each node of the tree. In Section 3, we introduce a new permutation to reduce the size of such intervals and thus limit the number of operations, first using geometric considerations for a regular nested dissection ordering, then with a purely algebraic approach that can be applied for the general case. We call it the flat tree algorithm because of the analogy with the ordering that one would obtain when "flattening" the tree. In Section 4, an original blocking algorithm is then introduced to further improve the flat tree ordering. It aims at defining a limited number of blocks of right-hand sides to minimize the number of operations while preserving parallelism. Section 5 gives experimental results on a set of systems coming from two geophysics applications relying on Helmholtz or Maxwell equations. Section 6 discusses adaptations of the nested dissection algorithm to further decrease computation and Section 7 shows why this work has a broader scope than solving Equation (2) and presents possible applications.

1. Nested dissection, sparse direct solvers, and triangular solve. Ordering the variables of A according to nested dissection [START_REF] George | Nested dissection of a regular finite-element mesh[END_REF] limits the operation count and storage of sparse direct methods. We introduce a 3 × 3 × 3 domain in Figure 1(a). It is first divided by a 3 × 3 constant-x plane separator u 0 and each subdomain is then divided recursively. By ordering the separators after the subdomains, large empty blocks appear in Figure 1(b). The elimination tree [START_REF] Liu | The role of elimination trees in sparse factorization[END_REF] represents the dependencies between computations. It can be compressed thanks to the use of supernodes, leading to the separator tree of Figure 1(c) when choosing supernodes identical to separators1 . The order of the tree nodes (u 111 , u 112 , u 11 , u 121 , u 122 , u 12 , . . . , u 0 ), partially represented on the right of the matrix, is here a postordering: nodes in any subtree are processed consecutively. For u ∈ T composed of α u variables, the diagonal block associated with u in the factors matrix is formed of two lower and upper triangular matrices L 11 and U 11 of order α u . The β u nonzero off-diagonal rows (resp. columns) below (resp. to the right of) L 11 (resp. U 11 ) can be compacted into a contiguous matrix L 21 of size β u × α u (resp. U 12 of size α u × β u ). For the example of u 1 in Figure 1, α u1 = 3 (corresponding to variables 7, 8, 9), and β u1 = 9 corresponding to the off-diagonal rows 19, . . . , 27. Starting from y ← b and from the bottom of the tree, at each node u, the active components of y can be gathered into two dense vectors y 1 of size α u and y 2 of size β u to perform the operations . Finding these reachable vertices can be done using the transitive reduction of G(L T ), which is a tree (the elimination tree) when L results from the factorization of a matrix with symmetric (or symmetrized) structure. Since we work with a tree T with possibly more than one variable to eliminate at each node, let us define V b as the set of nodes in T including at least one element of struct(b). The structure of L -1 b is obtained by following paths from the nodes of V b up to the root. The tree consisting of these paths is the pruned tree of b, and we denote it by T p (b). The number of operations ∆ from Equation (4) now depends on b:

× × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × ×× × × × × × × × × × × × × ×× × × × ×× × × × × × × × × × × × × ×××× f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f u0 u2 u211 u22 . . .
(5) ∆(b) = u∈Tp(b) δ u .
We now consider the multiple RHS case of Equation [START_REF] Amestoy | On computing inverse entries of a sparse matrix in an out-of-core environment[END_REF], where RHS columns have different structures and we denote by B i the columns of B, for 1 ≤ i ≤ m. Rather than solving m systems each with a different pruned tree T p (B i ), we favor matrix-matrix computations by considering V B = 1≤i≤m V Bi , the union of all nodes in T with at least one nonzero from matrix B, and the pruned tree T p (B) = 1≤i≤m T p (B i ) containing all nodes in T reachable from nodes in V B . The triangular and update operations (3) become

Y 1 ← L -1 11 Y 1 and Y 2 ← Y 2 -L 21 Y 1 , leading to: (6) ∆(B) = m × u∈Tp(B) δ u .
An example is given in Figure 2, where B = [{B Illustration of multiple RHS and tree pruning. × corresponds to an initial nonzero in B and f to "fill-in" that appears in L -1 B, represented in terms of original variables and tree nodes. Gray parts of L -1 B (resp. of L -1 B i ) are the ones involving computations when RHS columns are processed in one shot (resp. one by one). tree pruning, or vertical sparsity, but perform extra operations by considering globally T p (B) instead of each individual pruned tree T p (B i ). Processing B by smaller blocks of columns would further reduce the number of operations at the cost of more traversals of the tree and a smaller arithmetic intensity, with a minimal number of operations ∆ min (B) = i=1,m ∆(B i ) reached when B is processed column by column, as in Figure 2(a)(right). We note that performing this minimal number of operations while traversing the tree only once (and thus accessing the L factor only once) may require performing complex and costly data manipulations at each node u with copies and indirections to work only on the nonzero entries of L -1 B at u. We now present a simpler approach which exploits the notion of intervals of columns at each node u ∈ T p (B). This approach to exploit what we call horizontal sparsity in B was introduced in another context [START_REF] Amestoy | Parallel computation of entries of A -1[END_REF].

× × × × × f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f u 0 u 2 u 211 u 22 . . . u 1 u 111 u 12 . . . × f f f × f f × f f f × f f f × f f f u 0 u 2 u
× f f f 2 × f f 3 × f f f 4 × f f f 5 × f f f (a) Structures of L -1 B, L -1 B 1 , . . . , L -1 B 5 .
Given a matrix B, we associate to a node u ∈ T p (B) its set of active columns ( 7)

Z u = {j ∈ {1, . . . , m} | u ∈ T p (B j )} .
The interval min(Z u ), max(Z u ) includes all active columns, and its length is

θ(Z u ) = max(Z u ) -min(Z u ) + 1.
Z u is sometimes defined for an ordered or partially ordered subset R of the columns of B, in which case we use the notation

Z u | R and θ(Z u | R ).
For u in T p (B), Z u is non-empty and θ(Z u ) is different from 0. As illustrated in Figure 3, the idea is then to perform the operations (3) on the θ(Z u ) contiguous columns min(Z u ), max(Z u ) instead of the m columns of B, leading to

(8) ∆(B) = u∈Tp(B) δ u × θ(Z u ). × f f f × f f × f f f × f f f × f f f u 0 u 2 u

Permuting RHS columns.

It is clear from Figure 3 that when exploiting horizontal sparsity thanks to column intervals, the number of operations to solve (2) depends on the order of the columns in B. We express the corresponding minimization problem as: [START_REF] George | Nested dissection of a regular finite-element mesh[END_REF] Find a permutation σ of {1, . . . , m} that minimizes ∆(B, σ) andθ(σ(Z u )) is the length of the permuted interval min(σ(Z u )), max(σ(Z u )) .

= u∈ Tp(B) δ u × θ(σ(Z u )), where σ(Z u ) = {σ(i) | i ∈ Z u } ,
Rather than trying to solve the global problem with linear optimization techniques, we will propose a cheap heuristic based on the tree structure. We first define the notion of node optimality. Definition 1. Given a node u in T p (B), and a permutation σ of {1, . . . , m}, we say that we have node optimality at u, or that σ is u-optimal, if and only if θ(σ(Z u )) = #Z u , where #Z u is the cardinality of Z u . Said differently, σ(Z u ) is a set of contiguous elements. θ(σ(Z u ))-#Z u , the number of columns (or padded zeros) on which extra computation is performed, is 0 if σ is u-optimal. Using the RHS of Figure 3(a), the identity permutation is u 0 -optimal because #Z u0 = #{1, 2, 3, 4, 5} = 5 = θ(Z u0 ), but not u 1 -and u 2 -optimal, because θ(Z u1 ) -#Z u1 = 2 and θ(Z u2 )-#Z u2 = 1, respectively. Our aim is thus to find a permutation σ that reduces the difference θ(σ(Z u )) -#Z u . After describing a permutation based on a postordering of T p (B), we present our new heuristic, which targets node optimality in priority at the nodes near the top of the tree.

(B j ) ∈ V Bj (see below). B is said to be postordered if and only if: ∀j 1 , j 2 , 1 ≤ j 1 < j 2 ≤ m, we have either u(B j1 ) = u(B j2 ), or u(B j1 ) appears before u(B j2
) in the postordering. In other words, the order of the columns B j is compatible with the order of their representative nodes u(B j ).

The postordering has been applied [START_REF] Amestoy | On computing inverse entries of a sparse matrix in an out-of-core environment[END_REF][START_REF] Tz | Parallel triangular solution in the out-of-core multifrontal approach for solving large sparse linear systems[END_REF][START_REF] Yamazaki | On partitioning and reordering problems in a hierarchically parallel hybrid linear solver[END_REF] to build regular chunks of RHS columns with "nearby" pruned trees, thereby limiting the accesses to the factors or the amount of computation. It was also experimented together with node intervals [START_REF] Amestoy | Parallel computation of entries of A -1[END_REF] to RHS with a single nonzero per column. In 3). On the other hand, the postorder permutation, σ PO , reorders the columns of B so that the order of their representative nodes u 112 , u 12 , u 211 , u 212 , u 221 is compatible with the postordering. In this case, there are no gray empty cells (see Figure 4(a)) and ∆(B, σ PO ) = ∆ min (B). More generally, it can be shown that the postordering induces no extra computations for RHS with a single nonzero per column [START_REF] Amestoy | Parallel computation of entries of A -1[END_REF]. For applications with multiple nonzeros per RHS, each column B j may correspond to a set V Bj with more than one node, in which case we choose the node that appears first in the sequence of postordered nodes of the tree [START_REF] Amestoy | On exploiting sparsity of multiple right-hand sides in sparse direct solvers[END_REF]. In Figure 4(b), some columns of B have several nonzeros. For example,

V B1 = {u 111 , u 211 , u 0 } is × f f f × f f × f f f × f f f × f f f u 0 u 2 u 22 u 222 u 221 u 21 u 212 u 211 u 1 u 12 u 122 u 121 u 11 u 112 u 111
→ σPO × f f f × f f × f f f × f f f × f f f 5 2 4 1 3
(a) RHS with one nonzero per column.

× f f × f f × × f f f × × f f × f f f × f f × f f × × f f f 19-27 16-18
→ σPO × f f × f f × × f f f × × f f × f f f × f f × f f × × f f f 1 4 2 5 6 3 → σFT × f f f × f f f × f f × f f × × f f × f f × × f f f × × f f 4 2 1 5 6 3
(b) RHS with multiple nonzeros per column.

(σ PO (Z u211 )) -#Z u211 = 5 -2 = 3.
The quality of σ PO depends on the tree postordering. If u 111 and u 112 were exchanged in the original postordering, B 1 and B 4 would be swapped, further reducing ∆. One drawback of the postorder permutation is that, since the position of a column is based a single representative node, some information is unused. In Figure 4(b), we represented another permutation σ FT that yields ∆(B, σ FT ) = 1140. This more general and powerful heuristic is presented in the next subsection.

The flat tree permutation.

With the aim of satisfying node optimality (see Definition 1), we present another algorithm to compute the permutation σ by first illustrating its geometric properties and then extending it to rely only on algebraic properties.

Geometrical intuition.

For u ∈ T , the domain associated with u is defined by the subtree rooted at u and is denoted by T [u]. The set of variables in T [u] corresponds to a subdomain created during the nested dissection algorithm. As an example, the initial 2D domain in Figure 5(a) (left) is T [u 0 ] and its subdomains created by dividing it with u 0 are T [u 1 ] and T [u 2 ]. In the following, T [u] will equally refer to a subdomain or a subtree. We make the strong assumption here that the nonzeros in an RHS column correspond to geometrically contiguous nodes in a regular domain on which a perfect nested dissection has been performed. The flat tree algorithm relies on the evaluation of the position of each RHS column compared to separators. The name flat tree comes from the fact that, given a parent node with two child subtrees, the algorithm orders first RHS columns included in the left subtree, then RHS columns associated to the parent (because they intersect both subtrees), and finally, as in an inorder, RHS columns included in the right subtree. The RHS columns in Figure 5(a) may be identified by three different types noted a, b and c according to their positions according to the root separator u 0 . An RHS column is of type a (resp. c) when its nonzero structure is included in

T [u11] T [u12] u1 T [u21] T [u22] u2 u0 × f f × × × f × f f × f × f × × × × × × × × × f × f × × f f × × × f × f f (b) Flat tree step 2.
T [u 1 ] (resp. T [u 2 ]),
and b when it is divided by u 0 . First, we group the RHS according to their type (a, b, or c) with respect to u 0 which leads to the creation of submatrices/subsets of RHS columns noted a, b and c. Second, we make sure to place b between a and c. Since all RHS in a and b have at least one nonzero in T [u 1 ], u 1 belongs to the pruned tree of all of them, hence the dense area filled in gray in the RHS structure. The same is true for b and c and u 2 . By permuting B as [a,b,c] ([c,b,a] is also possible), a and b, and b and c, are contiguous. Thus, θ(Z u1 ) = #Z u1 , θ(Z u2 ) = #Z u2 and we have minimized [START_REF] Eisenstat | The theory of elimination trees for sparse unsymmetric matrices[END_REF] for u 1 and u 2 thanks to u 1 -and u 2 -optimality. The algorithm proceeds recursively on each submatrix created to obtain local node optimality. First, d, e, f (resp. j, k, l) form subsets of the RHS of a (resp. c) based on their position/type with respect to u 1 (resp. u 2 ), see Figure 5(b). Second, thanks to the perfect nested dissection assumption, u 1 and u 2 can be combined to form a single separator that subdivides the RHS of b into three subsets g, h and i. During this second step, B is permuted as [d,e,f ,g,h,i,j,k,l]. The algorithm stops when the tree is fully processed or the RHS sets contain a single RHS.

Algebraic approach.

Let us consider the columns of B as an initially unordered set of RHS columns that we denote R B = {B 1 , B 2 , . . . , B m }. R ⊂ R B is a subset of the columns of B and r ∈ R is a generic element of R (one of the columns B j ). A permuted submatrix of B can be expressed as an ordered sequence of RHS columns. For two subsets of columns R and R , [R, R ] denotes a sequence of RHS columns in which the RHS from the subset R are ordered before those from R , without the order within R and R to be necessarily defined. We found this framework of RHS sets and subsets better adapted to formalize our algebraic algorithm than matrix notation with complex index permutations. We now characterize the geometrical position of a RHS using the notion of pruned layer: for a given depth d in the tree, and for a given RHS r ∈ R B , we define the pruned layer L d (r) as the set of nodes at depth d in the pruned tree T p (r). In the example of Figure 5(a), L 1 (r) = {u 1 } for all r ∈ a, L 1 (r) = {u 2 } for all r ∈ c, and L 1 (r) = {u 1 , u 2 } for all r ∈ b. The notion of pruned layer formally identifies sets of RHS with common characteristics in the tree, without geometric information. This is formalized and generalized by Definition 3. Definition 3. Let R ⊂ R B be a set of RHS, and let U be a set of nodes at depth d of the tree

T . We define R[U ] = {r ∈ R | L d (r) = U }
as the subset of RHS with pruned layer U .

We have for example, see Figure 5:

R[{u 1 }] = a, R[{u 2 }] = c and R[{u 1 , u 2 }] = b at depth d = 1.
The algebraic recursive algorithm is depicted in Algorithm 1. Its arguments are R, a set of RHS and d, the current depth. Initially, d = 0 and R = R B = R[u 0 ]. At each recursion step, the algorithm builds the distinct pruned layers U i = L d+1 (r) for the RHS r in R. Then, instead of looking for a permutation σ to minimize u∈Tp(R B ) δ u ×θ(σ(Z u )), it orders the R[U i ] by considering the restriction of (9) to R and to nodes at depth d + 1 of T p (R). Furthermore, with the assumption that T is balanced, all nodes at a given level of T p (R) are of comparable size. δ u may thus be assumed constant per level and needs not be taken into account. The algorithm is a greedy topdown algorithm, where at each step a local optimization problem is solved. This way, priority is given to the top of the tree, which is in general more critical because factor matrices are larger.

Algorithm 1 Flat Tree

procedure Flattree(R, d) 1) Build the set of children C(R) 1.1) Identify the distinct pruned layers (pruned layer = set of nodes)

U ← ∅ for all r ∈ R do U ← U ∪ {L d+1 (r)} end for 1.2) C(R) = {R[U ] | U ∈ U} 2) Order children C(R) as [R[U 1 ], . . . , R[U #C(R) ]]: return [Flattree(R[U 1 ], d + 1),. . .,Flattree(R[U #C(R) ], d + 1)] end procedure
The recursive structure of the algorithm can be represented by a recursion tree T rec defined as follows: each node R of T rec represents a set of RHS, C(R) denotes the set of children of R and the root is

R B . By construction of Algorithm 1, C(R) is a partition of R, i.e., R = ˙ R ∈C(R) R (disjoint union). Note that all r ∈ R such that L d+1 (r) = ∅ belong to R[∅],
which is also included in C(R). In this special case, R[∅] can be added at either extremity of the current sequence without introducing extra computation and the recursion stops for those RHS, as will be illustrated in Figures 6 and7(a). With this construction, each leaf of T rec contains RHS with indistinguishable nonzero structures, and keeping them contiguous in the final permutation avoids introducing extra computations. Assuming that for each R ∈ T rec the children C(R) are ordered, this induces an ordering of all the leaves of the tree, which defines the final RHS sequence. We now explain how the set of children C(R) is built and ordered at each step:

1) Building the set of children. The set of children of R ∈ T rec is built by first identifying the pruned layers U of all RHS r ∈ R. The different pruned layers are stored in U and we have for example (Figure 5, first step of the algorithm), U = {{u 1 }, {u 2 }, {u 1 , u 2 }}. We define

C(R) = {R[U ] | U ∈ U} (Definition 3
), which forms a partition of R. One important property is that all r ∈ R[U ] have the same nonzero structure at the corresponding layer so that numbering them contiguously prevents the introduction of extra computation.

2) Ordering the children.

The children sequence [R[U 1 ], . . . , R[U #C(R) ]
] at depth d + 1 should minimize the size of the intervals for nodes u at depth d + 1 of T p (R). The order inside each R[U i ] does not impact the size of these intervals (it will only impact lower levels). For any node

u at depth d + 1 in T p (R), we have θ(Z u | R ) = max(Z u | R ) -min(Z u | R ) + 1 = imax(u) i=imin(u) #R[U i ],
where Z u | R is the set of permuted indices representing the active columns restricted to R, and

i min (u) = min{i ∈ {1, . . . , #C(R)} | u ∈ U i } (resp. i max (u) = max{i ∈ {1, . . . , #C(R)} | u ∈ U i })
is the first (resp. last) index i such that u ∈ U i . Finally, we minimize the local cost function (sum of the interval sizes for each node at depth d + 1): [START_REF] Gilbert | Predicting structure in sparse matrix computations[END_REF] cost

([R[U 1 ], . . . , R[U #C(R) ]]) = u∈Tp(R) depth(u)=d+1 imax(u) i=imin(u) #R[U i ] To build the ordered sequence [R[U 1 ], . . . , R[U #C(R) ]],
we use a greedy algorithm that starts with an empty sequence, and at each step k ∈ {1, . . . , #C(R)} inserts a RHS set R[U ] picked randomly in C(R) at the position that minimizes [START_REF] Gilbert | Predicting structure in sparse matrix computations[END_REF] on the current sequence. To do so, we simply start from one extremity of the sequence of size k -1 and compute [START_REF] Gilbert | Predicting structure in sparse matrix computations[END_REF] for the new sequence of size k for each possible position 0 . . . k; if several positions lead to the same minimal cost, the first one encountered is chosen. In case u-optimality is obtained for each node u considered, then the permutation is said to be perfect and the cost function is minimal, locally inducing no extra operations on those nodes. Figure 6 shows the recursive structure of the RHS sequence after applying the algorithm on a binary tree. We refer to this representation as the layered sequence. For simplicity, the notation for pruned layers has been reduced from, e.g., {u 1 } to u 1 , and from

{u 1 , u 2 } to u 1 u 2 . From the recursion tree point of view, R[u 1 ], R[u 1 u 2 ], R[u 2 ] are the children of R[u 0 ] in T rec , R[u 11 ], R[u 11 u 12 ], R[u 12 ] the ones of R[u 1 ], etc.
Although we still use a binary tree, we make no assumption on the RHS structure, the domain, or the ordering in the example of Figure 7. The set of pruned layers corresponding to

R[u 0 ] is U = {u 1 u 2 , u 1 , ∅, u 2 }, so that C(R[u 0 ]) = {R[u 1 u 2 ], R[u 1 ], R[∅], R[u 1 u 2 ]}.
As can be seen in the non-permuted RHS structure, R[∅] = B 4 at depth 1 induces extra operations at nodes descendants of u 0 , which disappear when placing R[∅] at one extremity of the sequence. We choose to place it

R[u 0 ] d = 0 R[u 1 ] R[u 1 u 2 ] R[u 2 ] d = 1 R[u 11 ] R[u 11 u 12 ] R[u 12 ] * R[u 21 ] R[u 21 u 22 ] R[u 22 ] d = 2
× × × f f f f × × × × × × × f f × f × f f × f × f f u 0 u 2 u 22 u 21 u 1 u 12 u 11 1 2 3 4 5 6 7 σ FT f f × × × × × × f × f f × × × f f × f × f f × f f × 2 3 6 1 7 5 4 (a) RHS structure. u 0 u 2 u 22 u 21 u 1 u 12 u 11 (b) Separator tree T . R[u 0 ] B 1 , . . . , B 7 R[u 1 u 2 ] B 1 B 3 B 6 B 7 R[u 11 u 22 ] B 7 R[u 21 u 22 ] B 1 R[u 12 u 21 ] B 6 R[u 12 ] B 3 R[u 1 ] B 2 R[u 2 ] B 5 R[∅] B 4
(c) Recursion tree Trec. last and obtain the sequence

[R[u 1 ], R[u 1 u 2 ], R[u 2 ], R[∅]].
A recursive call is done on the identified sets, as illustrated in Figure 7 5(b)), may still occur with the flat tree algorithm. This will also be illustrated in Section 5, where ∆(B, σ FT ) is 39% larger than ∆ min (B), in the worst case. A blocking algorithm is now introduced to further reduce ∆(B, σ FT ).

(c). Since R[u 1 ], R[u 2 ] and

4.

Toward a minimal number of operations using blocks.

Geometrical intuition.

In terms of operation count, optimality (∆ min (B)) is obtained when processing the columns of B one by one, using m blocks. However, this requires processing the tree m times and will typically lead to a poor arithmetic intensity (and likely a poor performance). On the other hand, the algorithms of Section 3 only use one block, which provides a higher arithmetic intensity but requires extra operations. In this section, our objective is to create a minimal number of (possibly large) blocks while reducing the number of extra operations by a given amount.

On the one hand, two RHS or sets of RHS included in two different domains exhibit interesting properties, as can be observed for sets a ∈ T [u 1 ] and c ∈ T [u 2 ] from Figure 5(a). No extra operations are introduced between them: ∆([a, c]) = ∆(a) + ∆(c). We say that a and c are independent sets and they can be associated together. On the other hand, a set of RHS intersecting a separator (such as set b) has zeros and nonzeros in rows common to their adjacent RHS sets (a and c) which will likely introduce extra computation. In Figure 5

(b), we have for example ∆([a, b]) = ∆([d, e, f , g, h, i]) > ∆([d, e, f ]) + ∆([g, h, i]) = ∆(a) + ∆(b) and ∆([a, b, c]) > ∆(a) + ∆(b) + ∆(c).

We say that b is a set of problematic RHS. Figure 4 For this, we formalize the notion of independence, which will be the basis for our blocking algorithm.

Algebraic formalization.

In this section, we assume the matrix B to be flat tree ordered and the recursion tree T rec to be built and ordered. Using the notations of Definition 3, we give an algebraic definition of the independence property: Definition 4. Let U 1 , U 2 be two sets of nodes at a given depth of a tree T , and let

R[U 1 ], R[U 2 ] be the corresponding sets of RHS. R[U 1 ], R[U 2 ] are said to be independent if and only if U 1 ∩U 2 = ∅.
With Definition 4, we are able to formally identify independent sets that can be associated together. BuildGroups(g d i , d + 1) solves the aforementioned coloring problem on these RHS sets and builds the k groups (g d+1 j+1 , . . . , g d+1 j+k ). In Figure 8(right), there is initially a single group g 0 1 = [R[u 0 ]] with one set of RHS, which may be expressed as the ordered sequence

Take for example a

= R[u 1 ] and c = R[u 2 ] (Figure 5(a)), R[u 1 ] and R[u 2 ] are independent and ∆([R[u 1 ], R[u 2 ]]) = ∆(R[u 1 ]) + ∆(R[u 2 ]). On the contrary, when R[U 1 ], . . . ,
R[u0] g 0 1 d = 0 R[u1] R[u2] R[u1u2] g 1 1 g 1 2 d = 1 R[u11] R[u12] R[u21] R[u22] g 2 1 R[u11u12] R[u21u22] g 2 2 * 1 * 2 g 2 3 * 3 g 2 4 d = 2
[R[u 1 ]R[u 1 u 2 ]R[u 2 ]]. g 0 1 does not satisfy the independence property at depth 1 because u 1 ∩ u 1 u 2 = ∅ or u 2 ∩ u 1 u 2 = ∅. BuildGroups(g 0 1 , 1) yields g 1 1 = [R[u 1 ], R[u 2 ]] and g 1 2 = [R[u 1 u 2 ]
]. The algorithm proceeds until a maximal depth

d max : (g 2 1 , g 2 2 ) = BuildGroups(g 1 1 , 2), (g 2 3 , g 2 4 ) = BuildGroups(g 1 2 ,
2), etc. To illustrate the interest of property (ii), let us take sets

d = R[u 11 ], f = R[u 12 ], j = R[u 21 ] and l = R[u 22 ] from Figure 5(b). One can see that ∆([d, f , j, l]) = ∆(d) + ∆(f ) + ∆(j) + ∆(l) < ∆([d, j, f , l]). Compared to [d, f , j, l]
which respects the global flat tree ordering and ensures u 1 -and u 2 -optimality, [d, j, f , l] does not and would thus increase θ(Z u1 ) and θ(Z u2 ). Furthermore, Algorithm 2 ensures the property (whose proof can be found in [START_REF] Amestoy | On exploiting sparsity of multiple right-hand sides in sparse direct solvers[END_REF]), that the independent sets of RHS grouped together do not introduce extra operations:

Property 5. For any group

g d = [R[U 1 ], . . . , R[U n ]] created through Algorithm 2 at depth d, we have ∆([R[U 1 ], . . . , R[U n ]]) = n i=1 ∆(R[U i ]).
Interestingly, Property 5 can be used to prove, in the case of a single nonzero per RHS, the optimality of the flat tree permutation. Corollary 6. Let R B be a set of RHS such that ∀r ∈ R B , #V r = 1. Then the flat tree permutation is optimal:

∆(R B ) = ∆ min (R B ).
Proof. Since ∀r ∈ R B , #V r = 1, T p (r) is a branch of T . As a consequence, any set of RHS R[U ] built through the flat tree algorithm is represented by a pruned layer U containing a single node u. At each step of the flat tree algorithm (Algorithm 1), the RHS sets identified are thus all independent from each other. When applying Algorithm 2, a unique group R B is then kept until the bottom of the tree. Blocking is thus not needed and Property 5 applies at each level of the recursion. ∆(R B ) is thus equal to the sum of the ∆(R[U ]) for all leaves R[U ] of the recursion tree

T rec . Since ∆(R[U ]) = ∆ min (R[U ]) on those leaves (all RHS in R[U ] involve the exact same nodes and operations), we conclude that ∆(R B ) = ∆ min (R B ).
This proof is independent of the ordering of the children at step 2 of Algorithm 1. Corollary 6 is thus more general: any top-down recursive ordering keeping together RHS with identical pruned layer at each step is optimal, as long as the pruned layers identified at each step are independent. Back to the BuildGroups function, the solution of the coloring problem may not be unique. Even on the simple example of Figure 8, there are several ways to define groups, as shown in Figure 9 for g 1 1 : both strategies satisfy the independence property and minimize the number of groups. The CritPathBuildGroups strategy tends to create a large group g 2 1 and a smaller one,

R[u 1 ] R[u 2 ] g 1 1 d = 1 R[u11] R[u12] R[u21] R[u22] g 2 1 R[u11u12] R[u21u22] g 2 2 d = 2 R[u 1 ] R[u 2 ] g 1 1 d = 1 R[u11] R[u12] R[u21u22] g 2 1 R[u21] R[u22] R[u11u12] g 2 2 d = 2
g 2 2 .
In each group the computations on the tree nodes are expected to be well balanced because all branches of the tree rooted at u 0 might be covered by the RHS (assuming thus a reasonably balanced RHS distribution over the tree). The choice of CritPathBuildGroups can be driven by tree parallelism considerations, namely, the limitation of the sum of the operation counts on the critical paths of all groups. The RegBuildGroups strategy tends to balance the sizes of the groups but may create more unbalance regarding the distribution of work over the tree.

We note that for a given depth, applying BuildGroups on all groups may not always be necessary, and that for a given group, enforcing the independence property may create more than two groups. In the next section we minimize the number of groups created using greedy heuristics.

4.3.

A greedy approach to minimize the number of groups. Compared to Algorithm 2, Algorithm 3 adds the group selection, limits the number of groups created from a given group to two, and stops depending on a given tolerance on the amount of operations.

First, instead of stepping into each group as in Algorithm 2, we select among the current groups the one responsible for most extra computation, that is, the one maximizing ∆(g) -∆ min (g). This implies that groups that are candidate for splitting might have been created at different depths and we use a superscript to indicate the depth d at which a group was split, as in the notation g d 0 . Second, instead of a coloring problem which creates as many groups as colors obtained, we look (procedure BuildMaxIndepSet) inside the RHS sets of g d 0 for a maximal group of independent sets at depth d+1, denoted g d+1 imax . The other sets are left in another group g d c , whose depth remains equal to d. g d c may thus consist of dependent sets that may be subdivided later if needed 2 . Rather than an exact algorithm to determine g d+1 imax , we use a greedy heuristic. Finally, we define µ 0 as the tolerance of extra operations authorized. With a typical value µ 0 = 1.01, the algorithm stops when the number of extra operations is within 1% of the minimal number of operations ∆ min , returning G as the final set of groups.

Algorithm 3 Blocking algorithm G ← {R B }, ∆ min ← ∆ min (R B ), ∆ ← ∆(R B ) while ∆/∆ min > µ 0 do Select g d 0 such that ∆(g d 0 ) -∆ min (g d 0 ) = max g∈G (∆(g) -∆ min (g)) Group selection (g d+1 imax , g d c ) ← BuildMaxIndepSet(g d 0 , d+1) G ← G ∪ {g d+1 imax , g d c } \ {g d 0 } ∆ ← ∆ -∆(g d 0 ) + ∆(g d+1 imax ) + ∆(g d c
) end while 5. Experimental results. In this section, we report on the impact of the proposed permutation and blocking algorithms on the forward substitution (Equation (2)), using a set of 3D regular finite difference problems coming from seismic and electromagnetism modeling [START_REF] Amestoy | Fast 3D frequency-domain full waveform inversion with a parallel Block Low-Rank multifrontal direct solver: application to OBC data from the North Sea[END_REF][START_REF] Shantsev | Large-scale 3D EM modeling with a Block Low-Rank multifrontal direct solver[END_REF], for which the solve phase is costly. The characteristics of the corresponding matrices and RHS are presented in Table 1. In both applications, the nonzeros of each RHS correspond to a small set of close points, near the top of the 3D grid corresponding to the physical domain, with some overlap between RHS. Except in Section 5.3, a geometric nested dissection (ND) algorithm is used to reorder the matrix.

Impact of the flat tree algorithm.

We first introduce the terminology used to denote the different strategies developed in this study and that impact the number of operations ∆. DEN represents the dense case, where no optimization is used to reduce ∆, and TP means tree pruning. When column intervals are exploited at each tree node, we denote by RAN, INI, PO and FT the random, initial (σ = id), postorder (σ PO ) and flat tree (σ FT ) permutations, respectively.

The improvements brought by the different strategies are presented in Table 2. Compared to the dense case, TP divides ∆ by at least a factor 2. When column intervals are exploited at each node, the large gap between RAN and INI shows that the original column order holds geometrical properties. FT behaves better than INI and PO and gets reasonably close to ∆ min . Overall, FT 2 In case g d c consists of independent sets and is selected, g d+1 c = g d c will only be subdivided at depth d + 2. provides a 13% gain on average over PO. However, the gain on ∆ decreases from 25% on the 10Hz problem to 1% on the H116 problem. This can be explained by the fact that B is denser for the seismic applications than for the electromagnetism applications (see Table 1). Indeed, the sparser B, the closer we are to a single nonzero per RHS in which case both FT and PO are optimal. Second, we evaluate the impact of exploiting RHS sparsity on tree parallelism. Table 3 gives the maximal theoretical speed-up S that can be reached using tree parallelism only (node parallelism is also needed, for example on the root). It is defined as S = ∆/∆ cp where ∆ cp is the number of operations on the critical path of the tree. We observe that, when sparsity is exploited, tree parallelism is significantly smaller than in the dense case. This is because the depth of the pruned tree T p (B) is similar to that of the original tree (some nonzeros of B generally appear in the leaves), while the tree effectively processed is pruned and thus the overall amount of operations is reduced. For the same reason, S is smaller for test cases where D(B) is small. For the 5Hz, 7Hz, and 10Hz problems which have more nonzeros per column of B, besides decreasing the operation count more than the other strategies, FT exhibits equivalent or even better tree parallelism than PO. For such matrices, where D(B) is large, FT balances the work on the tree better than PO and reduces the work on the critical path more than the total work. Overall, FT reduces the operation count better than any other strategy and has good parallel properties.

Impact of the blocking algorithm.

First, we show that the blocking algorithm decreases the operation count ∆ while creating a small number of groups. Second, we discuss parallel properties of the clustering strategies illustrated in Figure 9. In Table 4, we report the value of ∆ N G /∆ min as a function of the number of groups created. x means that the blocking algorithm stopped because the condition ∆ N G /∆ min ≤ µ 0 was reached, with µ 0 for Algorithm 3 set to 1.01. Computing from Table 4 the ratio of extra operations reduction 1 -∆ N G -∆min ∆1-∆min for N G groups created, we observe an average reduction of 74% of the extra operations when N G = 2, i.e., when only two groups are created. Table 4 also shows that ∆ N G reaches a value close to ∆ min very quickly. In Table 5, we report the sum of operation counts on the critical paths ∆ cp over all groups created using CritPathBuildGroup and RegBuildGroup strategies, when the number of groups created is three, leading to a value ∆ close to ∆ min , see column "NG=3" of Table 4. In this case, the total number of operations ∆ during the forward solution phase on all groups is equal whether we use CritPathBuildGroup or RegBuildGroup. Tree parallelism is thus a crucial discriminant between both strategies, and we indeed observe in Table 5 that CritPathBuildGroup effectively limits the length of critical paths over the three groups created, justifying its use.

Experiments with other orderings.

As mentioned earlier, several orderings may be used to order the unknowns of the original matrix, thanks to the algebraic nature of our flat tree and blocking algorithms. Although local ordering methods (AMD, AMF as provided by the MUMPS package3 ) are known not to be competitive with respect to algebraic nested dissectionbased approaches such as SCOTCH 4 or METIS 5 on large 3D problems, we include them in Table 6 in order to study how the flat tree and blocking algorithms behave in general situations.

First, an important aspect of using other orderings is that they often produce much more irregular trees, leading to a large number of pruned layers to sequence. The FT permutation reduces the operation count significantly with SCOTCH and METIS, for which we observe an average 31% and 26% reduction compared to the PO permutation. Gains are also obtained with AMD for most test cases. However, FT does not perform well with AMF. This can be explained by the fact that AMF produces too irregular trees which do not fit well with our design of the FT strategy.

Second, we evaluate the blocking algorithm (BLK) and compare it with a regular blocking algorithm (REG) based on the PO permutation, that divides the initial set of columns into regular chunks of columns. Table 6 shows that the number of groups required to reach ∆/∆ min ≤ 1.01 is much smaller for BLK than for REG in all cases. Our blocking algorithm is very efficient with most orderings except AMF, where the number of groups created is high (but still lower than REG). 5.4. Performance. We analyse in this section the time for the forward substitution and for the flat tree and blocking algorithms on a single Intel Xeon core @2.2GHz. A performance analysis in multithreaded or distributed environment for the largest problems is out of the scope of this study. In Table 7, we report the time of the forward substitution of the MUMPS solver6 and the percentage of time spent in BLAS operations, excluding the time for data manipulation and copies. Timings with the regular blocking REG to reach our target number of operations (∆ N G /∆ min ≤ 1.01), at the cost of a larger number of groups, N G 1, are also indicated. Table 7 shows that the time reduction is in agreement with the operations reduction reported in Table 2. One can also notice that when exploiting sparsity, the proportion of time spent in BLAS operations increases. This is due to the fact that the relative weight of the top of the tree (with larger fronts for which time is dominated by BLAS operations) is larger when sparsity is exploited. When targeting a reduction of the number of operations with a regular blocking (REG), the large number of blocks (last column of Table 7) makes the granularity of the BLAS operations too small to reach the performance of the BLK strategy. On the other hand, Table 8 shows that larger regular blocks improve the Gflop rate but are not competitive with respect to BLK because of the increase in the number of operations. We also observe that the best block size for REG is problem-dependent. Table 9 relates the execution time of the flat tree and blocking algorithms. The execution times are reasonable compared to the corresponding execution time of the forward substitution. Moreover, we observe that the time of the flat tree algorithm only slowly increases with the problem size. The time of the blocking algorithm, due to its limited number of iterations, is not critical.

6. Guided Nested Dissection. Given an ordering and a tree, one may think of moving the unknowns corresponding to RHS nonzeros to supernodes higher in the tree with on the one hand, a smaller pruned tree, but on the other hand, an increase in the factor size due to larger supernodes. Better, one may guide the ordering to include as many nonzeros of B as possible within separators during the top-down nested dissection and prune larger subtrees. This will however involve a significant extra cost for applications where each RHS contains several contiguous nodes in the grid, e.g., form a small parallelepiped. For such applications, the geometry of the RHS nonzeros could however be exploited. A first idea avoids problematic RHS by choosing separators that do not intersect RHS nonzeros. Although this idea could be tested by adding edges between RHS nonzeros before applying SCOTCH or METIS, this does not appear to be so useful in our applications, where we observed significant overlap between successive RHS. Another idea, when all RHS are localized in a specific area of the domain, is to shift the separators from the nested dissection to insulate the RHS in a small part of the domain. Such a modification of the ordering yields an unbalanced tree in which the RHS nonzeros appear at the smaller side of the tree, improving the efficiency of tree pruning and resulting in a reduction of ∆ min , and thus ∆. This so called guided nested dissection was implemented and tested on the set of test cases shown in Table 10, where we observe that the number of operations ∆ min is decreased, as expected. Since the factor size has also increased significantly, a trade-off may be needed to avoid increasing too much the cost of the factorization. 7. Applications and related problems. We describe applications where our contributions can be applied. When only part of the solution is needed, one can show that the approaches described in Section 2 can be applied to the backward substitution (U X = Y ), which involves similar mechanisms as the forward substitution [START_REF] Tz | Parallel triangular solution in the out-of-core multifrontal approach for solving large sparse linear systems[END_REF][START_REF] Yeung | Interactively cutting and constraining vertices in meshes using augmented matrices[END_REF]. The backward substitution traverses the tree nodes from top to bottom so that the interval mechanism is reversed, i.e., the interval from a parent includes the intervals from its children and the properties of local optimality are preserved. If the structure of the partial solution requested differs from the RHS structure, another call to the flat tree algorithm must then be performed to optimize the number of operations. Exploiting sparsity also in the backward step can for instance be useful in some augmented approaches [START_REF] Yeung | AMPS: an augmented matrix formulation for principal submatrix updates with application to power grids[END_REF] to deal with small matrix updates without complete refactoring, and in some 3D EM geophysics applications [START_REF] Shantsev | Large-scale 3D EM modeling with a Block Low-Rank multifrontal direct solver[END_REF]. Another application of this work is the computation of Schur complements, where instead of truncating a factorization of the whole system ( A C B D ), one exploits the factorization of A to use triangular solves with sparse RHS. Taking the symmetric case where C = B T , the Schur complement S can be written S = D -BA -1 B T = D -B(LL T ) -1 B T = D -(L -1 B T ) T (L -1 B T ), as in the PDSLin solver [START_REF] Yamazaki | On partitioning and reordering problems in a hierarchically parallel hybrid linear solver[END_REF]. Since B is sparse, B = L -1 B T can be computed thanks to the algorithms developed in this article before computing the sparse product B T B .

Finally, we comment on related problems and algorithms. We indicated that the blocking algorithm is closely related to graph algorithms like coloring and maximum independent set. Concerning the minimization problem (9) which we addressed with the flat tree algorithm, it can also be regarded globally: using the structure of L -1 B, the objective is to find a permutation of the columns that minimizes the sum of the intervals weighted with δ u . This interval minimization problem is similar to a sparse matrix profile reduction problem [START_REF] Berry | Sparse marix reordering schemes for browsing hypertext[END_REF][START_REF] Reid | Reducing the total bandwidth of a sparse unsymmetric matrix[END_REF] (and we thus suspect it to be NP-complete). As mentioned in the introduction, hypergraph models have been used in the context of blocking algorithms, with different constraints and objectives compared to ours [START_REF] Amestoy | On computing inverse entries of a sparse matrix in an out-of-core environment[END_REF][START_REF] Yamazaki | On partitioning and reordering problems in a hierarchically parallel hybrid linear solver[END_REF]. Modeling L -1 B as an hypergraph might lead to other heuristics than the flat tree algorithm using some variants of hypergraph partitioning, although dense parts in L -1 B might need special treatment. One advantage of our permutation and blocking algorithms is that, instead of tackling the problem globally, they decompose the problem into easier subproblems with low complexity by making use of the separator tree T , thereby exploiting the fact that L -1 B has a very special structure closely related to the tree. In the context of general unsymmetric matrices, the structure of the solution of the forward step is given by the set of reachable vertices in the elimination dag of L T [START_REF] Gilbert | Elimination structures for unsymmetric sparse LU factors[END_REF]. To make the elimination dag a tree to be used in our context, one could consider adding a limited number of entries in L. Similarly to the case of matrices with a symmetric or quasi-symmetric pattern for which the elimination tree of the matrix A + A T is used, one could add entries in L having a symmetric counterpart in U . Another possibility is to use the work presented in [START_REF] Eisenstat | The theory of elimination trees for sparse unsymmetric matrices[END_REF] that extends the notion of elimination tree to unsymmetric matrices by considering paths in the factor matrices to characterize the elimination tree. How useful this generalisation of the elimination tree can be in our context would deserve to be further studied.

Conclusion.

We introduced permutation and blocking algorithms to improve the tree pruning [START_REF] Tz | Parallel triangular solution in the out-of-core multifrontal approach for solving large sparse linear systems[END_REF] and the node interval [START_REF] Amestoy | Parallel computation of entries of A -1[END_REF] algorithms. A first contribution is the "flat tree" algorithm which permutes RHS to reduce the cost of the forward substitution. A second contribution is a blocking algorithm that further decreases this cost by adequately choosing groups of RHS that can be processed together. Although both algorithms are based on geometrical observations, they are designed with an algebraic approach, giving a general scope to this work. Notions of node optimality and RHS independence were introduced and formalized, together with theoretical properties to provide insight and to support the proposed algorithms. Experimental results on real test cases showed the effectiveness of both the flat tree and the blocking algorithms. Compared to a postorder-based permutation, the flat tree permutation showed an average (resp. maximum) gain of 13% (resp. 25%) on the total operation count with a nested dissection ordering, and interesting parallel properties. Moreover, results with the blocking algorithm validated our approach since only a handful of groups is created compared to several hundreds when using a regular blocking technique. Finally, sequential performance results confirmed the good potential of the proposed approaches.

3 × 3 ×

 33 3 regular mesh.

  Separator tree T .

Figure 1 .

 1 Figure 1. (a) A 3D mesh with a 7-point stencil. Mesh nodes are numbered according to a nested dissection ordering. (b) Corresponding matrix with initial nonzeros (×) in the lower triangular part of a symmetric matrix A and fill-in (f ) in the L factor. (c) Separator tree, also showing the sets of variables to be eliminated at each node.

( 3 ) y 1 ← L -1 11 y 1 and y 2 ← y 2 -L 21 y 1 .y 1 2 .

 3122112 and y 2 can then be scattered back into y, so that y 2 will be used at higher levels of T . When the root is processed, y contains the solution of Ly = b. Denoting by δ u = α u × (α u -1 + 2β u ) the number of arithmetic operations for (3), the number of operations to solve Ly = b is Exploitation of sparsity in right-hand sides. Consider a non-singular n × n matrix A with a nonzero diagonal, and its directed graph G(A), with an edge from vertex i to vertexj if a ij = 0.Given a vector b, let us define struct(b) = {i, b i = 0} as its nonzero structure, and closure A (b) as the smallest subset of vertices of G(A) including struct(b) without incoming edges. Gilbert [10, Theorem 5.1] characterizes the structure of the solution of Ax = b by the relation struct(A -1 b) ⊆ closure A (b), with equality in case there is no numerical cancellation. In our context of triangular systems, ignoring such cancellation, struct(L -1 b) = closure L (b) is also the set of vertices reachable from struct(b) in G(L T ), where edges have been reversed [11, Theorem 2.1]

  Pruned tree Tp(B) = Tp(B 1 ) ∪ . . . ∪ Tp(B 5 ).

Figure 2 .

 2 Figure2. Illustration of multiple RHS and tree pruning. × corresponds to an initial nonzero in B and f to "fill-in" that appears in L -1 B, represented in terms of original variables and tree nodes. Gray parts of L -1 B (resp. of L -1 B i ) are the ones involving computations when RHS columns are processed in one shot (resp. one by one).

  Pruned tree Tp(B) with intervals.

Figure 3 .

 3 Figure 3. Column intervals for the RHS of Figure 2: in gray (a) and above/below each node (b). Taking for instance u 21 , there are nonzeros in columns 1 and 4, so that Zu 21 = {1, 4}. Instead of performing the solve operations on m = 5 columns at u 21 , computation is limited to the θ(Zu 21 ) = 4 columns of interval 1, 4 (and to a single column at, e.g., node u 221 ). Overall, ∆(B) is reduced from 1320 to 948 (while ∆ min (B) = 744).

3. 1 .Definition 2 .

 12 The postorder permutation. In Figure1, the sequence [u 111 , u 112 , u 11 , u 121 , u 122 , u 12 , u 1 , u 211 , u 212 , u 21 , u 221 , u 222 , u 22 , u 2 , u 0 ] used to order the matrix follows a postordering. Consider a postordering of the tree nodes u ∈ T , and a RHS matrix B = [B j ] j=1...m where each column B j is represented by one of its associated nodes u

  Figure 3(a), B has a single nonzero per column. The initial natural order of the columns (INI) induces computation on explicit zeros represented by gray empty cells and we had ∆(B) = ∆(B, σ INI ) = 948 and ∆ min (B) = 744 (see Figure

Figure 4 .

 4 Figure 4. Illustration of the postordering of two RHS with one or more nonzeros per column.

Figure 5 .

 5 Figure 5. Flat tree geometrical illustration. In (a) and (b), the figure on the left represents different types of RHS, and the one on the right the permuted RHS matrix. × or f in a rectangle indicate the presence of nonzeros in the corresponding submatrix, parts of the matrix filled in gray are fully dense and blank parts only contain zeros.

Figure 6 .

 6 Figure 6. A layered sequence built by the flat tree algorithm on a binary tree. Sets R[∅] (not represented) could be added at either extremity of the concerned sequence (e.g., right after R[u 2 ] for a RHS included in u 0 ). With the strong assumptions of Figure 5, * = R[u 11 u 21 ], R[u 11 u 12 u 21 u 22 ], R[u 12 u 22 ]. Otherwise, * is more complex.

Figure 7 .

 7 Figure 7. Illustration of the algebraic flat tree algorithm on a set of 7 RHS.

  (right) gives another example where extracting the problematic RHS B 1 and B 5 from [B 4 , B 2 , B 1 , B 5 , B 6 , B 3 ] suppresses all extra operations: ∆([B 4 , B 2 , B 6 , B 3 ]) + ∆([B 1 , B 5 ]) = ∆ min = 1056. Situations where no assumption on the RHS structure is made are more complicated and require a general approach.

Figure 8 .

 8 Figure 8. A first version of the blocking algorithm (left). It is illustrated (right) on the layered sequence of Figure 6. With the geometric assumptions of Figure 5, * 1 = R[u 11 u 21 ], * 2 = R[u 12 u 22 ], and * 3 = R[u 11 u 12 u 21 u 22 ]. in Figure 8 traverses T rec from top to bottom. At each depth d, each intermediate group g d i verifies the following properties: (i) g d i can be represented by a sequence [R[U 1 ], . . . , R[U n ]], and (ii) the sequence respects the flat tree order of T rec . Then, BuildGroups(g d i , d + 1) first builds the sets of RHS at depth d + 1, which are exactly the children of the R[U j ] ∈ g d i in T rec . Second,BuildGroups(g di , d + 1) solves the aforementioned coloring problem on these RHS sets and builds the k groups (g d+1 j+1 , . . . , g d+1 j+k ). In Figure8(right), there is initially a single group g 0 1 = [R[u 0 ]] with one set of RHS, which may be expressed as the ordered sequence[R[u 1 ]R[u 1 u 2 ]R[u 2 ]]. g 0 1 does not satisfy the independence property at depth 1 because u 1 ∩ u 1 u 2 = ∅ or u 2 ∩ u 1 u 2 = ∅. BuildGroups(g 0 1 , 1) yields g 1 1 = [R[u 1 ], R[u 2 ]] and g 1 2 = [R[u 1 u 2 ]]. The algorithm proceeds until a maximal depth d max : (g 2 1 , g 2 2 ) = BuildGroups(g 1 1 , 2), (g 2 3 , g 2 4 ) = BuildGroups(g 1 2 , 2), etc. To illustrate the interest of property

Figure 9 .

 9 Figure 9. Two strategies to build groups: CritPathBuildGroups (left) and RegBuildGroups (right).

  11,1 }, {B 6,2 }, {B 13,3 },{B 10,4 }, {B 2,5 }], V B = {u 212 , u 12 , u 221 , u 211 , u 112 }, δ(u

	0 ) = 72, δ(u 1 ) = δ(u 2 ) = 60, δ(u 11 ) = 12, δ(u 112 ) = 6, etc.,
	∆(B) = 5 × 264 = 1320, and ∆(B 1 ) + ∆(B 2 ) + . . . + ∆(B 5 ) = 744. At this point, we exploit
	1
	5
	10
	15
	20
	25

  R[U n ] are not pairwise independent, we group together independent sets of RHS, while forming as few groups as possible. This problem is equivalent to a graph coloring problem, where R[U 1 ], . . . , R[U n ] are the vertices and an edge exists between R[U i ] and R[U j ] if and only if U i ∩ U j = ∅. Several heuristics exist for this problem, and each color will correspond to one group. The blocking algorithm as depicted

	Algorithm 2 Blocking algorithm
	for d = 0 to d max do
	j ← 0 /* #groups at depth d + 1 */
	for all groups g d i at detph d do (g d+1 j+1 . . . g d+1 j+k ) ← BuildGroups(g d i , d+1)
	/* k new groups have been created */
	j ← j + k
	end for
	end for

Table 1

 1 Characteristics of the n × n matrix A and n × m matrix B for different test cases. D(A) = nnz(A)/n and D(B) = nnz(B)/m represent the average column densities for A and B, respectively.

	application matrix n(×10 6 ) D(A) sym	m D(B)
	seismic modeling	5Hz 7Hz 10Hz	2.9 7.2 17.2	24 25 26	no no no	2302 2302 2302	567 486 486
		H0	.3	13	yes	8000	9.8
		H3	2.9	13	yes	8000	7.5
	electro-magnetism modeling	H17 H116 S3 S21	17.4 116.2 3.3 20.6	13 13 13 13	yes yes yes 12340 8000 8000 yes 12340	6 6 19.7 9.5
		S84	84.1	13	yes 12340	8.6
		D30	29.7	23	yes	3914	7.6

Table 2

 2 Number of operations (×1013 ) during the forward substitution (LY = B) according to the strategy used (ND ordering).

	∆	DEN	TP	RAN	INI	PO	FT	∆min
	5Hz	1.73	.74	.74	.44	.36	.28	.22
	7Hz	5.94	2.54	2.52	1.46	1.21	.92	.69
	10Hz	20.62	9.01	8.92	4.78	3.85	2.87	2.26
	H0	.39	.11	.11	.086	.070	.057	.050
	H3	7.19	3.33	3.31	2.48	1.47	1.26	.95
	H17	81.34	37.15	36.97	27.52	10.41	10.21	10.12
	H116 990.02 448.31 445.91 327.89 123.79 121.76 120.68
	S3	13.36	4.98	4.91	3.73	2.65	2.17	1.71
	S21 156.20	49.04	48.07	35.42	25.73	22.53	19.43
	S84 983.48 286.57 282.70 222.59 161.87 138.56 118.51
	D30	71.60	39.78	39.38	19.49	10.93	10.21	7.31

Table 3

 3 Theoretical tree parallelism according to the strategy used (ND ordering).

	S	DEN	TP RAN INI	PO	FT
	5Hz 8.60 3.91 3.88 3.11 2.39 2.54
	7Hz 8.92 3.97 3.94 3.02 2.25 2.48
	10Hz 9.10 4.04 4.02 2.96 2.30	2.30
	H0 5.88 2.11 2.11 1.75 1.51	1.45
	H3 5.99 3.22 3.21 2.47 2.02	2.11
	H17 6.32 3.34 3.32 2.54 2.00	1.97
	H116 7.92 3.63 3.61 2.75 2.05	2.02
	S3 6.12 2.84 2.83 2.18 1.73	1.61
	S21 6.30 2.56 2.46 1.85 1.49	1.47
	S84 8.01 2.41 2.38 1.90 1.53	1.52
	D30 8.50 4.73 4.70 2.86 2.05	2.56

Table 4

 4 Impact of the number of groups NG on the normalized operation count, until ∆ N G /∆ min becomes smaller than the tolerance µ 0 = 1.01 (ND ordering).

	∆ N G /∆ min	FT NG= 2 NG= 3 NG= 4 NG= 5
	5Hz 1.283	1.111	1.001	x	x
	7Hz 1.321	1.116	1.002	x	x
	10Hz 1.269	1.029	1.002	x	x
	H0 1.148	1.029	1.010	1.002	x
	H3 1.329	1.068	1.027	1.005	x
	H17 1.009	x	x	x	x
	H116 1.009	x	x	x	x
	S3 1.275	1.120	1.045	1.012	1.003
	S21 1.160	1.037	1.015	1.003	x
	S84 1.169	1.041	1.015	1.002	x
	D30 1.397	1.082	1.058	1.024	1.004

Table 5

 5 Sum of critical paths' operations (×1013 ) for two grouping strategies when three groups are created.

	g ∆ cp (g)	5Hz 7Hz 10Hz	H0 H3
	CritPathBuildGroup .092 .30 1.00 .037 .50
	RegBuildGroup	.12	.43	1.58	.044 .72

Table 6

 6 Operation count ∆(×10 13 ) for permutation strategies PO and FT, and number of groups NG required to reach ∆ N G ∆ min ≤ 1.01 for blocking strategies REG and BLK. Different orderings (AMD, AMF, SCOTCH, METIS) are used.

				AMD					AMF				SCOTCH					METIS
	σ	PO REG	FT BLK ∆min	PO REG	FT BLK ∆min	PO REG	FT BLK ∆min	PO REG	FT BLK ∆min
		∆ NG	∆ NG		∆ NG	∆ NG		∆ NG	∆ NG		∆ NG	∆ NG
	5Hz	1.44	53 1.36	4	1.25	.75	51	.87	7	.68	.47 328	.32	3	.25	.43 230	.30	3	.24
	7Hz	5.03	38 4.44	4	4.35	15.3	18	17.8	12	14.9 1.60 287 1.14	3	.86	1.42 230 1.08	3	.82
	10Hz 19.3 154	19.8	3	15.3	96.9 253	99.1	11	82.1 5.86 288 4.21	3	3.05	4.67 281 3.44	3	2.53
	H0	.54 533	.53	4	.47	.12 333	.12	5	.091 .073 499 .063	3	.055	.077 615 .067	3	.057
	H3	135 380	106	5	9.07	101	63	134	17	9.26 2.19 615 1.76	5	1.18	1.95 533 1.53	5	1.12
	H17	183 266	226	7	136	467 173	558	50	396 22.8 380 19.0	5	12.6 21.49 242 16.8	4	12.3
	H116 2244	1	2244	1	2244 39383	1 39384	1 39383	291 109	224	4	153	264	78	215	4	157
	S3	20.9 725 17.9	6	15.1	20.7 184	24.8	10	17.8	4.5 771 3.41	5	2.72	3.24 771 2.64	5	2.09
	S21	393 685	349	5	311	1141 493	1352	77	831 50.9 492 39.6	4	31.7	34.3 223 28.5	5	24.9
	S84 3025 352 2848	5	2501 38664 725 45346 213 30977	289 286	228	4	193	207 171	174	4	151
	D30	115 111	121	8	94.5	1015 139	1280	75	825 16.7 156 12.8	5	8.77	15.5 144 13.0	5	8.61

Table 7

 7 Time (seconds) of the forward substitution according to the strategy used and number of groups NG created with BLK and REG (ND ordering). The percentage of actual computation time (BLAS) is indicated in parentheses.

		DEN	TP	INI	PO	REG	FT	BLK	NG BLK NG REG
	5Hz	3132 (34)	561 (80)	305 (88)	246 (89)	429 (95)	190 (89)	148 (91)	3	328
	7Hz	9101 (40) 1727 (89)	951 (93)	784 (94) 1137 (97)	594 (94)	460 (95)	3	255
	H0	862 (58)	161 (82)	125 (85)	97 (89)	88 (96)	79 (89)	67 (92)	4	328
	H3	12419 (72) 4478 (92) 3274 (94) 1901 (96) 1709 (98) 1624 (96) 1234 (97)	4	306
	S3	26328 (64) 6839 (89) 5005 (92) 3429 (95) 3450 (99) 2804 (96) 2195 (97)	5	536

Table 8

 8 Operation count (×1013 ) and time (seconds) of REG for different block sizes.

	REG	∆	64	Ts	128 ∆	Ts	256 ∆	Ts	512 ∆	Ts	1024 ∆	Ts
	5Hz	.24		188	.25	179	.27	188	.27	188	.31	209
	7Hz	.74		558	.78	538	.86	575	.93	611	.98	643
	H0	.051		77 .051	74 .052	72 .054	74 .058		79
	H3	.98 1447 1.02 1405 1.04 1369 1.09 1407 1.16 1499
	S3	1.75 2533 1.78 2448 1.80 2358 1.86 2379 1.91 2424

Table 9

 9 Time (seconds) of the flat tree and blocking algorithms for different orderings.

		ND	AMD	AMF	SCOTCH	METIS
		FT BLK	FT BLK	FT BLK	FT BLK	FT BLK
	5Hz	.77 .11	.64 .19 1.35 .34	.79 .20 1.02 .13
	7Hz 1.30 .30 1.56 .62 5.68 3.9 1.85 .56 2.48 .34
	H0	.09 .04 1.72 .04	.13 .04	.09 .02	.12 .04
	H3	.52 .26 1.56 .29 1.15 3.0	.55 .24 1.15 .54
	S3	.88 .33 2.20 .41 2.51 1.3	.84 .36 1.75 .49

Table 10

 10 Number of operations ∆ min and factor size for original (ND) and guided (GND) nested dissection orderings.

	Matrices	5Hz		7Hz		10Hz	H0	H3
	Strategy	ND GND	ND GND	ND GND	ND GND	ND GND
	∆ min (×10 13 )	.22	.19	.69	.62	2.26 1.99 .050 .025	.95	.81
	factor size (×10 9 ) 3.72	5.18 12.8	19.7 44.8	73.4	.24	.37 4.50	5.57

In this example, identifying supernodes with separators leads to relaxed supernodes: the sparsity in the interaction between u 1 and u 0 (and

u 2 and u 0 ) is not exploited to benefit from larger blocks.

http://mumps.enseeiht.fr/

http://www.labri.fr/perso/pelegrin/scotch/

http://glaros.dtc.umn.edu/gkhome/metis/metis/overview

http://mumps-solver.org

Acknowledgements. We thank the referees for their very constructive remarks, EMGS and Seiscope for providing the test cases, and F.-H. Rouet for his comments on an early version of this paper. This work was partially supported by LABEX MILYON (ANR-10-LABX-0070).


M Korti

B.-E Benkelfat

A Alfalou

Image Encryption based on Optical Coherence Multiplexing Technique

Introduction

Optical image encryption has received a great deal of attention [1]. Most of the proposed methods are based on double-random phase encoding [START_REF]ebraic implementation of double-random[END_REF][START_REF]counting double-random-phase encoding for secure image[END_REF]. Recently, Alfalou et al. proposed a new method based on Stokes-Mueller formalism [START_REF] Alfalou | Dual encryption scheme of images using polarized light[END_REF]. Other methods using nonlinear optical techniques [5] and quantum imaging [START_REF]-phase photon-counting double-random[END_REF] have been reported. None of these methods has exploited the properties of coherence modulation of light for image encryption, although this approach has already been used for signal encryption processing by Wacogne et al. [START_REF] Wacogne | Enhanced security in a coherence modulation system using optical path difference corruption IEEE[END_REF].

In this paper, we propose a novel approach of optical image encryption based on modulation coherence of light that allows several signals to be multiplexed as an optical path-difference (OPD) larger than the coherence length of the light source. This method has readily been implemented in a wide range of applications, including optical sensing, acousto-optic signal processing, and real-time optical image addition and subtraction [START_REF] Elwardi | Coherence Modulation of Light for Mathematical Operations[END_REF]. The proposed encryption technique enables two levels of security to be carried out simultaneously, on the principle and on the experimental implementation.

Theory

The main advantage of the coherence modulation of light is that it permits the multiplexing of several signals in a single light beam which allows real-time all optical arithmetic operations between two 2-D signals [10]. Using this advantage, we will show how to encrypt a target image.

The coding system is schematically shown in Fig. 1. The images y x f , Each CM is composed of a birefringent plate (Q i ) and a spatial light modulator (SLM i ) which permits to encode every pixel (x,y) of the images y x f i , as an optical path delay given by y x Kf y x , , is the key-random-image, K is a constant related to the SLM, 01 is the OPD introduced by Q 1 which is randomly variable around 0 . 02 is the OPD introduced by Q 2 so that 02 01 is always a constant. At the output of system, the intensity is given by,

d y x y x P y x I out , 2 cos 1 , 2 cos 1 4 1 , 2 1 0 (2)
where is the wave number and 0 P is the power-spectrum of the source centered at 0 .

The decoding system consists in using an interferometer with an OPD 02 01 d , and the intensity of the encrypted image is then given by, To evaluate the encryption level of our method, we used a histogram analysis which describes the distribution of pixel values in both target and encrypted images. Then we have compared the two images using a phase-only-filter (POF). Finally, we have performed a correlation technique between adjacent pixels. All these methods have shown that the target image has successfully been encrypted.

y x f y x f K I y x I , , 1 8 ,

Conclusion

To conclude, we have proposed a novel technique for real-time optical image encryption based on coherence modulation of light. We have presented the theoretical study and numerical simulation then evaluated the encryption level of our method using several techniques. The experimental verification of the proposed technique is in progress.

  the center wavelength and is the spectral width of the source).

  Fig.1. Schematic of the encryption setup

Fig. 2 .

 2 Fig.2. shows the target (a), the key-random-image (b) and the encrypted image (c) which are (100,100) pixels 256-gray-level-amplitude images.

.

  The target image (a), key-random-image (b) and the encrypted image (c).