Pierre Jolivet
email: jolivet@ann.jussieu.fr

Frédéric Hecht
email: hecht@ann.jussieu.fr

Frédéric Nataf
email: nataf@ann.jussieu.fr

Christophe Prud'homme
email: prudhomme@unistra.fr

Scalable Domain Decomposition Preconditioners

Keywords: G.1.8 [Numerical analysis]: Partial Differential Equationsdomain decomposition methods, finite element methods; D.2.13 [Software engineering]: Reusable Software-Reusable libraries Algorithms, Performance Linear solvers, divide and conquer, scalability

Domain decomposition methods are, alongside multigrid methods, one of the dominant paradigms in contemporary largescale partial differential equation simulation. In this paper, a lightweight implementation of a theoretically and numerically scalable preconditioner is presented in the context of overlapping methods. The performance of this work is assessed by numerical simulations executed on thousands of cores, for solving various highly heterogeneous elliptic problems in both 2D and 3D with billions of degrees of freedom. Such problems arise in computational science and engineering, in solid and fluid mechanics. While focusing on overlapping domain decomposition methods might seem too restrictive, it will be shown how this work can be applied to a variety of other methods, such as non-overlapping methods and abstract deflation based preconditioners. It is also presented how multilevel preconditioners can be used to avoid communication during an iterative process such as a Krylov method.

INTRODUCTION

Discretizations of partial differential equations used to model physical phenomena typically lead to larger and larger systems that cannot be solved directly and require advanced preconditioning techniques to ensure a fast convergence of iterative methods. As the current trend in high performance computing evolve towards more and more concurrency, recent results using domain decomposition preconditioners [START_REF] Klawonn | Highly scalable parallel domain decomposition methods with an application to biomechanics[END_REF] and multigrid methods [START_REF] Baker | Scaling hypre's multigrid solvers to 100,000 cores[END_REF][START_REF] Sundar | Parallel geometric-algebraic multigrid on unstructured forests of octrees[END_REF] clearly show why these are the methods of choice for achieving high-throughput finite element simulations. Both are examples of the divide and conquer paradigm. On the one hand, for scalar equations such as the simulation of flows through porous media, multigrid methods are often the preferred methods. One the other hand, domain decomposition methods are widely used in solid mechanics, in multiphysics or for code coupling. While these methods can offer high efficiency for solving PDE in parallel, special care must be taken when solving highly heterogeneous and complex systems for ensuring convergence rates independent of the number of computing nodes and heterogeneities involved in the simulation, as shown in figure 1. There are two main families of domain decomposition algorithms: overlapping Schwarz methods and iterative substructuring methods. While this work is more focused towards the former, the key idea is the same for both methods: to obtain convergence rates which are independent of the number of subdomains, also referred to as substructures, an additional component to provide global transfer of numerical information across all subproblems is needed. In this work, the construction of a so-called coarse operator is presented. In section 2, some important notions about domain decomposition methods are recalled, in section 3, the construction of the coarse operator is presented. Numerical results are gathered in section 3. [START_REF] Bangerth | deal.II -a general-purpose object-oriented finite element library[END_REF], where problems of up to 22 billions unknowns in 2D, and 2 billions unknowns in 3D are solved on more than 16K threads. Finally, it is shown in section 3.5 how one can use coarse grid correction to avoid communication during the iterations of a Krylov method.

DOMAIN DECOMPOSITION PRECON-DITIONERS

Let Ω ⊂ R d (d = 2 or 3) be a domain whose associated mesh can be partitioned into N non-overlapping meshes #iterations Relative residual error "Basic" preconditioning "Advanced" preconditioning Figure 1: Convergence of the GMRES preconditioned by two domain decomposition methods on 16 subdomains. The "basic" method is oblivious to the heterogeneities in the domain, while the "advanced" method take them into account in the preconditioner.

{Ti} 1 i N using graph partitioners such as METIS [START_REF] Karypis | A fast and high quality multilevel scheme for partitioning irregular graphs[END_REF] or SCOTCH [START_REF] Chevalier | PT-Scotch: a tool for efficient parallel graph ordering[END_REF]. Let V be the finite element space spanned by the finite set of n basis functions {φi} 1 i n defined on Ω. Typical finite element discretizations of a symmetric, coercive bilinear form a : V × V → R yield the following system to solve:

Ax = b, (1)
where (Aij) 1 i,j n = a(φj, φi), and (bi) 1 i n = (f, φi), f being in the dual space V * . Now let {Vi} 1 i N be the local finite element spaces defined on the domains associated to each {Ωi} 1 i N where Ωi is the subdomain defined as the union of all mesh elements in Ti, for all 1 i N . If δ is a positive integer, the overlapping decomposition T δ i 1 i N is defined recursively as follow: T δ i is obtained by including all elements of T δ-1 i plus all adjacent elements of T δ-1 i . For δ = 0, T δ i = Ti. An example of such a construction is given figure 2, for δ = 2. Now consider the restrictions {Ri} 1 i N from V to V δ i 1 i N , the local finite element spaces on Ω δ i 1 i N , and a local partition of unity {Di}

1 i N such that N j=1 R T j DjRj = In×n . (2)
Algebraically speaking, if {ni} 1 i N denotes the number of degrees of freedom in each local finite element spaces, then Ri is a boolean matrix of size ni × n, and Di is a diagonal matrix of size ni × ni for all 1 i N . In our experiments, we use a pretty simple partition of unity, already used for example in [START_REF] Kimn | Restricted overlapping balancing domain decomposition methods and restricted coarse problems for the Helmholtz problem[END_REF]. Let χi be a continuous piecewise linear function (hence mesh elements vertices and degrees of freedom coincide) defined on Ω δ i as such:

χi = 1 on all nodes of T 0 i 1 - m δ on all nodes of T m i \ T m-1 i ∀m ∈ [1; δ]
The local partition of unity is defined as:

χi = χi N j=1 χj| V δ i ∩V δ j ,
so that the support of the non negative function χi is V δ i and

N i=1 χi = 1.
Using a linear interpolant from the finite element space of continuous piecewise linear function from Ω δ i to V δ i (which is typically of higher order), Di can be obtained from χi for all 1 i N . For more complex partitions of unity, see for example [START_REF] Cai | Restricted additive Schwarz preconditioners with harmonic overlap for symmetric positive definite linear systems[END_REF]. Using the partition of unity, a common one-level preconditioner for system (1) introduced in [START_REF] Cai | Restricted additive Schwarz preconditioner for general sparse linear systems[END_REF] is:

P -1 RAS = N i=1 R T i Di(RiAR T i ) -1 Ri . (3)
Equation ( 3) clearly shows the need of the globally assembled matrix, or to be more precise, the need of N assembled submatrices (commonly referred to as "Dirichlet" matrices). These are typically hard to assemble directly and independently with classical finite element packages, in comparison with unassembled matrices (or "Neumann" matrices). Two approaches can be considered to build each "Dirichlet" matrix.

1. Build A, then extract each assembled subproblem. While this is the natural approach, it usually require some communications to build a parallel structure capable of handling distributed degrees of freedom, with ghost elements.

2. Build the stiffness matrices A δ+1 i yielded by the discretization of a on V δ+1 i , then, remove the columns and rows associated to degrees of freedom lying on elements of T δ+1 i \ T δ i . This yields Ai := RiAR T i , the global assembled matrix A is never assembled.

The second approach does not require any additional parallel information or communication: there is no need for a global ordering, associated to a global partition of the degrees of freedom. Local matrices Ai are symmetric positive definite, as A is. Usually, one subdomain is mapped to a single MPI process. For that reason, a given MPI process has a natural access to ui := Riu for any function u ∈ V . In all the following, we will assume that a finite element function in V δ i (resp. V ) can be interpreted as vector of R n i (resp. R n ) for all 1 i N . In order to compute a global sparse matrix-vector product Ax, one has to notice that, thanks to the partition of unity and the duplication of unknowns on the overlap, it can be proven that:

RiAR T j Djxj = RiR T j RjARjDjxj , (4) so that
(Ax)i = RiAx = RiA N j=1 R T j DjRjx = Ri N j=1 R T j RjARjDjRjx = N j=1
RiR T j AjDjxj .

