Rocío Carratalá-Sáez

Mathieu Faverge

Grégoire Pichon

Guillaume Sylvand

Enrique S Quintana-Ortí

Tiled Algorithms for Efficient Task-Parallel H-Matrix Solvers

Keywords: Hierarchical matrices, LU factorization, linear systems, task-parallelism, multicore processors Matrices hiérarchiques, factorisation LU, systèmes linéaires, parallélisme à base de tâches, processeurs multicoeurs struct HCHAM_desc_s {

In this paper, we describe and evaluate an extension of the Chameleon library to operate with hierarchical matrices (H-Matrices) and hierarchical arithmetic (H-Arithmetic), producing efficient solvers for linear systems arising in Boundary Element Methods (BEM). Our approach builds upon an open-source H-Matrices library from Airbus, named Hmat-oss, that collects sequential numerical kernels for both hierarchical and low-rank structures; the tiled algorithms and task-parallel decompositions available in Chameleon for the solution of linear systems; and the StarPU runtime system to orchestrate an efficient task-parallel (multi-threaded) execution on a multicore architecture. Using an application producing matrices with features close to real industrial applications, we present shared-memory results that demonstrate a fair level of performance, close to (and sometimes better than) the one offered by a pure H-Matrix approach, as proposed by Airbus Hmat proprietary (and non open-source) library. Hence, this combination Chameleon + Hmat-oss proposes the most efficient fully open-source software stack to solve dense compressible linear systems on shared memory architectures (distributed memory is under development).

Introduction

Many real-life simulations require handling large collections of data together with a considerable number of computations. Solving linear systems of equations is often a key kernel in these applications, and usually represents one of the most expensive operations, both in terms of memory usage and computational cost. The matrices that realize these systems present considerable structural variations, covering a wide spectrum of different configurations, from dense to sparse, depending on the underlying simulation they represent. Thus, reducing the cost of this operation is interesting in order to 1) solve larger problems, and 2) reduce the time-to-solution as well as other associated costs (e.g., memory or energy).

Since their introduction [START_REF] Hackbusch | A sparse matrix arithmetic based on hmatrices. part i: Introduction to h-matrices[END_REF][START_REF] Grasedyck | Construction and arithmetics of h-matrices[END_REF], hierarchical Matrices (H-Matrices) have gained substantial momentum in Boundary Element Methods (BEM) and elliptic partial differential operators [START_REF] Hackbusch | Hierarchical Matrices: Algorithms and Analysis[END_REF]. In particular, thanks to the combination of dense and low-rank blocks in a nested structure, H-Matrices nowadays provide a powerful numerical tool to reduce storage costs and execution time significantly [START_REF] Grasedyck | Construction and arithmetics of h-matrices[END_REF]. For example, the LU Factorization of an n × n H-Matrix (H-LU) requires Θ(n k 2 log 2 n) floating-point operations (flops) in H-Arithmetic, where the parameter k can be tuned to control the accuracy of the approximation (that is, it establishes the local rank for the H-Matrix subblocks). In contrast, the same factorization costs Θ( 23 n 3 ) flops in the dense case. Improving the performance of H-Arithmetic operations is an active area of research that has recently produced a fair number of libraries as well as many interesting algorithmic developments (see Section 3). These research efforts are surely motivated by the relevance of the applications that can be efficiently tackled with H-Matrices, but also because of the complexity and benefits of H-Arithmetic. On the one hand, H-operations require dealing with both low-rank and dense blocks -which often implies re-compressing some of the intermediate results -while following a nested structure, and usually a recursive algorithm. On the other hand, in general, the definition and storage of H-Matrices implies complex data accesses. This fact promoted the appearance of alternative structures, such as Block Low-Rank (BLR) [START_REF] Amestoy | Improving multifrontal methods by means of block low-rank representations[END_REF][START_REF] Mary | Block Low-Rank multifrontal solvers: complexity, performance and scalability[END_REF][START_REF] Pichon | Sparse supernodal solver using block low-rank compression: Design, performance and analysis[END_REF] and lattice H-Matrices [START_REF] Ida | Lattice h-matrices on distributed-memory systems[END_REF][START_REF] Yamazaki | Distributed-memory lattice h-matrix factorization[END_REF], that trade off slightly higher time and memory costs in exchange for superior simplicity. One asset of these approaches it that it is easier to exploit parallelism, as they present more regular structures.

In this paper we combine existing efficient numerical kernels for hierarchical, low-rank and fullrank matrices, together with an efficient task-based implementation designed to solve dense linear systems. As we will show, this allows us to leverage modern programming models and runtime systems yielding to the single open-source solution that achieves fair parallel performance, while avoiding re-implementing pure H-Arithmetic. Concretely, our work integrates the following three components:

• For H-Arithmetic, we choose the Hmat-oss [START_REF]Hmat-oss library github repository[END_REF] kernels to operate with low-rank blocks, as they have been proved to offer fair efficiency in industrial applications [START_REF] Lizé | Résolution directe rapide pour les éléments finis de frontière en électromagnétisme et acoustique: H-Matrices[END_REF]. This is a library for H-Matrices maintained by Airbus, whose public version is sequential.

• In addition, for the task-based implementation of matrix operations, we leverage the Chameleon [START_REF] Agullo | Faster, cheaper, better -a hybridization methodology to develop linear algebra software for gpus[END_REF][START_REF] Agullo | Achieving high performance on supercomputers with a sequential task-based programming model[END_REF][START_REF]Chameleon software home page[END_REF] library, a dense linear algebra package that relies on the sequential task flow programming model to schedule tiled algorithms on top of multiple runtime systems such as OpenMP [START_REF]The openmp api specification for parallel programming home[END_REF], PaRSEC [START_REF] Hoque | Dynamic task discovery in parsec: A data-flow task-based runtime[END_REF], StarPU [START_REF]Starpu runtime system home page[END_REF][START_REF] Augonnet | Starpu: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF], or Quark [START_REF] Yarkhan | Quark users' guide[END_REF].

• Finally, we focus on the specific StarPU runtime system [START_REF] Augonnet | Starpu: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF] support of the Chameleon library to exploit the task-parallelism in our solution. StarPU is in charge of issuing the tasks to the systems cores while fulfilling inter-task dependencies.

Our strategy to re-utilitze structures that are similar to H-Matrices while reducing their complexity is close to what is referred to as "lattice H-Matrices" in [START_REF] Ida | Lattice h-matrices on distributed-memory systems[END_REF]. As opposed to that work, in our case, we need to split the initial matrix into regular tiles (i.e., tiles of the same size) to meet the requirements of Chameleon algorithms requirements, and each of these tiles are individually turned into H-Matrices. We will refer to our structure as "Tile H-Matrix".

The rest of the paper is structured as follows: in Section 2 we present in brief details about H-Matrices, the LU factorization, the H-LU algorithm, and its challenges; in Section 3 we present some related work, such as existing libraries, different compressed structures approaches, and runtime strategies to tackle H-Arithmetics. In Section 4 we describe H-Chameleon, giving details about the libraries that conform the basis of this work, our approach to represent the data, the original kernels from Chameleon and Hmat-oss we leverage, and how we construct the matrices. In Section 5 we describe the test case employed in the subsequently presented experiments, as well as the target platform in which we perform the executions. Finally, in Section 6 we expose the main remarks and conclusions we extract from this work.

Solution of H-Matrix Linear Systems

We open this section with a concise review of H-Matrices and H-Arithmetic. More details can be found in [START_REF] Hackbusch | A sparse matrix arithmetic based on hmatrices. part i: Introduction to h-matrices[END_REF][START_REF] Grasedyck | Construction and arithmetics of h-matrices[END_REF][START_REF] Hackbusch | Hierarchical Matrices: Algorithms and Analysis[END_REF]. Then, we detail the principle of the LU factorization for H-Matrices and how it can be parallelized with modern task based runtime systems.

