Sylvain Contassot-Vivier
email: sylvain.contassotvivier@loria.fr

Jean-François Couchot
email: jean-francois.couchot@univ-fcomte.fr

Pierre-Cyrille Héam
email: pierre-cyrille.heam@femto-st.fr

Gray Codes Generation Algorithm and Theoretical Evaluation of Random Walks in N-Cubes

Keywords: Gray Codes, Hamiltonian cycles, N-cube, mixing time, PRNG

In previous works, some of the authors have proposed a canonical form of Gray Codes (GCs) in Ncubes (hypercubes of dimension N). This form allowed them to draw an algorithm that theoretically provides exactly all the GCs for a given dimension N. In another work, we first have shown that any of these GC can be used to build the transition function of a Pseudorandom Number Generator (PRNG). Also, we have found a theoretical quadratic upper bound of the mixing time, i.e., the number of iterations that are required to provide a PRNG whose output is uniform. This article, extends these two previous works both practically and theoretically. On the one hand, another algorithm for generating GCs is proposed that provides an efficient generation of subsets of the entire set of GCs related to a given dimension N. This offers a large choice of GC to be used in the construction of Choatic Iterations based PRNGs (CI-PRNGs), leading to a large class of possible PRNGs. On the other hand, the mixing time has been theoretically shown to be in N log(N), which was anticipated in the previous article, but not proven.

Introduction

Random walk in the N-cube is the process of randomly crossing edges in the hypercube of N dimensions. This theoretical subject has received a lot of attention in recent years [START_REF] Lovász | Random walks on graphs[END_REF][START_REF] Diaconis | Asymptotic analysis of a random walk on a hypercube with many dimensions[END_REF][START_REF] Gregory | Random walk: a modern introduction[END_REF][START_REF] Scoppola | Exact solution for a class of random walk on the hypercube[END_REF][START_REF] David | Markov chains and mixing times[END_REF]. The article [START_REF] Contassot-Vivier | Random Walk in a N-Cube Without Hamiltonian Cycle to Chaotic Pseudorandom Number Generation: Theoretical and Practical Considerations[END_REF] is situated in this thematic field. Indeed, it focuses on a Random Walk in a N-cube where a balanced Hamiltonian cycle has been removed. It is worth recalling that a Hamiltonian cycle in the N-cube is also a Gray Code (GC). To compute such kind of cycle, an algebraic solution of an undeterministic approach [START_REF] Suparta | Totally balanced and exponentially balanced gray codes[END_REF] is presented. Solving the induced system of linear equations is an easy and practicable task. However, the number of solutions found by this method is infinitely smaller than what generally exists if the approach [START_REF] Suparta | Totally balanced and exponentially balanced gray codes[END_REF] was not used.

In another work [START_REF] Contassot | Canonical Form of Gray Codes in N-cubes[END_REF], we have proposed a canonical form of GCs in N-cubes. This form allows to draw an algorithm, based on Wild's scheme [START_REF] Wild | Generating all cycles, chordless cycles, and hamiltonian cycles with the principle of exclusion[END_REF], that theoretically provides exactly all the GCs for a given dimension N. However, due to the exponential growth of the number of GCs with the dimension of the N-cube, the number of produced GCs must be specified as a boundary argument. Moreover, we have observed that the time to generate a first GC may be quite large when the dimension of the N-cube increases.

Once a GC is produced, it can be used to build a random walk in the N-cube based post-processing of the output of any classical PRNG [START_REF] Couchot | Pseudorandom number generators with balanced gray codes[END_REF][START_REF] Contassot-Vivier | Random Walk in a N-Cube Without Hamiltonian Cycle to Chaotic Pseudorandom Number Generation: Theoretical and Practical Considerations[END_REF][START_REF] Bakiri | CIPRNG: A VLSI family of chaotic iterations post-processings for $\mathbb {F}_{2}$ -linear pseudorandom number generation based on zynq mpsoc[END_REF]. This post-processing has many advantages: compared to the original PRNG which does not have a chaotic behavior according to Devaney, it first adds this property [START_REF] Bahi | On the link between strongly connected iteration graphs and chaotic boolean discrete-time dynamical systems[END_REF]. It moreover preserves the cryptographical security of the embedded PRNG [START_REF] Guyeux | Efficient and cryptographically secure generation of chaotic pseudorandom numbers on GPU[END_REF], which may be required in a cryptographic context. Finally, it largely improves the statistical properties of the on board PRNG [START_REF] Bakiri | FPGA implementation of f2-linear pseudorandom number generators based on zynq mpsoc: A chaotic iterations post processing case study[END_REF] leading to a PRNG that successfully passes not only the classical statistical tests like [START_REF] Marsaglia | Diehard: a battery of tests of randomness[END_REF], but also the most difficult ones, namely the TestU01 [START_REF] Pierre | Testu01: A c library for empirical testing of random number generators[END_REF]. In particular, it has been already presented in ([2] Chap. 5, p. 84+) that such post-processing provides a new PRNG at minimal cost which is 2.5 times faster than the PRNGs which pass the whole TestU01 [START_REF] Pierre | Testu01: A c library for empirical testing of random number generators[END_REF]. Moreover, it is worth noticing that most of the classical PRNGs [START_REF] Ecuyer | Handbook of Computational Statistics, chapter Random Number Generation[END_REF][START_REF] Pierre | History of uniform random number generation[END_REF] do not pass the TestU01.

In such a context, the mixing time [START_REF] David | Markov chains and mixing times[END_REF] corresponds to the number of steps that are required to reach any vertex with the same probability according to a given approximation tolerance ε. In other words, it is the time until the set of vertices is ε-close to the uniform distribution. Such a theoretical model provides information on the statistical quality of the obtained PRNG. The lower the mixing time, the more uniform the PRNG output.

We first have shown that any of the generated GC can be used to build the transition function of a Pseudorandom Number Generator (PRNG) based on random walk in the N-cube. Then, we had found a theoretical quadratic upper bound of the mixing time.

This article is twofold as it extends these two previous works both practically and theoretically. In the first part (Section 2), we present another algorithm that can generate many GCs in a given N-cube more efficiently than our previous algorithm. Moreover, filtering rules can be applied during the generation in order to choose codes that are especially adapted to the use in PRNG based on random walks. This leads to a very large set of possible PRNGs. In the second part (Section 4), the mixing time has been theoretically shown to be in N. log 2 (N), which was anticipated, but not yet proven. Then, experimental results are given and discussed in Section 5.

Gray Codes Generation Algorithms

In this part, we propose another algorithm that efficiently generates GCs in a given N-cube, starting from a list of known GCs. The initial list may be reduced to a singleton.

Inverting Algorithm

Many approaches to generate GCs in a N-cube start from one or more known codes in dimensions smaller than N [START_REF] Robinson | Counting sequences[END_REF][START_REF] Girish | Balanced gray codes[END_REF][START_REF] Suparta | Totally balanced and exponentially balanced gray codes[END_REF][START_REF] Suparta | A construction of gray codes inducing complete graphs[END_REF][START_REF] Bykov | On locally balanced gray codes[END_REF]. However, such methods present either strong limitations over the type of cycles that can be generated, or large amount of computations that lead to non efficient generation algorithms.

The approach we propose in this paper is quite different as it builds new Gray codes that have the same dimension as the initial known code. The principle consists in transforming one GC into another non-isomorphic GC by an adequate operation. In the sequel, we call neighbors, two GCs that are directly linked by such an operation. Hence, the key idea of the approach is to explore the space of possible GCs by generating the neighbors of the initial code and then re-applying the same strategy to the obtained codes. Therefore, in a sense, we perform a transitive closure according to the considered code transformation as depicted in Figure 1.

... ... initial GC

Figure 1: Transitive exploration of the space of Gray Codes (GCs) from the initial black code. Each dot represents one GC. The color gets lighter at each additional transform.

GC Transformation Principle

The core issue in the proposed approach above is to find a code transformation that generates another valid Gray code that is not isomorphic to the initial code. In this context, we propose what we call the inversion transformation. It consists in finding a particular sub-sequence of the initial code that can be inverted to produce another valid Gray code. If C is a GC in dimension N whose sequence of vertices is (C0, ..., C 2 N -1 ), the problem is to find a sub-sequence (Ci, ..., Cj) in C such that the code (C0, ..., Ci-1, Cj, Cj-1, ..., Ci+1, Ci, Cj+1, ..., C 2 N -1 ) is also a valid GC.

We have identified conditions over the four vertices Ci-1, Ci, Cj and Cj+1 that fulfill that constraint:

1. the four vertices must be on the same face of the N-cube, i.e., the differences between the four vertices are contained in only two dimensions of the N-cube

2.

Ci and Cj must be diagonally opposed inside the face and, as a consequence, Ci-1 and Cj+1 are diagonally opposed too (this is equivalent to say that transitions (Ci-1, Ci) and (Cj, Cj+1) go in the same direction along their dimension)

Figure 2 illustrates the configuration providing a valid sub-sequence inversion. Those conditions comes from the fact that in order to obtain a valid GC, the exchanges of Ci and Cj in the sequence must satisfies the condition of Gray codes. This implies that the number of cycle Hamiltonian Face of the N-cube

C i-1 C j+1 C i C j C i-1 C i C j C j+1
Sub-seq inversion Figure 2: Inversion of sub-sequence (C i , ..., C j ) inside a face of the N-cube (dotted square).

differences between Ci-1 and Cj must be only one, and this is the same for Ci and Cj+1. Therefore, as Ci and Cj are at a Hamming distance of 1 from Ci-1 and they are different vertices (by construction of the initial GC), they must be exactly at a Hamming distance of two to each other. So, they are diagonally opposed on a same face of the N-cube and Ci-1 is also on that face (at one of the other two corners). The same reasoning applies for Ci, Cj and Cj+1, leading to the conclusion that Cj+1 must be at the last free corner of the same face. Therefore, it is diagonally opposed to Ci-1, leading to the configuration in Figure 2. By construction, it is sure that the obtained cycle is a valid Hamiltonian cycle (Gray code in the N-cube). Indeed, as the inversion operation is equivalent to change only some transitions of the input code, the output code contains exactly the same set of vertices as the input one. Therefore, the generated code satisfies the property of Hamiltonian cycle as it contains the total number of vertices of the considered graph. Moreover, as the Gray code property is independent of the cycle direction, the inverted sub-sequence (Cj, ..., Ci) is still a valid part of a Gray code. Also, as explained before, the two conditions enumerated above ensure that both transitions (Ci-1, Cj) and (Ci, Cj+1), connecting the inverted sub-sequence and the unchanged part of the initial code, take place along one dimension of the N-cube, thus satisfying the Gray code property. Therefore the generated code fulfills all the constraints of Gray codes in the N-cube. An example of code generation in dimension 4 is given in Figure 3. [START_REF] Contassot | Canonical Form of Gray Codes in N-cubes[END_REF][START_REF] Contassot-Vivier | Random Walk in a N-Cube Without Hamiltonian Cycle to Chaotic Pseudorandom Number Generation: Theoretical and Practical Considerations[END_REF][START_REF] Guyeux | Efficient and cryptographically secure generation of chaotic pseudorandom numbers on GPU[END_REF][START_REF] Gregory | Random walk: a modern introduction[END_REF] in the 4-cube.

