Hélène Coullon

Julien Bigot
email: julien.bigot@cea.fr

Christian Perez

The Multi-Stencil Language: orchestrating stencils with a mesh-agnostic DSL

Keywords: Domain Specic Language (DSL), stencil, numerical simulation, data parallelism, task parallelism, scheduling, MPI, OpenMP Langage de domaine sp¢ique, stencil, simulation num¯ique, parall¨isme de donnés, parall¨isme de tâches, ordonnancement, MPI, OpenMP

As the computation power of modern high performance architectures increases, their heterogeneity and complexity also become more important. One of the big challenges of exascale is to get programming models which gives access to high performance computing (HPC) to many scientists and not only to a few HPC specialists. One relevant solution to ease parallel programming for scientists is Domain Specic Language (DSL). However, one problem to avoid with DSLs is to not design a new DSL each time a new domain or a new problem has to be solved. This phenomenon happens for stencil-based numerical simulations, for which a large number of languages has been proposed without code reuse between them. The Multi-Stencil Language (MSL) presented in this paper is a language common to any kind of mesh used into a stencil-based numerical simulation. It is said that MSL is mesh-agnostic. Actually, from the description of a numerical simulation, MSL produces an empty parallel pattern, or skeleton, of the simulation which will be lled using other existing parallel languages and libraries. Thus, MSL, by nding a common language for dierent kinds of stencil-based simulation, facilitates code reuse. MSL is evaluated on a real case simulation which solves shallow-water equations. It is shown that MSL does not introduce overheads on data parallelism up to 16.384 cores, and that the hybrid parallelism (data and task) introduced improves performance of the simulation.

Introduction

As the computation power of modern high performance architectures increases, their heterogeneity and complexity also become more important. For example, the current world's top supercomputer Tianhe-2 1 is composed of multi-cores processors and accelerators, and is able to reach a theoretical peak performance around thirty peta oating point operations par seconds. However, to be able to use such a machine, multiple programming models, such as MPI (Message Passing Interface), OpenMP, CUDA etc., and multiple optimization techniques, such as cache or vectorization optimizations, have to be combined. Moreover, current architectures evolution lets think that heterogeneity and complexity in HPC will continue to grow in future.

One of the big challenges to solve to be able to use exascale computers is to propose programming models which gives access to high performance computing (HPC) to many scientists and not only to a few HPC specialists [START_REF]ETP4HPC Strategic Research Agenda Achieving HPC leadership in Europe[END_REF]. Actually, applications which run on supercomputers and which need such a computation power are physics, weather or genomic applications, which are not implemented by HPC specialists most of the time.

One possible runtime execution model for HPC is to propose dynamic scheduling of task graphs combined to message passing models [START_REF] Augonnet | StarPU: a unied platform for task scheduling on heterogeneous multicore architectures[END_REF][START_REF] Gautier | Xkaapi: A runtime system for dataow task programming on heterogeneous architectures[END_REF][START_REF] Wu | Hierarchical DAG Scheduling for Hybrid Distributed Systems[END_REF] (to be able to use more than one machine). Those models increase HPC code portability and reach an interesting performance onto heterogeneous architectures, which is interesting to reach exascale programming. At a higher abstraction level, general purpose parallel languages, such as OpenMP [START_REF] Dagum | Openmp: an industry standard api for shared-memory programming[END_REF] and OpenCL [START_REF] Stone | Opencl: A parallel programming standard for heterogeneous computing systems[END_REF] follow the direction of task graph scheduling proposed by those execution models. However, for non-experts end-users, general purpose languages still are dicult to use and to tune for a given application onto a given architecture. The current easiest (with the higher abstraction level), but still ecient, programming model for end-users is Domain Specic Language (DSL). Such a language is specic to the end-user domain and proposes a grammar which is easy to understand. The DSL compiler is able, because of the specic knowledge onto the domain targetted, to automatically apply parallelization and specic optimizations to produce a high performance back-end code. Thus, a DSL is able to split end-user concerns from HPC concerns which is relevant to reach exascale programming models.

However, many DSLs have been proposed for many domains yet, and, as far as we know, only a few of them have been able to reuse work already done by another language [START_REF] Sujeeth | Composition and reuse with compiled domain-specic languages[END_REF]. In other words, software engineering properties have to be integrated into DSL conception, such that a new DSL can be seen as a composition of parallelization, optimizations or even languages semantics already proposed by others.

For example, to numerically solve a set of partial dierential equations (PDEs), iterative methods are frequently used to approximate the exact solution through a discretized phenomena. A very well known and usual way to discretize PDEs is to transform them to explicit numerical schemes, also often called stencils. Many DSLs have been proposed for stencil computations [START_REF] Camier | Improving performance portability and exascale software productivity with the ∇ numerical programming language[END_REF][START_REF] Christen | PATUS: A Code Generation and Autotuning Framework for Parallel Iterative Stencil Computations on Modern Microarchitectures[END_REF][START_REF] Devito | Liszt: A domain specic language for building portable meshbased pde solvers[END_REF][START_REF] Ragan-Kelley | Halide: A language and compiler for optimizing parallelism, locality, and recomputation in image processing pipelines[END_REF][START_REF] Tang | The pochoir stencil compiler[END_REF], as it will be detailed in Section 2. Many of them use same kind of parallelization, data structures or optimizations, however each one has been built from scratch to deal with another additionnal specic case. We present the Multi-Stencil Language (MSL) DSL, also for stencil-based numerical simulations. MSL is a language with a light grammar to describe a numerical simulation without implementation details. From the description, the compiler has enough information to extract and build an empty parallel pattern of the simulation, which can be lled, in a second step, by implementation concerns. The parallel pattern generated by the language is able to use dierent existing languages and libraries as it is independent from implementation choices. Moreover, the parallelization performed by the language is large enough to be compatible with many architectures and back-end languages. Contributions presented in this paper are : the computational model of a multi-stencil program and its parallelization formalism; the MSL grammar and its compiler; a back-end implementation and its performance evaluation onto a real case numerical simulation up to 16.384 cores.

Section 2 introduces the related work on DSLs for stencils. Section 3 formally explains the targetted domain and its computational model. From this model can be extracted the grammar of MSL in Section 4. Section 5 shows how parallelism can be extracted from this light grammar of MSL. Sections 6 and 7 detail choices that have been done in this paper to evaluate MSL, and Section 8 shows performance results of the language. Finally, Section 9 concludes and proposes perspectives on this work. [START_REF] Bigot | A low level component model easing performance portability of hpc applications[END_REF] Related work

Many domain specic languages have been proposed for stencil computations. Each one has its own specicities and answers to a specic stencil case or to a specic additionnal optimization.

For example, Pochoir [START_REF] Tang | The pochoir stencil compiler[END_REF] works on cache optimization techniques for stencils applied onto Cartesian meshes. On the other hand, PATUS [START_REF] Christen | PATUS: A Code Generation and Autotuning Framework for Parallel Iterative Stencil Computations on Modern Microarchitectures[END_REF] proposes to add a parallelization strategy grammar to its stencil language to perform an auto-tuning parallelization strategy. Moreover, Halide [START_REF] Ragan-Kelley | Halide: A language and compiler for optimizing parallelism, locality, and recomputation in image processing pipelines[END_REF] proposes an optimization and parallelization of a pipeline of stencil codes, while ExaSlang [START_REF] Schmitt | Exaslang: A domain-specic language for highly scalable multigrid solvers[END_REF] is specic to multigrid numerical methods etc. The reason why each of those languages are implemented from scratch to answer its own particularities and its own optimizations is because each language is built and thought as a single block, which makes impossible or dicult code reuse from one language to another. In other words domain specic languages (for stencils or other domains) suer from a lack of software engineering properties, which could increase the productivity to build a new language, by codereuse, and also the maintainability of languages, with more separation of concerns. As far as we know, a single work proposes a framework to create DSL, which afterwards will be easier to compose and combine together [START_REF] Sujeeth | Composition and reuse with compiled domain-specic languages[END_REF], which is interesting for inter-domains applications. However, the point argued by this paper, and by the Multi-Stencil Language (MSL), is that DSL conception must be studied to maximize its usefulness for dierent types of applications, while keeping a good performance. This maximization is directly linked to the abstraction level proposed to the end-user, and also to separation of concerns and code reuse improvements.

