# InKS, a Programming Model to Decouple Algorithm from Optimization in HPC Codes

Ksander Ejjaouani, Olivier Aumage, Julien Bigot, Michel Mehrenberger, Hitoshi Murai, Masahiro Nakao, Mitsuhisa Sato

► To cite this version:

Ksander Ejjaouani, Olivier Aumage, Julien Bigot, Michel Mehrenberger, Hitoshi Murai, et al. InKS, a Programming Model to Decouple Algorithm from Optimization in HPC Codes. Journal of Supercomputing, 2019, 10.1007/s11227-019-02950-2. http://hal-02281963

HAL Id: hal-02281963

https://inria.hal.science/hal-02281963

Submitted on 9 Sep 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# InKS, a Programming Model to Decouple Algorithm from Optimization in HPC Codes

Ksander Ejjaouani · Olivier Aumage

Julien Bigot · Michel Méhrenberger

Hitoshi Murai · Masahiro Nakao

Mitsuhisa Sato

**Abstract** Existing programming models tend to tightly interleave algorithm and optimization in HPC simulation codes. This requires scientists to become experts in both the simulated domain and the optimization process and makes the code difficult to maintain or port to new architectures. In this paper, we propose the INKS programming model that decouples these concerns with two distinct languages: INKS<sub>pia</sub> to express the simulation algorithm and INKS<sub>pso</sub> for optimizations. We define INKS<sub>pia</sub> and evaluate the feasibility of defining INKS<sub>pso</sub> with three test-languages: INKS<sub>o/C++</sub>, INKS<sub>o/loop</sub>, INKS<sub>o/XMP</sub>. We evaluate the approach on synthetic benchmarks (NAS and heat equation) as well as on a more complex example (6D Vlasov-Poisson solver). Our evaluation demonstrates the soundness of the approach as it improves the separation of algorithmic and optimization concerns at no performance cost. We also identify a set of guidelines for the later full definition of the INKS<sub>pso</sub> language.

**Keywords** programming model · separation of concerns · HPC · DSL

Ksander Ejjaouani

Maison de la Simulation, CEA, CNRS, Univ. Paris-Sud, UVSQ, Université Paris-Saclay

Inria

Gif-sur-Yvette, France

E-mail: ksander.ejjaouani@inria.fr

Olivier Aumage

Inria, LaBri, Bordeaux, France

E-mail: olivier.aumage@inria.fr

Julien Bigot

Maison de la Simulation, CEA, CNRS, Univ. Paris-Sud, UVSQ, Université Paris-Saclay

Gif-sur-Yvette, France

E-mail: julien.bigot@cea.fr

Michel Mehrenberger

Université de Marseille, Marseille, France

E-mail: mehrenbe@math.unistra.fr

Hitoshi Murai · Masahiro Nakao · Mitsuhisa Sato

Riken CCS, Kobe, Japan

E-mail: {h-murai,masahiro.nakao,msato}@riken.jp
# 1 Introduction

It is more and more common to identify simulation as the third pillar of science together with theory and experimentation. Parallel computers provide the computing power required by the more demanding of these simulations. The complexity and heterogeneity of these architectures do however force scientists to write complex code (using vectorization, parallelization, accelerator specific languages, etc.). These optimizations heavily depend on the target machine, the code has to be adapted whenever it is ported to a new architecture.

As a result, scientists have to become experts in the art of computer optimizations in addition to their own domain of expertise. It is very difficult in practice to maintain a code targeting multiple distinct architectures. One fundamental cause for this situation is the tight interleaving of two distinct concerns imposed by most programming models. On the one hand, the algorithm comes from the expertise of the domain scientists and does not depend on the target architecture. On the other hand, optimizations form another domain of expertise and has to be adapted for a given architecture. Thereby, both algorithm and optimizations concerns are expressed within a single code. This mix impedes simulation codes maintainability and readability while hindering developers productivity.

Many approaches have been proposed to improve this situation in the form of libraries or languages [16, 17, 20, 18]. Approaches based on automated optimization processes typically isolate the algorithmic aspects well, but restrict their domain of applicability and/or the range of supported optimizations. Approaches based on optimization tools and libraries enable optimization specialists to express common optimizations efficiently but leave others mixed with the algorithm.

In this paper, we propose the Independent Kernel Scheduling (INKS) programming model to separate algorithm from optimization choices in HPC simulation codes. We define the INKS<sub>pia</sub> language used to express the algorithm of an application independently of its optimization. This separation aims to improve the readability and maintainability of codes while easing portability and new optimization expression. This approach is used for common optimizations while INKS<sub>o/C++</sub> for less common optimizations. Such a program can then be optimized using INKS<sub>o/XMP</sub> and INKS<sub>o/loop</sub>, two domain-specific languages (DSLs) which ask for optimization information only. While these DSLs target some common optimizations, INKS<sub>o/C++</sub> can be used for less common ones.

This paper makes the following contributions: 1) it defines the INKS programming model and its platform-independent algorithmic language INKS<sub>pia</sub>; 2) it proposes an implementation of INKS and tests the INKS<sub>pso</sub> approach with three optimization DSLs, INKS<sub>o/C++</sub>, INKS<sub>o/loop</sub> and INKS<sub>o/XMP</sub>; and 3) it evaluates the approach on the synthetic NAS parallel benchmarks [3] and on the 6D Vlasov-Poisson solving with a semi-Lagrangian method.

The remaining of the paper is organized as follows. Section 2 analyzes related works. Section 3 describes INKS and its implementation. Section 4 shows
the use of INKS on a 6D Vlasov-Poisson solver. Section 5 evaluates our approach. Section 6 concludes the paper.

## 2 Related Work

We now present approaches currently available to scientific programmers to help implementing simulation applications. A first widely used approach is based on imperative languages such as Fortran or C. Libraries like MPI extend this to distributed memory with message passing. Abstractions very close to the execution machine make fine-tuning possible to achieve good performance on any specific architecture. It does however require encoding complex optimizations directly in the code. As there is no language support to separate the algorithm and architecture-specific optimizations, tedious efforts have to be applied [13] to support performance portability. Algorithm and optimizations are instead often tightly bound together in codes.