A first consequence of the sub-sequence inversion is that the codes generated from one initial code have different dimension frequencies than their source. We remind the reader that the dimension frequencies of a GC are the numbers of times each dimension is traversed by the code. The mentioned consequence comes from the fact that when inverting the sequence (Ci, ..., Cj) in the initial code, the transitions (Ci-1, Ci) and (Cj, Cj+1) that are along the same dimension (opposite edges of the face), are respectively replaced by the transitions (Ci-1, Cj) and (Ci, Cj+1), that are along the other dimension of the face. This fact, that can be observed in Figure 2, implies that between the two cycles, the frequency of one dimension decreases by two while the frequency of another one increases by two. This is important when wanting to guide the generation towards cycles with specific properties, like the global balance of frequencies.

Another consequence is that the sub-sequence (Ci, ..., Cj) must contain at least five vertices to go from Ci to Cj by traversing one dimension at a time while avoiding the vertices Ci-1 and Cj+1. This reasoning applies also to the sub-sequence (Cj+1, ..., Ci-1) for the same reason. Such information can be exploited to optimize the generation algorithm.

Generation Algorithm from a Known GC

Algorithm 1 presents the sketch of the GC generation algorithm by sub-sequence inversion.

For each position in the initial sequence C, the algorithm takes the transition from this vertex as a reference and searches for any parallel transition in the sequence. By parallel transition, we mean a transition along the same dimension and in the same direction, whose vertices form a face of the N-cube together with the vertices of the reference transition.

In this algorithm, function dimOfTransition(x,y) returns the dimension index corresponding to the transition formed by the two neighbors x and y in the N-cube (complexity O (1)). Also, function bitInversion(x,i) inverses the bit i of vector x (O (1)). Function copyAndInvertSubSeq(C,a,b) returns a copy of cycle C in which the sub-sequence between indices a and b is inverted (O (|C|)). Finally, function exclusiveInsert(C,L) inserts cycle C in list L only if C is not yet in L, and it returns the Boolean true when the insertion is performed and false otherwise (O (log 2 (|L|)) thanks to the total order induced by the canonical form).

As depicted in Figure 4a, there are N -1 possible parallel transitions for a given transition in the N-cube. The dotted lines exhibit the dimensions that must be parsed to get the other neighbors of the reference transition source, which are not the reference transition destination. In this particular example, the algorithm would search for sub-sequences of the form (1, ..., 2, 3), (1, ..., 4, 5) and (1, ..., 8, 9). Each time one of these sub-sequences is found, a new GC can be produced. It is worth noticing that several sub-sequences may be valid for one reference transition, as shown in Figure 4b. In that example, the red transitions are valid parallel transitions of the black transition. Therefore, two distinct cycles can be generated by inverting either (1, ..., 2) or (1, ..., 4).

Algorithm 1: Function generateFrom(C). To be exhaustive, and to get all the possible inversions, the algorithm applies the search for parallel transitions successively to all the transitions in the initial sequence. Also, the constraint over the minimal length of the sub-sequence to be inverted and of the unchanged part of the cycle, mentioned in the previous section, has been taken into account in the inner part of our algorithm, when searching the considered parallel transition in the cycle (boundaries of variable k). Although this does not change the overall complexity of the algorithm, it saves computations and slightly reduces the execution time. For clarity sake, the algorithm is given in its simplest form, but additional features could be added, such as a filtering over the generated cycles according to some particular properties, or a more subtle selection and guiding of the sub-sequences searches that would avoid unnecessary computations when looking for specific GCs. Also, sorting the list lst according to the canonical forms of the generated cycles is recommended to perform faster exclusions or insertions.

Finally, we obtain an algorithm whose complexity to generate all the cycles that can be derived from one GC of length V is O V 2 .log 2 (V ) .This corresponds to the complexity O N.4 N for the corresponding N-cube of dimension N . Also, the memory complexity of the algorithm is O (c.V ), where c is the number of generated cycles (c = |lst|). This algorithm alone is not sufficient to generate as many cycles as wanted. In fact, it is trivial to deduce a coarse upper bound of the number of possible generations from one GC of length V , which is V.(log 2 (V ) -1). Thus, as explained in Section 2.1, we have designed a transitive closure algorithm that uses GC generation algorithm.

Transitive Closure Algorithm

Our transitive closure algorithm is an iterative algorithm that begins with a given set of known GCs (at least one). At the first step, each GC in the initial set is used as input to the GC generation algorithm described before. Then, the list of input GCs for the next step is rebuilt by including in it only the new cycles that have not yet been used as input GC. This clarification is important because one GC may be generated several times from different sources. So, it is necessary to maintain a list of all the GCs that have already been treated since the beginning of the process, in order to avoid duplicate work. At the next step, the GCs in the built list are used as input to the GC generation algorithm, and so on until no more GC is generated or a specific number of GCs has been reached. The overall process is detailed in Algorithm 2, in which the function generateFrom() is the one detailed in Algorithm 1.

Algorithm 2: Transitive closure generation algorithm. Concerning the stopping criteria of the iterative process, it must include the detection that no new code has been produced at the last iteration. Also, as soon as the dimension of the N-cube is greater than 5, the number of potential generations becomes very large and the stopping criteria must include a limitation of the number of generations. As for the previous algorithm, the lists lst and nl must be sorted according to the canonical forms in order to obtain fast searches and insertions.

The complexity of that algorithm is deduced by analyzing the nested loops. Although the number of iterations of the global while loop is variable, we can compute the total number of iterations of inner loops on i (line 9) and on j (line 11). This allows us to deduce an upper bound of the global complexity as the other actions performed in the while loop: the initialization of nl, the copy of nl into il and the stopping decision, take either constant or linear times (naive copy of il into nl, but optimization is possible) that are negligible compared to the total cost of the loop on i.

Concerning the loop on i, as the il list stores only the new cycles at each iteration of the closure process, we can deduce that each generated cycle is stored only once in this list, and thus is processed only once by the loop. Thus, the total number of iterations performed by this loop is equal to the total number of generated cycles (denoted by L). Inside this loop, the cost of generateFrom is N.4 N as seen in the previous sub-section. The number of cycles generated by generateFrom may vary from one cycle i to another. However, an amortized complexity analysis using α as the average number of generations per cycle allows us to deduce that the total number of iterations performed by the inner loop in j is α.L. This last loop is dominated by the exclusive insertions in the sorted list lst.

Here also, the size of that list increases during the whole process, from 0 to L elements. In order to compute the total cost of the successive calls to exclusiveInsert(), we define S as the sum of contributions of every call with a different size of list lst : S = L-1 i=0 log 2 (i) < L.log 2 (L). As there are α.L iterations instead of only L, this bound must be multiplied by α. So, the total cost of loop on j is O (α.L.log 2 (L)).

Finally, combining the costs of the L calls to generateFrom() and of the α.L iterations on j, we obtain a total cost of O L.(N.4 N + α.log 2 (L)) for the loop on i. This is clearly the dominant cost in the while loop, and thus the overall complexity of the transitive closure process. Also, the memory complexity is dominated by the list of generated cycles, leading to O (L.V ).

Finally, it is worth noticing that for a given initial GC, the final set of GCs obtained by the transitive closure based on sub-sequence inversions may not be equal to the entire set of all possible GCs in the corresponding N-cube, even when there is no restriction over the size of the final set. This comes from the fact that some GCs may have no valid parallel transition for any transition in the sequence, as show in Figure 5. It is worth noticing that the example on the left of this figure is a binary-reflected Gray code (BRGC) [START_REF] Marston | Explicit definition of the binary reflected gray codes[END_REF][START_REF] Bunder | On binary reflected gray codes and functions[END_REF]. By construction, BRGCs do not contain any valid parallel transition, leading to no possible sub-sequence inversion. Moreover, this property is not restricted to BRGCs, as depicted in the example on the right of the figure. Indeed, such GCs with no valid parallel transition correspond to singletons in the partition of the entire set of GCs, based on the sub-sequence inversion relation. So, since the partition is composed of several connected components, the number of produced GCs will depend on the connected component of the initial GC. However, we show in Section 5 that this phenomenon is very limited in practice as the probability to get such a singleton as initial GC is very low and decreases with the dimension of the N-cube. Moreover, in such case, it is not a problem to generate another initial GC by using any classical generation algorithm. Thus, this fact is not a significant obstacle to the generation of as many codes as wanted.

Using Gray Codes in Chaotic PRNGs

Full details on the use of Gray codes to build chaotic PRNGs can be found in [START_REF] Couchot | Pseudorandom number generators with balanced gray codes[END_REF][START_REF] Contassot-Vivier | Random Walk in a N-Cube Without Hamiltonian Cycle to Chaotic Pseudorandom Number Generation: Theoretical and Practical Considerations[END_REF]. For the sake of clarity, we briefly recall the main principles in this section. The chaotic PRNG can be seen as an iterative post-processing of the output of a classical PRNG of N bits (possibly with weak statistical properties) as shown in Algorithm 3. The required elements are the initial binary vector x 0 = (x 0 1 , x 0 2 , ..., x 0 N ), a number b of iterations and an iteration function f (x) = (f1(x), ..., fN (x)). This function is used in the chaotic mode F f (x, s), where s is the index of the only bit of x that is updated by F , so that we have:

x t+1 = F f (x t , st) = (x t
1 , x t 2 , ..., fs(x t ), ..., x t N ). The Gray code is used in the construction of function f : for a given dimension N , we start with the complete N-cube (Figure 6a), and we remove all the directed edges that belong to a given Gray code (Figure 6b) to obtain function f (Figure 6c). The removed edges are replaced by loop edges on every vertices.

Algorithm 3: Chaotic PRNG scheme.

1 Input: a function f , an iteration number b, an initial configuration x 0 (N bits) 2 Output: a configuration x (N bits) Applying the function F f on the pair (x, s) is equivalent to traverse the edge s from the node x in the modified N-cube presented in Figure 6c.

4 x ← x 0 5 for i=0 to b-1 do 6 s ← Random(N) 7 x ← F f (x, s) 8 endFor 9 return x

Reducing the Mixing Time Upper Bound of Boolean Functions in the N-Cube

As mentioned in the introduction and according to the PRNG presented in the previous section, the mixing time of a Boolean function in the N-cube provides important information about the minimal setting of the number of internal iterations of our chaotic PRNG (value b in Algorithm 3). In this section, a more accurate upper bound of this mixing time is proved by using Markov chains.

Theoretical Context