(

) 5
Applying RiR T j to a vector xj of V δ j for all 1 i = j N is equivalent to restricting xj to the degrees of freedom in V δ j that are duplicated within V δ i , sending the resulting vector to subdomain i, which then prolongates by 0 the received vector outside of the overlap V δ j ∩ V δ i . In the following, Oi will be the set of neighboring subdomains to i, i.e. j : j = i and V δ j ∩ V δ i = ∅ , and Oi = Oi ∪ {i}. There is no need to assemble each Ri, one only needs to know the action of {RiR T j } 1 i N,j∈O i .

Problem specification

It is well known that one-level domain decomposition methods as depicted in the introduction of this section, see (3), do suffer from poor conditioning when used with many subdomains, [START_REF] Quarteroni | Domain decomposition methods for partial differential equations[END_REF][START_REF] Smith | Domain decomposition: parallel multilevel methods for elliptic partial differential equations[END_REF][START_REF] Toselli | Domain decomposition methods -algorithms and theory[END_REF]. Indeed, if the subdomains are assumed to be of size O(H), then the condition number of the preconditioner grows as 1/H for overlapping methods. For a given global mesh, increasing the number of subdomains leads to decreasing H, hence increasing the number of iterations needed for the preconditioned iterative method to converge. To overcome this recurrent problem in overlapping and non-overlapping methods, one must introduce a so-called coarse operator. In this work, an already established coarse operator whose theoretical foundations are presented in [START_REF] Spillane | A robust two-level domain decomposition preconditioner for systems of PDEs[END_REF] is used. From a practical point of view, after building each local solver Ai, three dependent operators are needed :

1. a deflation matrix Z of size n × m, with m n, 2. a coarse operator E = Z T AZ of size m × m,

the actual preconditioner

P -1 A-DEF1 = P -1 RAS (I -AZE -1 Z T ) + ZE -1 Z T , (6)
thoroughly studied in [START_REF] Tang | Comparison of two-level preconditioners derived from deflation, domain decomposition and multigrid methods[END_REF].

The choice of using P -1

A-DEF1 instead of, for example,

P -1 A-DEF2 = (I -AZE -1 Z T )P -1 RAS + ZE -1 Z T , (7)
which is also studied in [START_REF] Tang | Comparison of two-level preconditioners derived from deflation, domain decomposition and multigrid methods[END_REF] is pretty simple. While both preconditioners have similar numerical properties, applying P -1 A-DEF1 to a vector u requires only one coarse problem solution (used in two different operations afterwards):

ZE -1 Z T u .
On the other hand, applying P -1

A-DEF2 requires two coarse problem solutions:

ZE -1 Z T u
and ZE -1 Z T P -1 RAS u . Because applying a coarse correction is the most communicationintensive operation when preconditioning an iterative method, as shown in section 3.2, it is best to compute only 1 correction per iteration for scalability purposes.

In overlapping and non-overlapping decomposition methods, the deflation matrix is usually defined as:

Z = R T 1 W1 R T 2 W2 • • • R T N WN ∈ R n × R N i=1 ν i
The question is then, how to choose all those Wi, for all 1 i N , to build a numerically scalable preconditioner ?

In [START_REF] Spillane | A robust two-level domain decomposition preconditioner for systems of PDEs[END_REF], they are built as:

Wi = DiΛi 1 DiΛi 2 • • • DiΛi ν i ∈ R n i × R ν i 1 i N (8)
A threshold criterion is used to select the νi eigenvectors {Λi j } 1 j ν i associated to the smallest eigenvalues in magnitude of the following local generalized eigenvalue problems:

A δ i Λi = λiDiR T i,0 Ri,0A δ i DiΛi (9)
where A δ i is the matrix yielded by the discretization of a on V δ i (unassembled submatrix), and Ri,0 is the restriction operator from V δ i to the overlap

V δ i ∩ ∪j∈O i V δ j . If Ω δ i ∩ ∂Ω = ∅, Ω δ
i is commonly referred to as a floating subdomain. Hence, the bilinear form a lacks essential boundary conditions on V δ i , so that A δ i is now symmetric positive indefinite (compared to Ai, the assembled submatrix, which is always symmetric positive definite under the assumption that A is). These independent and local eigenproblems are solved concurrently in order to find low eigenvalues that are in some sense close the lowest eigenvalues of the preconditioned system. It is well known that those are hurting the convergence rate of traditional Krylov methods, see [START_REF] Nicolaides | Deflation of conjugate gradients with applications to boundary value problems[END_REF][START_REF] Frank | On the construction of deflation-based preconditioners[END_REF] for further details. With this construction of the coarse operator, it can be proven that the condition number of the preconditioned system is now independent of the size of the subdomains, of the number of subdomains, and of the heterogeneities in the physical coefficients. Once again, interested readers are referred to [START_REF] Spillane | A robust two-level domain decomposition preconditioner for systems of PDEs[END_REF].

DESIGN OF THE COARSE OPERATOR

The goal of this section is to explain how a coarse operator built using deflation vectors can be efficiently assembled for large-scale simulations. This work is implemented in a light and versatile C++ framework that is not directly linked to domain decomposition methods, meaning that it is possible to use it to assemble coarse operator with other abstract deflation vectors, for example as defined in [START_REF] Grigori | A parallel two-level preconditioner for cosmic microwave background map-making[END_REF] for simulations in cosmology. When applicable, the MPI calls and operations related to linear algebra (either from dense BLAS [START_REF] Blackford | An updated set of basic linear algebra subprograms (BLAS)[END_REF] or sparse BLAS such as Intel Math Kernel Library or cuS-PARSE) are provided in typewriter font.

Assembling the coarse operator

Looking at the general formulation of E introduced in the previous section, two global sparse matrix-matrix products must be computed to assemble the coarse operator:

AZ then Z T (AZ) .
However, it is possible to exploit the sparsity pattern of Z to build E more efficiently than with these two consecutive computations. By construction, the deflation matrix Z is made of block of dense matrices of size ni × νi as displayed in figure 3. If there is more than one block for which there are non zeros in a given row of Z, it means that one is dealing with a duplicated unknown (those which are in the overlap). From the simple example in figure 3 with 4 subdomains, one can infer that O1 = {2}, O2 = {1, 3}, O3 = {2, 4}, and O4 = {3}. Just like the globally assembled matrix A, the global deflation matrix Z is never assembled, but its representation above is useful to understand the assembly of E. Instead, each subdomain has access to its local dense matrix Wi. The block (i, j) of E of size νi×νj is then equal to W T i RiAR T j Wj. From that algebraic definition, one can see that the sparsity pattern of E is linked to the connectivity between subdomains, because:

RiAR T j = 0n i ×n j ⇐⇒ V δ i ∩ V δ j = ∅ .
In the context of overlapping domain decomposition methods, using equations ( 4)- [START_REF] Chevalier | PT-Scotch: a tool for efficient parallel graph ordering[END_REF], the block (i, j) of E can be computed as

W T i RiAR T j Wj = W T i RiR T j RjAR T j Wj = WiRiR T j AjWj . (10)
As in equation ( 5), the previous equation shows how it is possible to take advantage of the duplicated unknowns on the overlap to compute a global product with only local computations and peer-to-peer transfers : the assembled matrices A and Z are not needed. Noticing that RiR T i = In i ×n i , the subdomain i needs to perform three tasks so that all blocks of E can be computed:

1. compute locally Ti = AiWi (csrmm) and Ei,i = W T i Ti (gemm), 2. send to each neighboring subdomain j ∈ Oi, Sj = RjR T i Ti, and receive from each neighboring subdomain Uj = RiR T j Tj,

compute locally

Ei,j = W T i Uj (gemm).
The cost of step 2 is approximately the same as one global sparse matrix-vector product. Using the same toy problem as in figure 3, the following figure is a representation of the sparsity pattern of E. Two colors are used to differentiate blocks that can be computed without any communication and those that cannot.

The three tasks to assemble all blocks of E can be customized to suit one's needs using C++ polymorphism. For example, in the context of non-overlapping methods, the sparsity pattern of E is typically more dense: a block (i, j) of E is not null if and only if j ∈ Oj, but also if there exists k ∈ Oj such that k ∈ Oi. This can be handled by our framework. Likewise, each local computations Ti as introduced in the first step, can be modified to involve more complex operations and communications patterns.