A second approach is offered by tools (libraries, frameworks or language extensions) that encode classical optimizations. *OpenMP* [5], *REPARA* [8] or *Kokkos* [4] support common shared-memory parallelization patterns. For example, Kokkos offers multidimensional arrays and iterators for which efficient memory mappings and iteration orders are selected independently. *UPC* [9] or XMP [17] support the partitioned global address space paradigm. For example, in XMP, directives describe array distribution and communications between nodes. These tools offer gains of productivity when the optimization patterns they offer fit the requirements. The separation of optimizations from the main code base also eases porting between architectures. Even if expressed more compactly optimizations do however remain mixed with the algorithm. For instance, in OpenMP or REPARA, parallel concerns are specified on top of an existing code which already carries optimization choices, such as loops order.

A third approach pushes this further with tools that automate the optimization process. For example, *PaRSEC* [12] or *StarPU* [1] support the many-tasks paradigm. In StarPU, the user expresses its code as a directed acyclic graph (DAG) of tasks with data dependencies that is automatically scheduled at runtime depending on the available resources. Other examples are *SkeTo* [21] or *Lift* [19] that offer algorithmic skeletons. Lift offers a limited set of parallel patterns whose combinations are automatically transformed by an optimizing compiler. Automating optimization improves productivity and clearly separate these optimizations which improves portability. The tools do however not cover the whole range of potential optimizations such as the choice of work granularity inside tasks in StarPU for example. The algorithm remains largely interleaved with optimization choices even with this approach.

A last approach is based on DSLs that restrict the developer to the expression of the algorithm only while optimizations are handled independently, such as *Pochoir* [20] or *PATUS* [6], DSLs for stencil problems. In Pochoir, the user specifies a stencil (computation kernel and access pattern), boundary
conditions and a space-time domain while all optimizations are handled by a
compiler. These approaches ensure a very good separation of concerns. The
narrower the target domain is, the more efficient domain and architecture-
specific optimizations are possible. However, it makes it less likely for the tool
to cover the needs of a whole application. On contrary, the wider the target
domain is, the less efficient optimizations are possible.

To summarize, one can consider a continuum of approaches from very gen-
eral approaches where the optimization process is manual to more and more
domain specific where the optimization process can be automated. The more
general approaches support a large range of optimizations and application do-
mains but yield high implementation costs and low separation of concerns and
portability. The more automated approaches reduce implementation costs and
offer good separation of concerns and portability but restrain the range of sup-
ported domains and optimizations. Ideally, one would like to combine all these
advantages: 1) the domain generality of imperative languages, 2) the ease
of optimization offered by dedicated tools and 3) the separation of concerns
and performance portability offered by DSLs with 4) the possibilities of fine
and manual optimizations offered by both imperative languages and dedicated
tools. The following section describes the INKS programming model that aims
to combine these approaches to offer such a solution.

3 The InKS programming model

This section first introduces the design of our INKS programming model, based
on the use of distinct languages to express the algorithm and optimization
choices separately. It then present a prototype implementation of the model
composed of INKS<sub>pia</sub>, the algorithm language, and two INKS<sub>pso</sub> optimization

languages. The simulation *algorithm* consists in the set of values computed, the formula used to produce them as well as the simulation inputs and outputs. *Optimization choices* include all that is not the algorithm: e.g. the computing unit selected for each computation, their ordering, the memory location for each value, etc. Multiple optimization choices can differ in performance but simulation results depend on the algorithm only. The INKS approach is summarized in Figure 1. The INKS<sub>pia</sub> language is used to express the algorithm with no concern for optimization choices. A compiler can automatically generate non-optimized choices from an INKS<sub>pia</sub> specification,

<figure><img src="image_1.png" /><figcaption>Fig. 1: The INKS model</figcaption></figure>

mostly for testing purposes. The INKS<sub>pso</sub> language is used to define
optimizations only, while other information is gathered from the INKS<sub>pia</sub> code.
Many versions of the optimization choices can be devised for a single algorithm; for example to optimize for multiple targets.

We now describe the $INKS_{pia}$ language, and proposes three preliminary DSLs to validate strategies for the design of the $INKS_{pso}$ language: $INKS_{o/XMP}$ handles domain decomposition, $INKS_{o/loop}$ focuses on efficient loop nest and $INKS_{o/C++}$ enables to use C++ and write arbitrary complex optimization.

## 3.1 The $INKS_{pia}$ language.

In $INKS_{pia}$ [2] values are stored in infinite multidimensional arrays based on dynamic single assignment (DSA, each coordinate can only be written once). Memory placement of each coordinate is left unspecified. Computations are specified by *kernel* procedures that **1)** take as parameters data arrays and integer coordinates; **2)** specify the coordinate they might read and will write in each array; and **3)** define either a C++ or INKS implementation. An INKS implementation defines kernels *validity domains*: coordinates where C++ kernels can generate values in arrays. Kernel execution order is left unspecified. The simulation entry point is a kernel marked *public*. Listing 1 presents a simple $INKS_{pia}$ code: a 1D 3-point stencil computation. The simulation consists of one 2D logical array, *Array* (line 12), two kernels, *stencil3* (line 1) and *boundary* (line 6), and is parameterized by two integers, *X* and *T* (line 11). Line 12 specifies that the simulation starts with a subset of *Array* (every values in the space dimension, at the first time-step) and expects, as output, another subset: every values in the space dimension at the last time-step. To do so, it can call the *stencil3* and *boundary* kernels with a specific set of values for their integer parameters *x* and *t* (validity domain) and with *Array* as their logical array parameter, as expressed on line 15 and 16.