It has been proven in [START_REF] Couchot | Pseudorandom number generators with balanced gray codes[END_REF] that walking in an N-cube where an Hamiltonian cycle has been removed induces a Markov M chain whose stationary distribution π is the uniform one. Let us first formalize the relation between mixing time and uniform distribution.

First of all, let be given two distributions µ and π on the same set Ω, the Total Variation distance µ -π TV is defined for by:

µ -π TV = max A⊂Ω |µ(A) -π(A)|.
Let then M (x, •) be the distribution induced by the x-th row of the Markov matrix M . If the Markov chain induced by M has a stationary distribution π, then we define

d(t) = max x∈Ω M t (x, •) -π TV.
Finally, let ε be a positive number, the mixing time with respect to ε is given by

tmix(ε) = min{t | d(t) ≤ ε}.
It defines the smallest iteration number that is sufficient to obtain a deviation lesser than ε to the stationary distribution. It has been proven in ([18] Equation (4.36)) that tmix(ε) ≤ log 2 (ε -1 ) tmix( 1 4 ). We are then left to evaluate tmix( 1 4 ). For technical issues, rather than studying the random walk induced by the PRNG, we investigate the associated lazzy Markov chain, that is at each step, with probability 1/2 the chain moves to the next element and with probability 1/2 nothing is done. The mixing time of a lazzy Markov chain is about twice the mixing time of the initial chain ([5] p. 12).

Let N be a positive integer and Ω = {0, 1} N . For X, Y ∈ Ω, we define the subset H(X, Y ) of {1, . . . , N } of indexes where X and Y differ:

H(X, Y ) = {i | X(i) = Y (i)}. We define the distance d on Ω by d(X, Y ) = |H(X, Y )|,
that is the number of components where X and Y differ.

We also consider a Hamiltonian path on the hypercube Ω. More precisely, let h be a bijective function from Ω to Ω inducing a cyclic permutation on Ω and such that for any X, d(X, h(X)) = 1.

Since d(X, h(X)) = 1, we denote by f (X) the single index such that X = h(X).

We consider the Markov Chain of matrix P on Ω defined by:

• P (X, X) = 1 2 + 1 2N , • P (X, Y ) = 1 2N if d(X, Y ) = 1 and Y = h(X), • P (X, Y ) = 0 otherwise.
The Markov Chain P can be defined by the following random mapping representation: Xn+1 = r(Xn, U, V ) where U is a uniform random variable on 1, . . . , N , V a uniform random variable on {0, 1} and r(X, U, V ) = X if U = f (X) and, otherwise, r(X, U, V ) is obtained from X by fixing its U -th index to V .

A Useful Random Process

The mixing time of the lazzy Markov chain will be studied in Section 4.3 using a classical coupling approach. In this section we propose a purely technical result useful for proving some properties on this coupling.

We consider a Markov chain (Zi) on {0, . . . , N }, of matrix K defined by:

• P(Zi+1 = 0 | Zi = 0) = 1, • P(Zi+1 = 0 | Zi = 1) = 1 2N , • P(Zi+1 = 2 | Zi = 1) = 1 N , • P(Zi+1 = 1 | Zi = 1) = 1 -3 2N , • P(Zi+1 = N | Zi = N ) = 1 N and P(Zi+1 = N -1 | Zi = N ) = 1 -1 N , • for every k ∈ {2, . . . , N -1}, P(Zi+1 = k -1 | Zi = k) = k-1 N and P(Zi+1 = k + 1 | Zi = k) = 1 N and P(Zi+1 = k | Zi = k) = 1 -k N , Let R d = min{k | Z k = 0 and Z0 = d} if it exists and ρ d = E[R d ]. Proposition 1. One has ρN ≤ 3N (1 + e) + N log(N -2).
Proof. According to Theorem 1.3.5 in [START_REF] Norris | Markov Chains[END_REF], the ρi are the solution of the following system:

       ρN = 1 + 1 N ρN + N -1 N ρN-1 ρi = 1 + 1 N ρi+1 + i-1 N ρi-1 + (1 -i N )ρi for 2 ≤ i ≤ N -1 ρ1 = 1 + 1 2N ρ0 + 1 N ρ2 + (1 -3 2N )ρ1 ρ0 = 0
This system is equivalent to:

       ρN = N N -1 + ρN-1 i N ρi = 1 + 1 N ρi+1 + i-1 N ρi-1 for 2 ≤ i ≤ N -1 3 2N ρ1 = 1 + 1 N ρ2 ρ0 = 0 Set, for 0 ≤ i ≤ N -1, αi = ρi+1 -ρi. Using the system, one has αN-1 = N N -1 , α0 = ρ1. Moreover, for 2 ≤ i ≤ N -1, substituting ρi+1 by αi + ρi in the equation i N ρi = 1 + 1 N ρi+1 + i-1 N ρi-1, one obtains that i -1 N ρi = 1 + 1 N αi + i -1 N ρi-1. Therefore, for 2 ≤ i ≤ N -1, i -1 N αi-1 = 1 + 1 N αi,
or, equivalently,

αi-1 = N i -1 + 1 i -1 αi. (1)
Since αN-1 = N N -1 and using Equation ( 1), one has

αN-2 = N N -2 + N (N -1)(N -2) and αN-3 = N N -3 + N (N -2)(N -3) + N (N -1)(N -2)(N -3)
.

By a direct induction, one obtains, for 1 ≤ k ≤ N -1:

α N -k = N k j=1 j-1 =0 1 N -k + . (2)
Consequently

α N -k ≤ N k j=1 j-1 =0 1 N -k = N k j=1 1 (N -k) j . It follows that, for 1 ≤ k ≤ N -2, α N -k ≤ N 1 N -k 1 - 1 N -k k 1 -1 N -k = N 1 N -k -1 1 - 1 N -k k ≤ N N -k -1 . (3)
may be observed, the executions times are quite similar. Times given below are averages of a series of executions. Two series of results are presented that exhibits the behavior of the performance our algorithm either in function of the dimension of the N-cube or in function of the number of GCs to generate.

In Table 1 are provided the execution times for generating 100,000 cycles in N-cubes of dimensions 5 to 10. Smaller dimensions have too few GCs to get representative times (only one GC in dimensions 2 and 3 and 9 in dimension 4). As expected, the execution time increases with the dimension. However, the progression is much slower than the theoretical complexity obtained in Section 2.4. This means that this upper bound of complexity is largely overestimated. Indeed, many iterations in the inner loop of the generation algorithm are not performed, since it stops as soon as a parallel transition has been found. This is not taken into account in the current estimation of the complexity but it will deserve a deeper study to get a more accurate expression. The execution times for generating different amounts of GCs in dimensions between 5 and 9 are given in Figure 7. We observe similar behaviors for all dimensions with a sharp increase of the execution time with the number of GCs. This mainly comes from the fact that in the closure algorithm, GCs may be generated as much times as their number of neighbors according to the sub-sequence inversion relation. Moreover, this number of duplicate generations increases with the size of the input list, which itself increases with the size of the generated set. This is why we observe such an exponential increase. A strategy to avoid such duplicates during the transitive closure will be investigated in future works. Nonetheless, the absolute performance of the current version is rather efficient when considering that the production of several thousands of GCs takes less than one minute. To conclude this part on the performance of the generation process, a first multi-threaded implementation has been tested, which provides gains between 16% and 90% depending on the number of generated cycles and the dimension of the N-cube. However, this first implementation is far from being optimal and additional development effort is planed as a future work.

Statistical Study of the Generation Algorithm

This second series of experiments aims at analyzing the impact of the presence of isolated GCs according to the sub-sequence inversion relation, as mentioned in Section 2.4. Table 2 gives the percentages of isolated GCs in dimensions 2 to 5 of N-cubes. It is not possible to get exact percentages for greater dimensions due to the exponential increase of the number of GCs. However, it is interesting to see that the ratio of isolated GCs drastically decreases when the dimension of the N-cube increases. It can be proved that this tendency continues with greater dimensions. This comes from the condition that a GC must satisfy to be isolated, i.e., to generate no other GC by the inversion algorithm. The condition is that it must not have any valid parallel transitions in its entire sequence.

However, for one transition, the number of possible parallel transitions is directly linked to the dimension of the N-cube. Moreover, the number of transitions in a sequence exponentially increases with the dimension of the N-cube. Therefore, the probability of having at least one valid parallel transition in the entire sequence exponentially increases with the dimension of the N-cube. Indeed, the constraint to satisfy becomes stronger and stronger as the dimension increases. Thus, it comes that the probability that a GC is isolated, i.e., is a singleton in the partition of the set of GCs, exponentially decreases when the dimension of the N-cube increases. Finally, this implies that the constraint over the choice of GCs as input to the closure algorithm is very weak, as most GCs will lead to the generation of other GCs. First of all, the mixing time is bounded by log(ε -1 ) ρN and thus, thanks to Proposition 1 by log(ε -1 ) 3N (1 + e) + N log(N -2).

For instance, for ε = 10 -6 , N = 5, or N = 8 these bounds are respectively equal to 319 and 550. In the same context, i.e., ε = 10 -6 , we practically observe that the largest mixing time for all the Markov chains when N = 5 is equal to 42. It can be thus notably deduced that theoretical bound provided in Proposition 1 is notably much larger than the practical one.

Conclusions

Practical and theoretical contributions have been presented in this paper, respectively related to the field of Gray codes generation and mixing time of Boolean functions in dimension N .

Concerning Gray codes generation, an algorithm based on sub-sequence inversion has been designed and implemented. It starts from one known GC and it produces new GCs of the same dimension. It has been explained why the number of generations strongly depends on the configuration of the input GC. Then, as the upper bound of the number of generations from a single GC is limited, a transitive closure algorithm has been proposed too, that can generate very large subsets of the entire set of GCs in a given N-cube. Such algorithms are essential to study the Gray Code properties that may present a particular interest for use in a chaotic PRNG based on random walks or jumps in a N-cube.

Experiments have shown that our algorithm can produce several thousands of GCs in time of the order of the minute. Also, it has been shown that for a given dimension of N-cube, the number of isolated GCs, i.e., the codes that do not allow any sub-sequence inversion, is a negligible part of the total number of GCs. Moreover, this ratio decreases when the dimension increases. Therefore, the probability that our algorithm does not produce any new GC from a given input, is very small.

Concerning the mixing time of Boolean functions in the N-cube, the complete proof of an upper bound in O (N.log 2 (N )) has been given. This is the most accurate upper bound that has been proved to this day. Such result is useful to limit significantly the number of iterations when computing mixing times of Boolean functions in the N-cube. This is of particular interest when looking for such functions that would be adequate for the use in a CI-PRNG, i.e., functions with minimal mixing time, among a large number of candidates.

In future works, we plan to address the potential algorithmic enhancements mentioned in the experiments, and to perform deep statistical analyzes of Gray codes according to the mixing times of their respective Boolean functions in the N-cube as well as other features. This will help to identify the particular properties of GCs that could be used to find the codes best suited to the use in PRNGs. The results presented in this paper directly contribute to make such statistical analyzes possible.