A brief introduction to H-Matrices

From a theoretical perspective, H-Matrices are built upon block cluster trees that represent a hierarchical partitioning of the column and row index sets of the original matrix. Mathematically, this is summarized in Definitions 1 and 2.

Definition 1 Let I be an index set with cardinality n = #I. The graph T I = (V, E), with vertices V and edges E, is a cluster tree over I if I is the root of T I and, for all v ∈ V , either v is a leaf of T I or v = ∪v ∈S(v) v , where S(v) denotes the set of sons of v. Then, considering the node b = p × q of the cluster tree T I , the block cluster tree T I×I over T I can be defined recursively for the node b, starting with the root I × I, as follows:

S(b) = ∅ if b is admissible, S(p) = ∅ or S(q) = ∅, S otherwise,
where S := { p × q : p ∈ S(p), q ∈ S(q) }.

The admissibility condition applied to decide whether to re-partition a block ensures that an admissible block (that is, one that is not re-partitioned) can be approximated by a low-rank block. Mathematically, H-Matrices are equivalent to a collection of low-rank blocks with a maximal rank k, as exposed in Definition 2.

Definition 2

The set of H-Matrices for a block cluster tree T I×I over a cluster tree T I is defined as:

H(T I×I , k) := { M ∈ R I×I | ∀p × q ∈ L(T I×I ) : rank(M | p×q ) ≤ k ∨ {p, q} ∩ L(T I ) = ∅ },
where L(T I ) is the set of leaves of T I and k ∈ N.

Inria

In practice, most of the operations involving H-Matrices need to truncate the rank of low-rank blocks via the Singular Value Decomposition (SVD) [START_REF] Golub | Matrix Computations[END_REF] (cf Section 2.2), in order to maintain the log-linear arithmetic complexity. Alternative solutions exist to reduce the cost of this operation by computing approximated low-rank compression such as Adaptive Cross Approximation (ACA) [START_REF] Rjasanow | Adaptive cross approximation of dense matrices[END_REF] or randomized techniques [START_REF] Martinsson | Randomized numerical linear algebra: Foundations & algorithms[END_REF].

LU factorization of general and H-Matrices

The LU factorization is an extremely useful decomposition for the solution of linear systems that, given a non singular matrix A ∈ R n×n (for the coefficients of the linear system), returns a unit lower triangular factor L ∈ R n×n , and an upper triangular factor U ∈ R n×n , such that A = LU .

There exists different variants of the algorithm to compute the LU factorization of a full-rank matrix. Over the years, the factorizations have evolved from blocked algorithms working by panels (LAPACK [START_REF]LAPACK Users guide[END_REF]), to the modern and now common tiled-algorithms, used in many libraries exploiting runtime systems to schedule the computation on the target architecture.

Considering a matrix A of n t × n t tiles/blocks with each block A ij ∈ R n b ×n b , the tiled LU factorization algorithm can be performed following the blocked Right-Looking (RL) variant [START_REF] Golub | Matrix Computations[END_REF], which is exposed in Algorithm 1. It consists of an outer loop that computes, at each iteration k, the LU factorization of the corresponding diagonal block A kk (xGETRF LAPACK operation); next it computes the k-th block column panel (from L) and the k-th block row panel (from U ) by solving several triangular systems (TRSMs); and then it updates (via the matrix-matrix product or GEMM) the remaining blocks of A with respect to the block panels of L and U calculated in this iteration. (L kk , U kk ) ← GETRF(A kk )

Algorithm 1 Tiled GETRF algorithm

3: for i = k + 1, k + 2, . . . , nt do 10:

for j = k + 1, k + 2, . . . ,
for j = k + 1, k + 2, . . . , nt do 11: Aij ← GEMM( L ik , U kj , Aij) 12:
end for 13:

end for 14: end for

The task-based version of Algorithm 1 identifies each GETRF, TRSM and GEMM kernel as a task, which can be enriched by keywords to specify the inputs and outputs of each kernel. With this additional information and the sequential execution of this algorithm, runtime systems are able to automatically infer data dependencies between the tasks and generate a direct acyclic graph (DAG) of the tasks in the algorithm. Each node refers to a task, and each edge to a data dependency. The DAG of a full-rank LU factorization for a tiled matrix of 3 × 3 is illustrated in Figure 1.

When moving to H-Matrices, the algorithm for the H-LU factorization can be directly constructed from Algorithm 1 by adding recursion into the kernels. The algorithm is applied on the first level of the hierarchy of the H-Matrix where the partitioning discussed previously defines the number of tiles. Then, the three kernels: GETRF, TRSM and GEMM are replaced by their H-Arithmetic versions, which recursively apply the tile algorithm for the partition of the lower level.

In the case of H-GETRF, this consists simply in recursively calling the Algorithm 1 in line 2, as long as the inout block is an H-Matrix. When reaching a leaf, the traditional LAPACK kernel is called.

=

C

A B The H-TRSM and H-GEMM are a little more complex, as they may involve matrices with different H-Structure. In the case of H-GEMM, with 3 matrices involved and 3 possible formats for each (low rank, full rank or subdivided), 27 different configurations exists. The simplest case occurs when all three matrices are subdivided: we then progress in the H-Matrix structure to the lower level. All the remaining cases have to be handled individually. Figure 2 illustrates the type of complex configuration that must be dealt with.

In this work, as detailed in section 4.1, we chose to rely on an existing implementation of the H-Matrix algorithm in order to skip this hassle: concretely, we picked the Hmat-oss [START_REF]Hmat-oss library github repository[END_REF] library, an open source sequential package for H-Arithmetic. Hence, in Algorithm 1, from Chameleon's point of view, the tiles manipulated become H-Matrices created by Hmat-oss, and the calls to LAPACK GEMM, TRSM and GETRF are replaced by calls to the Hmat-oss implementations of H-GEMM, H-TRSM and H-GETRF.

Inria 3 Related Work

Over the past few decades, several libraries have been made available for H-Matrices and H-Arithmetic. Starting from sequential approaches, shared memory or distributed memory versions have been proposed, as well as some solutions to exploit GPUs. Those works have demonstrated that implementing a fast and efficient library for general H-Matrix format is not an easy task, and as a consequence they often proposed solutions to simplify the problem.

We focus here on the libraries implemented with the task-based paradigm model on top of runtime supports. The key feature of the H-Matrices resides in the recursive structure of the matrix that both saves flops and memory. However, none of the modern runtime systems for task-based algorithms propose simple and easy to use features to describe recursive algorithms. Many works have thus explored different techniques to bear with this issue, starting by exploiting OpenMP task to increase the performance for multicore architectures [START_REF] Izadi | Hierarchical matrix techniques on massively parallel computers[END_REF][START_REF] Kriemann | H-lu factorization on many-core systems[END_REF][START_REF] Aliaga | Task-parallel lu factorization of hierarchical matrices using ompss[END_REF]. These solutions realized a bulk-synchronous parallelism that was limited by synchronizations at each level of the H-Structure. This was eventually extended to an MPI+OpenMP version in [START_REF] Ida | Parallel hierarchical matrices with adaptive cross approximation on symmetric multiprocessing clusters[END_REF].

