Xavier Lacoste

Mathieu Faverge
email: mathieu.faverge@inria.fr

Pierre Ramet

Samuel Thibault
email: samuel.thibault@inria.fr

George Bosilca
email: bosilca@icl.utk.edu

Taking advantage of hybrid systems for sparse direct solvers via task-based runtimes

The ongoing hardware evolution exhibits an escalation in the number, as well as in the heterogeneity, of computing resources. The pressure to maintain reasonable levels of performance and portability forces application developers to leave the traditional programming paradigms and explore alternative solutions. PASTIX is a parallel sparse direct solver, based on a dynamic scheduler for modern hierarchical manycore architectures. In this paper, we study the benefits and limits of replacing the highly specialized internal scheduler of the PASTIX solver with two generic runtime systems: PARSEC and STARPU. The tasks graph of the factorization step is made available to the two runtimes, providing them the opportunity to process and optimize its traversal in order to maximize the algorithm efficiency for the targeted hardware platform. A comparative study of the performance of the PASTIX solver on top of its native internal scheduler, PARSEC, and STARPU frameworks, on different execution environments, is performed. The analysis highlights that these generic task-based runtimes achieve comparable results to the application-optimized embedded scheduler on homogeneous platforms. Furthermore, they are able to significantly speed up the solver on heterogeneous environments by taking advantage of the accelerators while hiding the complexity of their efficient manipulation from the programmer.

I. INTRODUCTION

Emerging processor technologies put an emphasis on increasing the number of computing units instead of increasing their working frequencies. As a direct outcome of the physical multiplexing of hardware resources, complex memory hierarchies had to be instated to relax the memory bottleneck and ensure a decent rate of memory bandwidth for each resource. The memory becomes divided in several independent areas, capable of delivering data simultaneously through a complex and hierarchical topology, leading to the mainstream Non Uniform Memory Accesses (NUMA) machines we know today. Together with the availability of hardware accelerators, this trend profoundly altered the execution model of current and future platforms, progressing them toward a scale and a complexity unattained before. Furthermore, with the established integration of accelerators into modern architectures, such as GPUs or Intel Xeon Phis, high-end multi-core CPUs are consistently outperformed by these novel, more integrated, architectures both in terms of data processing rate and memory bandwidth. As a consequence, the working environment of today's application developers evolved toward a multi-level massively parallel environment, where computation becomes cheap but data movements expensive, driving sup the energetic cost and algorithmic overheads and complexities.

With the advent of APIs for GPU programming, such as CUDA or OpenCL, programming accelerators has been rapidly evolving in the past years, permanently bringing accelerators into the mainstream. Hence, GPUs are becoming a more and more attractive alternative to traditional CPUs, particularly for their more interesting cost-per-flop and watts-per-flop ratios. However, the availability of a particular programming API only partially addresses the development of hybrid algorithms capable of taking advantage of all computational resources available, including accelerators and CPUs. Extracting a satisfactory level of performance, out of such entangled architectures, remains a real challenge due to the lack of consistent programming models and tools to assess their performance. In order to efficiently exploit current and future architectures, algorithm developers are required to expose a large amount of parallelism, adapt their code to new architectures with different programming models, and finally, map it efficiently on the heterogeneous resources. This is a gargantuan task for most developers as they do not possess the architectural knowledge necessary to mold their algorithms on the hardware capabilities in order to achieve good efficiency, and/or do not want to spend new efforts with every new generation of hardware.

Solving large sparse linear systems of equations, Ax = b, is one of the most important and time-consuming parts in many scientific and engineering algorithms, a building block toward more complex scientific applications. A significant amount of research has been done on dense linear algebra, but sparse linear algebra on heterogeneous system is a work-in-progress. Multiple reasons warrant this divergence, including the intrinsic algorithmic complexity and the highly irregular nature of the resulting problem, both in terms of memory accesses and computational intensities. Combined with the heterogeneous features of current and future parallel architectures, this depicts an extremely complex software development field.

The PASTIX solver is a sparse direct solver that can solve symmetric definite, indefinite, and general problems using Cholesky, LDL T , and LU factorizations, respectively. The PASTIX implementation relies on a two-level approach using the POSIX Thread library within a node, and the Message Passing Interface (MPI) between different nodes. Historically, PASTIX scheduling strategy was based on a cost model of the tasks executed that defines the execution order used at runtime during the analyze phase. In order to complement the lack of precision of the cost model on hierarchical architectures, a dynamic scheduler based on a work-stealing strategy has been developed to reduce the idle times while preserving a good locality for data mapping [START_REF] Faverge | Fine grain scheduling for sparse solver on manycore architectures[END_REF]. More recently, the solver has been optimized to deal with the new hierarchical multi-core architectures [START_REF]Dynamic Scheduling for sparse direct Solver on NUMA architectures[END_REF], at the level of internal data structures of the solver, communication patterns, and scheduling strategies.

In this paper, we advance the state-of-the-art in supernodal solvers by migrating PASTIX toward a new programming paradigm, one with a promise of efficiently handling hybrid execution environments while abstracting the application from the hardware constraints. Many challenges had to be overcome, going from exposing the PASTIX algorithms using a task-based programming paradigm, to delivering a level of task parallelism, granularity, and implementation allowing the runtime to efficiently schedule the resulting, highly irregular tasks, in a way that minimizes the execution span. We exposed the original algorithm using the concept of tasks, a self-contained computational entity, linked to the other tasks by data dependencies. Specialized taskgraph description formats were used in accordance with the underlying runtime system (PARSEC or STARPU). We provided specialized GPU-aware versions for some of the most compute intensive tasks, providing the runtimes with the opportunity to unroll the graph of tasks on all available computing resources. The resulting software is, to the best of our knowledge, the first implementation of a sparse direct solver with a supernodal method supporting hybrid execution environments composed of multi-cores and multi-GPUs. Based on these elements, we pursue the evaluation of the usability and the appeal of using a taskbased runtime as a substrate for executing this particular type of algorithm, an extremely computationally challenging sparse direct solver. Furthermore, we take advantage of the integration of accelerators (GPUs in this context) with our supporting runtimes, to evaluate and understand the impact of this drastically novel portable way of writing efficient and perennial algorithms. Since the runtime system offers a uniform programming interface, dissociated from a specific set of hardware or low-level software entities, applications can take advantage of these uniform programming interfaces for ensuring their portability. Moreover, the exposed graph of tasks allows the runtime system to apply specialized optimization techniques and minimize the application's time to solution by strategically mapping the tasks onto computing resources by using state-of-the-art scheduling strategies.

The rest of the paper is organized as follows. We describe the supernodal method of the PASTIX solver in Section III, followed by a description of the runtimes used IV. Section V explains the implementation over the DAG schedulers with a preliminary study over multi-core architectures, followed by details on the extension to heterogeneous architectures. All choices are supported and validated by a set of experiments on a set of matrices with a wide range of characteristics. Finally, section VI concludes with some prospects of the current work.

II. RELATED WORK

The dense linear algebra community spent a great deal of effort to tackle the challenges raised by the sharp increase of the number of computational resources. Due to their heavy computational cost, most of their algorithms are relatively simple to handle. Avoiding common pitfalls such as the "fork-join" parallelism, and expertly selecting the blocking factor, provide an almost straightforward way to increase the parallelism and thus achieve better performance. Moreover, due to their natural load-balance, most of the algorithms can be approached hierarchically, first at the node level, and then at the computational resource level. In a shared memory context, one of the seminal papers [START_REF] Buttari | The impact of multicore on math software[END_REF] replaced the commonly used LAPACK layout with one based on tiles/blocks. Using basic operations on these tiles exposes the algorithms as a graph of tasks augmented with dependencies between them. In shared memory, this approach quickly generates a large number of ready tasks, while, in distributed memory, the dependencies allow the removal of hard synchronizations. This idea leads to the design of new algorithms for various algebraic operations [START_REF] Buttari | A class of parallel tiled linear algebra algorithms for multicore architectures[END_REF], now at the base of well-known software packages like PLASMA [START_REF] Agullo | Numerical linear algebra on emerging architectures: The PLASMA and MAGMA projects[END_REF].

This idea is recurrent in almost all novel approaches surrounding the many-core revolution, spreading outside the boundaries of dense linear algebra. Looking at the sparse linear algebra, the efforts were directed toward improving the behavior of the existing solvers by taking into account both task and data affinity and relying on a two-level hybrid parallelization approach, mixing multithreading and message passing. Numerous solvers are now able to efficiently exploit the capabilities of these new platforms [START_REF]Dynamic Scheduling for sparse direct Solver on NUMA architectures[END_REF], [START_REF] Schenk | Solving unsymmetric sparse systems of linear equations with pardiso[END_REF]. New solvers have also been designed and implemented to address these new platforms. For them the chosen approach follows the one for dense linear algebra, fine-grained parallelism, threadbased parallelization, and advanced data management to deal with complex memory hierarchies. Examples of this kind of solver are HSL [START_REF] Hogg | Design of a multicore sparse cholesky factorization using dags[END_REF] and SuperLU-MT [START_REF] Li | Evaluation of sparse LU factorization and triangular solution on multicore platforms[END_REF] for sparse LU or Cholesky factorizations and SPQR [START_REF] Davis | Algorithm 915, SuiteSparseQR: Multifrontal multithreaded rank-revealing sparse QR factorization[END_REF] and qr_mumps [START_REF] Buttari | Fine-grained multithreading for the multifrontal QR factorization of sparse matrices[END_REF] for sparse QR factorizations.

With the advent of accelerator-based platforms, a lot of attention has shifted toward extending the multi-core aware algorithms to fully exploit the huge potential of accelerators (mostly GPUs). The main challenges raised by these heterogeneous platforms are mostly related to task granularity and data management: although regular cores require fine granularity of data as well as computations, accelerators such as GPUs need coarse-grain tasks. This inevitably introduces the need for identifying the parts of the algorithm which are more suitable to be processed by accelerators. As for the multi-core case described above, the exploitation of this kind of platform was first considered in the context of dense linear algebra algorithms.