Figure 3 :

 3 Figure 3: Inversion of sub-sequence (11, ..., 14) (dotted line) on face[START_REF] Contassot | Canonical Form of Gray Codes in N-cubes[END_REF][START_REF] Contassot-Vivier | Random Walk in a N-Cube Without Hamiltonian Cycle to Chaotic Pseudorandom Number Generation: Theoretical and Practical Considerations[END_REF][START_REF] Guyeux | Efficient and cryptographically secure generation of chaotic pseudorandom numbers on GPU[END_REF][START_REF] Gregory | Random walk: a modern introduction[END_REF] in the 4-cube.

1 6 d 18 found ← true 19 /

 61819 Input: a Gray code C of dimension n 2 Output: a list lst of Gray codes deduced from C 4 V ← 2 n // Number of vertices in C 5 for i from 0 to V-1 do // Parsing of sequence C ← dimOfTransition(C i , C (i+1)%V ) // Dimension used by the transition 7 for j from 0 to n-1 do // Parsing of dimensions of the N -cube 8 if j = d then // For any other dimension than d 9 src ← bitInversion(C i , j) // Neighbor of C i along dimension j 10 dst ← bitInversion(src, d) // Neighbor of src along dimension d 11 // vertices C i , C i+1 , src and dst form a face of the N -cube 12 // We look for the presence of transition src → dst in C 13 found ← false // Transition src → dst not yet found in C 14 k ← 5 // Starting search offset from i 15 while k ≤ V-5 and ¬found do // Search of transition src → dst in C 16 p ← i + k // probing transition src → dst at offset k from i 17 if C p%V = src and C (p+1)%V = dst then // transition found! / Bulding of a new cycle from C with sub-seq (C i+1 , ..., Cp) inverted 20 newCycle ← copyAndInvertSubSeq(C, i+1, p) 21 // Insertion of the new cycle in lst only if not yet present 22 exclusiveInsert(newCycle, lst)

  Possible parallel transitions (blue) for a reference transition (black).

  Two valid parallel transitions (red) for the reference transition (black) in the blue cycle.

Figure 4 :

 4 Figure 4: Examples of possible and valid parallel transitions in the 4-cube.

Figure 5 :

 5 Figure 5: Cycles in the 4-cube with no valid parallel transition.

Figure 6 :

 6 Figure 6: Construction of function f in dimension 3.

  dim

Figure 7 :

 7 Figure 7: Exec times (s) in function of the number of generated GCs for different dimensions of N-cube.

Table 2 :

 2 Percentages of isolated GCs in the partition of the set of GCs. 5.3 Discussion over the Mixing Time Distributions in Dimensions 5 and 8

14 endif 15 endforeach 16 endforeach 17

1 Input: a set I of initial Gray codes of dimension n 2 Output: a list lst of Gray codes obtained from I by transitive closure 4 lst ← I // Initial GCs are put in the list to avoid re-process 5 il ← I // Input list for the first step 6 finished ← false // Process not finished at the beginning 7 while ¬finished do // Main loop of the transitive closure 8 nl ← ∅ // Init of the list of GCs created for the 1st time 9 foreach i in il do // Processing of every GC i in the input list 10 gen ← generateFrom(i) // Generation of new GCs from source i 11 foreach j in gen do // Checking for insertion of each generated GC: 12 if exclusiveInsert(j, lst) then // in lst if first creation (not yet in it) 13 insert(j, nl) // and also in nl in this case (next step) il ← nl // New GCs at this step become input GCs of next step 18 finished ← decision(lst) // Stopping decision according to stopping criteria 19 endwhile 20 return lst

Table 1 :

 1 Generation times of 100,000 Gray Codes (GCs) in function of the dimension of the N-cube.

	5	6	7	8	9	10
	time (s) 19.61 19.93 28.17 53.60 174.43 332.29

Acknowledgments: This work is partially funded by the Labex ACTION program (contract ANR-11-LABX-01-01).

Now, for k = N -1, Equation (2) provides:

, where e is the Euler number. Consequently and using Equation (3), we have

Since the N -2th harmonic number N -2 k=1 1 k is lower than 1 + log(N -2), we have

Moreover, by definition of the αi's, N -1 k=1 α k = ρN -ρ1. It follows that ρN ≤ ρ1 + N e + N (1 + log(N -2)). Since ρ2 -ρ1 = α1 ≤ N e and 3 2 ρ1 = N + ρ2, one has 1 2 ρ1 ≤ N e + N , proving that ρ1 ≤ 2N (1 + e). In conclusion,

which concludes the proof.

A Coupling for the Markov Chain P

We now define a coupling for the lazzy Markov chain. An upper bound of the mixing time is obtained by computing the maximal expected collapsing time of the coupling ( [START_REF] David | Markov chains and mixing times[END_REF] Chap. 5).

Let (Un) be a sequence of i.i.d of random variables on {1, . . . , N }, and (Vn) be a sequence of i.i.d of random variables on {0, 1}.

Let (Xn, Yn)n ≥ 0 be a random sequence in Ω × Ω defined by: Xn+1 = r(Xn, Un, Vn), and Yn+1 = Xn+1 if Xn = Yn and, otherwise, Yn+1 = r(Yn, Wn, Zn) where:

Proposition 2. The sequence (Xn, Yn)n ≥ 0 is a coupling for the Markov Chain P .

Proof. By construction, (Xn) is a Markov Chain with transition matrix P . For (Yn), it is clear that Wn follows the uniform distribution on {1, . . . , N }.

Lemma 1. Let (dn) be the sequence of integers defined by dn = d(Xn, Yn). For any n, one has |dn+1 -dn| ≤ 1. Moreover:

N and P(dn+1 = 2) ≤ P(dn+1 = 0). Proof. Let X = Xn and Y = Yn. Assume first that dn ≥ 2, several cases may arise:

, then dn+1 = dn with probability 1 2 and dn+1 = dn -1 with probability 1 2 . It follows that, in this case, the probability that dn+1 = dn + 1 is null. The probability that

, then dn+1 = dn with probability 1 2 and dn+1 = dn -1 with probability 1 2 . If Un = f (Y ), then dn+1 = dn + 1 with probability 1 2 and dn+1 = dn with probability 1 2 . It follows that the probability that dn+1 = dn+1 is 1 2N , and the probability that

). This is a dual case of the previous one with similar calculus, switching f (X) and f (Y ).

, then dn+1 = dn with probability 1 2 and dn+1 = dn + 1 with probability 1 2 . Similarly, If Un = f (Y ), then dn+1 = dn with probability 1 2 and dn+1 = dn + 1 with probability 1 2 . Consequently, the probability that

∈ H(X, Y ), then dn+1 = dn. It follows that the probability that dn+1 = dn + 1 is null and the probability that

The following cases arise:

Therefore, the probability that dn+1 = dn +1 is null and the probability that

If Un = f (X), then dn+1 = 0 with probability 1 2 and dn+1 = dn = 1 with probability 1 2 . If Un = f (Y ), then dn+1 = dn + 1 = 2 with probability 1 2 and dn+1 = dn = 1 with probability 1 2 . Therefore, the probability that dn+1 = 2 is 1 2N and the probability that 2 . Therefore, the probability that dn+1 = 2 is 1 N and the probability that dn+1 = 0 = 1 N .

Let T (X, Y ) the random variable of the first coupling time of (Xn, Yn), with X0 = X and

Experiments

In this section is presented a set of different experiments in relation to the GC generation algorithm and the mixing time.

Performance of Gray Codes Generation Algorithms

The first series of experiments is devoted to the performance of the GC generation algorithm. The execution times shown below have been obtained on a desktop machine with two Intel(R) Xeon(R) CPU E5-2640 v2 @ 2.00GHz containing eight cores each and 128GB RAM. The OS is Linux Debian 3.16.43 and the program is written in C/C++ and uses the OpenMP multi-threading library. Indeed, the implementation of our algorithm is multi-threaded and takes advantage of the available cores in the machine. This reduces the execution times but does not change the overall complexity of the algorithm. Also, different compilers have been tested (GCC, Intel, PGI) and although slight variations


Eduardo H M Cruz

Matthias Diener

Laércio Lima Pilla

Philippe O A Navaux

Laércio L Pilla

EagerMap: A Task Mapping Algorithm to Improve Communication and Load Balancing in Clusters of Multicore Systems

published or not. The documents may come

Introduction

Optimizing the execution of parallel applications has become an important research topic in recent years [START_REF] Ballard | Communication lower bounds and optimal algorithms for numerical linear algebra[END_REF][START_REF] Deveci | Fast and High Quality Topology-Aware Task Mapping[END_REF]. Modern parallel systems consist of many interconnected cluster nodes that themselves constitute a parallel shared-memory machine with a deep memory hierarchy. In such architectures, communication and load imbalance can have a higher impact on application performance than computation [START_REF] Shalf | Exascale Computing Technology Challenges[END_REF][START_REF] Dally | GPU Computing to Exascale and Beyond[END_REF]. A common way to reduce the impact of communication is to improve the locality of communication, by increasing the use of local interconnections and reducing the use of remote ones [START_REF] Unat | Trends in data locality abstractions for hpc systems[END_REF][START_REF] Shalf | Exascale Computing Technology Challenges[END_REF]. Load imbalance can be improved by using dedicated load balancing algorithms. These optimizations can be performed in parallel applications that use message passing libraries such as MPI [START_REF] Chen | MPIPP: An Automatic Profile-guided Parallel Process Placement Toolset for SMP Clusters and Multiclusters[END_REF], but also in applications that use the shared memory paradigm, such as OpenMP or UPC [START_REF] Azimi | Enhancing Operating System Support for Multicore Processors by Using Hardware Performance Monitoring[END_REF][START_REF] Anbar | Exploiting Hierarchical Locality in Deep Parallel Architectures[END_REF].

In both types of parallel applications, the communication performance between the tasks of a parallel application are influenced by the levels of the hierarchy [START_REF] Wang | Performance Analysis of Thread Mappings with a Holistic View of the Hardware Resources[END_REF]. Communication through a shared cache memory or intra-chip interconnection is faster than communication between processors due to the slower inter-chip interconnections [START_REF] Luo | Thread data sharing in cache: Theory and measurement[END_REF][START_REF] Cruz | Dynamic thread mapping of shared memory applications by exploiting cache coherence protocols[END_REF]. Likewise, communication within a machine is faster than the communication between nodes in a cluster or grid. The network speeds also vary, increasing the difference in the communication latencies and bandwidths. Clusters and grids impose an additional challenge, since each node can have a different machine configuration. Tasks that communicate intensely should be mapped to PUs close together in the hierarchy. To improve load balance, we need to analyze the load of each task of the parallel application and map them to the cores in such a way that the loads of all cores are similar. In this context, the mapping of tasks to processing units (PUs) plays a key role in the performance of parallel applications [START_REF] Zhai | Efficiently Acquiring Communication Traces for Large-Scale Parallel Applications[END_REF].