To alleviate the limitations of bulk-synchronous parallelism, [START_REF] Lizé | Résolution directe rapide pour les éléments finis de frontière en électromagnétisme et acoustique: H-Matrices[END_REF] proposed an implementation of the HMAT library on top of the StarPU runtime system by enumerating the required dependencies for each submitted task. This solution may end up with a large number of dependencies that the runtime system must deal with. In H2Lib [START_REF]H2lib library home page[END_REF] and HLibPro [28], the authors have been working on a semi-automatic DAG generation to generate this set of dependencies while being able to trim the DAG from unnecessary edges [START_REF] Boerm | Semi-automatic task graph construction for h-matrix arithmetic[END_REF]. In [START_REF] Augonnet | A hierarchical fast direct solver for distributed memory machines with manycore nodes[END_REF], the authors introduce fake dependencies while submitting the DAG similar to anchors they attach to only if needed. This helps them to recursively introduce the missing dependencies, while not overloading the runtime systems with too many of them. Similarly to this work, the authors proposed in [START_REF] Carratalá-Sáez | Exploiting nested task-parallelism in the h-lu factorization[END_REF] an implementation of H2lib on top of OmpSs, where a fair parallel performance is achieved thanks to novel features included in OmpSs-2 programming model, such as weak dependencies and early release. These features have been made available through the Nanos++ [START_REF] Perez | Improving the integration of task nesting and dependencies in openmp[END_REF] runtime, which extends OpenMP.

Other approaches have looked into different compression formats that could better fit the taskbased programming model by compromising lower memory and/or flops benefits in exchange for simpler programming and better efficiency. The most common one is the Block Low-Rank [START_REF] Amestoy | Improving multifrontal methods by means of block low-rank representations[END_REF] (BLR) format. This partitions the matrix into blocks which can be either in a low-rank representation or full-rank. This provides a regular tiling of the matrix easy to parallelize following all the previous works made in dense linear algebra. The Block Separable (BS) [START_REF] Cheng | On the compression of low rank matrices[END_REF][START_REF] Gillman | Fast direct solvers for elliptic partial differential equations[END_REF][START_REF] Gillman | A direct solver with o(n) complexity for integral equations on one-dimensional domains[END_REF] format is very close to BLR but, as opposed to it, it uses a weak admissibility condition such that all off-diagonal blocks are low-rank. The main advantage of those last two formats is that they remove the complex hierarchical block clustering of the H-Matrices. In [START_REF] Amestoy | Bridging the gap between flat and hierarchical low-rank matrix formats: the multilevel blr format[END_REF], a new intermediate format called multilevel BLR, was introduced to bridge the gap between H-Matrices and BLR, by reaching similar asymptotic costs while providing a regular structure easier to maintain and parallelize. However, there is still no efficient implementation of this format.

Finally, the Lattice H-Matrices, or Tile-H format, proposed in [START_REF] Ida | Lattice h-matrices on distributed-memory systems[END_REF][START_REF] Yamazaki | Distributed-memory lattice h-matrix factorization[END_REF], presents another alternative format to attain better asymptotic costs with respect to BLR, while presenting a similar regular structure. The idea relies in flattening the first levels of the clustering tree to obtain an n t × n t tile layout where each tile is itself an H-Matrix. In [START_REF] Yamazaki | Distributed-memory lattice h-matrix factorization[END_REF], they proposed an OpenMP + MPI implementation of the LU factorization with this format, developed in the HACApk library [START_REF]Hacapk library[END_REF]. The parallelization model follows the one developed in the dense linear algebra library Slate: OpenMP tasks working on each tile are submitted to perform both the computation and the MPI communications, and computations within each tiles are done with the sequential version of HACApk. This format has also been tested in libHLR [38] and in [START_REF] Augonnet | A hierarchical fast direct solver for distributed memory machines with manycore nodes[END_REF].

The main differences between our approach and the existing ones are that we 1) we focus on the tiling idea to be able to expose a higher concurrency degree as BRL/BS do, but without losing efficiency, because each tile contains an H-Matrix; 2) we leverage existing libraries (Chameleon and Hmat-oss) instead of reimplementing the whole H-Arithmetic, contrarily to what is done in the referenced Lattice H-Matrices approach; 3) as a result from using StarPU combined with Chameleon, we already have tradidional task based parallelism arising from blocks/tiles-based structures, so we are able to attain a fair parallel performance in this context with an already existing and efficient runtime.

H-Chameleon

The main goal of this work is to validat whether it is possible to ontain an efficient open-source H-Matrices library at a small developer cost. For that purpose, we leverage existing libraries, software, programming models and runtime systems to efficiently scale H-Matrices solvers, via the classic tile-based approach commonly used in dense linear algebra, in order to deliver good parallel performance. In this section, we expose the main details of the proposed solution.

The building blocks

The work of this paper relies on three building blokcs or components to provide this new H-Matrix library: the dense linear algebra library Chameleon, the runtime system StarPU, and the H-Matrices library Hmat-oss. Chameleon [START_REF] Agullo | Faster, cheaper, better -a hybridization methodology to develop linear algebra software for gpus[END_REF][START_REF] Agullo | Achieving high performance on supercomputers with a sequential task-based programming model[END_REF] is a dense linear algebra open source software written in C. It relies on the sequential task flow programming model supported by runtime systems such as OpenMP tasks, StarPU, PaRSEC-DTD, or Quark. All dense linear algebra algorithms are expressed as tile algorithms with sequential tasks that are submitted to the underlying runtime system. The runtime automatically infers the data dependencies in these algorithms thanks to keywords that specify the data accesses. Chameleon covers all BLAS subroutines with these algorithms, as well as one-sided factorization (Cholesky, LU, QR), and supports multiple runtime systems to schedule the tasks. StarPU is one of them, and the most integrated one into Chameleon, besides being the second building brick of our proposal.

StarPU provides tools to describe the pieces of data of the user, such that the data transfers from device to device or node to node, are transparent to the algorithm developer. This is a clue feature when developing H-Matrices algorithms, as it will suffice to provide the runtime with pack/unpack functions in order to transfer data (if needed).

Hmat-oss is the open source version of the Hmat library developed by Airbus. It is a sequential library written in C++. This library provides up-to-date implementations of H-Matrix operations and techniques, clustering algorithms and orderings, and real world applications samples. The non-public Hmat library will be employed as a performance reference or baseline in the experiments in Section 5.

Tile-H matrices

The baseline realization of Chameleon only supports dense matrices, which are stored employing a descriptor that specifies, among other information, the matrix dimensions, number and size of the tiles which define the matrix partitioning, pointers to data addresses in memory, and some control parameters to test certain features (e.g., whether there is an overall data pointer).

In order to accommodate Tile-H Matrices in our hierarchical extension of Chameleon, we have expanded the reference descriptor structure to become a collection of tiles. Each of these tiles is potentially an H-Matrix, a low-rank block, or a full-rank matrix. To do this change, Inria we enriched the Chameleon matrix descriptor (Structure 1) with an array for the new tile structures, and a helper function, get_blktile, to extract the correct tile pointer from the tile indices in the matrix. Additionally to the main data structure, the tile description has been modified from the simple data pointer used in runtime systems, such as OpenMP or Quark, to handle the dependencies toward a more complex structure able to handle different data formats. The new datatype, CHAM_tile_t, described in Structure 2, allows to simply store different matrix formats, defined by the format field, and a pointer to the matrix, mat, which can be either a full-rank matrix, or a more complex datatype, such as an H-Matrix coming from an external library. The data dependencies that were initially tracked down using the pointer to the data in the full-rank matrix are thus now followed by the pointer to this tile descriptor. This means that all the algorithms from the Chameleon library could work out of the box with an H-Matrix format, considering that there exist kernels to handle this type of matrices. Finally, the global representation of the matrix that links a Chameleon descriptor (CHAM_desc_t) with an Hmat-oss descriptor (hmat_matrix_t) in a unique element is done via an additional new data structure (Structure 3). In this new datatype, super represents the Chameleon library tile descriptor; clusters stores an array of the cluster trees created to partition the original data; admissibilityCondition contains this parameter value, necessary to determine whether a certain block is already admissible (and consequently converted to a low-rank block) or needs to be re-partitioned; hi is the interface defined in the Hmat-oss library to deal with H-Matrices, which contains general information (for example, the clustering algorithm needed to construct the H-Matrix; and data precision format); hmatrix contains the H-Matrix built employing Hmatoss construction kernels (which is the descriptor of the library whose content will be employed in the Hmat-oss kernels operations); and perm stores the permutation array.