MSL is itself a domain specic language for stencil-based simulations. It oers a way to give a mesh-agnostic description of a simulation, which could be common to dierent cases of simulations, and thus facilitates reuse of existing underlying data structures, optimizations or parallelization techniques. In other words, MSL extracts where parallelization and optimizations are needed into the simulation, by producing an empty parallel pattern, but not how to do it, which is left to existing languages (SkelGIS and OpenMP in this paper), and which is an implementation concern.

Liszt [START_REF] Devito | Liszt: A domain specic language for building portable meshbased pde solvers[END_REF] and Nabla [START_REF] Camier | Improving performance portability and exascale software productivity with the ∇ numerical programming language[END_REF], both oer languages for stencils applied onto any kind of mesh, from Cartesian to unstructured meshes. The mesh which is needed into the simulation can be built from a set of available symbols in the grammar of each language. Thus, those languages generalize the denition of a mesh, as it is proposed into the MSL formalism. However, two main dierences can be noticed.

First, the formalism is more exible in MSL. For example, a computation can be applied onto a subset of the space domain which is not possible with Liszt and Nabla. This functionnality is important in numerical simulations yet. For example, many computations performed in the simulation studied in this paper have to be applied onto a subpart of the overall space domain.

Second, in Liszt and Nabla the description of a simulation is not splitted from its implemen-Inria tation concerns. In other words, the topology of a mesh as well as numerical codes are given at the same time than the description of the dierent computations to apply into the simulation, while MSL splits those two dierent phases. Thus, MSL improves separation of concerns. Actually, a numerician could describes the dierent computations to perform into the simulation, while another one, later, could focuss on the implementation and the numerical code. This also facilitates reuse of other languages. Finally, the parallelization techniques proposed in this paper, take place at paradigm level, with data parallelism and hybrid (data plus task) parallelism. Thus the parallel pattern of the simulation does not need details onto parallel architectures (distributed or shared memories, with or without accelerators). This makes possible a large panel of backend architectures and languages. In other words MSL improves portability. This section gives a complete formal description of what we call a stencil program and its computations. The proposed language MSL (Multi-Stencil Language) does not depend on the type of space discretization used. Thus, some details given in the following sections, are useful to understand the domain but do not have to be kept in mind to understand the Multi-Stencil Language.

Time, mesh and data

Ω is the continuous space domain of a numerical simulation (typically R n ). A mesh M denes the discretization of the continuous space domain Ω of a set of PDEs and is dened as follows.

Denition 1 A mesh is a connected undirected graph M = (V, E), where V ⊂ Ω is the set of vertices and E ⊆ V 2 the set of edges. The set of edges E of a mesh M = (V, E) does not contain bridges. It is said that the mesh is applied onto Ω.

M = (V, E) applied onto Ω = R n is denoted dim(M) = n.
A mesh can be structured (as Cartesian or curvilinear meshes), unstructured, regular or irregular (without the same topology for each element) as illustrated in Figure 1. One can notice that more than one type of mesh is also possible inside a single simulation. For example, an hybrid mesh can be dened as an unstructured mesh composed itself of a Cartesian mesh inside each of its vertices. However in this paper, single mesh simulations are addressed.

Denitions (mesh)

An entity φ of a mesh M = (V, E) is dened as a subset of its vertices and edges, φ ⊂ V ∪E.

A group of mesh entities G ∈ P(V ∪ E) represents a set of entities of the same topology and the mesh of a group is denoted mesh(G), i.e.

G ∈ P(V ∪ E) ⇔ mesh(G) = (V, E)
The set of entities groups used in a simulation is denoted Φ.

For example, in a 2D Cartesian mesh, a group of entities where each entity is made of four vertices and four edges could form the cells. Another group of entities are the vertices dened as all singletons formed of a single vertex of V . Both groups, then, would be part of Φ. This example is illustrated in Figure 2b. Denition 3 (Denitions (time)) The nite sequence T : (t n ) n∈ 0,Tmax represents the discretization of the continuous time domain T = R. To each discrete time-step n ∈ 0, T max , it associates a time value t n ∈ T .

The time discretization can be as simple as a constant time-step with a xed number of steps. It can, however, also be dened by recursion with both the time-step and the number of steps depending on the data variables (see denitions below).

Denitions (quantity)

V = ∆ ∪ S is the set of variables or quantities.

∆ are the mesh variables that to each entity of a group and time-step associates a value δ : G × T → V δ where V δ is a value type.

The group of entities a variable is mapped on is denoted entity(δ) = G.

S are the scalar variables that to each time-step associates a value s : T → V δ where V δ is a value type. Amongst the scalar variables is one specic variable conv ∈ S, the convergence criteria whose value is 0 except at the last step where it is 1. Thus, ∀t ∈ 0, T max-1 , conv(t) = 1, and conv(T max ) = 1.

This section has presented the general formalism of meshes, their entities, groups of entities, time discretization, and nally quantities (mapped on mesh and scalars). While the details about mesh and their topologies were useful to specify the following concepts, the Multi-Stencil Language presented in this paper is mesh-agnostic. These details will therefore not be needed in the remaining of the paper.

Computations

Denitions

A computation domain D is a subpart of a mesh entities group, D ⊆ G ∈ Φ.

Inria

The set of computation domains of a numerical simulation is denoted D. N are the neighborhood functions n : G i → G m j which for a given entity φ ∈ G i returns a set of m entities in G j . One can notice that i = j is possible. Most of the time, such a neighborhood is called a stencil shape. Denition 4 A computation kernel k of a numerical simulation is dened as k = (S, R, (w, D), comp), where S ∈ S is the set of scalar to read, (w, D) ∈ ∆ × D is the single data written by the kernel:

w is the single quantity (variable) modied by the computation kernel, D is the computation domain on which w is computed, D ⊆ entity(w), R ∈ ∆ × N is the set of tuples (r, n) representing the data read where r is a quantity read by the kernel to compute w, n : entity(w) → entity(r) m is a neighborhood function that indicates whose values of r are read to compute w. comp is the numerical computation which returns a value from a set of n input values, comp : V n i → V j , where V i and V j are value types. Thus, comp represents the actual numerical expression which is computed by a kernel.

At each time-step, a set of computations is performed. During a computation kernel, it can be considered that a set of old states (t -1) of quantities are read (R), and that a new state (t) of a single quantity is written (w).

Such a denition of a computation kernel covers a large panel of dierent computations. For example, the four usual types of computations (stencil, local, boundary and reduction) performed into a simulation can be dened as follow :

A kernel computation k(S, R, (w, D), comp) is a stencil kernel if ∃(r, n) ∈ R such that n = identity .
A boundary kernel is a kernel k(S, R, (w, D), comp) where D is a specic computation domain at the border of entities, and which does not intersect with any other computation domain.

A kernel computation k(S, R, (w, D), comp) is a local kernel if ∀(r, n) ∈ R, n = identity.

A kernel computation k(S, R, (w, D), comp) is a reduction kernel if w is a scalar.

Since we only consider explicit numerical schemes in this paper, a kernel can not write the same quantity it reads, except locally, i.e. if ∃(w, n) ∈ R ⇒ n = identity A reduction is typically used to compute the convergence criteria of the time loop of the simulation.

Denition 5 The set of n ordered computation kernels of a numerical simulation is denoted

Γ = [k i ] 0≤i≤n-1 , such that ∀k i , k j ∈ Γ, if i < j, then k i is computed before k j . Denition 6 A multi-stencil program is dened by the octuplet MSP(M, Φ, D, N , ∆, S, T, Γ) Mesh Cells Edgex (a)
Mesh and two groups of mesh entities.

x,y x,y