Four main steps are required in the process of mapping tasks to the architecture. The first step is to detect the communication between the tasks. In message passage environments, the messages sent between the tasks need to be monitored. The second step is to detect the topology of the hierarchy, which is also highly dependent on the type of architecture. In shared memory architectures, tools such as hwloc [START_REF] Broquedis | hwloc: A Generic Framework for Managing Hardware Affinities in HPC Applications[END_REF] can be used, while most cluster environments have vendor-specific tools to manage topologies. The third step is to use a mapping algorithm to generate a mapping of tasks to PUs, combining communication, load, and topology information. The final step is to execute the application with the determined mapping, or migrate tasks to their assigned PUs.

This paper focuses on the third step: the mapping algorithm. The task mapping problem can be defined as follows [START_REF] Cruz | Dynamic thread mapping of shared memory applications by exploiting cache coherence protocols[END_REF]. Consider two graphs, one representing the parallel application, and one representing the parallel architecture. In the application graph, vertices represent tasks and edges represent the amount of communication between them. In the architecture graph, vertices represent the machine components, including the PUs, cache memories, NUMA nodes, network routers, switches, and others organized hierarchically, while edges represent the links' bandwidth and latency. The task mapping problem consists of finding a mapping of the tasks in the application graph to the PUs in the architecture graph, such that the total communication cost is minimized and the load on the cores is even.

The complexity of finding an optimal mapping is NP-Hard [START_REF] Bokhari | On the Mapping Problem[END_REF]. Due to the high number of tasks and PUs, finding an optimal mapping for an application is unfeasible. Heuristics are therefore employed to compute an approximation of the optimal mapping. However, current mapping algorithms still present a high execution time, since they were developed focusing on static mappings and are mostly based on complex analysis of the graphs. This reduces their applicability, especially for online mapping, since their overhead may harm performance.

In this paper, we propose EagerMap, an efficient algorithm to generate task mappings. It coarsens the application graph by grouping tasks that communicate intensely. This coarsening follows the topology of the architecture hierarchy. Based on observations of application behavior, we propose an efficient greedy strategy to generate each group. It achieves a high accuracy and is faster than other approaches of the state of the art. After the coarsening, we map the group graph to the architecture graph.

EagerMap was initially proposed in [START_REF] Cruz | An Efficient Algorithm for Communication-Based Task Mapping[END_REF]. The main extensions and improvements presented in this paper are the following.

The original algorithm supported symmetric tree topologies only, which limited its applicability mostly to shared memory machines. We extended it to also support clusters and grids, where the topology of the network that connects the cluster/grid is arbitrary and the topology internal to each compute node is a symmetric tree.

We also designed a parallel version of EagerMap, allowing it to calculate mappings faster.

We added load balancing support to EagerMap that can handle oversubscribed scenarios (more threads/processes than processing units).

Related Work

Previous studies evaluate the impact of task mapping considering communication [START_REF] Wang | Performance Analysis of Thread Mappings with a Holistic View of the Hardware Resources[END_REF], showing that it can influence several hardware resources. In shared memory environments, communication-based task mapping reduces execution time, cache misses and interconnection traffic [START_REF] Cruz | Hardwareassisted thread and data mapping in hierarchical multicore architectures[END_REF]. In the context of cluster and grid environments, mapping tasks that communicate to the same computing node reduces network traffic and execution time [START_REF] Brandfass | Rank reordering for MPI communication optimization[END_REF][START_REF] Ito | Automatically optimized core mapping to subdomains of domain decomposition method on multicore parallel environments[END_REF]. Communication costs can also be minimized in virtualized environments [START_REF] Sonnek | Starling: Minimizing Communication Overhead in Virtualized Computing Platforms Using Decentralized Affinity-Aware Migration[END_REF], demonstrating its importance for cloud computing.

Several mapping algorithms have been proposed to optimize communication and load balancing. Most traditional algorithms are based on graph partitioning, such as Scotch [START_REF] Pellegrini | Static Mapping by Dual Recursive Bipartitioning of Process and Architecture Graphs[END_REF]. The algorithm is based on the idea of divide and conquer, recursively allocating subsets of processes to subsets of PUs. The algorithm bipartitions an unprocessed set of PUs into two disjoint subsets, and calls a graph bipartitioning algorithm to split the subset of tasks associated with the PUs across the two subsets. As bipartitionings are performed recursively, the set sizes decrease, until all sets have only one element. A disadvantage of Scotch is that the topology description requires information about the communication latency between the PUs, and such values have a high influence in the final mapping. This presents a challenge for Scotch, since it is difficult to determine an accurate value for its parameters. On the other hand, EagerMap does not depend on such information, which makes it easier to use and can cause more reliable results. Other works that follow the same basic idea of Scotch are Zoltan [START_REF] Devine | Parallel hypergraph partitioning for scientific computing[END_REF], METIS [START_REF] Karypis | A Fast and High Quality Multilevel Scheme for Partitioning Irregular Graphs[END_REF][START_REF] Karypis | Metis -unstructured graph partitioning and sparse matrix ordering system[END_REF], and Chaco [START_REF] Hendrickson | The chaco users guide version 2.0[END_REF]. The work described in [START_REF] Glantz | Algorithms for Mapping Parallel Processes onto Grid and Torus Architectures[END_REF] also uses graph partitioning to group tasks, but it uses a greedy technique to map the groups to the topology. Some of these algorithms have been parallelized, such as PT-Scotch [START_REF] Chevalier | PT-Scotch: A Tool for Efficient Parallel Graph Ordering[END_REF] and Par-METIS [START_REF] Karypis | Parallel Multilevel K-way Partitioning Scheme for Irregular Graphs[END_REF].

Tree representations are used in TreeMatch [START_REF] Jeannot | Process Placement in Multicore Clusters: Algorithmic Issues and Practical Techniques[END_REF], which can lead to more optimized algorithms in architectures that can be represented in this way. However, the algorithm used in TreeMatch to group tasks has an exponential complexity because it generates all possible groups of tasks for each level of the memory hierarchy. This can lead to unreasonable mapping calculation time depending on the number of tasks. Another limitation of TreeMatch is that it is not able to calculate mapping for most clusters or grids, as their topology may not be represented by a tree. TreeMatch has also been ported to the Charm++ environment [START_REF] Jeannot | Topology and Affinity Aware Hierarchical and Distributed Load-balancing in Charm++[END_REF], where support for load balancing was included. [START_REF] Bak | Multi-level load balancing with an integrated runtime approach[END_REF] also proposes a load balancer for Charm++ with a focus on the integrated OpenMP runtime.

MPIPP [START_REF] Chen | MPIPP: An Automatic Profile-guided Parallel Process Placement Toolset for SMP Clusters and Multiclusters[END_REF] is a framework to find optimized mappings for MPI-based applications. MPIPP initially maps each task to a random PU. At each iteration, MPIPP selects pairs of tasks to exchange PUs to reduce communication cost as much as possible. The communication costs depend on chosen parameters, and can be related to the number of messages or amount of data transferred between the MPI processes. This process is repeated several times, improving the quality of the mapping on each iteration, and can go on until no gains are achieved. One of the main drawbacks of MPIPP is that it depends too much on the initial random mapping, which does not have any guarantees regarding the mapping quality. Due to this, MPIPP can generate bad mappings. Other works that are based on the concept of refining an initial mapping at each iteration are [START_REF] Boillat | A Fast Distributed Mapping Algorithm[END_REF] and [START_REF] Jovanovic | A heuristic algorithm for dynamic task scheduling in highly parallel computing systems[END_REF].

Our previous work [START_REF] Cruz | Dynamic thread mapping of shared memory applications by exploiting cache coherence protocols[END_REF] uses Edmonds' graph matching algorithm [START_REF] Edmonds | Maximum matching and a polyhedron with 0,1-vertices[END_REF] to calculate mappings, which solves the maximum weight perfect matching for complete weighted graphs to generate mappings. The idea of the proposal is to call the graph matching algorithm several times to generate the groups of threads that should be mapped together. Each time the matching algorithm is called, it generates pairs of threads that have a high degree of communication. Since it always generates pairs of threads, the solution is limited to environments where the number of tasks and PUs is a power of two. Another problem is that the proposed technique must be coded for a particular number of threads and architecture, which does not work automatically.

EagerMap: Greedy Hierarchical Mapping

EagerMap receives three pieces of input: a communication matrix containing the amount of communication between each pair of tasks, the load of each task and a description of the architecture hierarchy. It outputs which PU executes each task. To represent the architecture hierarchy, we use a tree, in which the vertices represent objects such as PUs and cache memories, and the edges represent the links between them. Our task grouping is performed with an efficient greedy strategy that is based on an analysis of the communication pattern of parallel applications. Fig. 1 depicts the different communication patterns of several parallel benchmark suites [START_REF] Bailey | The NAS Parallel Benchmarks[END_REF][START_REF] Jin | The OpenMP implementation of NAS Parallel Benchmarks and Its Performance[END_REF][START_REF] Luszczek | Introduction to the HPC Challenge Benchmark Suite[END_REF][START_REF] Bienia | The PARSEC Benchmark Suite: Characterization and Architectural Implications[END_REF], which we obtained using the methodology described in Section 4.2. We observe three essential characteristics in the communication behavior of the applications that need to be considered for an efficient mapping strategy: 1. There are two types of communication behavior: structured and unstructured communication. In applications with structured communication, each task communicates more with a subgroup of tasks, such that mapping these subgroups to PUs nearby in the hierarchy can improve performance. In Fig. 1, all applications except Vips show structured communication patterns. Our mapping algorithm is designed to handle structured communication patterns, be-cause in applications with unstructured communication, there may not be a task mapping that can improve performance. 2. In applications with structured communication patterns, the size of the subgroups with intense internal communication is usually small when compared to the total number of tasks in the parallel application. For instance, in the communication pattern of CG-MPI (Fig. 1b), subgroups of 8 tasks communicate intensely, out of 64 tasks in total.

3. The amount of communication within each subgroup is much higher than the amount of communication between different subgroups.

In this section, we describe EagerMap in detail, give an example of its operation and discuss its complexity. We first describe the algorithm only considering communication, and then we extend it to consider both communication and load.

Description of the EagerMap Algorithm

The algorithm requires two previously initialized variables: nLevels and execElInLevel. nLevels is the number of shared levels of the architecture hierarchy plus two. This addition is required to create a level to represent application tasks (level 0) and another level to represent the processing units (level 1). execElInLevel is a vector with nLevels positions. execElInLevel[0] is not used. execElInLevel [START_REF] Anbar | Exploiting Hierarchical Locality in Deep Parallel Architectures[END_REF] contains the number of processing units. For positions i, such that 1 < i < nLevels, the value is the number of hardware objects on the respective architecture hierarchy level. Hardware objects are cores, caches, processors and NUMA nodes, among others. For instance, execElInLevel [START_REF] Azimi | Enhancing Operating System Support for Multicore Processors by Using Hardware Performance Monitoring[END_REF] can be the number of cores, execElInLevel [START_REF] Bailey | The NAS Parallel Benchmarks[END_REF] the number of last level caches, and execElInLevel[nLevels-1] the number of NUMA nodes. Since private levels of the architecture hierarchy are not important for our mapping strategy, we only consider the shared levels when preparing both nLevels and execElInLevel.