Stopping at that point in the construction of E, solving systems involving the coarse operator would require either to factorize a relatively small matrix on a large number of MPI processes, or to call an iterative method with once again, a distribution of data that would lead to a too fine-grained granularity (blocks of rows of E are typically of size νi ranging from 1 to 30). In both cases, the increased communication overhead would lead to bad performances when solving multiple systems involving the coarse operator. Another approach would consist in replicating E over all subdomains, and then performing independent factorizations, but that is simply not feasible for large decompositions, with arbitrary numbers of deflation vectors. In the following subsections, a more efficient data distribution is explained. While it will imply more communications during the setup of the coarse operator, it will cover all possible issues stated in this paragraph, regarding communication overhead, memory consumption and fast computation of coarse solutions.

Master-slave approach

The idea is to use only a "small" group of processes that will be in charge of factorizing the coarse operator and that will afterwards be called for computing solutions of systems involving E -1 using a distributed sparse direct solver. This is inspired by the famous master-slave approach. For the rest of the paper, the following notations will be thoroughly used (for the sake of completeness, we provide their type in our implementation):

• unsigned short P: the number of masters, chosen at runtime by the user,

• MPI_Comm masterComm: a communicator between all masters, set to MPI_COMM_NULL on slaves, on which will be instantiated the distributed solver,

• MPI_Comm splitComm: a communicator between a master and its slaves in which the rank of the master is always 0, and the ranks of the slave follow the same order as in MPI_COMM_WORLD.

A representation of such communicators is given figure 5.

Prior to factorizing E, the first step is to assemble it in a distributed matrix on the masters. Each master will be in charge of assembling all the values of its slaves. It is assumed that the format in which the distributed matrix is stored is a simple global CSR or global COO -the standard format for most linear solvers available nowadays, meaning that for each non-zero value, one must know the absolute row and column indices of the given value in the global matrix E. For a process i, the global row indices IE i of all the blocks {Ei,j} j∈O i range from ri to ri + νi, where ri = i-1 j=1 νj, and the global column indices JE j of all the blocks {Ej,i} j∈O i range from ri to ri + νi. The simplest approach would then be to: The number of non-zero values for a process i is:

size(valE i ) = νi × νi
Blue values in figure 4 + νi ×

j∈O i νj
Red values in figure 4 , [START_REF] Frank | On the construction of deflation-based preconditioners[END_REF] meaning that prior to the three MPI_Gatherv in step 3, a call to MPI_Gatherv(Oi) on splitComm with rank 0 as root must be made to allocate the right buffers for the distributed format IE, JE, and valE on each master. While this approach is somehow natural for assembling the distributed matrix on the masters -because of the ordering of the rank of the slaves, calling MPI_Gatherv is similar to concatenating all local chunks of E -it implies a lot of (unnecessary) communications. In particular, why should slaves send to masters the global row and column indices ? Indeed, at the end of assembly, only the masters have access to the distributed coarse operator, so it is their responsability to compute the indices. The slaves should not have to store or compute anything related to the distributed format. The following approach has the advantage of transferring only what is needed from one slave i to its master: the array of double valE i . The indices will be computed after reception by each master, meaning that the memory overhead on the slaves is null (no integers are allocated for storing any index). The new workflow is now:

1. perform a neighborhood collective operation MPI_Ineighbor_alltoall(νi) 1 on the communicator to which the distributed graph topology information of the connectivity between subdomains is attached (MPI_Dist_graph_create_adjacent). Then allocate accordingly the buffers Sj and Uj for all j ∈ Oi

call MPI_Gather([νi, |Oi|]) (array of 2 integers) on split-

Comm with rank 0 as root so that masters can preallocate the distributed CSR (IE, JE, valE),

1 new to the MPI-3 standard 3. assemble locally {Ei,j} j∈O i as previously and send the values to the master. Prepend to the beginning of the message, the values of Oi, i.e. the final size of the message is |Oi| + [START_REF] Frank | On the construction of deflation-based preconditioners[END_REF].

Additionally, the masters must concatenate all νi gathered in step 2, using MPI_Allgatherv on commMaster to be able to compute all cumulative sums ri, for all i ∈ splitComm. This call is equivalent to the MPI_Allgather in step 1 of the "natural" algorithm, but this time it doesn't involve any slave. When a master receives a message from a slave i, it knows that the global row index ranges from ri to ri+νi, and because the first values of the received message are a copy of Oi, it can compute the correct global column indices for the neighbors of this slave. The complete algorithm for assembling the coarse operator is summarized in algorithms 1 (construction of all blocks of E) and 2 (distributed assembly on the masters).

Electing the masters

In the previous paragraph, the masters are defined in a rather abstract way, as the processes that have a rank equal to 0 in splitComm. We have two ways to define the aforementioned MPI communicator. The first is the natural distribution: the process are spread uniformly and contiguously into P groups, the masters are of rank i • N/P , for all 0 i P -1. The second distribution is a little more advanced and better suited for assembling symmetric coarse operators. In that case, one only needs to assemble the upper part of the distributed CSR or COO, so that only the following blocks are computed and assembled: Ei,j, ∀1 i N, ∀j ∈ Oi : j i .

Moreover, only the upper parts of the dense diagonal blocks Ei,i for all 1 i N are needed. To ensure load balancing between masters, the processes are now spread contiguously but non-uniformly, with masters of rank pi, where pi is defined by the following sequence to ensure heuristically that the number of values within each quadrilateral in figure 5 is the same: The procedure compute in algorithm 1 is in charge of returning a dense array of values corresponding to the block of E given in argument, while the procedure assemble in algorithm 2 is in charge of computing the indices of the values of the block of E given in argument and store them in the distributed matrix representation.

p0 = 0 pi = N -(pi-1 -N ) 2 -N 2 /P + 0.5 , 1 i P -1 4

Applying a coarse operator correction

Once E has been assembled, it is involved in the solution of the following problem for a given vector u: ZE -1 Z T u. It is assumed that the right-hand side Z T u and the solution E -1 Z T u of the coarse problem are kept distributed on mas-terComm at all time, as input and ouput of the distributed solver. The computation of the correction can obviously be broken down into three basic operations:

1. compute Z T u = w ∈ R N i=1 ν i .
Once again, the structure of Z makes it possible to compute this sparse matrix-vector product without any explicit global representation of Z. Indeed,

Z T u = N i=1 W T i Riu
This is evaluated by having all subdomains i compute locally wi = W T i ui (gemv), and calling MPI_Gather(v) on each splitComm with rank 0 as root to assemble w with all local contributions of size νi, for all i in splitComm, on the masters.

compute E

-1 w = y ∈ R N i=1 ν i .
This operation must be as fast and reliable as possible, since it is carried out by only few masters. The numerical factorization of E as computed in the previous paragraph during the assembly phase is reused for each forward elimination and back substitution.

compute Zy

= z ∈ R N i=1 n i .
This is the exact dual of step 1. This operation reads:

Zy = N i=1 R T i Wiyi
First, a call to MPI_Scatter(v) is made on each split-Comm with rank 0 as root so that each subdomain i can retrieve its yi of size νi from its master. Then they all compute zi = Wiyi (gemv). Finally,

(Zy)i = N j=1 RiR T j zj (12)
is computed using the same communication procedure as for the sparse matrix-vector product, see [START_REF] Blackford | An updated set of basic linear algebra subprograms (BLAS)[END_REF].

Using this construction, it is clear that using a non-uniform criterion νi for each subdomain leads to using MPI communications with varying counts of data from each process (hence the v in parentheses). Because these communications scale as O(N ), it is preferable to call prior to assembling the coarse operator MPI_Allreduce(νi,MPI_MAX). That way, it is possible to use MPI communications with equal counts of data, which typically scale as O(log(N )). Moreover, the theoretical estimate on the condition number of the preconditioner system remains valid. There are also multiple simplifications possible during the assembly of E, because using a uniform criterion means that any process i already knows the value of νj, for all j ∈ Oi without needing the call to MPI_Ineighbor_alltoall and such.

Cost analysis

As a consequence of the observation made in the previous paragraph, it is assumed in this section that the number of deflation vectors per subdomain is chosen uniformly and set to ν. This will also be the case in section 3.4 for our scaling experiments, and we now propose a comprehensive summary of the cost of our two-level preconditioner, compared to a simple one-level method. As a quick reminder, the cost of applying a one-level preconditioner such as P -1 RAS in (3) is the same as computing N independent solutions on each subdomain -application of (RiAR T i ) -1 , then performing one global matrix-vector product.