x y+1 x y-1 x-1 y x+1 y A B (b) A is computed with a 4- neighborhood stencil applied on B.
A is computed onto a computation domain which does not include all entities of the group. Example For example, in Figure 2b, assuming that the computation domain (full lines) is denoted dc1 and the stencil shape is n1, the stencil kernel can be dened as:

R : {(B, n1)}, w : A, D : dc1, comp : A(x, y) = B(x + 1, y) + B(x -1, y) + B(x, y + 1) + B(x, y -1).
On the other hand, in the example of Figure 2c, assuming the computation domain is dc2 and the stencil shape is n2, the stencil kernel is dened as:

R : {(C, n2), (A, identity)}, w : A, D : dc2, comp : A(x, y) = A(x, y) + C(x1, y1) + C(x1 + 1, y1).
A stencil program has been formally dened in this section. This formalism is used in the next Section to dene two parallelization techniques of a multi-stencil program.

One can note that all denitions given in this section are independent from the topology of the mesh. This property will be kept in the rest of this paper to propose the mesh-agnostic MSL language.

1 program : : = "mesh : " meshid 2 "mesh e n t i t i e s : " l i s t g r o u p 3 " computation domains : " listcompdom 4 " i n d e p e n d e n t : " l i s t i n d e 5

" s t e n c i l s h a p e s : " l i s t s t e n c i l 6

"mesh q u a n t i t i e s : " l i s t q u a n t i t i e s 7

" s c a l a r s " l i s t s c a l a r 8 l i s t l o o p 9 l i s t g r o u p : : = g r o u p i d " , " l i s t g r o u p | g r o u p i d listcompdom : : = compdom listcompdom | compdom compdom : : = compdomid " i n " g r o u p i d l i s t i n d e : : = i n d e l i s t i n d e | i n d e i n d e : : = compdomid "and" compdomid l i s t s t e n c i l : : = s t e n c i l l i s t s t e n c i l | s t e n c i l s t e n c i l : : = s t e n c i l i d " from " g r o u p i d " t o " g r o u p i d l i s t q u a n t i t i e s : : = q u a n t i t y l i s t q u a n t i t i e s | q u a n t i t y q u a n t i t y : : = g r o u p i d l i s t q u a n t i t y i d l i s t q u a n t i t y i d : : = q u a n t i t y i d " , " l i s t q u a n t i t y i d | q u a n t i t y i d l i s t s c a l a r : : = s c a l a r i d " , " l i s t s c a l a r | s c a l a r i d l i s t l o o p : : = l o o p l i s t l o o p | l o o p l o o p : : = " time : " i t e r a t i o n " c o m p u t a t i o n s : " l i s t c o m p i t e r a t i o n : : = num | s c a l a r l i s t c o m p : : = comp l i s t c o m p | comp comp : : = w r i t t e n "=" compid " ( " l i s t r e a d " ) " w r i t t e n : : = q u a n t i t y i d " [ " compdomid " ] " | s c a l a r l i s t r e a d : : = d a t a r e a d l i s t r e a d | d a t a r e a d d a t a r e a d : : = q u a n t i t y i d " [ " s t e n c i l i d " ] " | q u a n t i t y i d | s c a l a r The grammar of the Multi-Stencil Language is given in Figure 3 and an example is provided in Figure 4. A Multi-Stencil program is composed of eight parts.

1. The mesh keyword (Fig. 3, l.1) introduces an identier for M, the single mesh of the simulation. For example cart in Fig. 4, l.1. Since the language is independent of the mesh topology, this is not used by the compiler but is syntactic sugar for the user.

2. The mesh entities keyword (Fig. 3, l.2) introduces identiers for the groups of entities G ∈ Φ. For example cell or edgex in Fig. 4, l.2. Again, this is syntactic sugar only as the compiler does not rely on the mesh topology.

3. The computation domains keyword (Fig. 3, l.3) introduces identiers for the computation domains D ∈ D. For example d1 and d2 in Fig. 4, l.4-5. For reference, each domain is associated to a group of entities (Fig. 3, l.12) such as cell for d1 in Fig. 4, l.4.

4. The independent keyword (Fig. 3, l.4) oers a way to declare that computation domains do not overlap, such as d1 and d2 in Fig. 4, l.7. This is used by the compiler to compute dependencies between computations.

[ d1 ] = k0 ( tau ,A) C [ d2 ] = k1 (B [ nec ] ) D[ d1 ] = k2 (C) E [ d1 ] = k3 (C) F [ d1 ] = k4 (D, C [ nce ] ) G[ d1 ] = k0 (mu, tau , E) H[ d2 ] = k6 (F) I [ d1 ] = k7 (G,H) J [ d1 ] = k8 (mu, I [ ncc ] )
Figure 4: Example of program using the Multi-Stencil Language.

5. The stencil shapes keyword (Fig. 3, l.5) introduces identiers for n ∈ N , the stencil shapes. For each shape, the source and destination group of entities (Fig. 3, l.16) are specied. For example nec in Fig. 4, l.11 associates cell entities to each edgex entity.

6. The mesh quantities keyword (Fig. 3, l.6) introduces identiers for δ ∈ ∆, the quantities with the group of entities they are mapped on (Fig. 3, l.16). For example the quantities C and H are mapped on edgex entities.

7. The scalars keyword (Fig. 3, l.7) introduces identiers for s ∈ S, the scalars. For example mu and tau in Fig. 4, l.15.

8. Finally, the last part (Fig. 3, l.8) introduces the computations loops. Each loop is made of two parts:

the time keyword (Fig. 3, l.22) introduces conv, the convergence criteria either as a number of iterations or as 0, 1 valued scalar (Fig. 3, l.24). For example, 500 iterations are specied in Fig. 4, l.16, the computations keyword (Fig. 3, l.23) introduces identiers for each computation k ∈ Γ = (S, R, (w, D), comp). Each computation (Fig. 3, l.26) species:

the quantity w written and its domain D, for example in Fig. 4, l.22, kernel k4

computes the quantity F on domain d1, Inria the read scalars S and mesh quantities with their associated stencil shape R, for example in Fig. 4, l.16, k4 reads C with the shape nce and D with the default identity shape, it uses no scalar.

One can notice that in the example of Figure 4, there are no (reduction) kernel associated to the scalars mu and tau. In this case, those scalars are in fact constants. One can also notice that the computation to execute for each kernel is not specied. This is indeed not handled by MSL, which is mesh-agnostic, and which only generates a skeleton of the application (a parallel pattern of the application). Actually, a numerical computation is naturally mesh-dependent and is left to other languages. The numerical computation is specied elsewhere by referencing the identier chosen here.

This section has introduced the grammar of MSL. The next section will show how parallelization can be extracted from this simple language and how an empty parallel skeleton (or pattern) of the application can be generated.

Parallelization of Multi-Stencil Programs

As previously explained, in a computation k(S, R, (w, D), comp), comp is not handled by MSL. As a result, in the rest of this paper, and to simplify notations, we denote the same computation k(S, R, (w, D)).

Data parallelism

In a data parallelization technique, the idea is to split quantities on which the program is computed into balanced sub-parts, one for each available resource. The same sequential program can afterwards be applied on each sub-part simultaneously, with some additioinal synchronizations between resources to update the data not computed locally, and thus to guarantee a correct result.

More formally, the data parallelization of a multi-stencil program

MSP(M, Φ, D, N , ∆, S, T, Γ)
consists in, rst, a partitioning of the mesh M in p balanced sub-meshes (for p resources) M = {M 0 , . . . , M p-1 }. This step can be performed by an external graph partitionner [START_REF] Coullon | Algorithmic skeleton library for scientic simulations: Skelgis[END_REF][START_REF] Lachat | PaMPA: Parallel Mesh Partitioning and Adaptation[END_REF][START_REF] Pellegrini | Scotch: A software package for static mapping by dual recursive bipartitioning of process and architecture graphs[END_REF] and is not adressed by this paper.

As entities and quantities are mapped onto the mesh, the set of groups of mesh entities and the set of quantities ∆ are partitionned the same way than the mesh:

Φ = {Φ 0 , . . . , Φ p-1 }, ∆ = {∆ 0 , . . . , ∆ p-1 }.
The second step of the parallelization is to identify in Γ the needed synchronizations between resources to update data, and thus to build a new ordered list of computations Γ sync . Denition 7 For n the number of computations in Γ, and for i, j such that i < j < n, a synchronization is needed between k i and k j , denoted k i ≺≺ ≺ k j , if ∃(r j , n j ) ∈ R j such that w i = r j and n j = identity (k j is a stencil computation). The quantity to synchronize is {w i }.

Actually, a synchronization is needed by the quantity read by a stencil computation (not local), if this quantity has been modied before, which means that it has been written before. This synchronization is needed because a neighborhood function n ∈ N of a stencil computation involves values computed on dierent resources. However, as a multi-stencil program is an iterative program, computations which happen after k j at the time iteration t have also been computed before k j at the previous time iteration t -1. For this reason another case of synchronization has to be dened. Denition 8 For n the number of computations in Γ and j < n, if ∃(r j , n j ) ∈ R j such that n j = identity and for all i < j, k i ≺≺ ≺ k j , a synchronization is needed between k l and k j , where j < l < n, denoted k l ≺≺ ≺ k j , if w l = r j . The quantity to synchronize is {w l }.

Denition 9 A synchronization between two computations k i ≺≺ ≺ k j is dened as a specic com- putation

k sync i,j (S, R, (w, D)),
where

S = ∅, R = {(r, n)} = {(w i , n j ∈ N }, (w, D) = (w i , φ∈Dj n j (φ))).
In other words, w i has to be synchronized for the neighborhood n j for all entities of D j .

Denition 10 If k i ≺≺ ≺ k j , k j is replaced by the list

[k sync i,j , k j ]
When data parallelism is applied, the other type of computation which is responsible for additional synchronizations is the reduction. Actually, the reduction is rst applied locally on each subset of entities, on each resource. Thus, p (number of resources) scalar values are obtained. For this reason, to perform the nal reduction, a set of synchronizations are needed to get the nal reduced scalar. As most parallelism libraries (MPI, OpenMP) already propose a reduction synchronization with its own optimizations, we simply choose to replace the reduction computation by itself anotated by red.

Denition 11 A reduction kernel k j (S j , R j , (w j , D j )), where w is a scalar, is replaced by k red j (S j , R j , (w j , D j )).

One can notice that both types of synchronizations are performed by all resources.

Denition 12 The concatenation of two ordered lists of respectively n and m computations

l 1 = [k i ] 0≤i≤n-1 and l 2 = [k i ] 0≤i≤m-1 is denoted l 1 • l 2
and is equal to a new ordered list

l 3 = [k 0 , . . . , k n-1 , k 0 , . . . , k m-1 ].
Denition 13 From the ordered list of computation Γ, a new synchronized ordered list Γ sync is obtained from the call Γ sync = F sync (Γ, 0), where F sync is the recursive function dened in Algorithm 1.

Algorithm 1 follows previous denitions to build a new ordered list which includes synchronizations. In this algorithm, lines 7 to 19 apply Denition [START_REF] Coullon | Implementation and performance analysis of skelgis for network mesh-based simulations[END_REF], lines 20 to 29 apply Denition [START_REF] Dagum | Openmp: an industry standard api for shared-memory programming[END_REF], and nally lines 34 and 35 apply Denition [START_REF] Ferrari | A new two-dimensional shallow water model including pressure eects and slow varying bottom topography[END_REF]. Finally, line 39 of the algorithm is the recursive call.

The nal step of this parallelization is to run Γ sync on each resource. Thus, for each resource

0 ≤ r ≤ p -1 the multi-stencil program MSP r (M r , Φ r , D r , N , ∆ r , S, T, Γ sync ), (1)
is performed.

Inria Algorithm 1 F sync recursive function 1: procedure F sync (Γ,j)

2:

k j = Γ[j] 3: list = [] 4: if j = |Γ| then 5:
return list

6:
else if ∃(r j , n j ) ∈ R j such that n j = identity then 7:

for all (r j , n j ) ∈ R j such that n j = identity do 8:

found = false 9:

for 0 ≤ i < j do 10:

k i = Γ[i] 11:
if k i ≺≺ ≺ k j then (w, D) = (w i , φ∈Dj n j (φ))) (w, D) = (w i , φ∈Dj n j (φ))) Example Figure 4 gives an example of a MSP program. From this example, the following ordered list of computation kernels can be extracted:

k i = Γ[i]
Γ = [k 0 , k 1 , k 2 , k 3 , k 4 , k 0 , k 6 , k 7 , k 8 ]
From this ordered list of computation kernels Γ, and from the rest of the multi-stencil program, synchronizations can be automatically detected from the call to F sync (Γ, 0) to get the synchronized ordered list of kernels:

Γ sync = [k 0 , k sync 0;1 , k 1 , k 2 , k 3 , k sync 1;4 , k 4 , k 0 , k 6 , k 7 , k sync 7;8 , k 8 ], (2)
where

k sync 0;1 = (∅, {(B, nce)}, (B, ∪ φ∈D1 nce(φ))), (3a)
k sync 1;4 = (∅, {(C, nec)}, (C, ∪ φ∈D4 nec(φ))), (3b)
k sync 7;8 = (∅, {(I, ncc)}, (I, ∪ φ∈D8 ncc(φ))). (3c)

Hybrid parallelism

A task parallelization technique is a technique to transform a program as a dependency graph of dierent tasks. A dependency graph exhibits parallel tasks, or on the contrary sequential execution of tasks. Such a dependency graph can directly be given to a dynamic scheduler, or can statically be scheduled. In this paper, we introduce task parallelism by building the dependency graph between kernels of the sequential list Γ sync . Thus, as Γ sync takes into account data parallelism, we introduce hybrid parallelism.

Denition 14 For two computations k i and k j , with i < j, it is said that k j is dependant from k i with a read after write dependency, denoted k i ≺ raw k j , if ∃(r j , n j ) ∈ R j such that w i = r j . In this case, k i has to be computed before k j .

Denition 15 For two computations k i and k j , with i < j, it is said that k j is dependant from k i with a write after write dependency, denoted k i ≺ waw k j , if w i = w j and D i ∩ D j = ∅. In this case, k i also has to be computed before k j .

Denition 16 For two computations k i and k j , with i < j, it is said that k j is dependant from k i with a write after read dependency, denoted

k i ≺ war k j , if ∃(r i , n i ) ∈ R i such that w j = r i .
In this case, k i also has to be computed before k j is started so that values read by k i are relevant.

Those denitions are known as data hazards classication. However, a specic condition on the computation domain, due to the specic domain of multi-stencils, is introduced for the write after write case.

Denition 17 A directed acyclic graph (DAG) G(V, A) is a graph where the edges are directed from a source to a destination vertex, and where, by following edges direction, no cycle can be found from a vertex u to itself. A directed edge is called an arc, and for two vertices v, u ∈ V an arc from u to v is denoted ( u, v) ∈ A.

From an ordered list of computations Γ sync , a directed dependency graph Γ dep (V, A) can be built nding all pairs of computations k i and k j , with i < j, such that k i ≺ raw k j or k i ≺ waw k j or k i ≺ war k j .

Denition 18 For two directed graphs G(V, A) and G (V , A ), the union (V, A) ∪ (V , A ) is dened as the union of each set (V ∪ V , A ∪ A ).

Denition 19 From the synchronized ordered list of computation kernels Γ sync , the dependency graph of the computations Γ dep (V, A) is obtained from the call F dep (Γ sync , 0), where F dep is the recursive function dened in Algorithm 2.

Algorithm 2 F dep recursive function

1: procedure F dep (Γ sync ,j) 2: k j = Γ sync [j] 3: if j = |Γ sync | then 4:
return ({}, {})