Top Level Algorithm

The top level mapping algorithm is shown in Algorithm 1. It calculates the mapping for each level on the architecture hierarchy. The groups variable represents the groups of elements for the level being processed. The previousGroups variable represents the groups of elements of the previous level. First, it initializes groups with the application tasks (loop in line 2). Afterwards, it iterates over all levels on the architecture hierarchy, in line 8. After generating the groups of tasks for a level (line 10) (explained in Section 3.1.2), the algorithm generates a new communication matrix (line 12, discussed in Section 3.1.3). This step is necessary since we consider each group of tasks as the base element for mapping Algorithm 2: GenerateGroupsF orLevel: Generates the groups for a level of the architecture hierarchy. on the next hierarchy level.

Input: commMatrix[][],
The groups variable implicitly generates a tree of groups. Level 0 represents the tasks. Level 1 represents groups of tasks. Level 2 represents groups of groups of tasks. In other words, on each level a new application graph is generated by coarsening the previous level. After the loop in line 8 finishes, the groups variable represents the hierarchy level nLevels -1 and contains nElements elements. We set up rootGroup to point to these elements of the highest level (for loop in line 17). Finally, the algorithm maps the tree that represents the tasks, rootGroup, to the tree that represents the architecture topology, archT opologyRoot. This procedure is explained in Section 3.1.4.

Generating the Groups for a Level of the Architecture Hierarchy

The GenerateGroupsF orLevel algorithm, described in Algorithm 2, handles the creation of all groups for a given level of the architecture hierarchy. It expects that the levels of hierarchy up to the previous processed level to be already grouped in previousGroups. The maximum number of groups is the number of hardware objects of that level. The selection of which elements belong to each group is performed by GenerateGroup.

The GenerateGroup algorithm, in Algorithm 3, groups elements that present a large amount of communication among themselves. For the grouping, the algorithm adopts a greedy strategy as follows. Each iteration of the loop in line 2 adds one element to the group. The added element, expressed by the Algorithm 3: GenerateGroup: Generates one group of elements that communicate. winner variable, is the one that presents the largest amount of communication relative to the elements already in the group. The chosen variable is used to avoid selecting the same element more than once. GenerateGroup can be parallelized in the loop of line 4, where each thread would compute its local winner in parallel, as we will explain in Section 3.5.

Input: commMatrix[][],

Computing the Communication Matrix for the next Level of the Architecture Hierarchy

RecreateM atrix, described in Algorithm 4, regenerates the communication matrix to be used for the next level of the architecture hierarchy. The new communication matrix has an order of nGroups. It contains the amount of communication between the groups. It is calculated by summing up the amount of communication between the elements of different groups.

Mapping the Group Tree to the Architecture Topology Tree

The algorithm to map the group tree to the architecture topology tree is M apGroupsT oT opology, detailed in Algorithm 5. It performs a recursion over the levels of the architecture hierarchy. The recursion stop condition happens when it reaches the lowest level of the architecture topology, the processing unit (PU) (line 2). If the stop condition is not fulfilled, the algorithm already knows that the maximum number of groups per level never exceeds the number of hardware objects of that level, as explained in GenerateGroupsF orLevel.

Therefore, if the level of the architecture hierarchy in the recursion is shared (line 6), the algorithm only assigns one hardware object of the following level to each group and recursively calls itself for each element of the following level.

Otherwise, the level is private and is not considered for mapping (line 10).

Mapping Example in a Multi-level Hierarchy

For a better understanding of how our mapping algorithm works, we present an example of mapping 16 tasks in a cluster of 2 machines, each consisting of 8 cores, with L3 caches shared by 2 cores and 2 processors per machine. Fig. 2a illustrates the hierarchy of the machine, which needs to be transformed in the topology tree shown in Fig. 2c. The topology tree can be generated either automatically, using tools such as hwloc [START_REF] Broquedis | hwloc: A Generic Framework for Managing Hardware Affinities in HPC Applications[END_REF], or manually by the programmer. Private levels, with arity 1, such as L1 and L2 caches, as well as the processor, are not relevant. The global variable nLevels has the value of 4, since there are 4 levels in the hierarchy. Therefore, execElInLevel has 4 positions. Position 0 represents the tasks and is not used. Position 1 represents the PUs (the cores in this case) and its value is 8. Position 2 represents the L3 caches and its value is 4. Finally, position 3 represents the machines and its value is 2.

Regarding Algorithm 1 (M apAlgorithm), the loop in line 8 is executed 3 times. In the first iteration, it generates the groups of tasks that execute in each core. Since the number of tasks is 16 and the number of cores is 8, it generates 8 groups of 2 tasks each. The communication matrix used in the first iteration is shown in Fig. 3a. We demonstrate how the first group is generated in Fig. 4 (GenerateGroup, Algorithm 3). When i is 0, the winners vector is empty and task 0 is selected to be the first one of the group. When i is 1, the task that presents the most communication to the tasks in the winners vector is task 1. Therefore, task 1 is selected as winner and enters the group. Algorithm GenerateGroup then returns a group formed by tasks [0, 1], implicitly stored in the group.elements vector. GenerateGroupsF orLevel (Algorithm 2) repeats the same procedure until all groups of the level are formed. Then, the communication matrix for the next level, in Fig. 3b, is calculated by RecreateM atrix (Algorithm 4). The second iteration of the loop in line 8 of M apAlgorithm behaves similarly, but selects which groups of tasks from the previous iteration share each L3 cache. Since the number of groups was 8 and there are 4 L3 caches, it generates 4 new groups with 2 elements from the previous level in each.

In the third iteration of the loop in line 8 of M apAlgorithm, it selects which groups of groups of tasks share each machine. The communication matrix illustrated in Fig. 3c is used. Since the number of groups generated in the previous iteration was 4 and there are 2 machines, it generates 2 new groups of 2 elements each. To complete the group tree, the rootGroup variable points to these 2 groups. Afterwards, in line 19 of M apAlgorithm, the task group tree shown in Fig. 2b is finished.

M apGroupsT oT opology (Algorithm 5) maps the task group vertices to ver-tices of the architecture hierarchy. The mapping begins from the root node. Tasks [0 -7] are mapped to Machine 0. Tasks [0 -3] are mapped to L3 cache 0. Tasks [0 -1] are mapped to PU 0 (core 0). The recursion continues until all groups of tasks are mapped to a PU (core). In this example, tasks x and x + 1, where x is even, will be mapped to core x /2.

Complexity of EagerMap

For the analysis of the complexity of EagerMap, we first introduce the variables used in the equations. E is the number of elements to be mapped in the current level of the architecture hierarchy (tasks or groups of tasks), which is different for each level. G is the number of elements per group. P is the number of processing units. N is the number of tasks to be mapped. L is the number of levels on the architecture hierarchy. To make it easier to calculate the complexity, we consider that P ≤ N . The complexity of GenerateGroup is:

G i=1 E j=1 i = G i=1 E • i ≤ E • G 2 (1)
The complexity of GenerateGroupsF orLevel is shown in Equation 2. The number of groups is E G . Also, E is an upper bound for G.

E/G i=1 GenerateGroup ≤ E/G i=1 E • G 2 ≤ E G • E • G 2 ≤ E 3 (2)
The complexity of RecreateM atrix is O(E 2 ).

The complexity of M apGroupsT oT opology is the same as performing a depth-first search, O(V + C), where V is the number of vertices and C is the number of edges. Since our groups variable implicitly forms a tree, we know that C = V -1. The last level of this tree has N vertices, and the penultimate level has P vertices. The number of vertices of the other levels is the number of the previous level, divided by the number of sharers. The number of sharers is greater than 1 since we only keep track of shared levels. Therefore, V ≤ N + P + 2 • P = O(N ). With this, the complexity of M apGroupsT oT opology is O(N ).

The complexity of the top level algorithm, M apAlgorithm, depends on all previous algorithms, as shown in Equation 3. To calculate the complexity, we have to take into account that the value of E changes on each level of the architecture hierarchy.

L i=0

GenerateGroupsF orLevel + RecreateM atrix +M apGroupsT oT opology

(3) = L i=0 O E 3 i + O E 2 i + O(N ) (4)
We can rewrite Equation 4 as Equation 5 by considering that the algorithm iterates L + 1 times (L levels of the architecture topology plus one level for the Algorithm 6: LbGenerateGroupsF orLevel: Generates the groups for a level of the architecture hierarchy. tasks) and that only shared topology levels are represented, which means that, in the worst case scenario, the topology will be a binary tree with P leaves and a maximum number of vertices equal to 2P -1.

Input: commMatrix[][],
≤ L i=1 O P 2 i-1 3 + O P 2 i-1 2 + O(N 3 ) + O(N 2 ) + O(N ) (5) ≤ 2 O(P 3 ) + O(P 2 ) + O(N 3 ) + O(N 2 ) + O(N ) = 3 O(N 3 ) + O(N 2 ) + O(N ) = O N 3 (6)
With the simplifications shown in Equation 6, we find that EagerMap has a polynomial complexity of O(N 3 ).

Considering the Task Load

In a parallel application, it is very common to have its tasks use the CPU for different amounts of time and thereby have different computational power requirements. We call the CPU power a task requires to complete its execution the task load. In order to consider the task load when mapping tasks, and thereby balancing the load, the algorithms presented in the previous sections need slight modifications. Most modifications are done in the algorithms GenerateGroupsF orLevel (Algorithm 2) and GenerateGroup (Algorithm 3). In these algorithms, the total number of elements is related to the number of Algorithm 7: LbGenerateGroup: Generates one group of elements that communicate. tasks, and they try to keep the number of tasks per group as close as possible to nElements/nGroups. To consider the task load, we modify the concept of an element to represent the task load. In this way, instead of the task, the algorithm uses the task load. And instead of the total number of elements, the algorithm uses the total task load. The idea is to group tasks that communicate in the same groups, but limit the group sizes by the load, not number of tasks.

Input: commMatrix[][],
The modified algorithms to consider the load can be found in Algorithms 6 (LbGenerateGroupsF orLevel) and 7 (LbGenerateGroup). The key idea is that the maximum load of a group (loadT hreshold) is recalculated on every iteration of the loop (line 9 in Algorithm 6). We do this to get a more even balance of the load for the next groups, since the actual load of the group rarely is exactly totalLoad/avlGroups due to the different loads between the tasks. If this is not done, the loads of the last groups are usually very low. Regarding Algorithm 7, the main difference is that the group size is limited by the average load loadT hreshold, in line 4, while in the non load balancing version of the algorithm the group size was limited to the average number of tasks per group.