Memory footprint. Because it is clear that no global structure is needed for the construction and the use of Z and E in our implementation, the memory footprint is straightforward to evaluate. For slaves, a dense array of size ν × ni is needed to store Wi. For masters, the same array must be stored, plus the factors of E -1 (whose size depends on the distributed solver used).

Messages. During the construction of E, each process i must send and receive a single message from each neighboring subdomains j ∈ Oi. The size of theses messages are ν × size of overlap between V δ i and V δ j , see Sj and Uj in algorithm 1. Then each slave i must send a single message of double, of size |Oi| + |Oi| × ν × ν to its master, see msg-ToMaster in algorithm 2, and each master must reciprocally receive these message for all slaves in splitComm. During the solution, at each iteration, a MPI_Gather and a MPI_Scatter are called on each splitComm. The messages passed at the end of the correction step in equation ( 12) are of the exact same sizes as for the global matrix-vector product for the one-level method, c.f. ( 5).

Computational intensity.

The major contribution is the solution of each local and independant eigenvalue problem. Afterwards, for the setup of E, the BLAS computations stated in section 3.1 are negligible with respect to the time spent transferring messages, and at each iteration the costly operation in section 3.2 is the solution of the coarse correction -application of E -1 using a distributed direct solver. Note that for big subdomains, this cost is negligible compared to the one of computing each local solution for the one-level preconditioner P -1 RAS , c.f. (3).

Numerical results

Results in this section were obtained on Curie, a Tier-0 system for PRACE 2 composed of 5040 nodes made of 2 eightcore Intel Sandy Bridge processors clocked at 2.7 GHz. The interconnect is an InfiniBand QDR full fat tree and the MPI implementation used was BullxMPI version 1.1.16.5. Intel compilers and Math Kernel Library in their version 13.1.0.146 were used for all binaries and shared libraries, and as the linear algebra backend for both dense and sparse computation in our framework. LU or LDL T decomposition of the local problems are computed with either instances of MUMPS [START_REF] Amestoy | A fully asynchronous multifrontal solver using distributed dynamic scheduling[END_REF][START_REF] Amestoy | Hybrid scheduling for the parallel solution of linear systems[END_REF] or PaStiX [START_REF] Hénon | PaStiX: a high performance parallel direct solver for sparse symmetric positive definite systems[END_REF] on MPI_COMM_SELF, or preferably using Intel MKL PARDISO, University of Lugano PAR-DISO [START_REF] Schenk | Solving unsymmetric sparse systems of linear equations with PARDISO[END_REF][START_REF] Schenk | On fast factorization pivoting methods for sparse symmetric indefinite systems[END_REF], or WSMP [START_REF] Gupta | WSMP: Watson sparse matrix package -part I: Direct solution of symmetric systems[END_REF][START_REF] Gupta | WSMP: Watson sparse matrix packagepart II: Direct solution of general systems[END_REF] which are by design more suited for multithreaded computations. LU or LDL T decomposition of the coarse operator is computed using MUMPS, PaStiX, or PWSMP. Eigenvalue problems to compute the 2 Partnership for Advanced Computing in Europe deflation vectors as defined in ( 9) are solved using ARPACK [START_REF] Lehoucq | ARPACK users' guide: solution of large-scale eigenvalue problems with implicitly restarted Arnoldi methods[END_REF]. Finite element matrices are obtained from FreeFem++, but the design of our domain decomposition framework makes it simple to port to other Domain-Specific (Embedded) Language such as Feel++ [START_REF] Prud'homme | Feel++: A computational framework for Galerkin methods and advanced numerical methods[END_REF], FEniCS [START_REF] Logg | Automated solution of differential equations by the finite element method[END_REF] or finite element libraries like deal.II [START_REF] Bangerth | deal.II -a general-purpose object-oriented finite element library[END_REF], GetFem++ and such. We display the speedup and efficiency in terms of number of MPI processes. In these experiments, each MPI process is assigned a single subdomain, and 2 OpenMP threads. The GMRES is stopped when a relative 10 -6 decrease of the residual is reached.

First, the system of linear elasticity with highly heterogeneous elastic moduli is solved with a minimal geometric overlap of one mesh element. Its variational formulation reads:

a(u, v) = Ω λ∇ • u∇ • v + 2µε(u) : ε(v) + Ω f • v + ∂Ω g • v
where

• λ and µ are the Lamé parameters such that µ = E 2(1 + ν) • ε is the linearized strain tensor, f are the body forces (in this case, only the gravity), and g are the surface force (in this case, a vertical loading is imposed on some parts of the geometries).

and λ = Eν (1 + ν)(1 -2ν) (E
Such an equation typically arises in computational solid mechanics, for modeling small deformations of bodies. In 2D, we use piecewise cubic basis functions (∼ 33 nnz per row).

The system is of constant size equal to approximately 2 billions unknowns. In 3D, piecewise quadratic basis functions are used (∼ 83 nnz per row). The system is of constant size equal to approximately 300 millions unknowns. Both geometries are displayed in figure 6 and were meshed by Gmsh [START_REF] Geuzaine | Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities[END_REF] and partitioned with METIS. After the partitioning step, each local mesh is refined concurrently by splitting each triangle or tetrahedron into multiple smaller elements. This means that we start the simulation with a relatively "coarse" global mesh (26 million triangles in 2D, 10 million tetrahedra in 3D), which is then refined in parallel (thrice in 2D, twice in 3D). We get a nice speedup from 1024 × 2 = 2048 to 8192 × 2 = 16384 threads as shown in figure 8. According to the table in figure 8, the costly operations in the construction of the preconditioner are the solution of each local eigenvalue problem (9) (column deflation), and the factorization of each local solver Ai (column factorization). In 3D, the complexity of such operations typically grows superlinearly with respect to the number of unknowns. That explains why we can achieve superlinear speedup. At peak performance, on 16384 threads, the speedup relative to the runtime on 2048 threads equals 530.56 51.76 ≈ 10. In 2D, the computation costs are lower, and tend to scale better with the number of unknowns, which makes it harder to achieve high speedup for larger number of subdomains. At peak performance, on 16384 threads, the speedup relative to the runtime on 2048 threads equals 213.20 34.54 ≈ 6. In both cases, the solution of the eigenproblem is the limiting factor for achieving better speedups. This can be explained by the fact that the Arnoli method, which ARPACK is based on, tends to perform better for larger ratios ni νi , but these values decrease as the subdomains get smaller. We will get back to this issue in the conclusion. For larger but fewer subdomains, the time to compute the solution (column solution), i.e. the time for the GMRES to converge, is almost equal to the forward eliminations and back substitutions in the subdomains times the number of iterations. When the decompositions become bigger, subdomains are smaller, hence each local solution is computed faster and global communications have to be taken into account. For a comparison similar to the one in figure 1 in this bigger configuration, we display in figure 7 the convergence histogram of a simple one-level method versus this two-level method. One can easily understand that, while the cost of building the preconditioner cannot be neglected, it is necessary to ensure the convergence of the Krylov method.

Moving on to the weak scaling properties of our framework, the problem we now solve is a scalar equation of diffusivity with highly heterogeneous coefficients (varying from 1 to 3 • 10 6 as displayed in figure 9) on RAS and P -1 A-DEF1 for the problem of linear elasticity in 2D using 1024 subdomains. Timings for the setup and solution phases using P -1

A-DEF1 are available in figure 8. Using P -1 RAS , the convergence is not reached after 600 seconds. row). Its variational formulation reads:

a(u, v) = Ω κ∇u • ∇v + Ω f • v
where f is a source term. Such an equation typically arises for modeling flows in porous media, or in computational fluid dynamics. No change is needed in our framework, since all operations are algebraic. The only work needed outside of our framework is changing the mesh used for computations, as well as the variational formulation of the problem in the FreeFem++ DSL. On average, there is a constant number of degrees of freedom per subdomain equal to 280K in 3D, and near 2.7 millions in 2D. As for the strong scaling experiment, after building and partitioning a global "coarse" mesh (with few millions of elements), each local mesh is refined independently to ensure a constant size system per subdomain as the decomposition becomes bigger. The efficiency remains near the 90% mark, thanks to almost no variability in the factorization of the local problems and the construction of the deflation vectors. In 3D, the initial problem of 74 millions unknowns is solved in 200 seconds on 512 threads. Using 16384 threads, the problem is now made of approximately ≈ 96%. At such scales, the most penalizing step in the algorithm is the construction of the coarse operator, specially in 3D, with a non negligible increase in the time spent to assemble E. Finally, we present in this last paragraph the performances of our framework to assemble and factorize the coarse operator E for all the previous simulations. Tables in figures 8 and 10 already included these timings in their last column total (> factorization + deflation + solution), but for more in depth analysis, they are reported next separately.