Moreover, one constant becomes clear: a need for a portable layer that will insulate the algorithms and their developers from the rapid hardware changes. Recently, this portability layer appeared under the denomination of a taskbased runtime. The algorithms are described as tasks with data dependencies in-between, and the runtime systems are used to manage the tasks dynamically and schedule them on all available resources. These runtimes can be generic, like the two runtimes used in the context of this study (STARPU [START_REF] Augonnet | StarPU: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF] or PARSEC [START_REF] Bosilca | DAGuE: A generic distributed DAG engine for High Performance Computing[END_REF]), or more specialized like QUARK [START_REF] Yarkhan | Dynamic task execution on shared and distributed memory architectures[END_REF]. These efforts resulted in the design of the DPLASMA library [START_REF] Bosilca | Flexible development of dense linear algebra algorithms on massively parallel architectures with DPLASMA[END_REF] on top of PARSEC and the adaptation of the existing FLAME library [START_REF] Igual | The FLAME approach: From dense linear algebra algorithms to highperformance multi-accelerator implementations[END_REF]. On the sparse direct methods front, preliminary work has resulted in mono-GPU implementations based on offloading parts of the computations to the GPU [START_REF] George | Multifrontal Factorization of Sparse SPD Matrices on GPUs[END_REF]- [START_REF] Yu | A CPU-GPU Hybrid Approach for the Unsymmetric Multifrontal Method[END_REF]. Due to its very good data locality, the multifrontal method is the main target of these approaches. The main idea is to treat some parts of the task dependency graph entirely on the GPU. Therefore, the main originality of these efforts is in the methods and algorithms used to decide whether or not a task can be processed on a GPU. In most cases, this was achieved through a threshold based criterion on the size of the computational tasks.

Many initiatives have emerged in previous years to develop efficient runtime systems for modern heterogeneous platforms. Most of these runtime systems use a task-based paradigm to express concurrency and dependencies by employing a task dependency graph to represent the application to be executed. Without going into details, the main differences between these approaches are related to their representation of the graph of tasks, whether they manage data movements between computational resources, the extent they focus on task scheduling, and their capabilities to handle distributed platforms.

III. SUPERNODAL FACTORIZATION

Sparse direct solvers are algorithms that address sparse matrices, mostly filled with zeroes. In order to reduce the Figure 1: Decomposition of the task applied while processing one panel number of operations, they consider only non-zeroes of the matrix A. During factorization, new non-zero entries -called fill-in -appear in the factorized matrix and lead to more computation and memory consumption. One of the main objectives of those solvers is to keep the fill-in to its minimum to limit the memory consumption. The first step of any sparse direct solver is the computation of a nested dissection of the problem that results in a permutation of the unknowns of A. This process is a graph partitioning algorithm applied to the connectivity graph associated with the matrix. The computed permutation reduces the fill-in that the factorization process will generate, and the elimination tree [START_REF] Liu | The role of elimination trees in sparse factorization[END_REF] is generated out of the separators discovered during the graph partitioning process. Basically, each node of the tree represents the set of unknowns that belongs to a separator, and edges are connections between those separators in the original graph. The edges of the tree connect a node to almost all nodes in the path that connect it to the root of the tree. They represent contributions from one node to another during the factorization. The second step of sparse solvers is the analysis stage which predicts the non-zero pattern of the factorized matrix through a symbolic factorization. The resulting pattern is grouped in blocks of non-zero elements to allow for more efficient BLAS calls. Blocks can be enlarged, if extra fill-in is allowed, for better performance, or split to generate more parallelism.

Once the elimination tree and the symbolic factorization are computed, two different methods can be applied: multifrontal [START_REF] Duff | The multifrontal solution of indefinite sparse symmetric linear[END_REF] or supernodal [START_REF] Ashcraft | Progress in sparse matrix methods for large linear systems on vector supercomputers[END_REF]. The PASTIX solver uses a supernodal method. Each node of the elimination tree, or supernode, represents a subset of contiguous columns, also called a panel, in the matrix. To each node of the elimination tree, we associate a task, called 1D task, that performs three steps associated with the panel A, as shown in the Figure 1: 1) Factorization of the diagonal block, 2) Triangular solve on the off-diagonal blocks in the panel (TRSM), and

3) For each off-diagonal block A i , apply the associated update to the facing panel C i (GEMM) -we call C i the facing panel with the diagonal block owning the same rows as the off-diagonal block A i . Figure 1 represents a lower triangular matrix used in case of symmetric problems with Cholesky factorization, but is also valid for non-symmetric cases with PASTIX. In a general manner, PASTIX works on the matrix A+A T , which produces a symmetric pattern. In non-symmetric cases, the steps 2 and 3 are then duplicated for the L and U matrices of the LU factorization. Besides, the PASTIX solver doesn't perform dynamic pivoting, as opposed to SuperLU [START_REF] Demmel | A supernodal approach to sparse partial pivoting[END_REF], which allows the factorized matrix structure to be fully known at the analysis step. In the rest of the paper, we will discuss only the Cholesky implementation. LDL T and LU factorizations follow the same method.

The PASTIX solver relies on a cost model of this 1D task to compute a static scheduling. This static scheduling associates ready tasks with the first available resources among the computational units. The complexity of such an algorithm depends on the number of tasks and resources. Hence, the 1D task is kept as a large single task to lower the complexity of the analysis part. However, it is obvious that more parallelism could be extracted from those tasks, but would increase the analysis step complexity.

First, the triangular solves, applied on off-diagonal blocks of each panel, are independent computations that depend only on the diagonal block factorization. Thus, each panel is stored as a single tall and skinny matrix, such that the TRSM granularity can be decided at runtime and is independent of the data storage. At lower levels of the elimination tree, the small block granularity might induce a large overhead if they are considered as independent tasks. On the contrary, at higher levels, the larger supernodes (Order of N 2 3 for a 3D problem of N unknowns, or sqrt(N ) for a 2D problem) might be split to create more parallelism with low overhead. That is why supernodes of the higher levels are split vertically prior to the factorization to limit the task granularity and create more parallelism. In this study, to compare to the existing PASTIX solver, we keep all TRSM operations on a single panel grouped together as a single operation.

Second, the same remark as before applies to the update tasks with a higher order of magnitude as before. Each couple of off-diagonal blocks (A i , A j ), with i < j in a panel, generates an independent update to the trailing submatrix formed by their outer product. To adapt to the small granularity of off-diagonal blocks in sparse solvers, those updates are grouped together. Two variants exists. Leftlooking: all tasks contributing to a single panel are associated in a single task, they have a lot of input edges and only one in-out data. Right-looking: all updates generated by a single panel are directly applied to the multiple destination panels. This solution has a single input data, and many panels are accessed as in-out. PASTIX uses the right-looking variant. Therefore, the nature of the supernodal algorithm is in itself a DAG of tasks dependent on the structure of the factorized matrix, but independent of the numerical content thanks to the static pivoting strategy. However, since we consider a data as a panel, and both of the targeted runtime systems take a fixed number of dependencies per tasks, one update task will be generated per couple of panels instead of one per panel, as in PASTIX.

IV. RUNTIMES

In our exploratory approach toward moving to a generic scheduler for PASTIX, we considered two different runtimes: STARPU and PARSEC. Both runtimes have been proven mature enough in the context of dense linear algebra, while providing two orthogonal approaches to task-based systems.

The PARSEC [START_REF] Bosilca | DAGuE: A generic distributed DAG engine for High Performance Computing[END_REF] distributed runtime system is a generic data-flow engine supporting a task-based implementation targeting hybrid systems. Domain specific languages are available to expose a user-friendly interface to developers and allow them to describe their algorithm using highlevel concepts. This programming paradigm constructs an abridged representation of the tasks and their dependencies as a graph of tasks -a structure agnostic to algorithmic subtleties, where all intrinsic knowledge about the complexity of the underlying algorithm is extricated, and the only constraints remaining are annotated dependencies between tasks [START_REF] Cosnard | Compact dag representation and its symbolic scheduling[END_REF]. This symbolic representation, augmented with a specific data distribution, is then mapped on a particular execution environment. The runtime supports the usage of different types of accelerators, GPUs and Intel Xeon Phi, in addition to distributed multi-core processors. Data are transferred between computational resources based on coherence protocols and computational needs, with emphasis on minimizing the unnecessary transfers. The resulting tasks are dynamically scheduled on the available resources following a data reuse policy mixed with different criteria for adaptive scheduling. The entire runtime targets very fine grain tasks (order of magnitude under ten microseconds), with a flexible scheduling and adaptive policies to mitigate the effect of system noise and take advantage of the algorithmic-inherent parallelism to minimize the execution span.

The experiment presented in this paper takes advantage of a specialized domain specific language of PARSEC, designed for affine loops-based programming [START_REF] Bosilca | Flexible development of dense linear algebra algorithms on massively parallel architectures with DPLASMA[END_REF]. This specialized interface allows for a drastic reduction in the memory used by the runtime, as tasks do not exist until they are ready to be executed, and the concise representation of the task-graph allows for an easy and stateless exploration of the graph. In exchange for the memory saving, generating a task requires some extra computation, and lies in the critical path of the algorithm. The need for a window of visible tasks is then pointless, the runtime can explore the graph dynamically based on the ongoing state of the execution. STARPU [START_REF] Augonnet | StarPU: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF] is a runtime system aiming to allow programmers to exploit the computing power of clusters of hybrid systems composed of CPUs and various accelerators (GPUs, Intel Xeon Phi, etc) while relieving them from the need to specially adapt their programs to the target machine and processing units. The STARPU runtime supports a taskbased programming model, where applications submit computational tasks, with CPU and/or accelerator implementations, and STARPU schedules these tasks and associated data transfers on available CPUs and accelerators. The data that a task manipulates is automatically transferred among accelerators and the main memory in an optimized way (minimized data transfers, data prefetch, communication overlapped with computations, etc.), so that programmers are relieved of scheduling issues and technical details associated with these transfers. STARPU takes particular care of scheduling tasks efficiently, by establishing performance models of the tasks through on-line measurements, and then using well-known scheduling algorithms from the literature. In addition, it allows scheduling experts, such as compilers or computational library developers, to implement custom scheduling policies in a portable fashion.

The differences between the two runtimes can be classified into two groups: conceptual and practical differences. At the conceptual level the main differences between PARSEC and STARPU are the tasks submission process, the centralized scheduling, and the data movement strategy. PARSEC uses its own parameterized language to describe the DAG in comparison with the simple sequential submission loops typically used with STARPU. Therefore, STARPU relies on a centralized strategy that analyzes, at runtime, the dependencies between tasks and schedules these tasks on the available resources. On the contrary, through compiletime information, each computational unit of PARSEC immediately releases the dependencies of the completed task solely using the local knowledge of the DAG. At last, while PARSEC uses an opportunistic approach, the STARPU scheduling strategy exploits cost models of the computation and data movements to schedule tasks to the right resource (CPU or GPU) in order to minimize overall execution time. However, it does not have a data-reuse policy on CPUshared memory systems, resulting in lower efficiency when no GPUs are used, compared to the data-reuse heuristic of PARSEC. At the practical level, PARSEC supports multiple streams to manage the CUDA devices, allowing partial overlap between computing tasks, maximizing the occupancy of the GPU. On the other hand, STARPU allows data transfers directly between GPUs without going through central memory, potentially increasing the bandwidth of data transfers when data is needed by multiple GPUs.