Parallel EagerMap

The function GenerateGroup is the only function in the critical path of Ea-gerMap that can be parallelized without changing the serial version behavior. RecreateM atrix and M apGroupsT oT opology could also be parallelized. However, the overhead of parallelizing them surpasses the benefits. The parallel version of EagerMap consists of parallelizing the loop of line 4 of Algorithm 3. This loop is responsible for finding the next task that will belong to a group. The idea of the parallel algorithm is that we can search for the next task in parallel, where each thread evaluates a subset of the total tasks and finds its local task that maximizes the group communication.

The parallel algorithm is shown in Algorithm 8. We express the parallelization in a notation similar to the OpenMP standard. In the beginning of the parallel block, in line 4, we initialize the vector threadResult containing the local task that maximizes the group communication. The local task is searched in parallel in the for loop in line 10. After that, the task that maximizes communication is selected from the local task found by each thread, in line 24. Although there are no critical regions, we do not expect a linear speedup because the parallel region has a short duration. The speedup should be higher with a larger number of tasks.

The function calls to GenerateGroup related to the first hierarchy level are the most expensive procedures, since they are processing the grouping of tasks. Thereby, for the first hierarchy level, we use the parallel version. For the subsequent levels, the cost of GenerateGroup falls dramatically, such that is better to use the sequential version.

Running EagerMap in Non-Tree Topologies

The previous algorithms presented in this section work only in symmetric tree topologies. This makes them applicable only to shared memory machines, or clusters where all machines have the same topology. In this section, we explain how to extend EagerMap to work on any kind of cluster or grid. The idea is to make a pre-computation of which tasks will run on each machine of the cluster, and then call the previously presented algorithms for each machine separately. We consider that the communication within each machine is much faster than the communication between different machines.

The algorithms that perform this pre-computation are very similar to the others and are also based on a greedy grouping strategy. The first thing the algorithm does is to calculate the number of tasks that each machine will execute. For that, it counts the total number of available PUs in all machines, and spreads the tasks in such a way that every machine will have the closest possible number of tasks per PU. For instance, if the number of tasks is equal to the total number of PUs, each machine will have 1 task per PU. After that, it generates the group of tasks that each machine will run in a greedy way similarly to Algorithm 3, where the next task to be added to a group is the one that maximizes the communication to the elements already within the group.

The final step, as already explained, is to call EagerMap separately for each group. In this way, we can handle any kind of cluster/grid regardless of the topology or if the machines have different configurations. The load balance can be taken into account in this step in a very similar way to what is explained in Section 3.4. The division of the tasks that will run on each group can be parallelized as explained in Section 3.5.

Evaluation Methodology of EagerMap

In this section, we show how we evaluate our proposal. We discuss the benchmarks and architecture used in our evaluation, how we obtain the communication matrices, and other mapping strategies used for comparison.

Benchmarks

For the evaluation of the mapping algorithms, we use applications from the MPI implementation of the NAS parallel benchmarks (NPB) [START_REF] Bailey | The NAS Parallel Benchmarks[END_REF], the OpenMP NPB implementation [START_REF] Jin | The OpenMP implementation of NAS Parallel Benchmarks and Its Performance[END_REF], the High Performance Computing Challenge (HPCC) benchmark [START_REF] Luszczek | Introduction to the HPC Challenge Benchmark Suite[END_REF] and the PARSEC benchmark suite [START_REF] Bienia | The PARSEC Benchmark Suite: Characterization and Architectural Implications[END_REF]. The NAS benchmarks were executed with input size B, HPCC with an input matrix with 4000 2 elements and PARSEC was executed with the native input size.

Most experiments are based on applications with a static communication behavior. HPCC, which has several communicative phases, was used to show that efficient online algorithms, such as EagerMap, are important for applications with a dynamic behavior. All applications were configured to create a number of threads equivalent to the number of PUs of the machine.

Generating the Communication Matrices

For the MPI-based benchmarks, we used the eztrace framework [START_REF] Trahay | EZTrace: a generic framework for performance analysis[END_REF] to trace all MPI messages sent by tasks and built a communication matrix based on the number of messages sent between tasks. We also generated the communication matrices using the number of bytes and, although the value of each individual cell was different, the overall pattern was the same. For the HPCC benchmark, which has a different communication pattern for each step, we also use eztrace to generate the matrices statically, since the communication pattern does not change between executions. For the benchmarks that use shared memory for implicit communication using memory accesses, we built a memory tracer based on the Pin binary analysis tool [START_REF] Luk | Pin: Building Customized Program Analysis Tools with Dynamic Instrumentation[END_REF], similarly to [START_REF] Barrow-Williams | A Communication Characterisation of Splash-2 and Parsec[END_REF]. This tool traces all memory accesses from the tasks. We build a communication matrix by comparing memory accesses to the same memory address from different tasks and increment the matrix on every match.

Generating the Task Load

To measure the task load, we count the amount of instructions executed by each thread, which better represents the load than the traditional CPU time [START_REF] Cruz | Improving communication and load balancing with thread mapping in manycore systems[END_REF]. Therefore, the load between the tasks are imbalanced when the amount of instructions executed by each thread is different.

Hardware Architecture

The hardware architecture used in our experiments regarding shared memory consists of 4 Intel Xeon X7550 processors for a total of 64 PUs. In this architecture, there are three possibilities for communication between tasks. Tasks running on the same core can communicate through the fast L1 or L2 caches and have the highest communication performance. Tasks that run on different cores have to communicate through the slower L3 cache, but can still benefit from the fast intra-chip interconnection. When tasks communicate across physical processors, they need to use the slow inter-chip interconnection.

For the cluster experiments, we used OpenMPI 1.6.5 and two types of machines: (A) 2 Intel Xeon E5530, with a total of 16 PUs, and (B) 2 Intel Xeon CPU E5-2640v2, with a total of 32 PUs. We used only half of the PUs of the machines because it presented better performance than using all PUs due to SMT and poor support from OpenMPI. We evaluated [START_REF] Barrow-Williams | A Communication Characterisation of Splash-2 and Parsec[END_REF] To evaluate the load-aware EagerMap, we used a Xeon Phi Coprocessor model 3120P, with the Knights Corner architecture. We used this architecture because, due to the large number of cores (57 cores, each 4-way SMT), load balancing plays a key role in its performance.

Comparison

We compare EagerMap to (i) a Random mapping, which is an average result of 30 different random mappings; (ii) TreeMatch [START_REF] Jeannot | Near-optimal placement of MPI processes on hierarchical NUMA architectures[END_REF][START_REF] Jeannot | Process Placement in Multicore Clusters: Algorithmic Issues and Practical Techniques[END_REF]; (iii) Scotch [START_REF] Pellegrini | Static Mapping by Dual Recursive Bipartitioning of Process and Architecture Graphs[END_REF]; and (iv) MPIPP [START_REF] Chen | MPIPP: An Automatic Profile-guided Parallel Process Placement Toolset for SMP Clusters and Multiclusters[END_REF]. Most performance results are normalized to the performance of the default operating system mapping, the Linux Completely Fair Scheduler [START_REF] Jones | Inside the linux 2.6 completely fair scheduler[END_REF], which focuses on keeping the load balanced between the tasks and cores.

Results

We evaluate the performance and quality of our algorithm, as well as the performance improvements obtained when mapping the tasks of the applications.

Performance of the Mapping Algorithms

Fig. 5 shows the execution time of the four mapping algorithms for all benchmarks. For 128 tasks, EagerMap is about 10 times faster than Scotch, 1000 times faster than TreeMatch, and 100,000 times faster than MPIPP. TreeMatch has an exponential complexity, so a much higher execution time is expected for it. Due to this, the difference in execution time between TreeMatch and EagerMap increases together with the number of tasks. MPIPP is much slower than the other mechanisms because it needs to perform several iterations to refine its initial random mapping. It did not finish executing when the number of tasks was higher than 128.

Quality of the Mapping

The quality of the calculated mapping determines the benefits that can be achieved. Quality is measured by the amount of locality achieved. We use Equation 7 to calculate the quality, which is provided in the source code of TreeMatch. In this equation, n is the number of tasks, M [i][j] is the amount of communication between tasks i and j, map[x] is the processing unit (PU) mapped, and lat[a][b] is the latency of the PUs in the hierarchy. We calculated the latencies using LMbench [START_REF] Mcvoy | Lmbench: Portable Tools for Performance Analysis[END_REF].

M appingQuality = n-1 i=1 n j=i+1 M [i][j] lat map[i] map[j] (7)
Fig. 6 presents the mapping quality results for the communication matrices previously illustrated in Fig. 1. Applications with more structured communication patterns presented the highest improvements, as expected. CG-MPI, LU-MPI and HPCC (Phase 7) presented the best improvements because their communications can be easily optimized by mapping neighboring tasks together. In BT-MPI and MG-MPI, near and distant tasks communicate, presenting more challenges for the mapping algorithm. The reason is that if neighboring tasks are mapped together, the communication between distant tasks does not improve. Likewise, mapping distant tasks together does not improve communication between neighboring tasks. In applications with less structured patterns, such as Streamcluster, a lower improvement is expected because the ratio of communication between tasks that communicate more and tasks that communicate less is lower. Vips is the only application with unstructured communication, such that no task mapping was able to improve its communication.

The quality obtained with MPIPP is similar to the random mapping, since MPIPP is based on refining an initial random mapping. Although MPIPP can work with few tasks, when the number of tasks increases, the possibilities of permutation increase exponentially. This makes it more difficult to find new combinations to improve the initial solution. EagerMap, TreeMatch and Scotch presented similar results for all applications. This result demonstrates that EagerMap is able to achieve results as good as more complex algorithms due to the characteristics of the communication patterns we observed in Section 3.

Speedup of the Parallel Version of EagerMap

The speedup of the parallel version of EagerMap when mapping different numbers of tasks can be found in Fig. 7. We analyze the speedup relative to the number of tasks because it is the parameter with most influence in execution time. The task tree structure has little impact on the speedup, since the parallel algorithm is applied only to the first levels of the hierarchy (usually only to the first level), as explained in Section 3.5. As we explain in Section 3.5, we do not expect a linear speedup because the duration of the parallel phase is short. The parallel phase is restarted several times during the mapping calculation, imposing overhead in the parallelization. On the other hand, the impact of the overhead is decreased when the number of tasks to be mapped increases, as can be seen in the figure. This happens because, with more tasks, the duration of the parallel phase increases. With 131072 tasks, the speedup using 32 threads was 13.7, and should be higher for more tasks.

It is important to note that, when the number of tasks is large, the usage of the parallel version of EagerMap is encouraged. This happens because the execution time grows fast with the increase of the number of tasks. For instance, with 32768 tasks, the single threaded execution required 43 minutes. With 131072 tasks, the single threaded execution time grew to 25 hours. Using the parallel version with 32 threads, the execution time was reduced to 5.8 minutes and 1.8 hours, respectively.