A $INKS_{pia}$ code specifies a parameterized task graph (PTG) [7]. This graph is encoded using the Polyhedron Model [10]. It provides a compact representation of static control programs. This representation covers a large range of problems but imposes a few limitations. Mostly, all the problem parameters must be known at launch time, it does for example not support adaptive mesh or time-steps. It is still possible to express these concerns outside $INKS_{pia}$ and call the INKS implementation multiple times with different parameters.

Once the algorithm is specified in $INKS_{pia}$, the goal is to write the optimization choices meaning choosing a memory layout and a scheduling of the kernels. Two approaches exist to produce these choices: the automatic compiler or $INKS_{pso}$. All approaches rely on the *Kernel* file. This file is generated by the $INKS_{pia}$ compiler which translates the $INKS_{pia}$ kernels to C++ functions. The $INKS_{pia}$ compiler can produce generic optimization choices with a valid but non-optimized computations scheduling and memory allocations to execute them. Scheduling and memory layouts are computed using the *Integer Set Library* [22] and recent works on modular mapping in the polyhedron model [14]. Arbitrarily complex versions of optimization choices can also be written manually in plain C++. These functions can be called from any ex-
```c
1 kernel stencil3(x, t) : (
2   double A(2) {in: (x-1:x+1, t-1) | out: (x, t)}
3 )
4 #CODE(C) A(x, t+1) = 0.5*A(x, t-1) + 0.25*(A(x-1, t-1)+A(x+1, t-1)); #END
5
6 kernel boundary(x, t) : (
7   double A(2) {in: (x, t-1) | out (x, t)}
8 )
9 #CODE(C) A(x,t) = A(x, t-1); #END
10
11 public kernel inks_stencil(X, T) : (
12   double Array(2) {in: (0:X, 0) | out: (0:X, T-1)}
13 )
14 #CODE(INKS)
15 stencil3 (1:X-1, 1:T) : (Array),
16 boundary {(0, 1:T), (X-1, 1:T)} : (Array)
17 #END
```

<caption>Listing 1: 1D stencil computation on a 2D domain in INKS<sub>pia</sub></caption>

isting code whose language supports the C calling convention. However, that approach requires information present in INKS<sub>pia</sub> to be repeated. The INKS<sub>pso</sub> DSL thus interface the optimization process with INKS<sub>pia</sub>, offering the optimization specialists to specify optimization only.

## 3.2 INKS<sub>pso</sub> DSLs implementation

Implementing a complete optimization DSL for the INKS model is a long-term objective. In this study, we present the two DSLs INKS<sub>o/XMP</sub> and INKS<sub>o/loop</sub>. Both DSLs solely describe optimizations, and get missing information from INKS<sub>pia</sub> code.

INKS<sub>o/XMP</sub> (illustrated in Listing 2) handles distributed memory domain decomposition by combining C and directives based on XMP and adapted for INKS. The compiler replaces these directives by C and XMP code. The inks decompose directive supports static or dynamic allocation of logical arrays described in the algorithm. The domain size is extracted from INKS<sub>pia</sub> source and the user only has to specify its mapping onto memory. As in XMP, INKS<sub>o/XMP</sub> supports domain decomposition mapped onto an XMP topology. In INKS<sub>pia</sub> code, there are no concerns for memory optimization such as dimension ordering or memory reuse. Therefore, INKS<sub>o/XMP</sub> supports dimension reordering and folding which consists in reusing the same memory address for subsequent indices in a given dimension. The exchange directive supports halo exchanges. The user specifies which dimension should be exchanged and which computational kernel will be executed after the exchange. From these informations and the INKS<sub>pia</sub> kernel specification, the INKS<sub>o/XMP</sub> compiler then computes the halo size. While XMP requires halo values to be stored contiguously with the domain, INKS<sub>o/XMP</sub> supports a dynamic halo extension where halo values are stored in dedicated, dynamically allocated buffers to reduce memory footprint.

INKS<sub>o/loop</sub> (illustrated in Listing 3) offers to specify manually loop nests for which the compiler generates plain C++ loops. Plain C++ is usable in combination with INKS<sub>o/loop</sub>. The loop keyword introduces a nest optimization.
```c
1 /***InKSo/XMP specification*/
2 #pragma xmp nodes p[nNodes]
3 void inks_stencil(T& Input, size_t X, size_t T){
4 #pragma inks decompose % Array // allocation of "Array" algorithmic array
5 (2*2, 1) // dimension reordering, time dimension is fold by 2 and not distributed
6 with t onto p // block decomposition mapped on the XMP topology
7
8 for(int t=1; t<T; t++){
9 boundarychristain(X, 0, t);
10 for(int x=1; x<X-1; x++)
11 stencil3christain(X, x, t);
12 boundary(X, X-1, t);
13 }
14 }
15
16 /***XMP + C result*/
17 #pragma xmp nodes p[nNodes]
18 #pragma xmp template t[:][::] // xmp 1d logical array
19 #pragma xmp distribute t[*][block] onto p
20 void inks_stencil(T& Input, size_t X, size_t T){
21 #pragma xmp align Array[t][x] with t[t][x]
22 #pragma xmp templatehil fix t[2][X]
23 Array = (double(*)[X]) xmp_malloc(xmp_desc_of(f6d), 2, X);
24 for(int t=1; t<T; t++){
25 boundary(X, 0, t);
26 for(int x=1; x<X-1; x++)
27 stencil3(X, x, t);
28 boundary(X, X-1, t);
29 }
30 }
```

Listing 2: INKS<sub>o</sub>/XMP optimization choices of the INKS<sub>pia</sub> code presented in Figure 1 and the C + XMP result

```c
1 /*** stencil.iloop InKSo/Loop optimization choices file */*
2 loop stencil3_loop(t, X) : stencil3 { // set "t" value
3 // "Set" not specified -> loop bounds are computed, with a fixed "t"
4 Order: x; // order of the loop
5 }
6
7 /*** C++ result */
8 void stencil3_loop(T& Array, size_t t, size_t X){
9 for (int x=0; x<X; x++)
10 stencil3(X, x, t);
11 }
12
13 /*** InKSo/Loop result used */
14 void inks_stencil(T& Input, size_t X, size_t T){
15 for(int t=1; t<T; t++){
16 boundary(X, 0, t);
17 stencil3_loop(X, t, X);
18 boundary(X, X-1, t);
19 }
20 }
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
```