5: else if j < |Γ sync | then 6: G = ({}, {}) 7:
for 0 ≤ i < j do 8: end if 15: end procedure This constructive function is possible because the input is an ordered list. Actually, if k i ≺ k j then i < j. As a result, k i is already in V when the arc ( k i , k j ) is built.

k i = Γ sync [i] 9: if k i ≺ raw k j or k i ≺ waw k j or k i ≺ war k j then 10: G = G ∪ (k j , {( k i , k j })
One can notice that Γ dep is the dependency graph of the computations of a multi-stencil program, but it only takes into account a single time iteration. A complete dependency graph of the simulation could be built. This is a possible extension of this work. Proposition 20 The directed graph Γ dep is an acyclic graph.

As a result of the hybrid parallelization, each resource 0 ≤ r ≤ p -1 perform a multi-stencil program, dened by

MSP r (M r , Φ r , D r , N , ∆ r , T, Γ dep ).
The set of computations Γ dep is a dependency graph between computation kernels k i of Γ and synchronizations of kernels added into Γ sync . Γ dep can be built from the call to

F dep (F sync (Γ, 0), 0).
Example Figure 4 gives an example of MSP program. From Γ sync that has been built in Equation ( 2), the dependency DAG can be built. For example, as k 0 computes B and k 1 reads B, k 0 and k 1 becomes vertices of Γ dep , and an arc ( k 0 , k 1 ) is added to Γ dep . The overall Γ dep built from the call to F dep (Γ sync , 0) is drawn in Figure 5. In this section we detail a static scheduling of Γ dep by using minimal series-parallel directed acyclic graphs. Such a static scheduling may not be the most ecient one, but it oers a simple fork/join task model which make possible the design of a performance model. Moreover, such a scheduling oers a simple way to propose a fusion optimization.

k 0 k sync 0;1 k 1 k 2 k sync 1;4 k 3 k 4 k 5 k 6 k 7 k sync 7;8 k 8

Series-Parallel graph

In 1982, Valdes & Al [START_REF] Valdes | The recognition of series parallel digraphs[END_REF] have dened the class of Minimal Series-Parallel DAGs (MSPD). Such a graph can be decomposed as a serie-parallel tree, denoted T SP , where each leaf is a vertex of the MSPD DAG it represents, and whose internal nodes are labelled S or P to indicate respectively the series or parallel composition of sub-trees. Such a tree can be considered as a fork-join model and as a static scheduling.

Valdes & Al [START_REF] Valdes | The recognition of series parallel digraphs[END_REF] have identied a forbidden shape, or subgraph, called N , such that the following property is veried : Theorem 21 The transitive reduction of a DAG G is MSPD if and only if it does not contain N as an induced subgraph.

To remove the forbidden N-shapes from the transitive reduction of Γ dep = (V, E), we have chosen to apply an over-constraint with the relation k 0 ≺ k 3 , such that a complete bipartite graph is created for the sub-dag, and can be translated to a series-parallel decomposition, as illustrated in Figure ??.

k 0 k 1 k 2 k 3
Figure 6: Over-constraint on the forbidden N shape.

After these over-constraints are applied, Γ dep is MSPD. Valdes & Al [START_REF] Valdes | The recognition of series parallel digraphs[END_REF] have proposed a linear algorithm to know if a DAG is MSPD and, if it is, to decompose it to its associated binary decomposition tree. As a result, the binary tree decomposition algorithm of Valdes & Al can be applied on Γ dep to get the T SP static scheduling of the multi-stencil program.

Example The Serie-Parallel tree decomposition of the example given in Figure 4, which is built from the dependency graph of Figure 5 is given in Figure 7. In this subsection are introduced two performance model, one for the data parallelization technique, and one for the hybrid data and task parallelization technique, both previously explained.

The performance model of a data parallelization technique is inspired from the Bulk Synchronous Parallel model. Actually data parallelization technique consider that each process has its own part of the domain. Thus the performance model reveals that the computation time is the sum of the reference sequential time divided by the number of processes, and of the time spent in communications between processes. Thus, for T SEQ the sequential reference time, P the total number of processes, T COM the communications time, the total computation time is

T = T SEQ P + T COM . (4)
Thus, when the number of processes P increase in data parallelization, the performance model limit is T COM lim

P →+∞ T = T COM . (5)
As a result, the critical point for performance is when T COM ≥ T SEQ P , which happens naturally in data parallelization as T COM will increase with the number of processes, and

T SEQ P
decrease with the number of processes.

This limitation is always true, but can be delayed by dierent strategies. First, it is possible to perform communications through the network while computations are performed simultaneously. Second, it is possible to introduce another kind of parallelization, task parallelization. Thus, for the same total number of processes, only a part of them are used for data parallelization, and the rest are used for task parallelism. As a result, T SEQ P will continue to decrease but T COM will increase later. This second strategy is the one studied in the following hybrid performance model.

For an hybrid (data and task) parallelization technique, and for P data the number of processes used for data parallelization, P task the number of processes used for task parallelization, such that P = P data × P task is the total number of processes used, T task the overhead time due to task parallelization technique, and F task the task parallelization degree of the application, the total computation time is

T = T SEQ P data × F task + T COM + T task (6)
The time overhead due to task parallelization can be represented as the time spent to create a pool of threads and the time spent to synchronize those threads. Thus, for T cr the total time to create the pool of threads (may happened more than once), T sync the total time spent to synchronize threads, the overhead is

T task = T cr + T sync .
The task parallelization degree of the application F task is the limitation of a task parallelization technique. As explained before, a task parallelization technique is based on the dependency graph of the application. Thus, this dependency graph must expose enough parallelism for the number of available threads. For this performance model we consider that

F task = P task .
However, the upper bound of F task is constrained by the dependency graph of the application and by the time spent in each task.

As a result when P data is small a data parallelization technique may be more ecient, while an hybrid parallelization could be interesting at some point to improve performance. The question asked here is when is it interesting to use hybrid parallelization.

To answer this question lets consider the two parallelization techniques, data only and hybrid. We denote P data1 the total number of processes entirely used by the data only parallelization, P data2 the number of processes used for data parallelization in the hybrid parallelization, and P task the number of processes used for task parallelization in the hybrid parallelization, such that P data1 = P data2 × P task .

We search the point where the data parallelization is less ecient than the hybrid parallelization. Thus,

T SEQ P data1 + T COM 1 ≥ T SEQ P data2 × P task + T COM 2 + T task .
This happens when

T COM 1 ≥ T COM 2 + T task (7)
This performance model will be validated and will help to explain results the Section 8. Inria

Fusion optimization

Using MSL, it is possible to ask for data parallelization of the application, or for an hybrid parallelization. Even though the MSL language is not dedicated to produce very optimized stencil codes, but to produces the parallel pattern of the application, building the T SP tree make available an easy optimization when the data parallelization technique is the only one used. This optimization consists in proposing a valid merge of some computation kernels inside a single space loop. As a result, the user can use this valid fusion of kernels or not when implementing those. Those fusions can be computed from the canonical form of the T SP tree decomposition. The canonical form consists in recursively merging successive S vertices or successive P vertices of T SP .

The fusion function F f us is described in Algorithm 3, where parent(k) returns the parent vertex of k in the tree, and where k f us i;j represents the fusion of k i and k j keeping the sequential order i; j if i is computed before j in T SP . Finally, type(k) returns comp if the kernel is a computation kernel, and sync otherwise.

Algorithm 3 F f us

1: procedure F f us (T SP (V, E)) 2: for (k i , k j ) ∈ V 2 do 3: if parent(k i )==parent(k j ) then 4: if type(k i ) = type(k j ) = comp then 5: if parent(k i )==S then 6: if D i == D j then 7:
propose the fusion k f us i;j end for 17: end procedure We are not arguing that such a simple fusion algorithm could be as good as complex cache optimization techniques which can be found in stencil DSLs for example [START_REF] Tang | The pochoir stencil compiler[END_REF]. However, this fusion takes place at a dierent level and can bring performance improvments as it will be illustrated in Section 8. This fusion algorithm relies on very simple statements: Two successive computation kernels k i and k j which are under the same parent vertex S in TSP are, by construction, data dependant. As a result, what is written by the rst one is read by the second one. Thus, at least one data is common to those computations (the one written by k i ). Thus, if the computation domains verify D i = D j , the fusion of k i and k j will decrease cache misses. Two successive computation kernels k i and k j which are under the same parent vertex P in TSP are not, by construction, data dependant. However, if the computation domains verify D i = D j , and if R i ∩ R j = ∅ cache misses could also be decreased by the fusion k f us i;j .

7 Implementation

The MSL compiler takes a le written in the grammar described in Section 4 as input and generates an application skeleton. In this skeleton, the user still has to ll the functions corresponding to the various computation kernels. The overall behavior of the compiler is as follow:

1. it parses the input le and generates Γ, the list of computation kernels, 2. from Γ, it builds Γ sync , the list including synchronizations for data parallelism using Algorithm F sync introduced in Section 5, 3. from Γ sync , it builds Γ dep , the DAG supporting hybrid parallelism using Algorithm F dep introduced in Section 5, 4. it then removes the N-Shapes from Γ dep to get a MSPD graph, and generates its serieparallel binary tree decomposition T SP , 5. it performs the fusion of kernels in T SP if required, 6. nally it dumps an application structure matching T SP .

The abstract denition of MSL enables to use a very wide range of languages for the dump. In this paper, we use a compiler that targets the Low Level Components [START_REF] Bigot | A low level component model easing performance portability of hpc applications[END_REF] (L 2 C), a component model for C++ supporting high-performance applications. Similar to classes in object-oriented models, components specify the services they provide, but in addition they also specify the services they require. This enables to build applications by assembling components in a second phase through the matching of requirements of some components with the services provided by others. In the case of MSL, this assembly of components oers a way to clearly separate in distinct components the features that are implemented: by reusing the same code more than once, such as for the time iteration logic, by generating some code that depends on T SP and thus from the input le, such as for the task parallelism logic, or by letting the user ll-in an empty skeleton component, such as for the computation kernels.

By enforcing the use of well dened interfaces in the generated code, L 2 C makes easier the substitution of an implementation by another, for example to reconsider an implementation choice. It also simplies the addition of new features to MSL or the coupling of MSL with another domain-specic language, for example.

The two main other choices in the code generation concern the technologies used for data and task parallelizations. For the data-parallelization, we rely on SkelGIS, a C++ embedded DSL [START_REF] Coullon | The SIPSim implicit parallelism model and the SkelGIS library[END_REF] which oers distributed data structures (meshes) for numerical simulations, and programming interfaces to easily write codes using them. SkelGIS is implemented over MPI [START_REF] Graham | The MPI 2.2 Standard and the Emerging MPI 3 Standard[END_REF] and can therefore be used on distributed memory architectures such as clusters. Two kind of meshes are oered by SkelGIS: a) a distributed two dimensional Cartesian mesh [START_REF] Coullon | Algorithmic skeleton library for scientic simulations: Skelgis[END_REF], and b) a distributed graph of Cartesian meshes (hybrid mesh) [START_REF] Coullon | Implementation and performance analysis of skelgis for network mesh-based simulations[END_REF]. The use of dierent back-ends for data-parallelism, such as for example the distributed unstructured mesh proposed by PaMPA [START_REF] Lachat | PaMPA: Parallel Mesh Partitioning and Adaptation[END_REF], will be the subject of future work.