The figure 11 includes all timing relative to algorithms 1 and 2 described in section 3.1.1: the construction of the communicators, the assembly of E, and its numerical factorization. The most consuming part of the algorithm is the actual transfer and the assembly by the masters. Especially in 3D, when the coarse operator is becoming less and less sparse (directly linked with the average value of |Oi|), it is likely to become a problem for even larger decomposition. Note that the MPI implementation used for these experiments is not thread compliant (and in particular, does not support a level of thread support equal to MPI_THREAD_MULTIPLE), meaning that some "unnecessary" #pragma omp critical had to be used during the assembly by the masters. At these scales, another problem is the factorization of E. Indeed, increasing the number of masters P does not always have a beneficial effect for this concern, because distributed solvers have difficulties scaling beyond ∼ 128 processes.

Communication avoiding multilevel preconditioners

Classical Krylov methods require global synchronizations during the orthogonalization and the normalization of the basis vectors at each iteration that can hamper parallel scalability. Efforts have been made to reduce the number of such global operations [START_REF] Sturler | Reducing the effect of global communication in GMRES(m) and CG on parallel distributed memory computers[END_REF][START_REF] Chronopoulos | s-step iterative methods for symmetric linear systems[END_REF], and recent alternatives, such as the p 1 -GMRES introduced in [START_REF] Ghysels | Hiding global communication latency in the GMRES algorithm on massively parallel machines[END_REF], predict considerable improvement in speedups because they can overlap computations of sparse matrix-vector products and global synchronizations.

For further references in this paper, it is best to quote the actual computational loop of the p 1 -GMRES.

Algorithm Main loop of the p 1 -GMRES introduced in [13] 1: for i = 0, . . . , m do 2:

w ← Avi 3: if i > 1 then 4: vi-1 ← vi-1/hi-1,i-2 zi ← zi/hi-1,i-2 5: w ← w/hi-1,i-2 hi-1,i-1 ← hi-1,i-1/h 2 i-1,i-2 6: hj-1,i-1 ← hj-1,i-1/h 2 i-1,i-2 , j = 0, . . . , i -2 7: end if 8: zi+1 ← w -i-1 j=0 hj,i-1zj+1 9: if i > 0 then 10: vi ← zi -i-1 j=0 hj,i-1vj hi,i-1 ← ||vi||2 11: end if 12:
hj,i ← zi+1, vj , j = 0, . . . , i 13: end for When using a multigrid or domain decomposition based preconditioner, line 2 becomes w ← P -1

A-DEF1 Avi and it gets harder to hide the latency of the global reductions line 10 and 12 because applying the preconditioner might require a global synchronization (induced by the solution of a coarser problem).

Applying the preconditioner P -1

A-DEF1 built in this section requires one MPI_Gather(v) on each splitComm, followed by a coarse solution, followed by one MPI_Scatter(v) on each splitComm with rank 0 as root. While this shows that there is no global communication, data dependencies introduce a "virtual" global synchronization (using non-blocking collective operations is not possible). Moreover, the two global communications on each splitComm involve mainly short messages if the number of deflation vectors is low (typically νi < 30 for each subdomain), meaning that passing extra data is unlikely to hurt too much the execution time of these operations. This is why we propose a fused p 1 -GMRES, where:

• no global reduction is performed line 10 and 12 (only local computations),

• during the call to MPI_Gather(v) at the i-th iteration of the GMRES line 2 on each splitComm with rank 0 as root, each process also sends to its master local values of hi-1,i-2 and hj,i-1 for j = 0, . . . , i -1 (these are the non-reduced values of the previous iteration), • prior to calling the coarse solver, each master reduces the values received after the call to MPI_Gather(v) and then launches a MPI_Iallreduce on masterComm.

It is then time for them to compute the solution of the coarse problem, after which they wait for the nonblocking reduction to finish,

• eventually, the reduced values of hi-1,i-2 and hj,i-1 for j = 0, . . . , i -1 are scattered back with the coarse solution using a single MPI_Scatter(v) on each split-Comm with rank 0 as root.

Each iteration of this fused p 1 -GMRES can be seen as a coarse operator correction, without a single additional global communication or synchronization and only a single additional call to MPI_Iallreduce between the masters on masterComm, whose size is greatly smaller than the number of subdomains N . The communications induced by the MPI_Iallreduce can be overlapped with the computation of the coarse solution. While this approach looks promising, our test cases are too small and the numbers of iterations are relatively small, so that both pipelined GMRES are performing approximately the same as the reference GMRES. Such a direction of research is however becoming essential as the concurrency of modern supercomputers is increasing tremendously, and minimizing synchronizations is becoming essential for large decompositions, as shown by theoretical models in [START_REF] Ghysels | Hiding global communication latency in the GMRES algorithm on massively parallel machines[END_REF] for the unpreconditioned GMRES.

CONCLUSION AND OUTLOOK

In this paper, we assess the efficiency of our implementation of a new adaptative two-level preconditioner suited for various problems such as linear elasticity or Darcy's law with high-contrast coefficients. The preconditioner was recently theoretically introduced in [START_REF] Spillane | A robust two-level domain decomposition preconditioner for systems of PDEs[END_REF] and it provides a sound background for evaluating the performance of our portable framework. Using state of the art multithreaded direct linear solvers and eigensolvers and distributed direct linear solvers, we show experimentally that our approach is well suited for large-scale simulations in both 2D and 3D with high-order finite element methods, on up to 16384 threads, enabling simulations for two classes of elliptic problems of more than 22 billion unkowns in 2D and 2 billion unkowns in 3D.

For building such a preconditioner, the deflation vectors used in our coarse operator are currently determined a priori by solving independently local eigenvalue problems. It is also possible to retrieve them a posteriori during the convergence of the iterative method, using for example approximations of the Ritz vectors. Using new theoretical results, we are currently evaluating a new a posteriori construction that doesn't require the solution of the aforementioned eigenvalue problems. This will hopefully improve our results for more fine-grained granularity in our strong scaling experiments. Finally, we are currently investing other types of systems for which domain decomposition methods have been proven to be efficient, namely nonlinear problems in computational solid mechanics, using non-overlapping methods. Thanks to the versatile design of our framework, and the flexibility of current finite element Domain-Specific (Embedded) Languages such as FreeFem++, new experiments should be soon possible.

Figure 2 :

 2 Figure 2: Decomposition of the SC conference logo (no copyrights infringement intended) into three color-coded subdomains. On the left, the decomposition is non-overlapping, δ = 0. On the right, two consecutive extensions are performed, δ = 2, and represented in black overlay, to yield Ω 2 i 1 i 3 .

Figure 3 :

 3 Figure 3: Sparsity pattern of the deflation matrix Z (grey blocks represent nonzero entries) with 4 subdomains.

Figure 4 :

 4 Figure 4: Sparsity pattern of the coarse operator E. Blue blocks involve only local computations (step 1) while red blocks also involve peer-topeer transfers (steps 2 and 3).

Figure 5 :

 5 Figure 5: Distribution of E when build with 16 subdomains using 4 masters. Each color represents a different splitComm, each number represents the rank of the master (in MPI_COMM_WORLD) of a given split-Comm. On the right, the number of values per split-Comm is roughly the same if the values below the diagonal are dropped (symmetric coarse operator).

  being Young's modulus and ν Poisson's ratio). They are chosen to vary between two sets of values, (E1, ν1) = (2 • 10 11 , 0.25), and (E2, ν2) = (10 7 , 0.45).

Figure 6 :

 6 Figure 6: Tripod used for our 3D and cantilever used for our 2D strong scaling experiments. Black and light grey are used to represent the variations of the Young's modulus (200 GPa and 0.01 GPa) and Poisson's ratio (0.25 and 0.45).