Clustering algorithm

In order to use Tile-H matrices in the Chameleon library, we need an adapted clustering tree. Indeed, the Chameleon library works exclusively on regular tile sizes, with the exception of the padding row and column. Thus, it is not sufficient to flatten the first levels of the clustering tree as it is done in [START_REF] Yamazaki | Distributed-memory lattice h-matrix factorization[END_REF]. We extended the Hmat-oss library with a recursive tile clustering algorithm, named "NTilesRecursive", that recursively divides a given cluster tree (CT) into clusters that respect the regular partitioning into tiles of size N B.

This process is illustrated in Algorithm 2, and the parameters and functions employed are described next:

• CT is the cluster tree to partition;

• axis is the main axis of the current slice, which is exploited by the geometric clustering techniques to split along the largest dimension;

• NB is the desired tile size;

• offset represents the coordinates (values) of the first value in the current cluster tree;

• size is the size of the current cluster tree;

• the function slice returns a portion of the current cluster according to the given offset and size;

• the function largestDimension returns the largest dimension in the current cluster;

• and the function sortByDimension orders the current cluster of unknowns accordingly to the given dimension;

At each level, this function performs a pseudo-bisection aligned with the tile size along the largest dimension and returns the concatenation of the recursive call to each subset of unknowns. This provides a regular clustering of the unknowns that matches both the constant size requirement of the Chameleon library and the Tile-H format. A median bisection algorithm is then called within each cluster to refine the clustering of each tile.

Algorithm 2 NTilesRecursive( CT, n b , offset, size, axis)

1: n t = size N B 2: if n t == 1 then 3:
return CT 4: end if 5: dim =getLargestDimension(CT, axis) 6: sortByDimension(CT, dim)

7: offset L = offset 8: size L = N B * nt 2 9: offset R = offset + size L 10: size R = size -size L 11: CT L = slice( CT, offset, offset + size L ) 12: L = NTilesRecursive( CT L , N B, offset L , size L , dim) 13: CT R = slice( CT, offset + size L , size R ) 14: R = NTilesRecursive( CT R , N B, offset R , size R , dim) 15: return (L, R) Inria

Hmat-oss kernels

In order to implement a hierarchical LU factorization, we have leveraged a number of sequential numerical kernels from the Chameleon and Hmat-oss libraries, which provide the necessary operations to factorize our Tile-H matrix. In order to do these changes, we modified the main kernels of the Chameleon libraries. These kernels, as the data dependencies tracking system, integrated pointers to full-rank matrices. To limit the impact on the library, we just introduced an intermediate layer to enable the switch between full-rank and H kernels. Thus, the task insertion functions are not modified by this change, and the CHAM_tile_t datatype helps us to switch from one kernel type to another, thanks to the format field.

On the Hmat-oss library side, we provide a similar interface to BLAS for GETRF, TRSM, and GEMM operations, and an intermediate internal layer which allows us to switch from this H-Matrix interface to the more classical BLAS interface.

Numerical Experiments

In this section, we first present the environment from which the data comes from, and a description of the PlaFRIM platform used for the experiments. Afterwards, we analyze the parallel performance of H-Chameleon in a multicore system.

Experimental context

The test case used is the application Test_FEMBEM [39] that generates a real or complex matrix, which has features similar to real industrial applications in aeronautics. For any number of unknowns n, we create a cloud of points (x i ) 1≤i≤n located on the surface of a cylinder of chosen height and width, as illustrated on Figure 3. These points are equally spaced in both directions on the cylinder surface. Then, we define the interaction kernel between two points x i and x j separated by a distance d = |x i -x j | as K(d) = exp(ikd)/d in the complex case, and K(d) = 1/d in the real case. In the complex case, k plays the role of a wave number, it is chosen with the "rule of thumb" of having 10 points per wavelength, which is a rule commonly used in the wave propagation community. The singularity at d = 0 is simply removed by setting d equal to half the mesh step in that case.

In the real case, the rank of the H-Matrix blocks will be mostly independent of their sizes, and therefore most of the data (in terms of storage) will be located near the diagonal of the matrix. On the other hand, in the complex case, the rank will grow with the size of the blocks, and the data will be much more evenly distributed in the matrix. Hence, the amount of storage and work is a lot more important in the complex case, and the work distribution is much more difficult too.

Experiments platform

All our experiments have been performed on the PlaFRIM [40] test bench, and more specifically on the bora cluster. Each node is equipped with two Intel Xeon Skylake Gold 6240 of 18cores, running at 2.60 GHz, and equipped with 192 GB of memory. The application is compiled with GCC 9.2.0, and Intel MKL 2019 is used for the BLAS and LAPACK kernels. The 1.3.0 version of StarPU is used, and our proposal is built on revisions 33aa719 of Hmat-oss [START_REF]Hmat-oss library github repository[END_REF] and 08cf0cd1 of Chameleon [START_REF]Chameleon software home page[END_REF].

Experiments

We next present the experiments performed to evaluate H-Chameleon.

Our first experiment pursues to highlight that our implementation, though simpler than constructing a classical H-Matrix, still enables a good compression ratio. To this end, Figure 4 shows a comparison between Hmat-oss (dashed lines) and H-Chameleon (full lines) compression ratios for double (left plot) and complex (right plot) precisions, employing different matrix dimensions, from 10K to 200K, and various tile sizes, from 500 to 10K. Th results show that the difference is negligible in all cases, so that we can affirm that the clustering with fixed tile sizes does not impact the compression ratio on the studied test case, and can even provide better results than the classical median bisection used in Hmat (-oss).

Second, as precision is the bargaining chip in H-Scenarios to permit time and memory savings, it is also necessary to control that the proposed clustering does not impact the numerical accuracy of the H-Chameleon H-LU factorization operation. With the purpose of evidencing that precision differences are negligible, Figure 5 presents forward error measurements: ||x -x 0 || f /||x|| f , for different H-LU executions with the same matrix configurations employed in the previous experiment on the compression ratio. Note that the accuracy parameter is set at 10 -4 , both in Hmat and H-Chameleon. The largest observable differences are around 1.5e -4, which means we stay in the same magnitude order.

Figures 6 and7 offer multicore parallel performance comparisons (employing up to 36 threads), with various matrix dimensions from 10K to 200K, both in real and complex double precision scenarios. These figures study different scheduling strategies proposed by the StarPU runtime system, while comparing their performance to that provided by the StarPU based implementation of the Hmat library that deals with all the fine grain dependencies.

Three scheduling strategies are studied:

• The work stealing strategy (ws) uses a queue per worker and schedules the tasks on the worker which released them by default. Whenever a worker becomes idle, it steals a task from the most loaded worker.

• The locality work stealing (lws) strategy similarly uses a queue per worker that is now sorted by the priorities assigned to the tasks. New ready tasks are scheduled on the worker Tile size (NB) Compression ratio wrt to full-rank N q q q q q q q q 10000 20000 30000 which released them by default. Whenever a worker becomes idle, it steals a task from neighbor workers respecting the priority order.

• The priorities based (prio) approach uses a central task queue in which ready tasks are sorted by decreasing priority. All threads try to pull work out of this central queue.

Note that for our proposed implementation we never go over 35 worker threads to keep a dedicated core to the task submission. Experiments have shown that this was more efficient than oversubscribing with 36 worker threads.