For the task parallelism, we use OpenMP [START_REF] Dagum | Openmp: an industry standard api for shared-memory programming[END_REF]. OpenMP targets shared-memory platforms only. However, as the level of parallelism introduced by task parallelism technique is low compared to data parallelism, shared-memory, or intra-node parallelism is a good architecture choice Inria for task parallelism.While the version 4 of OpenMP has introduced explicit support for tasks, our implementation only requires version 3 whose fork-join model is well suited for the static scheduling introduced in Section 6. The use of dynamic schedulers such as provided by libgomp2 , StarPU [START_REF] Augonnet | StarPU: a unied platform for task scheduling on heterogeneous multicore architectures[END_REF], or XKaapi [START_REF] Gautier | Xkaapi: A runtime system for dataow task programming on heterogeneous architectures[END_REF] to directly execute the DAG Γ dep will be the subject of some future work.

As a result, the MSL compiler generates an hybrid code which uses both SkelGIS and OpenMP. It also generates empty components where the user must provide local sequential implementations of the kernels using well dened interfaces of the component to iterate over the local data subdomain. These implementations can of course be written by hand, but ideally one should also be able to generate them using a stencil compiler such as Pochoir or PATUS thus enabling the combinations of approaches used by those compilers. The following section will evaluate this compiler as well as the code it generates.

Evaluation

This section evaluates the MSL compiler and the code it generates. We rst evaluate the generated code in the case where only data parallelization and fusion are applied and compare this to a full SkelGIS program. Then, we analyze a case taking advantage of the full the hybrid parallelization combining data and task parallelism and validate the performance model proposed in Section 6.1. All evaluations presented in this section are based on a real case study of the shallow-Water Equations as solved in the FullSWOF2D 3 [START_REF] Coullon | The SIPSim implicit parallelism model and the SkelGIS library[END_REF][START_REF] Ferrari | A new two-dimensional shallow water model including pressure eects and slow varying bottom topography[END_REF] code from the MAPMO laboratory, University of Orléans. We have developed a MSL version of FullSWOF2D that contains 3 mesh entities, 7 computation domains, 48 data and 98 computations (32 stencil kernels and 66 local kernels).

Compiler

Table 1 illustrates the execution time of each step of the MSL compiler for the FullSWOF2D example. This has been computed on a laptop with a dual-core Intel Core i5 1.4 GHz, and 8 GB of DDR3. While the overall time of 4.6 seconds remains reasonable, one can notice that the computation of the T SP tree is by far the longest step. As a matter of fact, the complexity of the algorithm for N-shapes removal is O(n 3 ). The replacement of the static scheduling by a direct scheduling of the DAG by dedicated tools should solve this in the future.

Step In this evaluation, we disable task-parallelism (x its degree to 1) to focus on data-parallelism.

In the implementation evaluated in this paper, the data parallelism is handled using Skel-GIS that has already been evaluated in comparison with a native MPI implementation for the FullSWOF2D example [START_REF] Coullon | The SIPSim implicit parallelism model and the SkelGIS library[END_REF]. In addition to SkelGIS, MSL automatically inserts the required MPI synchronizations and identies possible fusion of kernels to reduce cache misses.

Cluster TGCC Curie Thin Nodes This section compares the performance of the code produced by MSL for FullSWOF2D with a plain SkelGIS program where synchronizations and fusion choices have been done by the developer without any automatic support. This enables to evaluate both the choices made by the compiler as well as the potential overheads of using L 2 C [2] that is not used in the plain SkelGIS version. The evaluations have been performed on the Curie supercomputer (TGCC, France) described in Table 2. Each evaluation has been performed ten times and the median is presented in results.

Weak scaling Figures 8 and9 show a weak scaling for a 400 × 400 points by core and a 800 × 800 points by core from 16 cores to 16.384 cores. Minimum and maximum values are also shown as error bars in gures.

One can notice that the MSL code produces a slightly better execution time on the 400 × 400 domain and a small overhead on the 800 × 800 domain. This might be explained by slightly dierent compiler optimization ags since L 2 C requires compilation as dynamic libraries with the -fpic compilation ag while the SkelGIS version does not. This ag can have slight positive or negative eect on code performance depending on the situation and might be responsible for the observed dierence. In any case, these dierence are barely noticeable and often included in the error bar of the measure and both version globally perform similarly. The MPI synchronizations performed by MSL perform as good as manually specied ones, and L 2 C does not introduce any noticeable overhead.

Strong scaling Figure 10 shows the number of iteration per second for a 10k×10k global domain size from 16 to 16.384 cores. The ideal strong scaling is also illustrated.