V. SUPERNODAL FACTORIZATION OVER DAG

SCHEDULERS

Similarly to dense linear algebra, sparse direct factorization relies on three types of operations: the factorization of the diagonal block (POTRF), the solve on off-diagonal blocks belonging to the same panel (TRSM), and the trailing panels updates (GEMM). Whereas the task dependency graph from a dense Cholesky factorization [START_REF] Buttari | A class of parallel tiled linear algebra algorithms for multicore architectures[END_REF] is extremely regular, the DAG describing the supernodal method contains rather small tasks with variable granularity and less uniform ranges of execution space. This lack of uniformity makes the DAG resulting from a sparse supernodal factorization complex, accruing the difficulty to efficiently schedule the resulting tasks on homogeneous and heterogeneous computing resources.

The current scheduling scheme of PASTIX exploits a 1Dblock distribution, where a task assembles a set of operations together, including the tasks factorizing one panel (POTRF and TRSM) and all updates generated by this factorization. However, increasing the granularity of a task in such a way limits the potential parallelism, and has a growing potential of bounding the efficiency of the algorithm when using many-core architectures. To improve the efficiency of the sparse factorization on a multi-core implementation, we introduced a way of controlling the granularity of the BLAS operations. This functionality dynamically splits update tasks, so that the critical path of the algorithm can be reduced. In this paper, for both the PARSEC and STARPU runtimes, we split PASTIX tasks into two sub-sets of tasks:

• the diagonal block factorization and off-diagonal blocks updates, performed on one panel;

• the updates from off-diagonal blocks of the panel to one other panel of the trailing sub-matrix. Hence, the number of tasks is bound by the number of blocks in the symbolic structure of the factorized matrix.

Moreover, when taking into account heterogeneous architectures in the experiments, a finer control of the granularity of the computational tasks is needed. Some references for benchmarking dense linear algebra kernels are described in [START_REF] Volkov | Benchmarking GPUs to tune dense linear algebra[END_REF] and show that efficiency could be obtained on GPU devices only on relatively large blocks -a limited number of such blocks can be found on a supernodal factorization only on top of the elimination tree. Similarly, the amalgamation algorithm [START_REF] Hénon | On finding approximate supernodes for an efficient ILU(k) factorization[END_REF], reused from the implementation of an incomplete factorization, is a crucial step to obtain larger supernodes and efficiency on GPU devices. The default parameter for amalgamation has been slightly increased to allow up to 12% more fill-in to build larger blocks while maintaining a decent level of parallelism.

In the remaining of the paper, we present the extensions to the solver to support heterogeneous many-core architectures. These extensions were validated through experiments conducted on Mirage nodes from the PLAFRIM cluster at

A. Multi-core Architectures

As mentioned earlier, the PASTIX solver has already been optimized for distributed clusters of NUMA nodes. We use the current state-of-the-art PASTIX scheduler as a basis, and compare the results obtained using the STARPU and PARSEC runtimes from there. Figure 2 reports the results from a strong scaling experiment, where the number of computing resources varies from 1 to 12 cores, and where each group represents a particular matrix. Empty bars correspond to the PASTIX original scheduler, shaded bars correspond to STARPU, and filled bars correspond to PARSEC. The figure is in Flop/s, and a higher value on the Y-axis represents a more efficient implementation. Overall, this experiment shows that on a shared memory architecture the performance obtained with any of the above-mentioned approaches are comparable, the differences remaining minimal on the target architecture.

We can also see that, in most cases, the PARSEC implementation is more efficient than STARPU, especially when the number of cores increases. STARPU shows an overhead on multi-core experiments attributed to its lack of cache reuse policy compared to PARSEC and the PASTIX internal scheduler. A careful observation highlights the fact that both runtimes obtain lower performance compared with PASTIX for LDL T on both PmlDF and Serena matrices. Due to its single task per node scheme, PASTIX stores the DL T matrix in a temporary buffer which allows the update kernels to call a simple GEMM operation. On the contrary, both STARPU and PARSEC implementations are using a less efficient kernel that performs the full LDL T operation at each update. Indeed, due to the extended set of tasks, the life span of the temporary buffer could cause large memory overhead. In conclusion, using these generic runtimes shows similar performance and scalability to the PASTIX internal solution on the majority of test cases, while providing a suitable level of performance and a desirable portability, allowing for a smooth transition toward more complex heterogeneous architectures.

B. Heterogeneous Architectures Implementation

While obtaining an efficient implementation was one of the goals of this experiment, it was not the major one.

The ultimate goal was to develop a portable software environment allowing for an even transition to accelerators, a software platform where the code is factorized as much as possible, and where the human cost of adapting the sparse solver to current and future hierarchical complex heterogeneous architectures remains consistently low. Building upon the efficient supernodal implementation on top of DAG based runtimes, we can more easily exploit heterogeneous architectures. The GEMM updates are the most computeintensive part of the matrix factorization, and it is important that these tasks are offloaded to the GPU. We decide not to offload the tasks that factorize and update the panel to the GPU due to the limited computational load, in direct relationship with the small width of the panels. It is common in dense linear algebra to use the accelerators for the update part of a factorization while the CPUs factorize the panel; so from this perspective our approach is conventional. However, such an approach combined with look-ahead techniques gives really good performance for a low programming effort on the accelerators [START_REF] Yamazaki | One-sided Dense Matrix Factorizations on a Multicore with Multiple GPU Accelerators[END_REF]. The same solution is applied in this study, since the panels are split during the analysis step to fit the classic look-ahead parameters.

It is a known fact that the update is the most compute intensive task during a factorization. Therefore, generally speaking, it is paramount to obtain good efficiency on the update operation in order to ensure a reasonable level of performance for the entire factorization. Due to the embarrassingly parallel architecture of the GPUs and to the extra cost of moving the data back and forth between the main memory and the GPUs, it is of greatest importance to maintain this property on the GPU. As presented in Figure 1, the update task used in the PASTIX solver groups together all outer products that are applied to a same panel. On the CPU side, this GEMM operation is split in two steps due to the gaps in the destination panel: the outer product is computed in a contiguous temporary buffer, and upon completion, the result is dispatched on the destination panel. This solution has been chosen to exploit the performance of vendor provided BLAS libraries in exchange for constant memory overhead per working thread.

For the GPU implementation, the requirements for an efficient kernel are different. First, a GPU has significantly less memory compared with what is available to a traditional processor, usually in the range of 3 to 6 GB. This forces us to carefully restrict the amount of extra memory needed during the update, making the temporary buffer used in the CPU version unsuitable. Second, the uneven nature of sparse irregular matrices might limit the number of active computing units per task. As a result, only a partial number of the available warps on the GPU might be active, leading to a deficient occupancy. Thus, we need the capability to submit multiple concurrent updates in order to provide the GPU driver with the opportunity to overlap warps between different tasks to increase the occupancy, and thus the overall efficiency.

Many CUDA implementations of the dense GEMM kernel are available to the scientific community. The most widespread implementation is provided by Nvidia itself in the CUBLAS library [START_REF] Nvidia | Cublas library[END_REF]. This implementation is extremely efficient since CUDA 4.2 allows for calls on multiple streams, but is not open source. Volkov developed an implementation for the first generation of CUDA enabled devices [START_REF] Volkov | Benchmarking GPUs to tune dense linear algebra[END_REF] in real single precision. In [START_REF] Tan | Fast implementation of dgemm on fermi gpu[END_REF], authors propose an assembly code of the DGEMM kernel that provides a 20% improvement on CUBLAS 3.2 implementation. The MAGMA library proposed a first implementation of the DGEMM kernel [START_REF] Nath | An improved MAGMA GEMM for Fermi graphics processing units[END_REF] for the Nvidia Fermi GPUs. Later, an auto-tuned framework, called ASTRA, was presented in [START_REF] Kurzak | Autotuning GEMM Kernels for the Fermi GPU[END_REF] and included into the MAGMA library. This implementation, similar to the ATLAS library for CPUs, is a highly configurable skeleton with a set of scripts to tune the parameters for each precision.

As our update operation is applied on a sparse representation of the panel and matrices, we cannot exploit an efficient vendor-provided GEMM kernel. We need to develop our own, starting from a dense version and altering the algorithm to fit our needs. Due to the source code availability, the coverage of the four floating point precisions, and it's tuning capabilities, we decided to use the ASTRA-based version for our sparse implementation. As explained in [START_REF] Kurzak | Autotuning GEMM Kernels for the Fermi GPU[END_REF] the matrix-matrix operation is performed in two steps in this kernel. Each block of threads computes the outer-product tmp = AB into the GPU shared memory, and then the addition C = βC + αtmp is computed. To be able to compute directly into C, the result of the update from one panel to another, we extended the kernel to provide the structure of each panel. This allows the kernel to compute the correct position directly into C during the sum step. This introduces a loss in the memory coalescence and deteriorates the update parts, however it prevents the requirement of an extra buffer on the GPU for each offloaded kernel.

One problem in the best parameters used in the MAGMA library for the ASTRA kernel is that it has been determined that using textures gives the best performance for the update kernel. The function cudaBindTexture and cudaUnbindTexture are not compatible with concurrent kernel calls on different streams. Therefore, the textures have been disabled in the kernel, reducing the performance of the kernel by about 5% on large square matrices.