Figure 7 :

 7 Figure 7: Convergence of the GMRES(40) preconditioned by P -1RAS and P -1 A-DEF1 for the problem of linear elasticity in 2D using 1024 subdomains. Timings for the setup and solution phases using P -1A-DEF1 are available in figure8. Using P -1 RAS , the convergence is not reached after 600 seconds.

Figure 8 :

 8 Figure 8: Strong scaling experiments.

Figure 9 :

 9 Figure 9: Diffusivity κ used for our 2D weak scaling experiment with channels and inclusions. Black and light grey are used to represent the variations of κ from 1 to 3 • 10 6 .

Figure 10 :

 10 Figure 10: Weak scaling experiments.

Figure 11 :

 11 Figure 11: Timings for assembling coarse operator. Results are gathered two-by-two, the first column is for the diffusivity problem, the second is for the elasticity problem.

  1. call MPI_Allgather(νi) on MPI_COMM_WORLD to be able to compute the cumulative sums {rj} j∈O i and to allocate the buffers Sj and Uj for all j ∈ Oi, 2. assemble locally {Ei,j} j∈O i as previously and store the values in CSR (or COO): (IE i , JE i , valE i ), 3. call MPI_Gatherv(IE i ), MPI_Gatherv(JE i ), and MPI_Gatherv(valE i ) on splitComm with rank 0 as root.

Algorithm 1

 1 Schematic construction of Ei,j, for all 1 i N and j ∈ Oi

		MPI_Ineighbor_alltoall(νi)
	2: MPI_Gather([νi, |Oi|], splitComm, 0)
		compute(Ti)	Ti = AiWi
	4: for j ∈ Oi do	After completion of line 1
		MPI_Isend(Sj, j, MPI_COMM_WORLD)	Sj = RjR T i Ti
	6:	MPI_Irecv(Uj, j, MPI_COMM_WORLD, rq[j])
		end for
	8: compute(Ei,i)	Diagonal block
		for j ∈ Oi do
	10:	MPI_Waitany(rq, &index)
		compute(E i,index )	Off-diagonal block
	12: end for
	Algorithm 2 Schematic assembly of E on the masters
		buildComm(P, splitComm, masterComm)
	14: if masterComm != MPI_COMM_NULL then	Master
		MPI_Allgatherv(νi)
	16:	Now receive all messages from the slaves
		for j = 1, . . . , MPI_Comm_size(splitComm) do
	18:	MPI_Irecv(msgFromSlave[j], j -1, splitComm,
		rq[j])
		end for
	20:	assemble(Ei,i)
		for k = 1, . . . , |Oi| do
	22:	assemble(E i,k )
		end for
	24:	Blocks local to the masters have been assembled
		for j = 1, . . . , MPI_Comm_size(splitComm) do
	26:	MPI_Waitany(rq, &index)
		assemble(E index,index );
	28:	for k = 1, . . . , |O index | do
		assemble(E index,msgFromSlave[index][k] )
	30:	end for
		end for
	32:	Blocks from the slaves have been assembled
		numericalFactorization(E)
	34: else	Slave
		msgToMaster = Oi
	36:	concatenate(msgToMaster, Ei,i)
		for k ∈ Oi do
	38:	concatenate(msgToMaster, E i,k )
		end for
	40:	MPI_Isend(msgToMaster, 0, splitComm)
		Send the index of the neighbors as well as the local
		rows of E computed in algorithm 1
	42: end if

  2.3 billions unknowns, and it is solved in 215 seconds, which yields an efficiency of 200 • 2.3 • 10 9 215 • 7.4 • 10 7 • 32 ≈ 90%. In 2D, the initial problem of 695 millions unknowns is solved in 175 seconds on 512 threads. Using 16384 threads, the problem is now made of approximately 22.3 billions unknowns, and it is solved in 187 seconds, which yields an efficiency of 175 • 2.2 • 10 10 187 • 6.7 • 10 8 • 32

  .2 13.7 143.0 • 10 6 85.3 • 10 6 10.75 s 8.04 s 8192 64 48 131072 131072 14.7 14.6 285.4 • 10 6 196.7 • 10 6 19.92 s 17.30 s

	3D
	N 256 512 1024 8 4 6 2048 12 12 40960 P dim(E) 5120 10240 8 20480 22528 13.0 12.0 |Oi| (average) 11.5 12.4 40960 13.8 12.9 4096 18 22 73728 73728 142D	nnz(E -1 ) 2.6 • 10 6 7.7 • 10 6 24.5 • 10 6 14.6 • 10 6 74.6 • 10 6 34.5 • 10 6	Time 2.78 s 3.35 s 4.42 s 11.25 s 6.91 s 5.68 s
	N 256 512 1024 10 8 P 2 4 2048 14 12 38912 dim(E) 5376 10240 20480 24576 5.7 |Oi| (average) 5.5 5.6 5.5 40960 5.8 5.7 4096 22 18 81920 73728 5.9 5.8 8192 36 36 163840 122880 5.9 5.8	nnz(E -1 ) 959.7 • 10 3 2.1 • 10 6 5.0 • 10 6 5.5 • 10 6 10.8 • 10 6 9.3 • 10 6 10.05 s Time 9.39 s 9.96 s 9.92 s 10.14 s 6.20 s 28.4 • 10 6 5.10 s 17.7 • 10 6 10.87 s 73.3 • 10 6 30.4 • 10 6 13.27 s 6.96 s

ACKNOWLEDGMENTS

This work has been supported in part by ANR through COSINUS program (project PETALh no. ANR-10-COSI-0013 and project HAMM no. ANR-10-COSI-0009). It was granted access to the HPC resources of TGCC@CEA made available within the Distributed European Computing Initiative by the PRACE-2IP, receiving funding from the European Community's Seventh Framework Programme (FP7/2007-2013) under grant agreement RI-283493. The first author would like to thank A. Gupta and O. Schenk for providing him suitable licenses of their respective software (WSMP and PARDISO).


Etienne Hendrickx
email: etienne.hendrickx@univ-brest.fr

Mathieu Paquier

Vincent Koehl

SHOULD A MOVIE HAVE TWO DIFFERENT SOUNDTRACKS FOR ITS STEREOSCOPIC AND NON-STEREOSCOPIC VERSIONS? A STUDY ON THE FRONT/REAR BALANCE

Keywords: Cinema, stereoscopy, sound mixing, balance, surround, 3D

Few psychoacoustic studies have been made on the influence of stereoscopy on the sound mixing of movies. Yet very different opinions can be found among scientific, esthetical or technical communities. Some argue that sound needs to be mixed differently for stereoscopic movies, whereas others pretend that image has actually caught up with sound, that was already "three-dimensional" and should not therefore be affected by stereoscopy. In the present experiment, expert subjects were asked to achieve surround sound ambiance mixings for eleven short sequences presented in both stereoscopic and nonstereoscopic versions. The results suggest that the influence of stereoscopy on the front/rear balance strongly depends on the content of the sequence and only appears in a few specific situations.

INTRODUCTION 1.Technical Considerations

In Hollywood, very different approaches can be found: for some sound engineers, stereoscopic images completely change our auditory perception and two different soundtracks are to be mixed for the 2D and 3D versions of a movie. Others argue that this influence is weak, even negligible.

Michael Semanick1 affirms that there were many differences between the stereoscopic and the nonstereoscopic mixes of Tim Burton's "Alice in Wonderland": For the 3D version, they added some surround effects in areas and pulled music out into surround a bit more. They also pushed the reverberation and the surrounds of the backgrounds further, and panned a bit more the dialog [START_REF] Gambier | Nouvelle approche pour la bande sonore d'un film en relief[END_REF].

Paul Martin Smith, editor of Eric Brevig's "Journey to the center of the earth", also defends the idea that the sound mixer should always work with the image projected in stereoscopy, because it influences the way one places sound sources in the space [START_REF] Krohn | Entretien avec Paul Martin Smith[END_REF].

On the other hand, the mixing crew of Martin Scorsese's "Hugo Cabret" considers the 3D actual revival as a phenomenon concerning image rather than sound [START_REF] Coleman | The sound of Hugo[END_REF]. They argue that sound was already "in 3D with the surround speakers" and should not be too much affected by stereoscopy. James Cameron and its sound crew seem to share the same opinion. On "Avatar", they worked with the 2D version, and occasionally checked if their mixing worked with the stereoscopic version. According to them, they hardly made any modifications [START_REF] Coleman | The sound of Avatar[END_REF].