In the case of the H-Chameleon implementation, all tile sizes presented in compression and accuracy curves have been tested and we chose the best one for each dimension and precision.

In most of the test cases, we observe that H-Chameleon presents slower execution times when using 1, 2, or 3 threads. The tile size being optimized for the 36 threads case induces an overhead of memory and required flops, which impacts the executions with a low number of threads. However, when a larger number of threads is used, this is compensated by the parallelism it provides, and so it enables a good scalability of the library. Hmat is not impacted by the tile size and manages to get better performance for the execution on the small numbers of threads.

The comparison of the real and complex double precision results shows that the Hmat performs better on the complex cases, while H-Chameleon has a better scaling on the real ones. This can be explain by the difference in the number of operations of the two test cases, due to the arithmetic and to their configuration, as previously explained in Section 5.1. In the complexarithmetic scenario, the cost of the kernels is high enough to cover the overhead of handling the large number of dependencies, which benefits Hmat. However, in the real case, the cost of handling all fine grain dependencies become too important with respect to the computational tasks, and therefore H-Chameleon outperforms Hmat.

The comparative study of the StarPU scheduling strategies in these figures reflects that, in general, the three variants deliver similar execution times. However, the strategies based on priorities provide higher performance, and the simple priority strategy turns to be the best in most of the cases, except the smaller dimensions. In the real double precision cases with N = 10K and N = 20K, the priority scheduler does not provide the fastest solution, as the computational tasks are too small and the idle threads create a contention in the single global task queue of this scheduler.

Concluding Remarks

In this paper, we have proposed an extension of the Chameleon library that takes advantage of H-Matrices and H-Arithmetic to accelerate the execution time and reduce the memory footprint of the LU factorization. More precisely, our approach takes advantage of the sequential kernels of Hmat-oss to perform H-Arithmetic, and the task-based approach of Chameleon to manage parallelism. The original large matrix is split into a set of tiles, and each tile can be represented either as a dense or an H-Matrix. Then, a runtime system, as StarPU does in our experiments, schedules the tasks on the system, and manages the parallelism, following the approaches developed in Chameleon for the dense case.

We have conducted experiments on a shared memory machine for large real-life cases, and our results have demonstrated that this approach is competitive with the proprietary Hmat library. Thus, it provides one of the first open-source library that is able to reach a good level of performance using H-Matrices. A main asset of the approach is that it will directly benefit from all improvements on the runtime side, as it is now integrated in the Chameleon library.

For future work, we plan to study the behavior of this approach for the distributed case, where the main challenge is to correctly handle communications, when the size of the structures, Inria depending on the ranks of matrices, cannot be known statically. It is particularly challenging to study such an approach with the use of dynamic runtime systems. The distributed H-Matrices implementations are also known to be largely unbalanced, and this work will provide a large test suite to work on data distribution and distributed load-balancing algorithms.

  Require: A ∈ R n×n with n = nt * n b 1: for k = 1, 2, . . . , nt do 2:

Figure 1 :

 1 Figure 1: DAG of the tasks that conform a full-rank LU factorization (following Algorithm 1) for a tiled matrix of 3 × 3 tiles.

Figure 2 :

 2 Figure 2: Example of H-GEMM configuration. The dark gray part of C is updated by contributions coming from both light gray and dark grays parts of B and C, with some of these blocks being full rank and others being low rank.