Listing 3: INKS<sub>o</sub>/loop optimization choices of a loop from the INKS<sub>pia</sub> code presented in Figure 1 with its C++ result and use

tion with a name, the list of parameters from the algorithm on which the loop bounds depend and a reference to the optimized kernel. Loop bounds can be automatically extracted from INKS<sub>pia</sub>, but the Set keyword makes it possible to restrict these bounds. The Order keyword specifies the iteration order on the dimensions named according to the INKS<sub>pia</sub> code. The Block keyword enables the user to implement blocking. It takes as parameters the size of block for the loops starting from the innermost one. If there are less block
$$ 
\left\{
\begin{array}{l}
\displaystyle \frac{\partial f(t, x, v)}{\partial t} + v \cdot \nabla_x f(t, x, v) - E(t, x) \cdot \nabla_v f(t, x, v) = 0 \\
- \Delta \phi(t, x) = 1 - \rho(t, x) \\
E(t, x) = - \nabla \phi(t, x) \\
\rho(t, x) = \int f(t, x, v) dv
\end{array}
\right.
\qquad (1)
$$

<figure><img src="image_2.png" /><figcaption>Fig. 2: The 6D Vlasov-Poisson algorithm</figcaption></figure>

sizes than loops, the remaining loops are not blocked. The Buffer keyword supports copying data in a local buffer before computation and back after to ensure data continuity and improve vectorization. The compiler uses data dependencies from the INKS_pia code to check the validity of the loops order and generate vectorization directives whenever possible.

## 4 The 6D Vlasov/Poisson problem

The 6D Vlasov-Poisson equation, presented in (1), describes the movement of particles in a plasma and the resulting electric field. We study its resolution for a single species on a 6D Cartesian mesh with periodic boundary conditions. We solve the Poisson part using a fast Fourier transform (FFT) and rely on a Strang splitting (order 2 in time) for the Vlasov part. This leads to 6 1D advections: 3 in space dimensions ($x_1, x_2, x_3$) and 3 in velocity dimensions ($v_1, v_2, v_3$). Each 1D advection relies on a Lagrange interpolation of degree 4. In the space dimensions, we use a semi-Lagrangian approach where the stencil is not applied around the destination point but at the foot of characteristics, only known at runtime, as described in more details in [18].

The main unknown is $f$ (f6D in the code), the distribution function of particles in 6D phase space. Due to the Strang splitting, a first half time-step of advections is required after f6D initialization but before the main time-loop. These advections need the electric field $E$ as input. $E$ is obtained through the FFT-based Poisson solver that in turn needs the charge density $\rho$ as input. $\rho$ is computed by a reduction of f6D. The main time-loop is composed of 4 steps: advections in space dimensions, computation of the charge density (reduction) and electric field (Poisson solver) and advections in velocity dimensions. The algorithm of the simulation is presented in Figure 2.

The remaining of the section presents two optimizations implemented in the Selalib [16] version of the 6D Vlasov-Poisson problem.

The 6D nature of f6D requires a lot of memory, but the regularity of the problem means it can be distributed in blocks with good load-balancing.
<figure><img src="image_3.png" /><figcaption>Fig. 3: Dynamic halo exchange representation on a 2D domain</figcaption></figure>

Halos are required to hold values of neighbors for the advections. Connected halo zones would increase the number of points in all dimensions and consume too much memory. Split advections mean that halos are required in a single dimension at a time though. We therefore use dynamic halos composed of two buffers, one for each boundary of the advected dimension (denoted "right" and "left"). Figure 3 shows this optimization on a 2D domain. Listing 4 shows the INKS<sub>o</sub>/XMP implementation of this strategy on 6D Vlasov-Poisson.

```c
1 /** InKS<sub>o</sub>/XMP code for the f6d domain decomposition and dynamic halo exchange */*
2 #pragma xmp nodes p6d[pV3] [pV2] [pV1] [pZ] [pY] [pX]
3 #pragma inks decompose dynamic
4 % f6d      // dynamic allocation of f6d algorithmic array
5 (8:, 7:, 6, 5, 4, 3, 2, 1) // dimension reordering, dim 7 and 8 are folded
6 with t6d onto p6d // block decomposition mapped on the XMP topology
7 // Dynamic halo exchange on the 4th dimension, halo sizes are computed automatically
8 // R and L are now allocated buffers and contain the halo values
9 #pragma inks exchange periodic f6d(4, advection4) to R and L
10 foo(R, L);
11 /** Generated C+XMP code for the f6d domain decomposition */*
12 double (f6d*)[] [] [] []; // need to declare f6d global, valid in xmp
13 #pragma xmp nodes p6d[pV3] [pV2] [pV1] [pZ] [pY] [pX] // xmp 6d cartesian node topology
14 #pragma xmp template t6d[:]:[:]:[:]:[:]:// xmp 6d logical array
15 #pragma xmp distribute t6d [block][block][block][block][block][block] onto p6d
16 // map element of f6d to element of t6d
17 #pragma xmp align f6d[n][m][l][k][j][i] with t6d[n][m][l][k][j][i]
18 #pragma xmp template_fix t6d[N][M][L][K][J][I]
19 f6d = (double(*)[M][L][K][J][I]) xmp_malloc(xmp_desc_of(f6d), N, M, L, K, J, I);
```

Listing 4: 6D decomp. & halo exchange in INKS<sub>o</sub>/XMP and resulting code

Advections account for the main computational cost of the problem, up to 95% of the sequential execution time. Six loops surround the stencil computation of each advection and in a naive version, the use of a modulo to handle periodicity and application along non-contiguous dimensions slow down the computation. To enable vectorization and improve cache use, we copy f6D elements into contiguous buffers along with the *left* and *right* halos. Advections are applied on these buffers before copying them back into f6D. Blocking further improves performance by copying multiple elements at a time. Listing 5 corresponds to the INKS<sub>o</sub>/loop implementation of these optimizations in a sequential version and presents the generated code. The Poisson solver
1 /*** advec.iks InKSpia algorithmic file ***/
2 kernel advection3 (i, j, k, l, m, n, t, K, step) : (
3   double f6d {in: (i, j, 0:K, l, m, n, t, step-1) |out: (i, j, k, l, m, n, t, step)},
4   int disp {in: (n)},
5   double coef {in: (n, 0:4)}
6 )
7 #CODE (C)
8 f6d(i, j, k, l, m, n, t, step) =
9   coef(n, 0) * f6d(i, j, k-1+disp(n), l, m, n, t, step-1)
10  + coef(n, 1) * f6d(i, j, k+0+disp(n), l, m, n, t, step-1)
11  + coef(n, 2) * f6d(i, j, k+1+disp(n), l, m, n, t, step-1)
12  + coef(n, 3) * f6d(i, j, k+2+disp(n), l, m, n, t, step-1);
13 #END
14 public kernel main_code(t, I, J, K, L, M, N, Niter) : (
15   double coef(2),
16   int disp(1),
17   double f6d(6) {in: (i, j, 0:K, l, m, n, t, step-1) |out: (i, j, k, l, m, n, t, step)}
18 )
19 #CODE(INKS)
20 advection3 (0:I, 0:J, 0:K, 0:L, 0:M, 0:N, 1:Niter, K, 2) : (f6d, disp, coef)
21 /* ... */
22 #END
23 /*** advec.iloop InKSo/Loop optimization choices file ***/
24 loop advection3_loops(t, I, J, K, L, M, N, Niter) : advection3 { // set "t" value
25   // "Set" not specified -> loop bounds are computed, with a fixed "t"
26   Order: n, m, l, k, i, j; // order of the loop
27   Block: 16; // blocking on the inner dimension j
28   Buffer: f6d(3); // copy the third dimension of f6d to a 1d buffer
29 }
30 }