Figure 3 shows the study we made on the GEMM kernel and the impact of the modifications we did on the ASTRA kernel. These experiments are done on a single GPU of the Mirage cluster. The experiments consist of computing a representative matrix-matrix multiplication of what is typically encountered during sparse factorization. Each point is the average performance of 100 calls to the kernel that computes: C = C -AB T , with A, B, and C, matrices respectively of dimension M -by-N , K-by-N , and M -by-N . B is taken as the first block of K rows of A as it is the case in Cholesky factorization. The plain lines are the performance of the CUBLAS library with 1 stream (red), 2 streams (green), and 3 streams (red). The black line represents the peak performance obtained by the CUBLAS library on square matrices. This peak is never reached with the particular configuration case studied here. The dashed lines are the performance of the ASTRA library in the same configurations. We observe that this implementation already looses 50GFlop/s, around 15%, against the CUBLAS library, and that might be caused by the parameters chosen by the auto-tuning framework which has been run only on square matrices. Finally, the dotted lines illustrate the performance of the modified ASTRA kernel to include the gaps into the C matrix. For the experiments, C is a panel twice as tall as A in which blocks are randomly generated with average size of 200 rows. Blocks in A are also randomly generated with the constraint that the rows interval of a block of A is included in the rows interval of one block of C, and no overlap is made between two blocks of A. We observe a direct relationship between the height of the panel and the performance of the kernel: the taller the panel, the lower the performance of the kernel. The memory loaded to do the outer product is still the same as for the ASTRA curves, but memory loaded for the C matrix grows twice as fast without increasing the number of Flop to perform. The ratio Flop per memory access is dropping and explains the decreasing performance. However, when the factorization progresses and moves up the elimination trees, nodes get larger and the real number of blocks encountered is smaller than the one used in this experiment to illustrate worst cases.

Without regard to the kernel choice, it is important to notice how the multiple streams can have a large impact on the average performance of the kernel. For this comparison, the 100 calls made in the experiments are distributed in a round-robin manner over the available streams. One stream always gives the worst performance. Adding a second stream increases the performance of all implementations and especially for small cases when matrices are too small to feed all resources of the GPU. The third one is an improvement for matrices with M smaller than 1000, and is similar to two streams over 1000.

This kernel is the one we provide to both runtimes to offload computations on GPUs in case of Cholesky and LU factorizations. An extension of the kernel has been made to handle the LDL T factorization that takes an extra parameter to the diagonal matrix D and computes: C = C -LDL T . This modified version decreases the performance by 5%.

C. Heterogeneous experiments

Figure 4 presents the performance obtained on our set of matrices on the Mirage platform by enabling the GPUs in addition to all available cores. The PASTIX run is shown as a reference. STARPU runs are empty bars, PARSEC runs with 1 stream are shaded and PARSEC runs with 3 streams are fully colored. This experiment shows that we can efficiently use the additional computational power provided by the GPUs using the generic runtimes. In its current implementation, STARPU has either GPU or CPU worker threads. A GPU worker will execute only GPU tasks. Hence, when a GPU is used, a CPU worker is removed. With PARSEC, no thread is dedicated to a GPU, and they all might execute CPU tasks as well as GPU tasks. The first computational threads that submit a GPU task takes the management of the GPU until no GPU work remains in the pipeline. Both runtimes manage to get similar performance and satisfying scalability over the 3 GPUs. In only two cases, MHD and pmlDF, STARPU outperforms PARSEC results with 3 streams. This experimentation also reveals that, as was expected, the computation takes advantage of the multiple streams that are available through PARSEC. Indeed, the tasks generated by a sparse factorization are rather small and won't use the entire GPU. This PARSEC feature compensates for the prefetch strategy of STARPU that gave it the advantage when compared to the one stream results. One can notice the poor performance obtained on the afshell test case: in this case, the amount of Flop produced is too small to efficiently benefit from the GPUs.

VI. CONCLUSION

In this paper, we have presented a new implementation of a sparse direct solver with a supernodal method using a taskbased programming paradigm. The programming paradigm shift insulates the solver from the underlying hardware. The runtime takes advantage of the parallelism exposed via the graph of tasks to maximize the efficiency on a particular platform, without the involvement of the application developer. In addition to the alteration of the mathematical algorithm to adapt the solver to the task-based concept, and to provide an efficient memory-constraint sparse GEMM for the GPU, contributions to both runtimes (PARSEC and STARPU) were made such that they could efficiently support tasks with irregular duration, and minimize the non-regular data movements to, and from, the devices. While the current status of this development is already significant in itself, the existence of the conceptual task-based algorithm opened an astonishing new perspective for the seamless integration of any type of accelerator. Providing computational kernels adapted to specialized architectures has become the only obstruction to a portable, efficient, and generic sparse direct solver exploiting these devices. In the context of this study, developing efficient and specialized kernels for GPUs allowed a swift integration on hybrid platforms. Globally, our experimental results corroborate the fact that the portability and efficiency of the proposed approach are indeed available, elevating this approach to a suitable programming model for applications on hybrid environments.