Performance Improvements in a Shared-Memory Environment

We executed applications using the mapping obtained with the algorithms in the machine with Intel Xeon X7550 processors. The execution time results for the MPI and OpenMP NAS Benchmarks are shown in Fig. 8. For these applications, since their communication pattern is stable, we calculated the mapping statically. As a case study for online mapping, we used the HPCC Benchmark, shown in Fig. 9, since it contains 16 phases with different communication behaviors. For each phase, we call the mapping algorithm and migrate tasks accordingly. Execution time results correspond to the average of 30 executions, and are normalized to the execution time of the original policy of the operating system. The confidence interval is smaller than 1% for all algorithms.

In general, applications with more structured communication patterns present higher performance improvements. For instance, CG-MPI's performance was significantly improved when compared to the OS mapping, since its communication pattern shows high potential for mapping, as discussed in Section 5.2. On average, EagerMap improved performance by 10.3%, while TreeMatch, Scotch and MPIPP showed improvements of 9.6%, 8.7% and 8.4%, respectively. The results of the HPCC Benchmark (Fig. 9) show a larger improvement by EagerMap than related work. Also, the lower overhead of our algorithm does not harm the application performance as much as TreeMatch and MPIPP.

Performance Improvements from the Load-Aware Algorithm

The performance improvements from the load-aware EagerMap in Xeon Phi is shown in Fig. 10. On average, EagerMap improved the performance by 44.9%, while Scotch, Scatter and Compact improved by an average of 17.8%, 19.3% and 9.6%, respectively. Scatter and Compact do not perform any analysis on the communication and load, following only pre-defined rules. This shows that EagerMap, by actually considering the communication and load aspects of the applications, can provide a better performance improvement. On the other hand, Scotch also takes into account both the load and communication to calculate the mapping. This indicates that EagerMap can make a better use of this information than Scotch.

EagerMap

Performance Improvements in Clusters

The results of the cluster experiments with applications from the MPI NAS benchmarks can be found in Fig. 11. The results are normalized to the random mapping. The mapping was performed only statically because OpenMPI does not support task migration between nodes.

EagerMap achieves performance results similar to Scotch. However, as shown in Section 5.1, the time EagerMap requires to calculate the mapping is about 10 times faster than Scotch. In the experiments with 32 tasks, in the 4×8, 2×16 and Mixed-32 clusters, the application CG was the one that achieved the best results, where the performance was improved by up to 394% in Mixed-32. In 4×16, with 64 tasks, LU was the application with the best improvements, of 286%. BT and SP use less cores than the total in all clusters except 4×16 and Mixed-64 (the applications use a quadratic number of tasks), and Scotch presents better results. In 4×16 and Mixed-64, EagerMap was better. This shows that EagerMap is better at optimizing for communication, while Scotch is better at balancing the load between the machines. Regarding MG, the Interleave policy achieves good results because, coincidentally, its mapping is very close to the best mapping possible for MG.

Conclusions

The mapping of tasks to PUs in parallel architectures influences the communication performance and load balance. The communication performance can be improved by mapping tasks that communicate to PUs nearby in the memory hierarchy or to the same node in clusters or grids, making use of faster interconnections. The load balance can be improved by mapping the tasks in such a way that the load is evenly distributed among the PUs. The mapping algorithm selects which PU will execute each task and plays a key role in these types of mapping. When the mapping algorithm considers both metrics, more challenges arise, since these metrics can lead to contradictory decisions.

In this paper, we proposed EagerMap, that originally worked only in sharedmemory architectures and only considered the communication pattern. We added support for load balancing, mapping in clusters and grids, and proposed a parallel version to accelerate its execution. In contrast to previous work, it adopts a more efficient method to select which tasks should be mapped together, based on an analysis of the communication pattern of parallel applications. We performed experiments with a large set of benchmarks with different communication characteristics. Results show that EagerMap calculated better task mappings than the state of the art, with a drastically lower overhead and better scaling, which makes EagerMap more suitable for online mapping.

As future work, we intend to study EagerMap to allow it to execute in a cluster in a distributed way.

EagerMap is licensed under the GPL and is available at http://github.com/ehmcruz/eagermap.

Figure

  Figure 1: Example communication matrices for parallel applications consisting of 64 tasks. Axes represent task IDs. Cells show the amount of communication between tasks. Darker cells indicate more communication. The communication matrices were generated as explained in Section 4.2.

  Task group tree. (c) Architecture hierarchy tree.

Figure 2 :

 2 Figure 2: Mapping 16 tasks in an architecture with 8 PUs, L3 caches shared by 2 PUs, 2 processors per machine, and 2 machines.

Figure 3 :

 3 Figure 3: Communication matrices used in our example. A box is drawn around the groups of tasks generated by the algorithm.

Figure 4 :

 4 Figure 4: Example of the GenerateGroup algorithm, when generating the first group in our example.

Figure 5 :

 5 Figure 5: Execution time (in ms) of the mapping algorithms, for different numbers of tasks to be mapped.

Figure 7 :

 7 Figure 7: Speedup of the parallel version of EagerMap with different number of tasks to be mapped, varying the number of execution threads.

Figure 8 :

 8 Figure 8: Performance improvement compared to the operating system mapping of the NAS-MPI and NAS-OMP benchmarks in a shared memory environment.

Figure 10 :

 10 Figure 10: Performance improvements from the load-aware EagerMap in Xeon Phi for the NAS OpenMP benchmarks compared to the operating system mapping.

Figure 11 :

 11 Figure 11: Performance improvement compared to the operating system mapping of the NAS-MPI benchmarks in the clusters.

  The top level algorithm of EagerMap.

	Algorithm 1: M apAlgorithm: Input: commMatrixInit[][], nTasks
		Output: map[]
		LocalData: nElements, i, nGroups, rootGroup, commMatrix[][], groups[],
		previousGroups[]
		GlobalData: nLevels, execElInLevel[], hardwareTopologyRoot
	1 begin
	2	for i←0 ; i<nTasks ; i←i+1 do
	3	groups[i].id ← i;
	4	groups[i].nElements ← 0;
	5	end
	6	nElements ← nTasks;
	7	commMatrix ← commMatrixInit;
	8	for i←1 ; i<nLevels ; i←i+1 do
	9	previousGroups ← groups;
		/* GenerateGroupsForLevel is implemented in Algorithm 2.	*/
		1: Example communication matrices for parallel applications consisting
	of 64 tasks. Axes represent task IDs. Cells show the amount of communication
	between tasks. Darker cells indicate more communication. The communication
	matrices were generated as explained in Section 4.2.

10 [nGroups, groups] ← GenerateGroupsForLevel(commMatrix, nElements, i, previousGroups, execElInLevel[i]); 11 if i < nLevels-1 then /* RecreateMatrix is implemented in Algorithm 4. */ 12 commMatrix ← RecreateMatrix(commMatrix, groups, nGroups); 13 end 14 nElements ← nGroups; 15 end 16 rootGroup.nElements ← nElements; 17 for i←1 ; i<nElements ; i←i+1 do 18 rootGroup.elements[i] ← groups[i]; 19 end /* MapGroupsToTopology is implemented in Algorithm 5. */ 20 MapGroupsToTopology(archTopologyRoot, rootGroup, map); 21 return map; 22 end

Algorithm 4 :

 4 RecreateM atrix: Calculates the communication matrix for the next level.

		Machine		Machine
		Processor	Processor	Processor	Processor
		Core Core	Core Core	Core Core	Core Core
	Input: commMatrix, groups[], nGroups L1/L2 L1/L2 L1/L2 L1/L2	L1/L2 L1/L2	L1/L2 L1/L2
	Output: newCommMatrix[][] LocalData: i, j, k, z, w L3	L3	L3	L3
	begin
		for i←0 ; i<nGroups-1 ; i←i+1 do for j←i+1 ; j<nGroups ; j←j+1 do Interconnection	Interconnection
	4	w ← 0;
	5	for k←0 ; k<groups[i].nElements; k←k+1 do
	6	for z←0 ; z<groups[j].nElements; z←z+1 do (a) Topology used in the example.
		root
		0,1,2,3 4,5,6,7
	0,1,
		end
	10	else
	11	MapGroupsToTopology(hardwareObj.linked[0], group, map);
	12	end
	13 end

7 w ← w + commMatrix[ groups[i].elements[k].id ][ groups[j].elements[z].id ]; 8 end 9 end 10 newCommMatrix[i][j] ← w; 11 newCommMatrix[j][i] ← w; 12 end 13 end 14 return newCommMatrix ; 15 end Algorithm 5: M apGroupsT oT opology: Maps the group tree to the hardware topology tree. Input: hardwareObj, group, map[] LocalData: i begin if hardwareObj.type = ProcessingUnit then for i←0 ; i<group.nElements ; i←i+1 do 4 map[ group.elements[i].id ] ← hardwareObj.id; end else if hardwareObj.nSharers > 1 then for i←0; i<group.nElements; i←i+1 do 8 MapGroupsToTopology(hardwareObj.linked[i], group.elements[i], map);

  Algorithm 8: P arallelGenerateGroup: Generates one group of elements that communicate.

	Input: commMatrix[][], totalElements, groupElements, chosen[], previousGroups[]
	Output: group
	LocalData: i, j, wMax, nt, winners[], winner, threadResult[]
	1 begin
	2	for i←0 ; i<groupElements ; i←i+1 do
	3	wMax ← -1;
	4	#pragma omp parallel
		ThreadData: id, j, w, k
		end
		end
		end
		end
		wMax ← -1;
		for j←0 ; j<nt ; j←j+1 do
		if threadResult[j].wMax>wMax then
		wMax ← threadResult[j].wMax;
		winner ← threadResult[j].winner;
		end
		end
		chosen[winner] ← 1;
		winners[i] ← winner;
		group.elements[i] ← previousGroups[winner];
		end
		return group
	35 end

5 id ← getThreadID(); 6 #pragma omp master 7 nt ← getNumberOfThreads(); 8 threadResult[id].wMax ← -1; 9 #pragma omp for for j←0 ; j<totalElements ; j←j+1 do if chosen[j]=0 then w ← 0; for k←0 ; k<i ; k←k+1 do 14 w ← w + commMatrix[j][winners[k]]; end if w > threadResult[id].wMax then 17 threadResult[id].wMax ← w; 18 threadResult[id].winner ← j;

  Mixed-32 2 machines of type A and 1 of type B, 32 tasks in total Mixed-64 4 machines of type A and 2 of type B, 64 tasks in total

		different clusters with
	these machines:
	4×8	4 machines of type A, 32 tasks in total
	2×16	2 machines of type B, 32 tasks in total
	4×16	4 machines of type B, 64 tasks in total
	8×16	8 machines of type B, 128 tasks in total

Acknowledgment

This research received funding from the EU H2020 Programme and from MCTI/RNP-Brazil under the HPC4E project, grant agreement n. o 689772. It was also supported by Intel.