<caption>Listing 5: A loop nest in $INKS_{pia}$ optimized in $INKS_{o/loop}$</caption>

relies on a remapping scheme, where the domain decomposition is modified
between each FFT execution.

5 Evaluation

This section evaluates the INKS model on the NAS benchmark, a simple
stencil code and the 6D Vlasov-Poisson problem. We have implemented the
algorithm of the following programs using INKS<sub>pia</sub>.

- 4 sequential NAS kernels (IS, FT, EP and MG), C++ version [11] as reference;

- finite difference 3D heat resolution (7-point stencil) ([15] as reference);

- 6D Vlasov-Poisson, using Fortran/MPI Selalib [16] as reference.

For the 3d heat equation solver, two strategies were implemented: one uses dou-
ble buffering (Heat/Buf) and the other implements a cache oblivious strategy
(Heat/Obl). To evaluate the INKS programming model on Vasov-Poisson 6D,
we have conduct three experiments. The first one is based on optimization
choices written in plain C++ and cover the whole simulation. Then, we eval-
uate INKS<sub>o/XMP</sub> and INKS<sub>o/loop</sub> on Vlasov-Poisson separately as they target
different optimizations and are not usable together currently. A first experi-
ment focuses on the sequential aspects with the intra-node optimization of the
v<sub>1</sub> advection using either INKS<sub>o/loop</sub> or INKS<sub>o/C++</sub>. A second experiment fo-
cuses on the parallel aspects with the charge density computation, the Poisson
solver and a halo exchange optimized either with C/XMP or with INKS<sub>o/XMP</sub>.
<table><caption>Table 1: Summary of the INKS<sub>o/C++</sub>, INKS<sub>o/loop</sub> and INKS<sub>o/XMP</sub> optimization choices implemented based on six INKS<sub>pia</sub> algorithm: the 3D heat solver, 4 NAS kernels (EP, FT, MG, IS) and the Vlasov-Poisson 6D solver.</caption><thead><tr><th>Algorithm</th><th>INKS<sub>o/C++</sub></th><th>INKS<sub>o/loop</sub></th><th>INKS<sub>o/XMP</sub></th></tr></thead><tbody><tr><td>Heat</td><td>Buf<br>Obl</td><td colspan="2">3 versions: Same as reference, 2D and 3D cache blocking</td></tr><tr><td>NAS</td><td>EP<br>FT<br>MG<br>IS</td><td colspan="2">Red</td></tr><tr><td>VP6D</td><td>Initialization<br>Advections<br>Reduction<br>Poisson solver</td><td colspan="2">adv_v1</td></tr></tbody></table>