Listing 1 :

 1 typedef struct c h a m e l e o n _ d e s c _ s { .CHAM_desc_t datatype modifications to handle more generic tile formats.

Listing 2 :

 2 typedef struct c h a m e l e o n _ t i l e _ s { int8_t format ; i n t m , n , ld ; void * mat ; } CHAM_tile_t ; CHAM_tile_t data structure to accommodate any format of tiles in the Chameleon library.

Figure 3 :

 3 Figure3: Illustration of the test case used for the experiments. On the left, the mesh of the cylinder with the distribution of the unknowns on the surface for 10k points. On the middle, the associated compressed real matrix in the HMAT format (classical H-Matrix). On the right, the associated compressed real matrix in the proposed fixed-sized Lattice or Tile-H based matrix format. In the matrices, low-rank blocks are represented in green (with a number specifying the rank), while dense blocks are coloured in red.

Figure 4 :

 4 Figure 4: Comparison of the compression ratio between Hmat-oss original clustering algorithm (dashed line) and H-Chameleon (full lines), in function of the tile size, for double precision (left) and complex precision (right).

Figure 5 :

 5 Figure 5: Comparison of the precision of the solver (based on forward error) between Hmat-oss original clustering algorithm (dashed line) and H-Chameleon ntiles clustering algorithm (full lines) in function of the tile size, for double precision (left) and complex precision (right).

Figure 6 :

 6 Figure 6: Comparison of the multicore parallel executions between H-Chameleon and Hmatoss LU factorization, employing up to 36 threads. Results are shown for real double (top) and complex double (bottom) precision, and for small matrix dimensions: 10K (with N B = 250 for d, N B = 500 for z), 20K (with N B = 500), 40K (with N B = 1000).

Figure 7 :

 7 Figure 7: Comparison of the multicore parallel executions between H-Chameleon and Hmatoss LU factorization, employing up to 36 threads. Results are shown for real double (top) and complex double (bottom) precision, and for larger matrix dimensions: 80K (with N B = 1000 for d, N B = 2000 for z), 100K (with N B = 1000 for d, N B = 2000 for z), 200K (with N B = 2000 for d, N B = 4000 for z).

RR n°9327

avenue de la Vieille Tour 33405 Talence Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgments

Experiments presented in this paper were carried out using the PlaFRIM experimental testbed, supported by Inria, CNRS (LABRI and IMB), Université de Bordeaux, Bordeaux INP and Conseil Régional d'Aquitaine (see https://www.plafrim.fr/). R. Carratalá-Sáez and E. S. Quintana-Ortí were supported by project TIN2017-82972-R from the Ministerio de Ciencia, Innovación y Universidades, Spain. R. Carratalá-Sáez was supported by the FPU program of the Ministerio de Educación, Cultura y Deporte, and a HiPEAC collaboration Grant (2019).

This work is also supported by the Agence Nationale de la Recherche, under grant ANR-19-CE46-0009.


Volker Zorngiebel
email: zorngiebel@isl.tm.fr

Sigo Scharnholz

Emil Spahn

Pierre Brosselard

N Arssi

Jean-Pierre Chante

Dominique Planson

B Spangenberg

H Kurz

J.-P Chante

Fabrication and Characterisation of High Voltage SiC-Thyristors

Keywords: SiC, high power devices, thyristors, GTO, EGR, MESA, JTE

come

Introduction

The physical properties, above all its high electric field, good thermal conductivity and high temperature capability make SiC to be the top choice for high power semiconducting switches [START_REF] Chow | SiC and GaN bipolar power devices[END_REF]. These properties could be used to design higher voltage devices with higher current densities in order to reduce the number of devices used in a switching system. On the basis of numerical simulations published by N. Arssi et al [START_REF] Arssi | Study based on the Numerical Simulation of a 5 kV Asymmetrical 4H-SiC Thyristor for High Power Pulses Application[END_REF] we designed a thyristor in reference to the studies for future pulsed power applications as reported by S. Scharnholz et al. [START_REF] Scharnholz | High-Power MOSFETs and Fast-Switching Thyristors Utilized as Opening Switches for Inductive Storage Systems[END_REF] and E. Spahn et al. [START_REF] Spahn | Modular Pulsed Power Supply Units for active Protection Systems[END_REF]. In order to define our devices, we looked about what is the state of the art for high power silicon thyristors.

Fabrication

Device Layout. Starting the simulation with a fixed point for the breakover voltage the structure for the thyristor was defined as in Fig. 1 to sustain a VBO of 5 kV. The first simulations were done without including the periphery into the calculations. We only looked to the pn-junction which has to block the voltage. So, with an 35 µm thick low doped epitaxial blocking layer, different kinds of periphery protection were calculated like MESA, EGR and JTE. The first device we realized had epitaxial guard rings (EGR's). In Fig. 2 shows a picture of the device displaying the layout of the structure with interdigitated gate and anode areas. The structures realized on the masks vary from 170 to 330 µm in length and 35 to 50 µm in width Moreover, the number of fingers was varied and so the total anode area ranges from 0.7x10 -3 to 5.7x10 -3 cm 2 . The sizes of the gate structures are varying in the same range as the anode. The number of guard rings was set to five with a width of 12 µm and a distance of 2 µm. The guard rings are terminated by a trench with a width of 45 µm. These parameters were also defined by the simulation. For the MESA protected devices we used the same variations for the size of the structures. The distance from the edge of the MESA to the anode structure was set to 13 µm as a minimum. Process Flow. For the realization of the devices we used three different sets of masks, one for each kind of termination. In principle the process flow of all three terminations has the same sequence. The first step was to etch the trench for the gate contacts through the top epi-layer. The mask was a Ti mask or alternatively a Ti/Cr-mask for the ECR RIE etching. With this step the structure of the anode and also the alignment marks for the following steps were defined. The second step was to etch the termination area. We used the same etching parameters as before and the same mask metal. Fabricating the devices with EGRs we tried to stop etching accurately on the p drift layer. As contrast to this for the MESA and for the JTE devices we etched deep into the p layer. In the finished devices we obtained an over etching for the EGRs of 150 nm and a MESA of 1,3 µm into the p layer. The simulation shows that the breakover voltage is enhanced with increasing the MESA depth. Limited by the technical equipment we are not able to etch more than 1,3 µm deep in one step. Also as a result from the etching the sidewall angle of the structures vary from 93° to 98°. The third step was to define the metal contacts for anode and gate after deposition of the passivation oxide. The oxide was deposited in two different ways. For the first samples we deposited 800 nm oxide by PECVD and for the second series a 1600 nm thick oxide was deposited by CVD at 900°C. The oxide was annealed in two steps with a final dry reoxidation as described from Scharnholz et al. [START_REF] Scharnholz | Utilization of RPECVD oxide in U-MOS structures and different MOSFETs on SiC[END_REF]. The metal contacts for the anode consisted of a Al/Ti double layer. For the gate and the cathode we used a Ti/Ni double layer. The contacts where annealed at 1000°C for 180 s in an RTA furnace.

For the JTE there was an additional step for the implantation which required one more mask layer to define the implantation area in a deposited protection oxide.

Up to now we realised thyristors with an EGR and a MESA termination. The devices terminated by a JTE and a MESA with larger etching depth are still in process.

Characterisation

Blocking-state. The results shown in Fig. 3 are measured on a sample of the first series with devices having a MESA protection. In this particular case the etched MESA was 1,3 µm in depth.

The maximum breakover voltage VBO measured on this sample was almost 1500 V. The measurement was done in oil to prevent a breakdown in air close to the surface of the devices. Because of the technological process flow, there was a small crack between the metal contacts and the passivation oxide. This was changed for the second series of samples. Manuscript for Proc. 4th European Conf. on SiC and Related Materials (ECSCRM), Linköping, Sweden, 1-5 Sep. 2002 N + = 7 10 18 cm -3 In a second case we measured the breakover voltage VBO on samples with EGRs. This time the measurement was done with the samples put into perfluorinated fluid Galden™ HT270. As to be seen in Fig. 4 the breakover voltage amounts up to about 1900 V for these devices. The measurement was stopped when the leakage current reached 1 µA to prevent from a premature destruction. Switching. Applying a gate current pulse with a length of 1 µs an amplitude of at least 27 mA was necessary to switch on the devices. Fig. 5 we show the gate pulse and the corresponding signals of anode-cathode voltage and the cathode current. Using the same pulse we could switch on the devices with anode-cathode voltages ranging from 45 V to 300 V. In all these cases the cathode current was limited to 100mA with series resistances ranging from 0,4 to 3 kW. The switch-off of the devices was investigated by using a gate pulse of 1µs and a variable amplitude. With an anode cathode voltage of 100 V we switched off cathode currents of more than 60 mA. We also investigated the influence of the duration of the gate pulse on the switch-off behaviour. Therefore we vary the pulse length from 1µs up to 1 ms. To switch off a device with a anode cathode voltage of 100 V and a cathode current of 100 mA we need a gate pulse amplitude of 61,5 mA for the short pulse and 43,5 mA for the long pulse.

Conclusion

Based on simulation results thyristors on SiC with different terminations were fabricated. We characterised the devices according to the breakover voltage as well as the switch-on and -off behaviour. The breakover voltage attains values up to 1900 V for devices terminated by epitaxial guard rings. Investigations concerning the switching behaviour showed that the devices could be turned on and off by short gate pulses. Besides this the experimental results show a good correlation to the simulations. Based on this further simulations were performed showing that high power thyristors are feasible with the chosen design and process flow. Therefore future work will focus on the realisation of such high power devices on SiC.

Fig. 1 :Fig. 2 :

 12 Fig.1: Half cell of the thyristor with the doping Fig. 2: Device with EGR termination (5 rings Level and thickness given by Creeä of 12µm width and a distance of 2 µm)

Fig. 3 : 4 :

 34 Fig. 3: Simulation of the dependence of the breakover Fig. 4: Measurement of the breakover voltage with IG = 0mA voltage from MESA depth on different devices at room temperature

Fig. 5 :

 5 Fig. 5: Switching on by different anode cathode voltage and constant current


domain: spi phys

text:
Hubert Guillaume
email: guillaume.hubert@onera.fr

Ricaud Philippe
email: philippe.ricaud@meteo.fr

Four Years of Measurements of the Energy Spectrum of Cosmic-Ray-induced-neutrons on the Concordia Antarctic Station Taking into Account Environmental and Systematic Effects

Neutron spectrometers are operated simultaneously in the framework of a remoted network. The CHINSTRAP, supported by IPEV, aims at recording neutron spectra at the Concordia station (+3233m) in Antarctica since December 2105. The neutron spectrometer measures the neutron spectrum over a wide energy range from meV up to tens of GeV with a short time resolution. Several parameters can influence the measurement: the atmospheric pressure, the hydrometric environment close to the instrument (snowfall) and the atmospheric water vapor. The atmospheric water vapor dominates the ground cosmic ray variations within the first few hundred meters above the ground. The HAMSTRAD polar project is dedicated to measure of the trends in water vapor and temperature profiles from the lower part of the troposphere to the lower part of the stratosphere, using a radiometer. This paper proposes to present data analyses, including short and long-term analyses (first four years) and taking into account environmental and systematic effects. In particular, study focusing on Forbush Decreases and Ground Level Enhancement occurred in September 2017 is proposed.

Introduction

The CHINSTRAP (Continuous High-altitude Investigation of the Neutron Spectra for Terrestrial Radiation Antarctic Project), supported by the French Polar Agency (IPEV), aims at recording neutron spectra at the Concordia station (+ 3233m) since December 2105. Several parameters can influence the measurement: the atmospheric pressure, the hydrometric environment close to the instrument (snowfall) and the atmospheric water vapor. The atmospheric water vapor dominates the ground cosmic ray variations within the first few hundred meters above the ground. Moreover, water vapor contents were recorded in Concordia station in the framework of the HAMSTRAD (H 2 O Antarctica Microwave Stratospheric and Tropospheric Radiometers) project. Thus, two contributions were considered to correct count rate data, the air pressure and the atmospheric water vapor.

The four years of measurements performed in Concordia have just been completed. This paper proposes to present data measurements (count rates, energy spectrum, atmospheric water vapor impacts), including short and long-term analyses and taking into account environmental and systematic effects. In particular, study focusing on solar event occurred in September 2017 is proposed.

Description of the CHINSTRAP and HAMSTRAD projects

The aim of the HAMSTRAD project [START_REF] Ricaud | HAMSTRAD-Tropo, A 183-GHz Radiometer Dedicated to Sound Tropospheric Water Vapor Over Concordia station, Antarctica[END_REF] [START_REF] Ricaud | Review of tropospheric temperature, absolute humidity and integrated water vapour from the HAMSTRAD radiometer installed at Dome C, Antarctica[END_REF] is to measure the trends in water vapour and temperature profiles from the lower part of the troposphere to the lower part of the stratosphere and their links with climate change. The HAMSTRAD radiometer is a genuine state-of-the-art microwave instrument dedicated for the detection of 1) the 60-GHz oxygen line to measure tropospheric temperature profile, and 2) the 183-GHz water vapor (H 2 O) line to get tropospheric H 2 O. The radiometer has been deployed in the Concordia polar station in January 2010 and is working since then. Liquid Nitrogen calibration is performed one or twice per year (summer/winter periods). The instrument is fully automated and data are sent to France on a daily basis. The instrument measures absolute humidity and temperature from the surface to about 10 km altitude, together with integrated water vapour (IWV), with a time resolution of about 7 minutes. The Liquid Water Path (LWP) parameter is retrieved since 2013, and is sensitive to the presence of supercooled liquid water clouds above the station.