First, one can notice that the strong scaling evaluated for the code generated by MSL is close to the ideal speedup up to 16.384 cores, which is a very good result. Moreover, no overheads are introduced by MSL which shows that automatic synchronization detections and automatic fusion detections are the same one that the one written manually into the SkelGIS code of FullSWOF2D. Finally, no overheads are introduced by components of L 2 C. A small behavior dierence can be noticed with 2 9 = 512 cores, however this variation is no longer observed with 1024 cores.

Fusion Finally, to evaluate the data parallelization technique automatically introduced by MSL, the fusion optimization is evaluated. Figure 11 shows the number of iterations per second as a function of the number of cores, for FullSWOF2D with and without fusion optimization and onto a 500 × 500 domain. As explained in Section 6.3, the MSL fusion happens at a high level and is most of the time done naturally by a computer scientist. However, for a non computer scientist which write its numerical codes, an automatic proposition of such fusions makes the implementation easier. Moreover, one can notice that the performance is clearly improved (around 40%) by this fusion.

Hybrid parallelism

In this evaluation, we introduce task parallelism to evaluate the hybrid parallelization oered by MSL that mixes data parallelization and task parallelization. In the implementation evaluated in this paper and in addition to SkelGIS, the hybrid parallelization relies on OpenMP. Table 3: Parallelism level (number of parallel tasks) and the number of times this level appears.

The series-parallel tree decomposition T SP of this simulation, extracted by the MSL compiler, is composed of 17 sequence nodes and 18 parallel nodes. Table 3 represents the number of time a given level of parallelism, i.e., the number of tasks to perform concurrently, is observed in the nal tree. One can notice that the level of task parallelism extracted from the Shallow water equations is limited by two sequential parts in the application (level 1). Moreover, a level of 16 parallel tasks is reached two times, and ve times for the fourth level. This means that if two threads are dedicated to task parallelism, two parts of the code at least will not take advantage of this and that no part of the code would benet from more than 16 threads. The task parallelism alone is therefore clearly not enough to take advantage of even a single node of a cluster that typically supports today between 16 and 32 threads (16 for Curie).

On the other hand, Figure 12 illustrates limitations of data parallelization technique alone. While the computation time decreases linearly with the number of core used, the communication behavior is much more erratic. Between 2 and 16 cores, the communication happen inside a single node and the time is small and nearly constant. There is a small oscillation that might be explained by dierence between number of processors that are power of 4, where the communications are all the same, and those that are nota power of 4, where there is an imbalance. Starting with 32 cores, there are multiple nodes and the communication time is typically modeled as L + S/B where L is the latency, S the data size and B the bandwidth. This explains the decrease of time from 32 to 128 cores where the data sizes communicated by each process decreases. The increased observed after that might be due to the fact that with the increased number of processes the fat-tree becomes deeper and the latencies increase.

All in all, when the number of core used increases, the computation/communication ratio becomes poorer and poorer. As a result, the data parallelism alone fails to provide enough parallelism to leverage the whole machine and other sources of parallelism have to be found. The blue curve of Figure 13 shows the total execution time for the same example as in Figure 12 and as expected, the speedup bends down from 256 to 2048 cores.

In addition to the blue curve, Figure 13 shows strong scalings for the same example but using an hybrid parallelization. For example, the purple curve shows the parallelization which uses 8 cores, to perform the task parallelization, for each process used for data parallelization (i.e., MPI process). As a result, for example, when using 2 machines of the TGCC cluster, with a total of 32 cores, 4 cores are used for SkelGIS MPI processes, for data parallelization, and for each one 8 cores are used for task parallelization (4 × 8 = 32).

As a result, and as explained in Section 6.2, data is less divided into sub-domains and the eect which is observe onto the blue curve is delayed. This gure shows a comparison with 2, 4, 8 and 16 cores per MPI process for task parallelization.

From 2 to 8 cores, the improvement of the strong scaling is clear. However, reaching 16 cores, an important initial overhead appears; moreover, the curve bends down rapidly instead of improving the one with 8 cores for task parallelization. Two dierent phenomena happen.

First, thin nodes of the TGCC Curie are built with two NUMA nodes each of 8 cores. As a result, when increasing the number of OpenMP cores for task parallelization from 8 to 16 cores, an overhead is introduced by exchanges of data between memories of the two NUMA nodes. This phenomena is illutrated in Figure 14. In this gure, a dierent strategy is used to bind threads onto available cores (using OpenMP). This strategy, called spread (instead of close in Figure 13), binds threads on cores in order to spread as much as possible onto resources, which means that the two NUMA nodes are directly used whatever the number of cores kept for task parallelization. As a result, and as shown in the gure, using 2, 4 and 8 cores an initial overhead is introduced.

The second phenomena, which happens in Figure 13 with 16 cores for tasks parallelization, is due to the level of parallelization introduced by the task parallelization technique. Actually, as illsutrated in Table 3, only two steps of the T SP static scheduling generated by the MSL compiler can take advantage of 16 cores among a total of 18 steps. This phenomena has been explained iterations per second Ideal MSL data parallelization only MSL using 2 cores for tasks MSL using 4 cores for tasks MSL using 8 cores for tasks MSL using 16 cores for tasks A spread OpenMP clause is used to bind threads onto cores. in Section 6.2 by the variable F task and the fact that it is not always true that F task = P task . This explains why using 16 cores for task parallelization in Figure 14 is still less ecient than using 8 cores even if the two NUMA nodes are always used in this evaluation.

Finally, to completely validate the performance model introduced in Section 6.2, and to understand when the hybrid parallelization becomes more interesting than the data parallelization, Figure 15 represents T COM 1 and T COM 2 + T task of Equation [START_REF] Coullon | Implementation and performance analysis of skelgis for network mesh-based simulations[END_REF], for the best case (i.e., when 8 cores are used for task parallelization), and with a close OpenMP bind of threads onto cores. 4 perfectly matches results observed in Figure 13 for 8 cores used for task parallelization per core used for data parallelization. As a result, the hybrid parallelization becomes better with a total of 512 cores.

Conclusion

We have presented the Multi-Stencil Language (MSL), a domain specic language for stencilbased numerical simulations. Compared to other existing stencil DSLs, MSL proposes a higher abstraction level for end-users. MSL oers a way to describe a numerical simulation independently from the type of mesh used. As a result, the language can be used for a larger number of simulations. It helps to clearly separate the description of the simulation from implementation choices, which facilitates reuse of existing languages and libraries, for example languages for data structures optimizations and parallelization strategies. In this paper, the formal denition of a multi-stencil program has been given. This formalism helps to understand what is dened into the light grammar of MSL, and also how parallelization can be extracted from it. MSL produces an empty parallel pattern of the numerical simulation described. This pattern indicates where synchronizations are needed, but not how they are performed, which is dependent from the type of mesh and from implementation choices. This second step is left to other existing languages, which are chosen as a back-end of MSL. For example, in this paper the chosen back-end is a code using SkelGIS [57], a templated C++ language for distributed meshes (using MPI), and OpenMP [START_REF] Dagum | Openmp: an industry standard api for shared-memory programming[END_REF].

Experiments presented in this paper show that the back-end code produced by MSL, i.e., by a mesh-agnostic language, does not introduce overheads up to 16.384 cores, and that the hybrid (data and task) parallelism automatically introduced by the language, which combines SkelGIS and OpenMP as a back-end, improves performance compared to a data parallelization only. Thus, with the same code asked to the user, two dierent parallel versions of the simulation can be obtained. However, if performances are improved by the hybrid parallelization, this performance is also limited by the scheduling strategy used in this paper. Actually, as described in Section 6, a static series-parallel (fork-join) scheduling is performed in this work. This helps to dene a predictive performance model, which has been validated in experiments. However, getting rid of global synchronizations (join) could be an interesting improvement to get better performance, especially for irregular simulations, where tasks are not naturally balanced. Our immediate perspective, thus, is to use OpenMP dynamic schedulers [START_REF] Augonnet | StarPU: a unied platform for task scheduling on heterogeneous multicore architectures[END_REF][START_REF] Gautier | Xkaapi: A runtime system for dataow task programming on heterogeneous architectures[END_REF][START_REF] Wu | Hierarchical DAG Scheduling for Hybrid Distributed Systems[END_REF] as back-ends to study performance improvments. Using such dynamic schedulers also opens to more heterogeneous architectures, such as GPUs or many-cores accelerators.