Future work will concentrate on smoothing the runtime integration within the solver. First, in order to minimize the scheduler overhead, we plan to increase the granularity of the tasks at the bottom of the elimination tree. Merging leaves or subtrees together yields bigger, more computationally intensive tasks. Second, we will pursue the extension of this work in distributed heterogeneous environments. On such platforms, when a supernode updates another non-local supernode, the update blocks are stored in a local extramemory space (this is called "fan-in" approach [START_REF] Ashcraft | A fanin algorithm for distributed sparse numerical factorization[END_REF]). By locally accumulating the updates until the last updates to the supernode are available, we trade bandwidth for latency. The runtime will allow for studying dynamic algorithms, where the number of local accumulations has bounds discovered at runtime. Finally, the availability of extra computational resources highlights the potential to dynamically build or rebuild the supernodal structures according to the load on the cores and the GPUs.

  a f s h e ll 1 0 ( D , L U ) F il te r V 2 Z , L U ) F la n ( D ,

Figure 3 :

 3 Figure 3: Multi-stream performance comparison on the DGEMM kernel for three implementations: CUBLAS library, ASTRA framework, and the sparse adaptation of the ASTRA framework.

aFigure 4 :

 4 Figure 4: GPU scaling study: GFlop/s performance of the factorization step with the three schedulers on a set of 10 matrices. Experiments exploit twelve CPU cores and from zero to three additional GPUs.

Table I :

 I Matrix description (Z: double complex, D: double).

				PASTIX with 1 core	STARPU with 1 core	PARSEC with 1 core
				PASTIX with 3 cores	STARPU with 3 cores	PARSEC with 3 cores
				PASTIX with 6 cores	STARPU with 6 cores	PARSEC with 6 cores
				PASTIX with 9 cores	STARPU with 9 cores	PARSEC with 9 cores
				PASTIX with 12 cores	STARPU with 12 cores	PARSEC with 12 cores
	Matrix	Prec	Method	Size	nnz A	nnz L	TFlop
	Afshell10	D	LU	1.5e+6	27e+6	610e+6	0.12
	FilterV2	Z	LU	0.6e+6	12e+6	536e+6	3.6
	Flan	D	LL T	1.6e+6	59e+6	1712e+6	5.3
	Audi	D	LL T	0.9e+6	39e+6	1325e+6	6.5
	MHD	D	LU	0.5e+6	24e+6	1133e+6	6.6
	Geo1438	D	LL T	1.4e+6	32e+6	2768e+6	23
	Pmldf	Z	LDL T	1.0e+6	8e+6	1105e+6	28
	Hook	D	LU	1.5e+6	31e+6	4168e+6	35
	Serena	D	LDL T	1.4e+6	32e+6	3365e+6	47

Figure

2

: CPU scaling study: GFlop/s performance of the factorization step on a set of nine matrices with the three schedulers.

INRIA Bordeaux. A Mirage node is equipped with two hexacore Westmere Xeon X5650 (2.67 GHz), 32 GB of memory and 3 Tesla M2070 GPUs. PASTIX was built without MPI support using GCC 4.6.3, CUDA 4.2, Intel MKL 10.2.7.041, and SCOTCH 5.1.12b. Experiments were performed on a set of nine matrices, all part of the University of Florida sparse matrix collection

[START_REF] Davis | University of Florida sparse matrix collection[END_REF]

, and are described in Table

I

. These matrices represent different research fields and exhibit a wide range of properties (size, arithmetic, symmetry, definite problem, etc). The last column reports the number of floating point operations (Flop) required to factorize those matrices and used to compute the performance results shown in this section.

ACKNOWLEDGMENT

The work presented in this paper is supported by the Inria Associated Team MORSE, by the National Science Foundation under Grant Number CCF-1244905, and by the ANR ANEMOS ANR-11-MONU-002. The authors would like to thanks E. Agullo, G. Boutin, and A. Guermouche.


may come

Separability of multiple deep crack defects with an NDE Eddy Current System

Rimond Hamia, Christophe Cordier,


G T Swindles

F De Vleeschouwer

G Plunkett

Dating peat profiles using tephra: stratigraphy, geochemistry and chronology

Keywords: cryptotephra, tephrochronology, volcanic ash layers. _______________________________________________________________________________________

In this article we provide a brief overview of the protocols for dating peat profiles using tephrochronology. A standardised methodology for the detection, extraction and analysis of tephras is presented and the relevant problems and limitations are discussed.

INTRODUCTION

Tephrochronology (dating using volcanic ash layers) is a particularly useful tool for developing robust chronologies from peat sequences (e.g. [START_REF] Langdon | The climate of Scotland over the last 5000 years inferred from multiproxy peatland records: inter-site correlations and regional variability[END_REF]) and for examining intersite records across specific isochrons or 'snapshots in time', allowing the spatial synchroneity of environmental changes to be scrutinised [START_REF] Hall | Tephra-linked landscape history around 2310 BC of some sites in counties Antrim and Down[END_REF][START_REF] Hall | Assessing the impact of Icelandic volcanism on vegetation systems in the north of Ireland in the fifth and sixth millennia BC[END_REF][START_REF] Langdon | Snapshots in time: precise correlations of peat-based proxy climate records in Scotland using mid-Holocene tephras[END_REF]. Tephrochronology may also prove helpful during periods when the chronological precision of 14 C is limited [START_REF] Plunkett | New dates for first millennium BC tephra isochrones in Ireland[END_REF][START_REF] Swindles | A delayed climatic response solar forcing at 2800 cal. BP: multiproxy evidence from three Irish peatlands[END_REF].

Tephra layers are sometimes visible to the naked eye, for example in Iceland, New Zealand and Kamchatka [START_REF] Braitseva | Holocene key-marker tephra layers in Kamchatka, Russia[END_REF][START_REF] Larsen | Geochemistry of historical-age silicic tephras in Iceland[END_REF][START_REF] Lowe | Tephras and New Zealand archaeology[END_REF], but invisible (and often microscopic) cryptotephra layers (sensu [START_REF] Lowe | A summary of terminology used in tephra-related studies[END_REF]) have been found in many more areas of the world at a distance from volcanic sources (e.g. [START_REF] Persson | Tephrochronological investigations of peat deposits in Scandinavia and the Faroe Islands[END_REF][START_REF] Dugmore | Icelandic volcanic ash in Scotland[END_REF][START_REF] Dugmore | Seven tephra isochrones in Scotland[END_REF][START_REF] Hall | Late-Quaternary Icelandic tephras in Ireland and Great Britain: detection, characterization and usefulness[END_REF][START_REF] Van Den Bogaard | Two tephra layers bracketing late Holocene paleoecological changes in northern Germany[END_REF][START_REF] Zillén | Calendar year ages of three middle Holocene tephra layers identified in varved lake sediments in west central Sweden[END_REF][START_REF] Boygle | Towards a Holocene tephrochronology for Sweden: geochemistry and correlation with the North Atlantic tephra stratigraphy[END_REF][START_REF] Chambers | Tephrostratigraphy of An Loch Mor, Inis Oirr, western Ireland: implications for Holocene tephrochronology in the northeastern Atlantic region[END_REF][START_REF] Pilcher | A Holocene tephra record from the Lofoten Islands[END_REF][START_REF] Wastegård | Late Quaternary tephrochronology of Sweden: a review[END_REF][START_REF] Gehrels | A continuous 5300-yr Holocene cryptotephrostratigraphic record from northern New Zealand and implications for tephrochronology and volcanic hazard assessment[END_REF][START_REF] Hang | First discovery of cryptotephra in Holocene peat deposits of Estonia, eastern Baltic[END_REF][START_REF] Payne | Using cryptotephras to extend regional tephrochronologies: An example from southeast Alaska and implications for hazard assessment[END_REF]. Some historic tephras can be assigned a known calendar date such as the eruptions of Hekla in Iceland in 1947, 1510and 1104(cf. Swindles et al. 2008). The ages of other tephras have been constrained through high-precision radiocarbon techniques including wiggle-match dating of the deposits in which they are found [START_REF] Pilcher | Dates of Holocene Icelandic volcanic-eruptions from tephra layers in Irish peats[END_REF][START_REF] Plunkett | New dates for first millennium BC tephra isochrones in Ireland[END_REF][START_REF] Barber | Dating the Glen Garry tephra: a widespread late-Holocene marker horizon in the peatlands of northern Britain[END_REF][START_REF] Wastegård | Age, geochemistry and distribution of the mid-Holocene Hekla-S/Kebister tephra[END_REF] and Bayesian age modelling approaches [START_REF] Blockley | Improved age modelling and highprecision age estimates of late Quaternary tephras, for accurate palaeoclimate reconstruction[END_REF]. Tephrochronology is an extremely powerful tool for constraining past environmental events registered in peat sequences. However, this can be achieved only if accurate, quantitative distributional and geochemical data are obtained from the tephra layer in question. The chemical stability of tephras in peat has been recognised and repeatedly demonstrated by the successful correspondence of distal tephra geochemistries with material from source eruptions (e.g. [START_REF] Dugmore | Geochemical stability of finegrained silicic Holocene tephra in Iceland and Scotland[END_REF][START_REF] Wastegård | Late Quaternary tephrochronology of Sweden: a review[END_REF].

This paper describes the analytical procedures for dating peat profiles using tephrochronology. Tephras are firstly detected in peat ashes using light microscopy. Where tephra layers are encountered, an acid digestion technique is used to extract the shards from the peat matrix and the tephras are then mounted on slides for electron probe microanalysis. Identification of tephras is carried out through comparisons of geochemical data with results from previous studies.

FINDING TEPHRA LAYERS IN PEAT PROFILES

Detecting and isolating tephras

Several non-destructive approaches may be used as a first stage in the detection of microscopic tephra layers in peat profiles, such as magnetic susceptibility, spectrophotometry and X-ray fluorescence (cf. [START_REF] Gehrels | Towards rapid assay of cryptotephra in peat cores: Review and evaluation of various methods[END_REF]. The methods for analysis of tephra in peats that we present here are adapted from [START_REF] Dugmore | Geochemical stability of finegrained silicic Holocene tephra in Iceland and Scotland[END_REF] and [START_REF] Hall | Late-Quaternary Icelandic tephras in Ireland and Great Britain: detection, characterization and usefulness[END_REF]. For an initial, rapid assessment of the occurrence of tephra, peat sequences can be subdivided into contiguous 5 cm-long samples along the length of the core or monolith, then the samples placed into crucibles and burnt in a muffle furnace at 600°C for six hours [START_REF] Hall | Late-Quaternary Icelandic tephras in Ireland and Great Britain: detection, characterization and usefulness[END_REF]. The resulting ashes are suspended in 10% HCl for up to 24 hours, and then washed with deionised water to remove all traces of acid. The aqueous samples are centrifuged at 3,000 r.p.m. for 10 minutes (no brake set) to concentrate the tephra at the bottoms of the tubes. For samples containing a significant proportion of minerogenic material, such as those from minerotrophic peats, sieving at 10 µm in an ultrasonic bath for 5-10 minutes can be carried out, although there is a risk that this may damage some fine glass shards. This stage is not usually needed for ombrotrophic peats. Samples rich in biogenic silica (e.g. diatoms, phytoliths) can be treated with hot 5% NaOH or KOH [START_REF] Rose | Selective concentration and enumeration of tephra shards from lake sediment cores[END_REF] but treatment should be limited to one hour to prevent chemical alteration of the volcanic glass.

Slide preparation for petrography

Samples are added in drops (using a Pasteur pipette) to a glass microscope slide on a hotplate at a low heat setting, and left until all the water evaporates. A drop of Histomount (or similar) is added to the slide and a coverslip applied. The slides should be left to dry and must be stored horizontally until the mounting medium has hardened.

Petrography

Slides are scanned using a transmitted light microscope at a total magnification of 100-400x.

Plane-polarised light may be used to differentiate between tephra and other minerogenic material [START_REF] Hall | Late-Quaternary Icelandic tephras in Ireland and Great Britain: detection, characterization and usefulness[END_REF][START_REF] Chambers | Tephrostratigraphy of An Loch Mor, Inis Oirr, western Ireland: implications for Holocene tephrochronology in the northeastern Atlantic region[END_REF]. Tephra can be recognised on the basis of its isotropism, colour, shape and vesicularity [START_REF] Heiken | An Atlas of Volcanic Ash[END_REF][START_REF] Westgate | Correlation techniques in tephra studies[END_REF]. Individual shards should be counted in transects and their colour and morphological characteristics may be noted (Figure 1).

A quantitative estimation of the mineral vs. the glass shard content can be achieved by counting a minimum of 300 grains using the method of [START_REF] Van Harten | On the estimation of relative grain frequencies in heavy mineral slides[END_REF]. Glass shards can be described in terms of their size and shape, but also the size and shape of their gas inclusions [START_REF] Heiken | Morphology and petrography of volcanic ashes[END_REF][START_REF] Heiken | An Atlas of Volcanic Ash[END_REF][START_REF] Schmid | Descriptive nomenclature and classification of pyroclastic deposits and fragments: Recommendations of the IUGS Subcommission on the systematics of igneous rocks[END_REF]. Thin sections of the dense mineral residue can also be analysed to give a better idea of the relative proportions of the principal minerals found in the tephra. The minerals present in thin section can provide an insight into the formation of the tephra in question (type of magma, geotectonic context) (cf. [START_REF] Mackenzie | Atlas of Rock-Forming Minerals in Thin Section[END_REF].

Specific minerals may be helpful in distinguishing volcanic centres when glass shards are absent or in such a low concentration that they do not allow for a subsequent chemical investigation. For example, titanite is a typical heavy mineral from the Laacher See tephra [START_REF] Van Den Bogaard | Laacher See Tephra: A widespread isochronous late Quaternary tephra layer in Central and Northern Europe[END_REF], and has been used to identify the event when glass shards were absent from sediments [START_REF] Juvigné | Distribution of widespread postglacial tephra sheets from the Eifel and the Central-Massif in NE France and neighbouring regions[END_REF].

Quantitative analysis

Following the initial detection of tephra at coarse resolution, the above procedures should be carried out using contiguous samples of peat in the section(s) of the core where shard concentrations have been identified to develop a tephrostratigraphy, i.e. a stratigraphic outline of tephra horizons, for the profile. Peat volume can be measured by calibrated syringe or a similar procedure. This stage can be combined with loss-on-ignition analysis if desired. Tablets of Lycopodium clavatum spores (cf. Stockmarr 1971) (or similar) may be added to the ashes during suspension in HCl for semiquantitative determination of tephra shard concentration (cf. [START_REF] Gehrels | A continuous 5300-yr Holocene cryptotephrostratigraphic record from northern New Zealand and implications for tephrochronology and volcanic hazard assessment[END_REF]. If this procedure is not used, then each individual tephra shard per 1 cm 3 sample must be counted to determine concentration accurately. Tephra concentration can be presented as number of shards per cm 3 (Figure 2).

ANALYSIS OF TEPHRA SHARDS

Preparation for chemical analyses

Tephra extracted by the combustion method cannot be used for geochemical analysis because the procedure alters the chemistry of the glass, so an acid digestion technique is used [START_REF] Hall | Late-Quaternary Icelandic tephras in Ireland and Great Britain: detection, characterization and usefulness[END_REF]. A sample of at least 1-3 cm 3 of peat, representing the horizon of maximum tephra concentration, is placed in a beaker on a hotplate (medium heat) inside a fume cupboard, then up to 200 ml of concentrated H 2 SO 4 is added until water is driven off. A small amount of concentrated HNO 3 can then be added; the solution reacts vigorously and dark fumes of NO 2 are given off. Further small amounts of HNO 3 are added until the reaction stops. When digestion of the peat is complete, the end product should be a straw yellow coloured solution. The solution should be left to cool and then diluted in one litre of water before fine-sieving (e.g. 10 μm). The sample must be washed with copious amounts of deionised water (2-5 L) to ensure that all traces of acid are removed. The sample can then be pipetted into a centrifuge tube and centrifuged at 3,000 r.p.m. for ten minutes (no brake set) to concentrate the tephra, and the supernatant is decanted. In addition we recommend that, if organic or waxy deposits remain, the digestion should be repeated [START_REF] Hall | Late-Quaternary Icelandic tephras in Ireland and Great Britain: detection, characterization and usefulness[END_REF]. Hot dilute alkali Figure 2. Tephrostratigraphy of Dead Island Bog (Northern Ireland), from the Hekla 4 isochron [START_REF] Larsen | H4 and other acid Hekla tephra layers[END_REF] which has been 14 C wiggle-match dated to 2339-2279 cal. BC [START_REF] Pilcher | Dates of Holocene Icelandic volcanic-eruptions from tephra layers in Irish peats[END_REF] to the historically dated Hekla AD 1947 tephra [START_REF] Thorarinsson | The eruptions of Hekla in historical times, I: The eruption of Hekla 1947-1948[END_REF]) (after [START_REF] Swindles | A 4500-year proxy climate record from peatlands in the North of Ireland: the identification of widespread summer 'drought phases[END_REF]. Tephra counts are expressed as number of shards cm -3 . treatment may also be required if the sample is rich in biogenic silica. The size of the peat sample should be chosen to ensure that enough glass shards are concentrated to obtain sufficient analyses.

Slide preparation for geochemical analysis

As each electron microprobe laboratory may have specific requirements regarding the way in which samples for analysis are mounted, it is worth contacting the system administrator in advance of slide preparation. Typically, the surface of a glass slide is coarsely ground using 600 µm alumina grit or carborundum paper to improve sample cohesion (note that the slide must be of suitable dimensions to fit within the instrument that will be used to analyse the tephra geochemistry). The slide should be labelled using a diamond scriber on the un-ground side and a pencil on the ground side. Small square cells and labels for the tephra samples can also be drawn (using a pencil) on the ground side depending on sample size (up to six per slide). The samples are pipetted onto their individual cells and the water evaporated on a hotplate at a low heat setting. One drop of epoxy resin (e.g. Araldite or equivalent, following the laboratory instructions) is applied over the dry specimens (ensuring a clean spatula is used for each sample) and left to dry on a hotplate. The slide should be ground on alumina grit or carborundum paper using progressively finer grades to leave a layer of resin 30-50 µm thick. Grinding should continue until tephra shards are nearly exposed (which should be checked regularly using a reflected light microscope to prevent over-exposure and loss of shards). The slide should then be washed in deionised water in an ultrasonic bath and polished using fine grade abrasive paper or polishing paste (alumina or diamond) at 6 µm, 3 µm and 1 µm (washing in an ultrasonic tank before each stage) to achieve a fine polish on the slide surface. The slide should then be meticulously washed in deionised water and carbon coated under a vacuum following the procedures of the electron probe laboratory. Instead of the glass slide method, some laboratories use brass or epoxy resin disks in which the tephra is mounted before being ground and polished.

Geochemical analysis

It is crucial that geochemical analysis of tephra follows protocols that provide consistently reliable data and researchers should be aware that interlaboratory discrepancies have been noted [START_REF] Westgate | Correlation techniques in tephra studies[END_REF][START_REF] Hunt | Tephra geochemistry: a discussion of some persistent analytical problems[END_REF], 1996;[START_REF] Potts | G-Probe-1 -An international proficiency test for microprobe laboratories -Report on round 1 February 2002 (TB-1 basaltic glass)[END_REF][START_REF] Coulter | Testing the reliability of the JEOL FEGSEM 6500F electron microprobe for quantitative major element analysis of glass shards from rhyolitic tephra[END_REF] Tephra geochemistry should be determined using wavelength-dispersive electron probe microanalysis. Although settings of the instrument may vary between laboratories, a defocused beam of 5-10 µm should be used and Na should be counted first to guard against migration [START_REF] Hunt | Tephra geochemistry: a discussion of some persistent analytical problems[END_REF], 2001). Algorithms such as ZAF, PAP or X-PHI should be used to correct for absorbance and fluorescence effects [START_REF] Sweatman | Quantitative electron-probe microanalysis of rock-forming minerals[END_REF][START_REF] Pouchou | Quantitativeanalysis of homogeneous or stratified microvolumes applying the model PAP[END_REF][START_REF] Merlet | An accurate computer correction program for quantitative electron probe microanalysis[END_REF]. Primary and secondary standards such as pure synthetic oxides, simple silicates, andradite garnet, standard basalts and Lipari obsidian should be used for calibration and to check the accuracy of the analyses before and during analysis runs [START_REF] Hunt | Interlaboratory comparison of electron probe microanalysis of glass geochemistry[END_REF]. The main species that should be obtained are SiO 2 , TiO 2 , Al 2 O 3 , K 2 O, CaO, FeO (total) , Na 2 O, MgO and MnO. Other species, such as F and P 2 O 5 can be analysed if possible. Systems with linked energy dispersive systems (EDS) may be of particular use to help pinpoint very sparse tephra shards in thin section. The tephra geochemistry can be compared with previous analyses from online databases (e.g. Tephrabase, http://www.tephrabase.org/) or published works, and presented in biplots and/or ternary diagrams (Figure 3). Principal Components Analysis (PCA) may also be useful for tephra discrimination in certain circumstances [START_REF] Pollard | Some numerical considerations in the geochemical analysis of distal microtephra[END_REF]). However, it should be interpreted alongside the more traditional bivariate or trivariate major oxide diagrams [START_REF] Pearce | Comment on "Some numerical considerations in the geochemical analysis of distal microtephra[END_REF].

It is common practice to omit data with less than 95% of the analytical total, and this may be appropriate in many cases [START_REF] Hunt | Tephra geochemistry: a discussion of some persistent analytical problems[END_REF]. However, it is not always possible to obtain totals >95%, especially from tephras containing a large amount of water (either primary magmatic or through post-depositional hydration of the glass) [START_REF] Pearce | Comment on "Some numerical considerations in the geochemical analysis of distal microtephra[END_REF]. The number of analyses needed for the identification of the tephra is variable and dependent on the discretion of the researcher. It is often said that twenty or more is a suitable number, but this may not be attainable with very sparse layers. Sometimes distinct tephras can be identified with five (or even fewer) analyses, especially if other dating evidence is available. However, researchers should strive to get as many analyses as possible from a tephra layer. Normalisation of tephra geochemical data is common practice in the Southern Hemisphere and in North and South America, as this can remove the variable effects of hydration in different samples (e.g. [START_REF] Pearce | Comment on "Some numerical considerations in the geochemical analysis of distal microtephra[END_REF]. This practice is probably best avoided, however, as it masks the quality of the data [START_REF] Pollard | Some numerical considerations in the geochemical analysis of distal microtephra[END_REF]) and does not allow other researchers access to the raw data.

Density separation

In peats containing a large amount of minerogenic material, density separation procedures may be needed to extract and concentrate the shards (cf. [START_REF] Turney | Extraction of rhyolitic component of Vedde microtephra from minerogenic lake sediments[END_REF]. We have achieved good results using LST Fastfloat (http://www.polytungstate.co.uk/), which is a solution of low toxicity sodium heteropolytungstates in water. Densities of 2.3-2.5 g cm -3 are commonly used in tephrochronology, which can be checked using a hydrometer [START_REF] Turney | Extraction of rhyolitic component of Vedde microtephra from minerogenic lake sediments[END_REF]. Bromoform (density 2.8 g cm -3 ) can also be used for the separation of glass shards. This method is relatively fast and allows the separation of minerals with grain sizes below 300 μm, which is typical of tephra layers [START_REF] Juvigné | Scheidetrichtermethode oder Zentrifugaltrennung zur quantitativen Gewinnung von Schwermineralen (Separating funnel or centrifugal method for the quantitative separation of heavy minerals)[END_REF]). However, bromoform should be used with caution as it is toxic to humans and the environment.

POTENTIAL AND LIMITATIONS

Although a standard protocol for tephrochronology is presented here, there is a growing body of literature presenting experimental methods for the detection, extraction, analysis and dating of tephra (e.g. [START_REF] Caseldine | A rapid, non-destructive scanning method for detecting distal tephra layers in peats[END_REF][START_REF] Buck | Bayesian tools for tephrochronology[END_REF][START_REF] Blockley | A new and less destructive laboratory procedure for the physical separation of distal glass tephra shards from sediments[END_REF], Payne & Blackford 2005[START_REF] References Auclair | Luminescence investigation of loess and tephra from Halfway House section, Central Alaska[END_REF][START_REF] Gehrels | Towards rapid assay of cryptotephra in peat cores: Review and evaluation of various methods[END_REF][START_REF] De | Development and application of high-resolution petrography on resin-impregnated Holocene peat columns to detect and analyse tephras, cryptotephras, and other materials[END_REF]. Tephrochronology is firmly established as a useful chronological tool for peat-based studies, but a number of problems and limitations are apparent. One limitation is that tephra layers are not found everywhere in the world, or are sometimes characterised by very low shard concentrations, making identification difficult. A further problem is that tephra fallout is often not spatially coherent, leading to a distinct 'patchiness' of tephra layer distribution [START_REF] Dugmore | Long-distance marker horizons from small-scale eruptions: British tephra deposits from the AD 1510 eruption of Hekla, Iceland[END_REF]. Modelling approaches have illustrated that climate variability is an important factor in the atmospheric transport of Icelandic tephra in the subpolar North Atlantic [START_REF] Lacasse | Influence of climate variability on the atmospheric transport of Icelandic tephra in the subpolar North Atlantic[END_REF]. 'Taphonomic' problems may also occur in peats as the tephra may have been deposited on an uneven surface and re-working can occur [START_REF] Holmes | A Tephra-dated study of Vegetation and Climate Change in the Mid-Holocene of North West Europe[END_REF], Payne et al. 2005[START_REF] Payne | The formation of tephra layers in peatlands: an experimental approach[END_REF]. In volcanic regions it is often difficult to interpret the abundance of inter-mixed tephra layers. Also, the zone of highest concentration of glass may not necessarily represent the primary ash-fall layer but instead reflect the erosion of exposures or surfaces containing tephra adjacent to a peat bog.

Some studies have suggested that there is chemical alteration of tephra shards in the postdepositional environment due to cationic leaching from the matrix, or complete destruction of the silica network (e.g. [START_REF] Pollard | Chemical alteration of tephra in the depositional environment: theoretical stability modelling[END_REF]. From their experimental work, Wolff-Boenish et al. (2004a, 2004b) proposed a relationship between the lifetime of 1 mm-radius natural glass spheres and their silica content at 25°C and pH 4. The dissolution rate depends mainly on the silica content of the tephras: rhyolitic tephras are generally more stable than basaltic [START_REF] Oelkers | General kinetics description of multioxide silicate mineral and glass dissolution[END_REF][START_REF] Pollard | Chemical alteration of tephra in the depositional environment: theoretical stability modelling[END_REF], Wolff-Boenish et al. 2004a, 2004b). The dissolution processes are slow and a 1 mm thick layer of glass can take up to 4,000 years to dissolve under these conditions (Wolff-Boenish et al. 2004a, 2004b). On the other hand, [START_REF] Hodder | Dissolution and depletion of ferromagnesian minerals from Holocene tephra layers in an acid bog, New Zealand, and implications for tephra correlation[END_REF] reported how glass and ferromagnesian minerals in very acidic anoxic peat bogs could be partly or totally dissolved, sometimes leading to complete dissolution of certain minerals (e.g. biotite) in relatively short time periods.

These different observations illustrate that a systematic and thorough analysis of sediments, tephra and related alteration products may be needed when studying distal tephras in acidic peatland environments. Nevertheless, the successful linkage of distal tephras to source material in peatlands across Europe (e.g. [START_REF] Hall | Late-Quaternary Icelandic tephras in Ireland and Great Britain: detection, characterization and usefulness[END_REF][START_REF] Van Den Bogaard | Two tephra layers bracketing late Holocene paleoecological changes in northern Germany[END_REF][START_REF] Wastegård | Late Quaternary tephrochronology of Sweden: a review[END_REF] and North America (e.g. [START_REF] Payne | Using cryptotephras to extend regional tephrochronologies: An example from southeast Alaska and implications for hazard assessment[END_REF] demonstrates that chemical alteration is probably not a critical problem when working with peat profiles.

Several studies have also provided insights into the effects of bacterial activities on glass dissolution (e.g. [START_REF] Thorseth | Textural and chemical effects of bacterial activity on basaltic glass: an experimental approach[END_REF][START_REF] Staudigel | Biologically mediated dissolution of volcanic glass in seawater[END_REF][START_REF] Brehm | The role of microorganisms and biofilms in the breakdown and dissolution of quartz and glass[END_REF]. Experimentally, [START_REF] Thorseth | Textural and chemical effects of bacterial activity on basaltic glass: an experimental approach[END_REF] found that the dissolution rate can increase by a factor of ten after six months of glass alteration due to the formation of new bacterial communities. [START_REF] Birkefeld | In situ investigation of dissolution of heavy metal containing mineral particles in acidic forest soil[END_REF] also emphasised that dissolution rates are greater in the field than in the laboratory, due to extensive biological influence from fungi, bacteria and plant roots. However, although bacterial action could induce a rapid destruction of the glass walls, bacteria will also form and/or adsorb amorphous particles (from organic or mineral composition) at the grain surface, developing an alteration film that could decrease or stop the alteration of tephra. Such protective barrier films have been studied in the case of basaltic glass alteration and the use of this material as radioactive stockpile (e.g. [START_REF] Techer | Dissolution kinetics of basaltic glasses: control by solution chemistry and protective effect of the alteration film[END_REF].

There may also be problems in discriminating tephra layers using major element geochemistry. For example, the historical Hekla 1510 and 1947 tephras cannot be separated on the basis of their major element geochemistry, and other stratigraphic or chronological information is needed for their reliable identification [START_REF] Swindles | Constraining the age of spheroidal carbonaceous particle (SCP) stratigraphies in peats using tephrochronology[END_REF]. Trace and rare earth elemental composition of glass shards using laser ablation ICP-MS or ion probe analysis may be needed to discriminate between such tephras [START_REF] Pearce | Correlation and characterisation of individual glass shards from tephra deposits using trace element laser ablation ICP-MS analyses: current status and future potential[END_REF][START_REF] Coulter | Testing the reliability of the JEOL FEGSEM 6500F electron microprobe for quantitative major element analysis of glass shards from rhyolitic tephra[END_REF].

Figure 1 .

 1 Figure 1. Photograph showing tephra shards with characteristic fluted and vesicular morphologies. The shards are from a currently unidentified tephra layer from the Shetland Isles, UK.

Figure 3 .

 3 Figure 3. Ternary plot showing FeO (total) , CaO and K 2 O for the Hekla 1947, 1510 and Hekla 4 tephra layers encountered in two peatlands in the North of Ireland (after Swindles et al. 2010). It is evident that Hekla 1947 (circles) and 1510 (crosses) are indistinguishable on the grounds of their major element oxide geochemistry, whereas Hekla 4 (triangles) can be easily distinguished from these layers.

ACKNOWLEDGEMENTS

Many thanks are expressed to Maria Gehrels and one anonymous reviewer; also to Ian Lawson for comments on an early draft of the manuscript.


Rimond Hamia

Christophe Cordier

Sébastien Saez

Christophe Dolabdjian

Giant Magneto impedance Sensor for Non Destructive Evaluation Eddy Current System


Souleymen Sahnoun

Pierre Comon

Deterministic Blind Identification in Antenna Array Processing

The estimation of directions of arrival is formulated as the decomposition of a 3-way array into a sum of rankone terms. However, a low-rank tensor approximation does not always exist. We propose an optimization technique based on differentiable angular constraints on the factors, ensuring the existence of the low-rank tensor decomposition. The efficiency of the proposed algorithm is demonstrated via numerical simulations, and compared to Cramér-Rao bounds.

I. INTRODUCTION

Estimation of Directions of Arrival (DoA) is a central problem in antenna array processing, including in particular radar, sonar, or telecommunications [1]. Traditional approaches are based on low-rank approximation of the covariance matrix of observations, and on detecting points of minimal distance with the so-called array manifold [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF] [START_REF] Van Trees | Optimum Array Processing[END_REF]. These approaches hence assume that (i) the measurements are weakly stationary over sufficiently long observation lengths, (ii) the number of sources of interest is smaller than the number of sensors, and (iii) spatial responses of all sensors are known, and in particular their location (in other words, the sensor array needs to be calibrated).

In [START_REF] Sidiropoulos | Parallel factor analysis in sensor array processing[END_REF], a deterministic approach has been proposed, which permits not only to work with short data lengths (and hence less stationary sources) , but also to localize more sources than sensors present in the reference array. This approach is based on the same rotational invariance as exploited in [START_REF] Roy | ESPRIT -estimation of signal parameters via rotational invariance techniques[END_REF]. It consists in storing the measurements in a 3-way array, and to decompose it into a sum of rank-one terms. One very interesting by-product of [START_REF] Sidiropoulos | Parallel factor analysis in sensor array processing[END_REF] is that source copies are also delivered for free, without any further estimation stage.

In this paper we revisit the approach of [START_REF] Sidiropoulos | Parallel factor analysis in sensor array processing[END_REF], where an important fact has been neglected: a low-rank tensor approximation does not always exist, so that the latter approach is actually ill posed. This fact has been already pointed out in [START_REF] Lim | Blind multilinear identification[END_REF], and additional constraints have been suggested, which enjoy a reliable physical meaning and at the same time ensure existence of a solution. In the present contribution, we shed some light on conditioning and algorithmic issues.

In the next section, we introduce notation and formulate the DoA estimation as a low-rank decomposition of a 3-way array. Then we describe in Section III an optimization algorithm using a new differentiable constraint ensuring existence of the low-rank tensor approximation. Finally, we show the benefit of the proposed algorithm via numerical experiments, and compare results to Cramér-Rao bounds.

II. MODELING AND NOTATION

In the following, vectors will be denoted by bold lowercases, e.g. a, whereas matrices or higher-order arrays will be denoted by bold uppercases, e.g. A. Moreover, a r will denote the rth column of matrix A.

Suppose that R narrow-band radiating sources impinge on an array of sensors, formed of L subarrays of K sensors each. We make the far-field assumption, that is, we assume that sources are located sufficiently far from the array, compared to the array dimensions, so that waves can be considered as plane. The key assumption made in [START_REF] Roy | ESPRIT -estimation of signal parameters via rotational invariance techniques[END_REF], [START_REF] Sidiropoulos | Parallel factor analysis in sensor array processing[END_REF], [START_REF] Lim | Blind multilinear identification[END_REF] is that, taking one subarray as reference, every subarray can be deduced from the reference one by an unknown translation in space, defined by some vector δ ℓ of R 3 , 1 ≤ ℓ ≤ L, δ 1 def = 0. Denote ς r (t) the signal transmitted by the rth source, d r its DoA viewed from the array, and s(y, t) the signal measured at a point in space defined by its coordinates y (we consider complex envelopes about the central frequency). Then we have:

s(y, t) = R r=1 ς r (t) a r (y), a r (y) def = exp{ ω C y T d r } (1)
where ω is the central pulsation of the narrow-band waves, C the wave celerity, and  = √ -1. Because waves are plane and narrow-band, the signal measured at another point y + τ , deduced from y by a translation τ takes the form:

s(y, τ , t) = R r=1 a r (y) b r (τ ) ς r (t), b r (τ ) def = exp{ ω C τ T d r }
(2) In other words, function s(y, τ , t) decomposes into a sum of R simpler functions whose variables separate. Now, if we discretize the R 3 space with the above defined subarrays, and take M time samples, we end up with a multilinear relationship. In fact, let p k be the coordinate vector of the kth sensor of the reference subarray, and δ ℓ the translation defining the location of the ℓth subarray, 1 < ℓ ≤ L. Then signal (2) can be stored in a K × L × M three-way array, which follows the model below:

T kℓm = R r=1 λ r A kr B ℓr S mr (3)
where 3) is related to the Canonical Polyadic decomposition (CP) 1 , which consists of decomposing a tensor T into a sum of decomposable tensors. For the sake of convenience, equation ( 3) is rewritten in vector form as

A kr = 1 √ K exp  ω C p T k d r , B ℓr = 1 √ L exp  ω C δ T ℓ d r , S mr = ς r (t m )/ ς r , and λ r = √ KL ς r . Model (
t = R r=1 λ r a r ⊗ b r ⊗ s r (4)
where ⊗ denotes the Kronecker product, and t = vec{T} is a column vector of dimension KLM containing the entries of the 3-way array T.

III. EXISTENCE AND UNIQUENESS

The goal is to identify the directions of arrival (DoA), θ r , of the R impinging plane waves and to estimate corresponding transmitted source signals ς r (t m ) up to a scale factor, given the whole array T. Clearly, a sufficient condition is to be able to identify all parameters in the RHS of (4).

A. Low rank approximation

Actually, observations are corrupted by noise, so that (3-4) do not hold exactly. A natural idea is then to fit model (4) by minimizing the error

Υ(A, B, S; Λ) = t - R r=1 λ r a r ⊗ b r ⊗ s r 2 (5)
where Λ denotes the diagonal matrix containing the λ r 's. Minimizing error (5) means finding the best rank-R approximate of T and its CP decomposition. However, the infimum of Υ may not be reached; see e.g. [START_REF] Comon | Tensors, usefulness and unexpected properties[END_REF], [START_REF] Lim | Blind multilinear identification[END_REF] and references therein. The reason is that the set of rank-R tensors is not closed for R > 1. The idea we promote here is to impose an additional constraint that will ensure the existence of a minimum, as elaborated in the next section. In addition, from the physical point of view, one can make the following observations:

• sources that are totally correlated need to be localized separately only if they are sufficiently well angularly separated. In that case they correspond to multi-paths. • sources that are located in the same direction need to be estimated separately only if they are sufficiently decorrelated. • otherwise, one can assimilate highly correlated sources arriving from close directions to a single fat source, spread out in space. The purpose of the section is to formalize these constraints.

1 also sometimes called Candecomp/Parafac in Psychometry.

B. Coherences

As in the compressed sensing literature [START_REF] Gribonval | Sparse representations in unions of bases[END_REF], [START_REF] Donoho | Optimally sparse representation in general (nonorthogonal) dictionaries via l1 minimization[END_REF], we define the coherence of a set of unit norm vectors as the largest value of cross scalar products:

µ A = sup k =ℓ |a H k a ℓ | (6)
The coherence of matrix A is defined this way, if a k denote its (unit norm) columns. Coherences of matrices B and S are defined similarly, and denoted by µ B and µ S , respectively. Let G be the R × R Gram matrix defined by:

G pq = (a H p a q )(b H p b q )(s H p s q )
Then for given matrices A, B and S, the optimal value Λ o minimizing error Υ is obtained by cancelling the gradient of (5) w.r.t. Λ, which leads to the linear system:

G λ o = f , (7)
where λ o = diag(Λ o ) and vector f in the right hand side is defined by the contraction f r = ijk T ijk A * ir B * jr S * kr , 1 ≤ r ≤ R. Equation [START_REF] Comon | Tensors, usefulness and unexpected properties[END_REF] shows that coherences play a role in the conditioning of the minimization problem. Also note that only the product between coherences appears, and not coherences individually.

C. Existence

We are now in a position to state conditions of existence. It has been shown in [START_REF] Lim | Blind multilinear identification[END_REF] that if

µ A µ B µ S < 1 R -1 (8)
then the infimum of ( 5) is reached. This happens because error (5) becomes coercive as soon as ( 8) is satisfied. And it must then reach its minimum since it is continuous. We refer to [START_REF] Lim | Blind multilinear identification[END_REF] for a discussion on the physical signification of this constraint, which may significantly modify the CP decomposition. Constraint [START_REF] Gribonval | Sparse representations in unions of bases[END_REF] needs some care because it involves max operators, which are not differentiable. For this reason, we propose to use the fact that L ∞ norm can be bounded by L 2ρ norms, and approximated for large values of ρ:

||z|| ∞ = max k {z k } ≤ ||z|| 2ρ = ( k z 2ρ k ) 1/2ρ , ∀ρ ≥ 1 for z k ∈ R + .
Applying this inequality to z k ≡ a H p a q allows to bound coherences above by a differentiable quantity, so that another (somewhat more constraining) sufficient condition can be obtained. More precisely:

µ A ≤ µ(A, ρ) def = p<q |a H p a q | 2ρ 1/2ρ
We subsequently call C ρ the constraint obtained by replacing the max operators by the L 2ρ norms in constraint (8):

C ρ def = 1 -R + µ(A, ρ) -1 µ(B, ρ) -1 µ(S, ρ) -1 > 0 (9)

D. Uniqueness

There exist sufficient conditions ensuring that the solution of ( 4) is unique, which involve coherences [START_REF] Lim | Blind multilinear identification[END_REF]. However, the condition below is much less constraining [10, p.13]:

R ≤ M and R(R -1) ≤ K(K -1)L(L -1) 2 (10)
and guarantees that there exists almost surely a unique solution. Other sufficient conditions for generic uniqueness exist [START_REF] Catalisano | Higher secant varieties of the Segre varieties[END_REF], [START_REF] Abo | Induction for secant varieties of Segre varieties[END_REF], but are less attractive when one dimension (i.e. M ) is large.

IV. COMPUTER RESULTS

A. Optimization

The constrained optimization is carried out with the help of gradient descent type algorithms, which handle constraints in different manners. Denote for conciseness x = vec{[A T , B T , S T ]} and define the objective function:

F ρ (x; λ) = Υ(x; Λ) + η exp(-γ C ρ (x)) ( 11
)
where η is the penalty weight, γ is introduced to control the importance of penalty C ρ (x), and λ is defined in [START_REF] Comon | Tensors, usefulness and unexpected properties[END_REF] and depends on x and t. This leads to the algorithm below ALGORITHM 1) Choose R satisfying [START_REF] Lathauwer | A link between canonical decomposition in multilinear algebra and simultaneous matrix diagonalization[END_REF].

2) Initialize (A(0), B(0), S(0)) to matrices with unit-norm columns satisfying C ρ ≥ 0. 3) Compute G(0) and f (0), and solve G(0) λ(0) = f (0) for λ, according to (7) 4) For k ≥ 1 and subject to a stopping criterion, do a) Compute the descent direction as the gradient w.r.t.

x: f (k) for λ, according to [START_REF] Comon | Tensors, usefulness and unexpected properties[END_REF]. In the algorithm, η is decreased as the reconstruction error Υ(x; Λ) decreases, whereas γ is kept fixed.

d(k) = -∇F ρ (x(k -1); λ(k -1)) b) Compute a stepsize ℓ(k) c) Update x(k) = x(k -1) + ℓ(k) d(k) d)
We give now some gradient expressions 2 necessary to determine the descent direction d(k) when F ρ is used:

∂Υ ∂A = 2 AM A -2 N A where M A pq def = jk λ p B jp S kp S * kq B * jq λ * q , N A ip def = jk T ijk B * jp S * kp λ * p and ∂ exp(-γ C ρ ) ∂A = γ exp (γ C ρ ) L A ρ A (A H A) ⊡ Q A -I
2 Matrix gradients are written with the conventions described in [START_REF] Comon | Estimation multivariable complexe[END_REF], [START_REF] Hjørungnes | Complex-valued matrix differentiation: Techniques and key results[END_REF].

where ⊡ denotes the Hadamard (entry-wise) product,

L A ρ def = q<p |a H p a q | 2ρ -1 2ρ -1 µ(B, ρ) -1 µ(S, ρ) -1 ,
and

Q pq def = |a H q a p | 2ρ-2
. Keep in mind that expressions above hold true because matrix A has unit-norm columns. And expressions are similar for matrices B and S, which also have unit-norm columns.

B. Advantage of the constraint C ρ (x)

To see the interest of constraint C ρ (x) used in the optimization algorithm, Figure 3 sketches the evolution of the reconstruction error Υ(x; Λ), and C ρ (x) as a function of iterations. The figure shows that thanks to the constraint C ρ : (i) iterates are incited to remain/turn back to the feasible region (where existence is guaranteed), (ii) the optimization algorithm converges quickly because iterates are allowed to move away from the feasible region (depending on parameters η and γ).

C. Monte Carlo experiments

In this section, we asses the performances of the proposed algorithm using several sensor subarrays in two different configurations. The scenario on which the algorithm is tested can be of interest in various applications. For instance, consider sonar buoys left floating on the surface and equipped with a device permitting to maintain their orientation towards North. Each subarray is rigid and fixed on a buoy, but its relative location with respect to other subarrays is unknown.

The first configuration consists of two identical subarrays (L = 2) separated by a distance of 25λ (i.e. δ 2 = [0, 25λ, 0] T ), where λ = ω/2πC is the wavelength. The second configuration consists of three identical subarrays (L = 3) with the following translations: δ 2 = [0, 25λ, 0] T , δ 3 = [0, 37.5λ, 5λ] T . In both configurations, each subarray is an ULA array of 4-element with half-wavelength spacing. For the two configurations, we consider three equal-power narrowband source signals (R = 3) impinging respectively from θ 1 = 77 • , θ 2 = 55 • and θ 3 = 5 • (angles w.r.t endfire). In all experiments, γ = 5, M = 200 time samples are used, and 50 Monte-Carlo simulations are run for each SNR level. The additive noise is complex-valued circular Gaussian. Figures IV-A and IV-C report the total root-mean square error (RMSE) of the DoA using respectively the first (L = 2) and the second (L = 3) configuration. The deterministic Cramér-Rao bound (CRB) is reported as a benchmark; exact expressions of DoA CRB will be detailed in a full length paper. In fact, CRB can be found in [START_REF] Liu | Cramér-Rao lower bounds for low-rank decomposition of multidimensional arrays[END_REF] for factor matrices, e.g. A, and in [START_REF] Ottersten | Performance analysis of the total least squares ESPRIT algorithm[END_REF] for DoAs obtained with L = 2 subarrays, but DoA CRB are not available for L > 2.

Configuration 1. This experiment shows that: (i) the proposed CP algorithm exhibits the same performances as ESPRIT, which makes sense, (ii) MUSIC performs the best, but exploits more information, namely the exact knowledge of sensor locations, whereas this information is actually not available in the present scenario. Hence MUSIC performances just serve as a reference.

Configuration 2. This experiment shows that the proposed algorithm yields better results than EPSRIT. The reason is that ESPRIT uses at most two subarrays, whereas the proposed algorithm uses all of them. Again, MUSIC is reported just as a reference benchmark.

V. CONCLUDING REMARKS

The DoA estimation is formulated as a CP decomposition under the assumptions presented in § II on the sensor array formed of L identical subarrays. We proposed an optimization algorithm using a new differentiable constraint C ρ (x) ensuring existence of the low-rank tensor approximation. It was shown that thanks to C ρ (x) the proposed algorithm converges quickly and is prevented to leave for long the feasible region. DoA estimation results show that, as expected, the CP algorithm exhibits better results than reference DoA estimation methods when L > 2. Performances of signal estimation will be reported in a full length paper, as well as expressions of CRBs.

  Extract the 3 blocks from x(k): A(k), B(k) and S(k) e) Normalize the columns of A(k), B(k) and S(k) f) Compute G(k) and f (k), and solve G(k) λ(k) =

Fig. 1 .

 1 Fig. 1. Total DoA error versus SNR, with L = 2 subarrays.

Fig. 2 .

 2 Fig. 2. Total DoA error versus SNR, with L = 3 subarrays.

Fig. 3 .

 3 Fig. 3. Reconstruction error Υ(x; Λ) and Cρ(x) as a function of the number of iterations.

  Work funded by the European Research Council under the European Community's Seventh Framework Programme FP7/2007-2013 Grant Agreement no. 320594.