The aim of the CHINSTRAP project is to operate a neutron spectrometer at the Concordia station [START_REF] Hubert | Analyses of cosmic ray induced-neutron based on spectrometers operated simultaneously at mid-latitude and Antarctica high-altitude stations during quiet solar activity[END_REF]. The geomagnetic conditions (~ 0.001 GV) and the high altitude (3233 m above the sea level) of Concordia let us collect cosmic rays (CRs) data and observe solar particle events (SEP) such as Forbush Decreases (FD) [START_REF] Cane | Coronal Mass Ejections and Forbush decreases[END_REF] or Ground Level Enhancement (GLE). The neutron spectrometer is based on a Bonner multi-sphere system extended to high neutron energies (energy range from meV to GeV). As detailed in previous works, this system was composed of spherical 3He proportional counters surrounded with spherical PEHD (high density polyethylene) moderators with different thicknesses. Additionally, the spectrometer includes two PEHD spheres with inner tungsten and lead shells in order to increase the response to neutrons above 20 MeV. The output pulse signals are treated by a multi-channel analyzer, translated as spectrum or count. The response functions were used to convert the measured counting rates to neutron energy spectrum. Fig. 1 presents both instruments, which are operated in the physics shelter of the Concordia station.

Analyses of energy spectrum of CR induced neutrons in Dome C

Neutron detectors and spectrometers are usually used to monitor the dynamics of incoming cosmic-ray-induced neutrons of high energy (cascade neutron i.e. E > 20 MeV) under the assumption that local environmental conditions do not influence the neutron fluence rate. Count rates dynamics are mainly due to environmental and systematic effects including the atmospheric pressure variations. Indeed, the cosmic ray flux increases at higher altitude as air pressure and the shielding effect of the atmosphere decrease. More precisely, the vertical distribution of temperature, pressure, density, and composition of the atmosphere constitutes atmospheric structure. These quantities also vary with season and location in latitude and longitude, as well as from night to day; however, under the topic of atmospheric structure, the focus is on the average variations with height above sea level.

Other effect concerns the albedo neutron [START_REF] Hubert | Modelling of ground albedo neutrons to investigate seasonal cosmic-ray-induced neutron variations measured at high-altitude stations[END_REF] produced by the interaction of air-shower neutrons with the surface of the Earth. Hydrogen in soil, air and snow determines the amount of ground albedo neutrons in the sensitive energy range from 1 eV to 10 MeV. The local environment of the experimental shelter in Concordia is mainly composed by ice water. Precipitation levels are very low in time, providing stable conditions. Albedo contribution does not vary significantly in time and it is not an obstacle to exploit the full energy range of spectra.

PoS(ICRC2019)293

Environmental and systematic effects.

Fig. 2 presents the temporal evolution of the LWP and IWV measured by HAMSTRAD above the altitude of Concordia station from January 2016 to January 2019. In addition, the atmospheric pressure measured in Concordia during the same period using a pressure detector integrated in the neuron spectrometer. Among all contributions, the atmospheric pressure water vapor dominates the ground CR variations in Concordia [START_REF] Rosolem | The effect of atmospheric water vapor on neutron count in the cosmic-ray soil moisture observing system[END_REF]. Temporal evolution of the LWP (kg.m -2 ) and IWV (g.m -2 ) as measured by HAMSTRAD above the altitude of Concordia station from January 2016 to January 2019. Also shown are the atmospheric pressure measured in Concordia during the same period thanks to the pressure detector integrated in the neuron spectrometer. Thus, two contributions were considered to correct count rate data: the atmospheric air pressure and the water vapor. Thus, corrected count rate by air pressure C P,corr is given by the equation 1: 0 ( ) , .

P P P corr corr C C e β - - = ( 1
)
where β is the barometric coefficient (% hPa -1 ) empirically determined and C corr the corrected count rates. The β factor, equals to 7.1×10 -3 hPa -1 , was empirically determined in a previous work [START_REF] Leuthold | Long-term measurements of cosmic ray neutrons by means of a Bonner Spectrometer at Mountain altitudes -First results[END_REF]. As described in [START_REF] Rosolem | The effect of atmospheric water vapor on neutron count in the cosmic-ray soil moisture observing system[END_REF], the atmospheric humidity in the count rate data is taken into account using a single scaling factor C WV ,corr determined by simulations [START_REF] Rosolem | The effect of atmospheric water vapor on neutron count in the cosmic-ray soil moisture observing system[END_REF]:

( ) , 1 0.0143 W V corr C C IW V = × + × ∆ (2)
where ∆IWV corresponds to the difference of the integrated vapor density value with the reference value on the day of the calibration. Concordia environment is well known and stable over time (low precipitation levels). It is one of the driest and coldest places over the world with precipitable water or Integrated Water Vapor (IWV) less than 1 kg m -2 over the year reaching less than 0.1 kg m -2 in winter [START_REF] Rosolem | The effect of atmospheric water vapor on neutron count in the cosmic-ray soil moisture observing system[END_REF] with exceptionally wet and warm episodes (IWV >> 1 kg m -2 ) as e.g. in January 2017 and December 2019. The atmosphere is considered dry in the Dome C

PoS(ICRC2019)293

when the absolute humidity reaches 0.3 g.m -3 , which is very low compared to values observed in typical moist conditions (i.e. 20 g.m -3 at middle latitudes).

Data confirm that the Concordia environment can be therefore considered extremely dry. As shown in a previous work [START_REF] Rosolem | The effect of atmospheric water vapor on neutron count in the cosmic-ray soil moisture observing system[END_REF], for low atmosphere water vapor content and IWV (in the order of few tenths of g.m -3 or kg. m -2 for absolute humidity and IWV, respectively), environmental or systematic effects on fast neutrons are negligible (i.e. a few tenths of a percentage) while they are significant for low-energy neutrons (epithermal and thermal).

Modeling approach based on previous works