Finally, to validate the separation of concerns introduced by MSL, a second back-end for unstructured meshes could be interesting. We think about PaMPA [START_REF] Lachat | PaMPA: Parallel Mesh Partitioning and Adaptation[END_REF] as a distributed data structure backend, instead of SkelGIS. Thus using a single language, the MSL language, two dierent kinds of numerical simulations could be dened: one onto Cartesian meshes; and one on unstructured meshes. This work has partially been supported by the PIA ELCI project of the French FSN. This work was granted access to the HPC resources of TGCC under the allocations t2015067470 and x2016067617 made by GENCI.

Figure 1 :Denition 2

 12 Figure 1: From left to right, Cartesian, curvilinear and unstructured meshes.

  A is computed with a 2-neighborhood stencil applied on C. A is computed onto a computation domain includes all entities of the group.

Figure 2 :

 2 Figure 2: (a) a Cartesian mesh and two kind of groups of mesh entities, (b) an example of stencil kernel on cells, (c) an example of stencil kernel on two dierent entities of the mesh.

Figure 3 :

 3 Figure 3: Grammar of the Multi-Stencil Language.

mesh : c a r t mesh e n t i t i e s : c e l

  n c i l s h a p e s : ncc from c e l l t o c e l l nce from c e l l t o edgex nec from edgex t o c e l l mesh q u a n t i t i e s : c e l l A, B, D, E , F ,G, I , J edgex C,H s c a l a r s : mu, tau time : 500 c o m p u t a t i o n s : B

R

  = {(w i , n j )} 15:

R

  = {(w i , n j )} 25:

  G ∪ F dep (Γ sync , j + 1)14:

Figure 5 : 4 6

 54 Figure 5: Γ dep of the example of program of Figure 4

6 Figure 7 :

 67 Figure 7: Serie-Parallel tree decomposition of the example of program of Figure 4

Figure 8 :

 8 Figure 8: weak-scaling with 400 × 400 domain per core.

Figure 9 :

 9 Figure 9: weak-scaling with 800 × 800 domain per core.

Inria

Figure 10 :

 10 Figure 10: Strong scaling on a 10k × 10k domain.

Figure 11 :

 11 Figure 11: Strong scaling on a 500x500 domain, with and without the fusion optimization.

Inria

Figure 12 :

 12 Figure 12: Computation vs communication times in the data parallelization technique.

Figure 14 :

 14 Figure 14: Strong scaling comparisons between data parallelization and hybrid parallelization.

Figure 15 :

 15 Figure 15: Execution times (seconds) of T COM 1 and T COM 2 + T task for 8 cores for task parallelization. Verication of the Equation (7).

Figure 15

 15 Figure15and Table4perfectly matches results observed in Figure13for 8 cores used for task parallelization per core used for data parallelization. As a result, the hybrid parallelization becomes better with a total of 512 cores.

  == D j and R i ∩ R j = ∅ then

	8:	end if
	9:	else if parent(k i )==P then
	10: if D i 11: propose the fusion k f us i;j 12: end if
	13:	end if
	14:	end if
	15:	end if
	16:

Table 1 :

 1 Execution times of the MSL compiler

			Parser Γ sync Γ dep	T SP
		Time (ms)	1	2	4.2	3998.5
		%	0.022 0.043 0.09	86.6
	8.2	Data parallelism

Table 2 :

 2 Hardware conguration of TGCC Curie Thin nodes.

	Processor	2×SandyBridge
		(2.7 GHz)
	Cores/node	16
	RAM/node	64 GB
	RAM/core	4GB
	Compiler [-O3]	gcc 4.9.1
	MPI	Bullxmpi
	Network	fat-tree Inniband

Table 4

 4 gives time details.

	Inria

Table 4 :

 4 Execution times (seconds) of T COM 1 , T COM 2 and T task for 8 cores for task parallelization. Verication of the Equation[START_REF] Coullon | Implementation and performance analysis of skelgis for network mesh-based simulations[END_REF].

	2 -11	16	32	64	128	256	512	1024	2048
						cores

urlwww.top500.org RR n°8962

RR n°8962

The Multi-Stencil Language From the formalism detailed in the previous section, the Multi-Stencil Language and its grammar can already be given. This grammar is sucient to automatically extract parallelism from the program as will be explained in the next section. Inria

Inria

Coullon & Bigot & Perez

https://gcc.gnu.org/projects/gomp/

RESEARCH CENTRE GRENOBLE -RHÔNE-ALPES


C Bouveyron
email: charles.bouveyron@parisdescartes.fr

J Jacques
email: julien.jacques@univ-lyon2.fr

funFEM: an R package for functional data clustering

Keywords:

Introduction

A new model-based clustering algorithm for times series (or more generally functional data), called FunFEM, has been proposed in [START_REF] Bouveyron | The discriminative functional mixture model for the analysis of bike sharing systems[END_REF]. It is based on a functional mixture model which allows the clustering of the data in a discriminative functional subspace. This model presents the advantage to be parsimonious and to allow the visualization of the clustered systems. This paper presents the funFEM package for R which implements this new clustering algorithm.

The discriminative functional mixture model

Let {x 1 , ..., x n } be the observed curves we want to cluster. In practice, the functional expressions of these curves are not known and we only have access to the discrete observations x i (t is ) at a finite set of ordered times {t is : s = 1, . . . , m i }. The first step in functional data analysis usually consists in recovering the functional nature of data, and for this we assume that the observed curves can be decomposed in a finite basis of function

x i (t) = p j=1 γ ij ψ j (t).
The functional mixture model proposed in [START_REF] Bouveyron | The discriminative functional mixture model for the analysis of bike sharing systems[END_REF] assumes that the basis expansion coefficient γ i = (γ i1 , . . . , γ ip ) t of curve x i follows a mixture of Gaussians:

p(γ) = K k=1 π k φ(γ; U µ k , U t Σ k U + Ξ), (1)
where π k is the prior probability of the kth group, φ is the standard Gaussian density function, U is a p × d matrix mapping the coefficient γ into the discriminative subspace (of dimension d < p), µ k and Σ k are the mean vector and covariance matrice (for cluster k) of the mapping of γ into the discriminative subspace, and Ξ the noise covariance matrix.

Model inference is based on an EM-like algorithm, including an additional step between the traditional E and M steps in which the orientation matrix U is updated. This algorithm is implemented in the funFEM package for R, available on the CRAN. The main function, funFEM has only two mandatory arguments: the functional data, defined as a functional object of the fda package, and the number of clusters (or a vector of). The outputs of funFEM are: the posterior probabilities and the estimated clusters, the model parameter estimation and several model selection criteria. These latter can be used in order to choose the optimal number of clusters. The use of funFEM is now illustrated on the Velov dataset (available in the funFEM package), which contains one week of loading curves (i.e. the proportion of available bikes) of the bike sharing system stations of Lyon (called Velov). The number of Velov stations in Lyon is 345, and the curves are sampled approximatively each hour. Due to the periodic nature of the curves, a Fourier basis is considered. Below is given the R script used to launch the clustering analysis with funFEM. # Load the velov data and smoothing R> library(funFEM) R> data(velov) R> basis <-create.fourier.basis(c(0, 181), nbasis=25) R> fdobj <-smooth.basis(1:181,t(velov$data),basis)$fd # Clustrering with FunFEM R> res = funFEM(fdobj,K=6,model="AkjBk",init="kmeans",lambda=0,disp=TRUE) Figure 1 presents the resulting mean curves of each cluster as well as the geographical repartition of clusters. Such results typically allow us to analyse the bike sharing system of Lyon.