Scientific Considerations

The few studies that have been made so far on sound related to stereoscopic images have mainly focused on physical realism. They rely on the hypothesis that a sound reproduction closer to a real sound field (such as Wave Field Synthesis, Ambisonics, or binaural techniques) should increase even more the feeling of being part of the movie, in the same way that stereoscopic images increase the sense of presence of the audience, as shown by Ijsselsteijn et al. in [START_REF] Ijsselsteijn | Effects of stereoscopic presentation, image motion, and screen size on subjective and objective corroborative measures of presence[END_REF]. However, it has not been proved yet that stereophonic images increased the sense of presence because it was closer to human perception, and besides several studies suggest that a physically realistic reproduction of sound may not be appropriate for movies: André et al. [START_REF] André | Sound for 3D cinema and the sense of presence[END_REF] played back a stereoscopic movie with different sound conditions, from stereo, that has a low spatial audiovisual coherence, to WFS, that has a very high spatial audiovisual coherence. Only 12 out of 33 participants felt that the sound condition had had an impact on their sense of presence, and they reported the WFS soundtrack as providing the lower sense of presence, when it was supposed to be the more physically realistic. However, only one sequence was used for the study (the first three scenes of an animation stereoscopic movie) with a sound reproduction system that was not 360° but only frontal. Further psychological studies are therefore needed to confirm that a physically realistic reproduction of sound can truly increase or decrease the sense of presence [START_REF] André | Adding 3D sound to 3D cinema: Identification and evaluation of different reproduction techniques[END_REF]. Moreover, a recent study has shown that too much spatial information could reduce attention to narrated text [START_REF] Schmidt | Interactive 3D audio: enhancing awareness of details in immersive soundscapes? 133 rd AES Convention[END_REF]. Even if this study was carried out without images, it may be a clue to explain the results of André et al..

Esthetical Considerations

For Chion [START_REF] Michel Chion. L'audiovision | [END_REF], people are so accustomed to the conventions of cinema that they have become new references of reality. Rumsey [START_REF] Rumsey | Spatial Quality Evaluation for Reproduced Sound : Terminology, Meaning, and a Scene-Based Paradigm[END_REF] also hypothesizes that reproduced sound may have acquired its own standard of realism, different from natural listening.

The ventriloquism effect is a good example to support this idea: in theaters, the voices are most of the time reproduced on the central speaker. Yet the voices seem to match with the position of the actors, even when they are at the borders of the screen or off screen. Sometimes, breaking this "rule" and panning the dialog can even be perceived in a negative way [START_REF] Chion | Un art sonore, le cinema[END_REF].

However, It may sometimes be interesting to draw inspiration from reality and try to adapt some physical aspects of a real sound field to fit in a "traditional" reproduction system, like the 5.1 configuration. For example, Chion [START_REF] Michel Chion. L'audiovision | [END_REF] suggests that surround gives more credit to ambiance sounds because it is sensorily more convincing, as it reproduces the fact that sound comes not only from front but also from behind.

Nevertheless, the goal of this study is not to determine whether a more physically realistic reproduction of sound can increase the sense of presence, but rather to study the differences of auditory perception that may occur when watching a movie in 2D or in 3D. The present study focuses on the influence that stereoscopy may have on the mixing of backgrounds. Subjects were asked to adjust the front/rear balance for eleven stereoscopic sequences, along with their non-stereoscopic version. If the results show that background surrounds are significantly pushed further when watching the sequences in 3D, it could mean that a stereoscopic image, as it is closer to human perception, motivates the subject to have a sound reproduction closer to reality as well, with no discrimination between the front and the rear sounds (as opposed to the "traditional" sound in 2D cinema, in which most of the energy comes from the front speakers [START_REF] Jullier | Le son au cinéma[END_REF]). Either way, it would suggest that new "conventions" (or "standards of realism") are needed, and that two different soundtracks have to be mixed for movies that are to be either projected in 2D or in 3D.

EXPERIMENTAL SETUP

Material

Eleven stereoscopic sequences were used for the test, along with their non-stereoscopic version (see Tab. 1).

Seven sequences were specially shot for the test, while the remaining four were taken from "Tonnerre de Brest", a 3D documentary by Pierre Souchar. All the sequences were shot with a Panasonic AG-3DP1 camera, and were chosen so that they would offer a variety of dynamics (static shots, dollying in, hand-held shots, etc.), types (close shots, medium shots, long shots, etc.) and categories (sea, city, crowd, interior, forest).

Also Guichardan had noticed that surrounds could be more or less efficient depending on the dramatic content of the stimulus used [START_REF] Guichardan | Interactions between surround level and the immersive feeling in the multichannel movie experience. 111 th AES Convention[END_REF]. Though his experience only concerned 2D images, it was nonetheless decided to vary the "degree of dramaturgy" of our sequences, from sequences in which basically "nothing happens" (long shot of a quiet sea) to dialog or music scenes.

Recordings

It was also decided to use different recording and mixing techniques that are commonly used for professional shootings (see Tab.1).

3 sequences were recorded using Double-M/S (see Fig. 1), one of the most established surround recording techniques for applications such as documentary sound and radio drama [START_REF] Wittek | Double M/S -a surround technique put to test[END_REF]. It consists of two cardioid microphones (one facing forward, and one facing backward) and one bidirectionnal microphone (angled 90°), whose signals are then matriced to generate multi-channel sound tracks. Fig. 1. Double-M/S system. The three microphones are supposed to be coincident. 4 sequences were mixed using decorrelated stereophonic backgrounds (one reproduced on the front speakers and the other reproduced on the rear speakers), which is a very common way of building ambiance sounds. Those backgrounds were recorded with an ORTF pair, which consists of two cardioid microphones spaced with 17 cm in a base angle of 110°, and which has been proved to yield a fine distance discrimination compared to other traditional arrays such as XY or M/S [START_REF] Hugonnet | Comparative spatial transfer function of six different stereophonic systems[END_REF].

2 sequences were recorded using Double-ORTF (see Fig. 2), a system with four cardioids arranged in a perfect square shape, with an angle of 110 degrees between the two front microphones and between the two rear microphones. Many tests have pointed out the efficiency of this system for surround sound [START_REF] Kornacki | Problems Related to Surround Sound Production[END_REF].

• Fig. 2. Double-ORTF layout (from [START_REF] Czyżewski | Some rules and methods for Creation of Surround Sound[END_REF]).

All those recording setups were implemented using Schoeps CMC6 microphones with MK4 cardioid directivity capsules (except for the Double M/S system, which also included a Schoeps CMC6 microphone with a MK8 figureeight directivity capsule), connected to a Sonosax SX-R4 audio recorder. A Neumann KMR 81 shotgun microphone was also used for the monophonic recording of dialogs and footsteps from sequences 4, 5, and 6.

At last, the two sequences in the church were recorded using a Fukada Tree. It is a non-coincident multichannel array that was the preferred system in the comparative study of Kassier, probably due to its pleasant spatial impression [START_REF] Kassier | An informal comparison between surround-sound microphone techniques[END_REF]. Several configurations exist for the Fukada Tree, and it was decided to use the same array as Hiekkanen et al. in [START_REF] Hiekkanen | Reproduction of virtual reality with multichannel microphone techniques[END_REF] (see Fig. 3), which consists of three cardioid microphones forming a triangle for the front channels, and of two cardioids facing backward for the rear channels.

Fig. 3. Fukada Tree layout (from [19]).

The tree was implemented using five Schoeps CMC6 microphones with MK4 cardioid directivity capsules, connected to the inputs of a RME Fireface 800 interface.

Reproduction Setup

The impact of stereoscopy was studied on a "traditional" 5.1 mix. The listening test took place in the mixing auditorium of the Image & Sound department of the University of Brest. This room is especially designed for musical mixing and film postproduction. Five professional monitoring loudspeakers (PSI Audio 25-3) were arranged in 3/2 stereo configuration. The loudspeakers were fed by a RME Fireface 800 interface connected to a MacBook Pro computer, and the image was projected by an Epson EH-TW6000 projector, with Epson ELPGS01 3D active glasses.