Some methods were proposed to model the secondary atmospheric CRs, some of which are based on nuclear transport simulations. The ATMORAD approach [START_REF] Hubert | Radiation effects investigations based on atmospheric radiation model (ATMORAD) considering GEANT4 simulations of extensive air showers and solar modulation potential[END_REF] combines the GEANT4 simulations of extensive Air Showers to primary GCR spectra define by the Force-Field Approximation [START_REF] Gleeson | Solar Modulation of Galactic Cosmic Rays[END_REF] which provides a simple parametric formula of the differential spectrum. Thus, the whole energy spectrum for protons and α-particles can be described using the solar potential φ. This parameter is a formal spectral index whose physical interpretation is not straightforward, especially on short time scales and during periods of active Sun. Several strategies have been developed for the reconstruction of time series of the modulation potential φ(t), mainly based on neutron monitor data and the yield functions [START_REF]Cosmic Ray Station[END_REF].

The secondary CR spectra can be determined as function of the altitude, the longitude, the latitude and the sun activity through the solar potential. In this work, ATMORAD was used to reconstruct the CRs temporal series in Concordia using the solar modulation issued from [START_REF]Cosmic Ray Station[END_REF]. The corrected measurements by water vapor and atmospheric pressure effects is compared with the calculated neutron flux using a reference solar potential.

Long-term analyses

Fig. 3 presents the total neutron flux corrected by water vapor and atmospheric pressure effects from January 2016 to January 2019. Experimental data are issued from CHINSTRAP and additional data by calculations using ATMORAD calculations and considering solar potential extract from Oulu data [START_REF]Cosmic Ray Station[END_REF]. It is very interesting to note the relevance between corrected measurement and calculations, implying that coupled data issued from HAMSTRAD and CHINSTRAP allow for investigating the dynamic of neutron spectrum in the energy range from meV to GeV.

The variations are attributable to the solar activity. These trends are also observed on most NMs. The impact of the solar cycle (24 th ) is particularly identifiable on this result. Indeed, this increase reflects the transition from strong to low solar activity. Cycle 24 reached its maximum in April 2014. Reversed polarity polar active sunspot regions observed from December 2016 indicate that the transitional phase to solar cycle 25 is in process. The next solar minimum will be reached probably in late 2019 or 2020. Total corrected (by water vapor and atmospheric pressure) neutron flux from January 206 to January 2019, issued from measurements and obtained by calculations using ATMORAD calculations and considering solar potential extract from Oulu data.

The daily neutron lethargic spectra (i.e. E×ϕ(E) in unit cm -2 .s -1 ) for the first of January of the year 2016, 17, 18 and 19 are presented in Fig. 4. The relative proportions of the each energy domain (thermal, epithermal, evaporation and cascade) are relatively invariant (Table I). This remark concerns four days for which solar conditions are quiet. Total corrected (by water vapor and atmospheric pressure) neutron flux from January 2016 to January 2019, issued from measurements and obtained by calculations using ATMORAD calculations and considering solar potential extract from Oulu data.

Short-term analysis: example of September 2017 solar event

During the first days of September 2017, the Sun was very active inducing solar flares. After a series of major flares including 27 M-class and 4 X-class flares, the solar active region had produced an X9.3 flare on the 6th September. The solar activity decreased on the 9th and 10th September. Bright gamma-ray emission from the X8.2 class solar flare was recorded on the 10th September 2017 by the Large Area Telescope and Fermi Gamma-ray Space Telescope [START_REF] Atwood | The Large Area Telescope on the Fermi Gamma-Ray Space Telescope Mission[END_REF]. This flare was characterized by a SEP event with high-energy particles, producing nuclear

PoS(ICRC2019)293

cascades in the atmosphere and detected on the ground. This flare was the second GLE (#72) of the current solar cycle. Fig. 5 presents the temporal series relative to (a) the water vapor absolute humidity (ground, + 100, + 250 and +500 m), (b) the IWV parameter, (c) the uncorrected and corrected (pressure and pressure + water vapor) hourly count rate obtained for the 7"Tungsten detector and (d) the total corrected neutron flux during the observed solar event (6 -13 September 2017).

As stated above, the data show that the Concordia environment can be therefore considered extremely dry during the solar event. Two FDs can be identified in Fig. 2(d). The second FD is composed by diurnal variations due to the anisotropy and superposed on the recovery step. It is a typical multi-step FD [START_REF] Bhaskar | Role of solar wind speed and interplanetary magnetic field during two-step Forbush decreases caused by Interplanetary Coronal Mass Ejections[END_REF], the amplitudes of steps are in the order of 3% and 11.3 %, respectively. Complementary analyses [START_REF] Hubert | Analysis of the Forbush Decreases and Ground-Level Enhancement on September 2017 using Neutron spectrometers operated in Antarctic and mid-latitude stations[END_REF] investigate the flux variations for each neutron energy domain and during the solar event including FD and GLE. During FD, fluxes relative to the four energy domains decrease according to a similar factor. However, during GLE, the thermal, epithermal and evaporation neutron populations increase more than cascade neutrons. The thermal peak presents a change in amplitude resulting from a change in the evaporation neutron population. Evaporation neutrons are moderated in the environment producing thermal neutrons and the increasing thermal neutrons are directly impacted by the increasing of evaporation neutrons component.

Conclusion

First four year of measurements of the energy spectrum of CR on the Concordia Station in Antarctica are presented. The particularity of these analyses is that the data from the HAMSTRAD project were used to take into account the environmental and systematic effects. Thus, temporal series of integrated water vapor, liquid water path, atmospheric pressure, neutron fluxes and spectra were presented. Long-and short-term analyses were presented, in particular during the solar event of September 2017. Simulations allow to show that variations observed in the long-term are mainly attributable to the solar activity. The relative proportions of each energy domain (thermal, epithermal, evaporation and cascade) are relatively invariant except during GLE event for which thermal, epithermal and evaporation domains are the most impacted.

Fig. 1 .

 1 Fig. 1. (a) Concordia station, (b) view of the neutron spectrometer in the Physic shelter and (c) view of the HAMSTRAD radiometer outside the station, it is possible to distinguish a phenomenon of diamond dust producing a halo around the sun from ice crystals in suspension.

Fig. 2 .

 2 Fig. 2.Temporal evolution of the LWP (kg.m -2 ) and IWV (g.m -2 ) as measured by HAMSTRAD above the altitude of Concordia station from January 2016 to January 2019. Also shown are the atmospheric pressure measured in Concordia during the same period thanks to the pressure detector integrated in the neuron spectrometer.

Fig. 3 .

 3 Fig. 3.Total corrected (by water vapor and atmospheric pressure) neutron flux from January 206 to January 2019, issued from measurements and obtained by calculations using ATMORAD calculations and considering solar potential extract from Oulu data.

Fig. 4 .

 4 Fig. 4.Total corrected (by water vapor and atmospheric pressure) neutron flux from January 2016 to January 2019, issued from measurements and obtained by calculations using ATMORAD calculations and considering solar potential extract from Oulu data.

Fig. 5 .

 5 Fig. 5.(a) water vapor absolute humidity (ground, + 100, + 250 and +500 m), (b) the IWV parameter (water vapor integrated along the vertical from the surface to the tropopause), (c) the uncorrected/corrected hourly count rate (7"Tungsten sphere) recorded in Concordia during solar events and (d) the total neutron flux.

  Water vapor abolute humidityWater vapor (g/m3

Table I :

 I Total neutron flux and contributions in % of thermal (Th.), epithermal (Epi.), evaporation (Eva.) and cascade (Cas.) neutron, for the OI January of 2016, 17, 18 and 19.

Total flux Contribution (in %) /cm²/s Th. Epi. Eva. Cas. 01/01/2016

	0.165	17.6	19.5 34.3 28.6
	01/01/2017 0.175	18.3	19.1 34.7 28.0
	01/01/2018 0.172	17.4	19.4 33.9 29.3
	01/01/2019 0.178	17.6	19.2 33.9 29.4