<table><caption>Table 2: Execution time of the INKS<sub>o/C++</sub>, INKS<sub>o/loop</sub> and INKS<sub>o/XMP</sub> implementations of the sequential NAS benchmark, class B - Time/iteration of the 3D heat equation (7point stencil), size (1024³) - Time/iteration of the INKS<sub>o/C++</sub> and Fortran implementation of the sequential 6D Vlasov-Poisson, size (32⁶). Median and standard deviation of 10 executions - GNU Complexity score of the implementation.</caption><thead><tr><th rowspan="2">Benchmark</th><th colspan="3">Execution time (second)</th><th colspan="2">Complexity</th></tr><tr><th>Reference</th><th>InKS<sub>o/C++</sub></th><th>Rel. dev.</th><th>Ref.</th><th>InKS</th></tr></thead><tbody><tr><td>NAS/FT</td><td>49.16 (±0.16)</td><td>44.00 (±0.05)</td><td>11.66%</td><td>6</td><td>5</td></tr><tr><td>NAS/IS</td><td>1.91 (±0.00)</td><td>1.90 (±0.01)</td><td>0.53%</td><td>55</td><td>52</td></tr><tr><td>NAS/MG</td><td>4.61 (±0.02)</td><td>4.29 (±0.01)</td><td>7.46%</td><td>20</td><td>12</td></tr><tr><td>NAS/EP</td><td>53.52 (±0.02)</td><td>55.35 (±0.37)</td><td>-3.30%</td><td>19</td><td>19</td></tr><tr><td>Heat/Buf</td><td>2.76 (±2.78)</td><td>2.93 (±3.53)</td><td>-6.22%</td><td>5</td><td>3</td></tr><tr><td>Heat/Obl</td><td>1.20 (±1.16)</td><td>1.18 (±0.96)</td><td>1.81%</td><td>22</td><td>13</td></tr><tr><td>VP6D</td><td>25.92 (±0.21)</td><td>21.92 (±0.16)</td><td>18.28%</td><td>N/A</td><td>N/A</td></tr><tr><td>Benchmark</td><td>Reference (C++)</td><td>InKS<sub>o/loop</sub></td><td>Rel. dev.</td><td>Reference (Fortran)</td><td></td></tr><tr><td>NAS/IS</td><td>1.91 (±0.00)</td><td>1.91 (±0.01)</td><td>0.00%</td><td>N/A</td><td></td></tr><tr><td>NAS/MG</td><td>4.61 (±0.02)</td><td>4.32 (±0.02)</td><td>6.84%</td><td>N/A</td><td></td></tr><tr><td>VP6D</td><td>4.18 (±0.05)</td><td>4.23 (±0.06)</td><td>-1.10%</td><td>3.63 (±0.09)</td><td></td></tr><tr><td>Benchmark</td><td>Reference</td><td>InKS<sub>o/XMP</sub></td><td>Efficiency (# cores)</td><td></td><td></td></tr><tr><td>NAS/EP</td><td>53.52 (±0.02)</td><td>14.00 (±0.13)</td><td>70.5% (4 cores)</td><td></td><td></td></tr></tbody></table>

The table 1 summarizes the optimization choices we have implemented. All codes are compiled with Intel 18 compiler (-O3 -xHost), Intel MPI 2018 and executed on the Irene cluster (TGCC, France) equipped with 192 GB RAM, two Skylake 8168 CPUs per node and a EDR InfiniBand interconnect. The table 2 summarizes the optimization choices we have implemented.

The NAS CG kernel relies on indirections not expressible in the polyhedron model of INKS<sub>pia</sub>. Its implementation would thus have to rely on a large C++ kernel whose optimization would be mixed with the algorithm. For the same reason, the NAS FT kernels was only partially implemented. INKS<sub>pia</sub> can however be used to express all other NAS kernels, the 3D heat equation solver as well as the 6D Vlasov-Poisson algorithm. Even if not as expressive as C or Fortran, INKS<sub>pia</sub>, through the polyhedron model, handles static controls programs. This covers the needs of a wide range of simulation domains and offers abstractions close to the execution machine rather than from a specific simulation domain. More specifically, it supports programs expressible as Parameterized Task Graphs [7]: a directed acyclic graph of tasks (here kernels)
in a compact representation whose dependence and scheduling are parame-
terized by integers fixed at the INKS entry point execution. Among others,
it can express computations such as FFTs or stencils with input coordinates
unknown at compile-time, as in 6D Vlasov-Poisson.

INKS separates the specification of algorithm and optimization in distinct files. Multiple optimization strategies can be implemented for a single algorithm, as shown for the 3D heat equation where each relies on a specific memory layout and scheduling. Similarly, based on the INKS<sub>pia</sub> algorithm implemented for the IS, EP and MG NAS kernels, we have developed multiple optimization choices based on INKS<sub>o/C++</sub> and either INKS<sub>o/loop</sub> or INKS<sub>o/XMP</sub>. For more complex cases such as the 6D Vlasov-Poisson problem, we were able to develop four optimization choices based on one INKS<sub>pia</sub> algorithm: (1) INKS<sub>o/XMP</sub>, (2) C with XMP, (3) INKS<sub>o/loop</sub> and (4) INKS<sub>o/C++</sub>. This proves, to some extent, that the separation of concerns is respected.

Finding the right metric to evaluate the easiness of writing a code is a
difficult question. As illustrated in Listings 1 and 5 however, algorithm expres-
sion in INKS<sub>pia</sub> is close to the most naive C implementation where loops are
replaced by INKS validity domains with no worry for optimization. The poly-
hedron model used to represent the INKS<sub>pia</sub> code is compatible with a subset
of the C language. Therefore, we could have choose the C as the algorithm
language. However, the optimization process would have suffer from such a
choice. Indeed, while INKS<sub>pia</sub> is designed to remove optimizations through the
use of DSA and fine granularity kernels, C enables its users to reuse memory
and specify a total scheduling. Although it is possible to analyze the code and
to retrieve the DSA form and the partial order of the code using the polyhe-
dron model, it will be complicated for optimization specialists to find which
loops can be broke or which arrays can be expended from a C code.

Concerning the specification of optimization choices, it is close to their
expression in C++. Table 2 compares the GNU complexity score of INKS op-
timizations to the reference code. INKS scores are slightly better, which indi-
cates that our language is not more complex than C++. The difference comes
from the extraction of the computational kernels, placed in the algorithm,
which hides parts of the complexity. In addition, the use of INKS<sub>o/C++</sub> to
write optimizations let optimization specialists reuse their preexisting knowl-
edge of this language. Similarly writing the INKS<sub>o/C++</sub> version of optimization
choices for the INKS<sub>pia</sub> version of Vlasov-Poisson 6D is close to the expression
to the same optimization in Fortran, in the reference version. These consider-
ations should not hide the fact that some information has to be specified both
in the INKS<sub>pia</sub> and INKS<sub>o/C++</sub> files with this approach leading to more code
overall.

On contrary, INKS<sub>o/XMP</sub> and INKS<sub>o/loop</sub> enables the developer to spec-
ify optimization choices only while algorithmic information is extracted from
INKS<sub>pia</sub> code. This is illustrated by Listing 2 presenting the INKS<sub>o/XMP</sub> 6D
domain decomposition and the XMP result. Both are equivalent, but the
INKS<sub>o/XMP</sub> expects only optimization choices parameters. Hence, one can test
another memory layout, such as a different dimension ordering, by changing
only a few parameters, while multiple directives must be modified in XMP. Moreover, using the domain decomposition directive offered by INKS<sub>o/XMP</sub> and XMP code, one can derive a simple parallel code from a INKS<sub>pia</sub> algorithm, as we did with the 3D Heat solver. Note that since INKS<sub>o/XMP</sub> is a wrapper for XMP using INKS<sub>pia</sub> to retrieve some information, it is usable in any codes that can be expressed using INKS<sub>pia</sub> and optimized with XMP. Similarly, with INKS<sub>o/loop</sub> (Listing 5), developers can easily test different optimization choices that would be tedious in plain C++. This is what we did with the 3D Heat solver based on the double buffering technique. As shown on Figure 1, using INKS<sub>o/loop</sub>, we have implemented 3 versions of the loops: the same as reference, one based on a 2D cache-blocking and a last one using a 3D cache-blocking by modifying almost nothing in the INKS<sub>o/loop</sub> code. Since INKS<sub>o/XMP</sub> and INKS<sub>o/loop</sub> are respectively usable with C and C++, INKS does not restrict the expressible optimization choices: one can still implement optimizations not handled by our DSLs in C/C++. Moreover, operations such as halo size computation or vectorization capabilities detection are automatized using the INKS<sub>pia</sub> code. In summary, the approach enables optimization specialists to focus on their specialty which make the development easier.

Regarding performance, the INKS approach makes it possible to express
optimizations that do not change the algorithm. Optimizations of the four
NAS parallel benchmarks and 3D heat equation solver in INKS were trivial
to implement and their performance match or improve upon the reference as
presented in Table 2. Investigation have shown that Intel ICC 18 does not vec-
torize properly the reference versions of MG and FT. The use of the Intel ivdep
directive as done on the INKS versions leads to slightly better performance.

For the full Vlasov-Poisson 6D problem, INKS<sub>pia</sub> and C++ enable us to implement the same complex optimization strategies written in the reference version. Therefore, our implementation match the reference in terms of performance, as shown on Table 2. For the v<sub>1</sub> advection, both the INKS<sub>o/C++</sub> and INKS<sub>o/loop</sub> optimizations of the INKS code achieve performance similar to the reference as shown in Table 2. For the parallel aspects, the INKS<sub>o/XMP</sub> optimization offers performance similar to XMP as shown on Figure 4. MPI is faster on all cases compared to both XMP and INKS<sub>o/XMP</sub>. At the moment, it seems that XMP does not optimize local copies which slows down the Poisson solver. Besides, some XMP directives are based on MPI RMA which make the comparison with MPI Send/Receive complex. Still, MPI is much harder to program: more than 350 lines of MPI and Fortran are required to handle domain decomposition, remapping for FFT and halo exchange in Selalib while 50 lines in XMP and 15 in INKS<sub>o/XMP</sub>.

However, INKS<sub>pia</sub>, INKS<sub>o/loop</sub> and INKS<sub>o/XMP</sub> have limitations. As mentioned earlier, INKS<sub>pia</sub> can not be used to express non static control program, such as the CG or the full FT NAS kernel. It also makes more complex the optimization of arrays passed in parameters (as input or output). INKS<sub>o/loop</sub> offers very limited option, making it unusable with complex loops such as the one in the 3D heat solver optimized using the cache-oblivious strategy. Similarly, loop nests calling multiple INKS<sub>pia</sub> computational kernels can not be expressed
<figure><img src="image_4.png" /><figcaption>Fig. 4: Weak and strong scaling for 3 parts of the Vlasov-Poisson solver up to 64 nodes (1 process/node) on a 32⁶ grid divided among processes (strong scaling) or 16⁶ grid per process (weak scaling). Median of 10 executions.</figcaption></figure>

using INKS<sub>o/loop</sub>. Moreover, it can not respect the constraints imposed by a specific memory mapping and offers only a few optimization. Adding new optimization strategies would require to add new keywords, such as unrolling, loop fusion, etc. INKS<sub>o/XMP</sub> offers memory allocations but no controls on the access to this memory, letting optimization specialists in charge of the good use of this memory. Besides, it is usable only with XMP.

Although we want to address some issues in the INKS<sub>pia</sub> language, INKS<sub>o/loop</sub> and INKS<sub>o/XMP</sub> were preliminary tests before a real INKS<sub>pso</sub> implementation. The goal was to propose a way to express optimizations and retrieving the algorithmic information in INKS<sub>pia</sub> code. These two DSLs enables us to highlight a set of guidelines for the full definition of the INKS<sub>pso</sub> language. First, this definition must be based on the concepts that INKS<sub>pso</sub> must express, i.e. the optimizations related to memory (allocations and layouts) and to computations scheduling. With INKS<sub>o/XMP</sub> we tried some tests on the memory aspects, providing a directive to allocate a logical array and to reorder its dimensions. INKS<sub>o/loop</sub> focus on computations scheduling by adding strong constraint on the order of the computations. Secondly, in order to express these optimization concepts only, it must be bound to its algorithmic counterpart, express in INKS<sub>pia</sub>. In INKS<sub>o/XMP</sub>, the domain decompistion directive makes reference directly to logical arrays described in INKS<sub>pia</sub>, making possible the computation of the size of each dimension and the halo size depending of the data being accessed. In INKS<sub>o/loop</sub>, the *Order* keyword refers to the structuring variable of a computational kernel defined in INKS<sub>pia</sub>. Finally, we must work on the mean to express these optimization concepts with the most generality. For this part, we think that working with existing tools is probably mandatory for most complex strategy, such as XMP for the domain decomposition. Using another existing tool for computations scheduling, making it more general. Similarly, a declarative language such as INKS<sub>o/loop</sub> may not be the best strategy: to handle the interaction between memory and computation, an imperative language could be easier.
# 6 Conclusion and future works

In this paper, we have presented the INKS programming model to separate the algorithm and the optimization choices and its implementation supporting two DSLs: INKS<sub>o/loop</sub> for loop optimizations and INKS<sub>o/XMP</sub> for domain decomposition. We have evaluated INKS on synthetic benchmarks and on the Vlasov-Poisson solving. We have demonstrated its generality and advantages in terms of separation of concerns to improve maintainability and portability while offering performance on par with existing approaches.

While this paper demonstrates the interest of the INKS model, it still requires some work to further develop it. We will improve the optimization DSLs; base INKS<sub>o/loop</sub> on existing tools and ensure interactions with INKS<sub>o/XMP</sub>. This will be done within the scope of the development of a complete INKS<sub>pso</sub> DSL enabling its users to manage memory placement and computations scheduling. This planed DSL will truly separates algorithm from optimization choices in static programs while its compiler could offer static analysis of the program to ensure correctness. We also want to target different architectures to demonstrate the portability gains of the INKS model.

# References

1. Augonnet C, Thibault S, Namyst R, Wacrenier PA (2011) StarPU: A unified platform for task scheduling on heterogeneous multicore architectures. Concurr Comput : Pract Exper 23(2), DOI 10.1002/cpe.1631

2. Aumage O, Bigot J, Ejjaouani K, Mehrenberger M (2017) InKS, a programming model to decouple performance from semantics in simulation codes. Tech. rep., Inria

3. Bailey DH, Barszcz E, Barton JT, Browning DS, Carter RL, Dagum L, Fatoohi RA, Frederickson PO, Lasinski TA, Schreiber RS, et al. (1991) The NAS parallel benchmarks. Journal of Supercomputing Applications

4. Carter Edwards H, Trott CR, Sunderland D (2014) Kokkos. J Parallel Distrib Comput 74(12), DOI 10.1016/j.jpdc.2014.07.003

5. Chandra R, Dagum L, Kohr D, Maydan D, McDonald J, Menon R (2001) Parallel Programming in OpenMP. Morgan Kaufmann

6. Christen M, Schenk O, Burkhart H (2011) PATUS: A Code Generation and Autotuning Framework for Parallel Iterative Stencil Computations on Modern Microarchitectures. In: Parallel & Distributed Processing Symposium (IPDPS) 2011, IEEE, DOI 10.1109/ipdps.2011.70

7. Cosnard M, Jeannot E (1999) Compact dag representation and its dynamic scheduling. Journal of Parallel and Distributed Computing 58(3), DOI 10.1006/jpdc.1999.1566

8. Danelutto M, García J, Miguel Sanchez L, Sotomayor R, Torquati M (2016) Introducing parallelism by using repara c++11 attributes. pp 354–358, DOI 10.1109/PDP.2016.115
9. El-Ghazawi T, Carlson W, Sterling T, Yelick K (2005) UPC: Distributed Shared Memory Programming. Wiley

10. Feautrier P, Lengauer C (2011) Polyhedron Model, Springer. DOI 10.1007/978-0-387-09766-4_502

11. Griebler D, Loff J, Mencagli G, Danelutto M, Fernandes LG (2018) Efficient NAS benchmark kernels with c++ parallel programming. In: 2018 26th Euromicro International Conference on Parallel, Distributed and Network-based Processing (PDP), DOI 10.1109/PDP2018.2018.00120

12. Hoque R, Herault T, Bosilca G, Dongarra J (2017) Dynamic task discovery in PaRSEC: A data-flow task-based runtime. In: 8th Workshop on Latest Advances in Scalable Algorithms for Large-Scale Systems, ACM

13. Höhnerbach M, Ismail AE, Bientinesi P (2016) The vectorization of the Tersoff multi-body potential: An exercise in performance portability. In: International Conference for High Performance Computing, Networking, Storage and Analysis, IEEE

14. Isoard A (2016) Extending polyhedral techniques towards parallel specifications and approximations. PhD thesis, École doctorale en Informatique et Mathématiques de Lyon

15. Kamil S (2012) StencilProbe: A microbenchmark for stencil applications. Accessed: 2017-08-25

16. Kormann K, Reuter K, Rampp M, Sonnendrücker E (2016) Massively parallel semi-lagrangian solution of the 6d vlasov-poisson problem.

17. Lee J, Sato M (2010) Implementation and performance evaluation of Xcal-ableMP: A parallel programming language for distributed memory systems. In: International Conference on Parallel Processing Workshops

18. Mehrenberger M, Steiner C, Marradi L, Crouseilles N, Sonnendrucker E, Afeyan B (2013) Vlasov on GPU (VOG project). ESAIM: Proc 43, DOI 10.1051/proc/201343003

19. Steuwer M, Remmelg T, Dubach C (2017) LIFT: A functional data-parallel IR for high-performance gpu code generation. In: 2017 IEEE/ACM International Symposium on Code Generation and Optimization (CGO)

20. Tang Y, Chowdhury RA, Kuszmaul BC, Luk CK, Leiserson CE (2011) The Pochoir stencil compiler. In: 23rd Symposium on Parallelism in Algorithms and Architectures, ACM, SPAA '11, DOI 10.1145/1989493.1989508

21. Tanno H, Iwasaki H (2009) Parallel skeletons for variable-length lists in SkeTo skeleton library. In: Proceedings of the 15th International Euro-Par Conference on Parallel Processing, Springer, Euro-Par '09, DOI 10.1007/978-3-642-03869-3_63

22. Verdoolaege S (2010) isl: An integer set library for the polyhedral model. In: Fukuda K, Hoeven Jvd, Joswig M, Takayama N (eds) Mathematical Software – ICMS 2010, Springer