All the sequences had been previously edited and mixed in that same auditorium by the experimenters. For each sequence, a front/rear balance was fixed by the experimenters, one for the stereoscopic version, and one for the non-stereoscopic version. The average of the two balances was then calculated to define for each sequence an original front/rear balance (nominally 0 dB) that would be retained as a reference for later analysis. This "0 dBbalance" did not necessarily mean having as much energy coming from the front speakers than coming from the rear speakers, which would have been irrelevant for some sequences such as sequences 1 and 2, in which there is only reverberation in the rear speakers.

The gains of the loudspeakers were adjusted so that they would individually produce at the listening position a sound pressure level of 85 dB, C-weighted, when playing a pink noise at -20 dBFS RMS. This calibration is in accordance with the Dolby recommendations and was maintained during the whole process, from the original mixings to the final listening tests. For each sequence, a global level was set subjectively by the experimenters themselves, as it is often the case in subjective tests [START_REF]AES recommended practice for professional audiosubjective evaluation of loudspeakers[END_REF].

Subjects and protocol

The listeners were all paid volunteers from the Image & Sound course at the University of Brest. They were master's degree students, which means they already had a strong experience in critical listening and mixing, and could therefore be considered as "experts" [START_REF]Sensory analysis -General guidance for the selection, training and monitoring of assessors -Part 2: Experts[END_REF].

The eleven stereoscopic sequences, along with their nonstereoscopic version, were presented to all the subjects. The order in which the 22 stimuli occured was random and different for each subject. They were asked to keep their 3D glasses at all time, even for the 2D sequences (during which the same image was sent to the left eye and to the right eye), to avoid an influence of the loss of brightness (that can go from 40 to 70% with active glasses). Each subject completed the test twice, with a fifteen-minute break between the two sessions. The average length of a session was about forty minutes. No subjects reported that they had experienced visual fatigue or discomfort due to prolonged 3D viewing, nor that the test had been too long or too demanding.

They were asked the question: "What front/rear balance would you like for this sequence?". They had to answer the question by setting themselves their proper balance using a knob from a Digidesign Command-8 controller, a MIDI interface that could provide 128 different steps for the balance. The knob was "infinite" and the gradation around it had been hidden, so that the listener did not have any visual or tactile feedback that could have influenced him. By turning the knob clockwise, the subject increased the level of front speakers while decreasing the level of surround speakers, until having only sound coming from front. A simple "sine-cosine" pan law was chosen for the evolution of intensity, which sounded more natural and closer to the subjects's everyday experience of mixing [START_REF] Griesinger | Stereo and Surround Panning in Practice[END_REF].

Let n be the MIDI value set by the subject, and GF and GR be respectively the amplification gain (expressed in dB) for the front and the rear speakers: For example, if the subject set the knob at its middle course (n = 64), then . Neither the level of the front speakers nor the level of the rear speakers were modified, which meant the subject chose the balance as it was initially mixed by the experimenters. If the subjects turned the knob to its maximal position (n = 127), then and , with sound only coming from front. However, to avoid dealing with infinite numbers and forbidden values, the subject could only set n between 1 and 126.

Each sequence was about 30 seconds long and automatically looped. When the subject was satisfied with his balance, he had to press a "push" button to go to the next sequence. Each sequence was initially presented with a random front/rear balance. Playback and data capture from the knob were controlled by a software implemented in Max/MSP.

For sequences 4, 5, and 6, the subjects could not modify the level of the dialogs and the footsteps that were reproduced on the central speaker.

RESULTS

Statistical analysis was performed, using ΔG as variable, to determine if there were any significant differences of balance between "stereoscopic mix" and "non-stereoscopic mix". As the distributions were not normal, a nonparametric statistical test appropriate to repeated measures had to be used: the Wilcoxon signed-rank test [START_REF] Hollander | Nonparametric Statistical Methods[END_REF].

When analyzing the first session of the test (see Fig. 4), only the tenth extract (long shot of the bow of a ship, with two sailors maneuvering) was shown as having a significant difference between its stereoscopic and non-stereoscopic mixes (p = 0,0322), with medians equal to + 5.4 dB for the stereoscopic mix and + 7.0 dB for the non-stereoscopic mix. However, the difference was not significant anymore (p = 0.0645) if one did not take into account the results of subject 3 (+ 38.2 dB for the stereoscopic mix, which corresponds to the maximal value the knob can take, and + 8.9 dB for the non-stereoscopic mix). Besides, the results of subject 3 were far less extreme in the second session (+ 4.9 dB for the stereoscopic mix and -7.6 dB for the non-stereoscopic mix).

Once the first session was finished, the subject would take a fifteen-minute break and would then complete the test a second time (with the sequences in a different order). The analysis of the second session showed significant differences for three sequences (see Fig. 5), with the nonstereoscopic mixes always more frontal than the stereoscopic mixes. Sequence 2 (close shot on a cellist speaking in a church) was significantly different (p = 0,0098), with medians equal to -1.2 dB for the 3D version and + 4.1 dB for the 2D version. Sequence 4 (two shot of a dialog in a café, with dialogs reproduced on the central speaker) was significantly different (p = 0,0488), with medians equal to + 7.6 dB for the 3D version and + 11.1 dB for the 2D version. At last, sequence 9 (long shot of a quiet sea, from the bow of a ship) was significantly different (p = 0,0186), with medians equal to + 0.8 dB for the 3D version and + 4.1 dB for the 2D version. All the other sequences were not significantly different, including sequence 10 (p = 0,1309).

DISCUSSION

When there are significant differences, the 2D mix is always more frontal than the 3D mix, as supported by some sound engineers (see section 1.1.). However, it only concerns a few sequences: three significant differences for the second session, and one questionable difference for the first session.

Many subjects reported that they had trouble mixing complex sequences, such as sequence 7 (bridge with a lot of cars and a tram passing by) or sequence 11 (large crowd in a street). Throughout the two sessions, they would always notice new auditory or visual objects and would try to take them into account in their balance. Sequence 1 (the tracking shot in the church) was also problematic for the subjects, as they had trouble setting a balance that worked for the entire sequence. Most of them would have preferred to mix it dynamically, having more and more sound coming from front as the shot gets closer to the cellist. On the other hand, the three significant sequences of session 2 were very simple scenes: they were static shots, with few sound or visual cues, and very limited movements. More significant differences were found in the second session than in the first one. It suggests that the subjects may have to go through a phase of learning in order to give relevant results. Besides, many subjects reported that it took them several sequences, sometimes an entire session, to realize that some sequences were actually in 2D (they were probably fooled by the 3D glasses that they were wearing permanently). More sessions are therefore needed to assert the results of this study.

CONCLUSION

The results suggest that the influence of stereoscopy on the front/rear balance strongly depends on the content of the sequence and only appears in a few specific situations. This is in accordance with the statements of the sound engineers of "Avatar", who claimed they had hardly made any modifications for the stereoscopic mix [START_REF] Coleman | The sound of Avatar[END_REF]. Stereoscopy tended to have an influence for simple sequences. For complex scenes, the number of parameters taken into account by the subjects for their balances became too important and significantly reduced the influence of the visual rendering.

Tab. 1 .

 1 Bach's third cello suite in a very reverberant church Dollying in Fukada Tree 2 Cellist speaking, still in the same church Close shot Fukada Tree 3 Interior of a childcare center, with backgrounds of children playing from outside Long shot Decorrelated Stereo 4 Dialog in a café, with dialogs reproduced on the central speaker Two shot Decorrelated Stereo 5 Same shot that 4) but with the dialogs pushed further in the central speaker (+ 6 dB) Two shot Decorrelated Stereo 6 Two girls walking in a forest, shot from behind, with their footsteps reproduced on the central speaker The eleven sequences, with their content, shot type, recording system and snapshot.

Fig. 4 .

 4 Fig. 4. Boxplots of the stereoscopic vs. non-stereoscopic balances of session 1, for the eleven sequences.

Fig. 5 .

 5 Fig. 5. Boxplots of the stereoscopic vs. non-stereoscopic balances of session 2, for the eleven sequences.

Two-time Academy Award winner sound mixer for "The Lord of the Rings: The Return of the King" (2003) and "King Kong" (2005).

ACKNOWLEDGMENTS

The authors would like to thank all the many contributors to this work, such as director Pierre Souchar (3D Fovéa), Erwan le Morvan (University of Western Brittany) and the students from ISB who took part in the subjective experiments. This work was funded by the European crossborder cooperation programme INTERREG IV A France (Channel) -England, co-funded by the ERDF, in the context of the Cross Channel Film Lab project.