Alain Darte

Alexandre Isoard

Tomofumi Yuki

Extended Lattice-Based Memory Allocation

Keywords: Static optimization, polyhedral analysis, memory reuse and allocation, array contraction, integer lattices

This work extends lattice-based memory allocation, an earlier work on memory (array) reuse analysis. The main motivation is to handle in a better way the more general forms of specifications we see today, e.g., with loop tiling, pipelining, and other forms of parallelism available in explicitly parallel languages. Our extension has two complementary aspects. We show how to handle more general specifications where conflicting constraints (those that describe the array indices that cannot share the same location) are specified as a (non-convex) union of polyhedra. Unlike convex specifications, this also requires to be able to choose suitable directions (or basis) of array reuse. For that, we extend two dual approaches, previously proposed for a fixed basis, into optimization schemes to select suitable basis. Our final approach relies on a combination of the two, also revealing their links with, on one hand, the construction of multi-dimensional schedules for parallelism and tiling (but with a fundamental difference that we identify) and, on the other hand, the construction of universal reuse vectors (UOV), which was only used so far in a specific context, for schedule-independent mapping.

Technique étendue d'allocation mémoire basée sur les réseaux entiers

Résumé : Ce travail étend l'allocation mémoire basée sur les réseaux entiers précédemment proposée en analyse de réutilisation mémoire (de tableaux). La motivation principale est de traiter de meilleure façon les formes plus générales de spécifications rencontrées aujourd'hui, comportant du tuilage de boucles, du pipeline, et d'autres formes de parallélisme exprimées dans les langages à parallélisme explicite. Notre extension a deux aspects complémentaires. Nous montrons comment nous pouvons prendre en compte des spécifications plus générales où les contraintes de conflit (celles qui décrivent les indices de tableaux qui ne peuvent pas partager le même emplacement mémoire) sont spécifiées par une union (non-convexe) de polyèdres. Au contraire des spécifications convexes, ceci requiert d'être capable de choisir des directions (c'est-à-dire une base) adéquates de réutilisation des cases de tableaux. Pour cela, nous étendons deux approches duales, précédemment proposées pour une base fixée, en des schémas d'optimisation permettant de choisir des bases adaptées. Notre approche finale consiste en une combinaison des deux approches, révélant également des liens avec, d'une part, la construction d'ordonnancements multi-dimensionnels pour le parallélisme et le tuilage (avec une différence fondamentale que nous identifions) et, d'autre part, la construction de vecteurs de réutilisation universelle (UOV), qui étaient utilisés jusqu'à présent uniquement dans un contexte spécifique, celui des allocations valides pour tout ordonnancement.

Mots-clés : Optimisations statiques, analyse polyédrique, réutilisation et allocation mémoire, contraction de tableaux, réseaux entiers

Introduction

As the gap between memory performance and compute power keeps increasing, the importance of efficient memory usage is also increasing. This is even more emphasized when exploiting hierarchical memories and/or when accelerators such as GPUs or FPGAs are used as they are often limited by the on-chip memory capacity and/or the bandwidth between the host and the accelerator. The problem of efficient memory allocation is further complicated by the trade-off between parallelism and memory usage. Memory reuse is a standard technique for allocating scalar variables to registers. Memory reuse for arrays, in particular intra-array reuse, is used to reduce statically the memory footprint of data-intensive applications, after analyzing the liveness of the different elements of an array. The need for such array contraction is of course important for high-level program specifications (for example array languages) where the programmer expresses its applications in an abstract view of the storage locations, possibly even using arrays in single assignment, thus without paying too much attention to memory usage. But such a memory allocation technique is also required within compilers themselves as a complementary step to many code transformations, for the design of intermediate buffers introduced by the compiler and the management of local memories, and to reduce the effect of some previous array expansion phases.

In this paper, we extend a technique called lattice-based memory allocation [START_REF] Darte | Lattice-based memory allocation[END_REF] that was originally proposed as a generalization of different strategies based on affine mappings with foldings by modulo operations (called modular mappings), formalized with integer lattices. The original work was aimed at handling regular kernels executed by sequential and/or limited forms of parallel schedules where simple optimization strategies appeared to be sufficient. We extend this framework on two main aspects:

• conflict set (a relation to express array elements that may not be mapped to the same memory location) as union of polyhedra, where the initial work is limited to a convex polyhedron;

• optimized basis heuristics to choose the direction of the modular mapping allocation, or the basis of the corresponding lattice.

The key insight is that optimizing the dimensions of a multi-dimensional modular mapping successively, greedily minimizing the resulting array size for each dimension, may lead to worse overall allocations, a phenomenon that is exacerbated when the conflict set is not convex. By designing an optimization procedure to build short reuse vectors (that indicate array elements mapped to the same location), sharing similarities with UOV (universal occupancy vector) based approaches [START_REF] Strout | Schedule-independent storage mapping for loops[END_REF][START_REF] Thies | A step towards unifying schedule and storage optimization[END_REF][START_REF] Yuki | Memory allocations for tiled uniform dependence programs[END_REF], we address this issue and propose a combined heuristic to address this issue. Section 2 presents a simple example, in an informal way, to give the key intuition behind our approach and the reasons why previous approaches would fail to handle it efficiently. Section 3 gives some necessary background and formal definitions, defining more precisely the notions of conflict sets an modular mappings. Section 4 is the heart of our paper. It details the theory and algorithms we use to extend two previously-proposed dual approaches for memory mapping [START_REF] Darte | Lattice-based memory allocation[END_REF] to non-convex conflict sets and with optimized basis. It also shows how parametric solutions can be obtained and how the two heuristics can be combined. Section 5 illustrates on several examples our technique, showing how a library such as isl, for manipulating integer set relations, can be used to implement it. Finally, as our work shows multiple connections with a priori separated earlier work, we review in Section 6 some of them, detailing their links and differences with our technique. We conclude in Section 7.

Intuition of the Approach

In this section, we illustrate the key intuition behind our work using an example. There are many existing array mapping techniques that work well when the conflict set is described with a single polyhedron [START_REF] Somashekaracharya | Automatic intra-array storage optimization[END_REF][START_REF] Darte | Lattice-based memory allocation[END_REF][START_REF] Lefebvre | Automatic storage management for parallel programs[END_REF][START_REF] Quilleré | Optimizing memory usage in the polyhedral model[END_REF]. We are interested in more complex cases which involve unions of polyhedra. A common situation that gives rise to such a non convex domain is when loop tiling [START_REF] Xue | Loop Tiling for Parallelism[END_REF] is applied and we need to store the live-out values of a tile [START_REF] Somashekaracharya | Automatic intra-array storage optimization[END_REF]. Note that we present here only a simplified view of the problem to give the important intuition. The actual problem is finding an allocation, given a schedule or a set of possible schedules, which may include parallelism and tiling, and whose legality is formulated using a relation, computed thanks to liveness analysis and representing the possible overlaps between the live-ranges of array elements.

Let us consider a case where we seek an optimized allocation for a reverse-L shaped region depicted in Figure 1 (see also the figure caption for detailed explanations). The number of live-out values is 4N -4 if N is the length of the square edge (24 in the figure, with N = 7). One "good" affine allocation to maps all these values to different locations is to "project" along (1, 1) with an additional modulo factor of 2 (as depicted in Figure 1a), which corresponds to an array of size 2(2N -1), thus only 2 elements more than the "optimal", with the corresponding mapping (x, y) → (x -y mod (2N -1), y mod 2).

Existing techniques struggle to find this allocation for different reasons. One of the earliest, yet powerful, method for memory allocation by Lefebvre and Feautrier [START_REF] Lefebvre | Automatic storage management for parallel programs[END_REF] consists in choosing some successive modulo folding for each dimension, restricting to mappings along the canonical axes. The first modulo should be larger than the maximal distance between two points (here N -1 along the x axis), then the second modulo larger than the maximal distance between two points with the same value of x, which is also N -1. This results here in an allocation of size N 2 , with the corresponding mapping (x, y) → (x mod N, y mod N ).

Several other techniques have been presented that explore allocations using non-canonical projections or mappings [START_REF] Somashekaracharya | Automatic intra-array storage optimization[END_REF][START_REF] Darte | Lattice-based memory allocation[END_REF][START_REF] Quilleré | Optimizing memory usage in the polyhedral model[END_REF]. We leave the detailed discussions of these work to Section 6 and only give here a high-level description of the most recent work by Bhaskaracharya, Bondhugula, and Cohen [START_REF] Somashekaracharya | Automatic intra-array storage optimization[END_REF], which has the same objectives as ours. They proposed a method that combines ideas from multi-dimensional affine scheduling with memory allocation. The key idea is to interpret each dimension of a multi-dimensional affine mapping as a family of parallel hyperplanes that separate points from each other. If each parallel hyperplane contains at most one point of the domain to be stored, the corresponding affine function represents a legal memory allocation mapped to a one-dimensional array where the "latency" [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, I, one-dimensional time[END_REF] (width or maximal distance), plus 1, is the number of elements in the array. Using this formulation, their approach first tries to find such a family of hyperplanes that contain at most one point each. If such hyperplanes cannot be found, their technique finds a second family of hyperplanes focusing on points that lie in the same hyperplane of the first family, and so on. The dimensionality of the final array is the number of linearly-independent hyperplanes generated.

One possible allocation that can be found by their technique is illustrated in Figure 1b. Since the mapping they compute also depends on how the domain is decomposed into a union of polyhedra, their technique may find other (worse) allocations, as explained in Section 6.2. The key reason of this inefficiency lies in the objective functions used. Their primary objective tries to minimize the number of hyperplane families, i.e., the dimensionality of the mapped array, by maximizing at each stage the number of polyhedra whose points are all separated by hyperplanes. Their secondary objective is to minimize the number of elements in the corresponding dimension, i.e., to find hyperplanes yielding minimal width.

The primary objective can be satisfied by using the hyperplanes shown in Figure 1b, parallel to the vector (2, 3) (or similarly parallel to (3, 2)), thus The reverse-L shape of width 2 may be viewed as the live-out values of a square tile of size N , when the data dependences (not depicted here) across tiles have a uniform length of 2 in each dimension. The memory allocation characterized by the reuse vectors in Figure 1a leads to the optimal affine mapping (projection plus modulo 2) that we seek. Figure 1b illustrates another possible affine mapping (a projection) that may be found by other techniques, with N extra storage. orthogonal to (3, -2), with width 5(N -1). The corresponding mapping is (x, y) → 3x -2y mod (5N -4). For N = 7, this gives, as depicted, 31 different memory locations. However, if a one-dimensional allocation cannot be found, the greedy heuristic will minimize the number of elements in the current dimension, which, as explained in Section 6.2, may favor in this example hyperplanes parallel to the canonical axes (i.e., the shortest one-dimensional schedules), leading to the same mapping as Lefebvre and Feautrier, with size N 2 . The optimal solution is not obvious in this iterative formulation, since the "good" hyperplanes parallel to (1, 1) (i.e., orthogonal to (1, -1)) do not satisfy the primary objective. Indeed, some points are still mapped to a common location, and it is not optimal for this dimension because its corresponding width is 2(N -1) while N -1 is achievable.

In our work, we use a different formulation to overcome inefficiency in these cases. For our example, it can be intuitively explained as finding the shortest vector that points to somewhere outside of the domain of interest, which is the vector (2, 2). In general, we need multiple linearly-independent vectors, captured through lattices.

Background

We now introduce the two key concepts used in our paper, conflict set, which gives the constraints for valid mappings, and modular mappings, a particular form of functions used for intra-array reuse.

Conflict Set

Lattice-based memory allocation [START_REF] Darte | Lattice-based memory allocation[END_REF][START_REF] Darte | Lattice-based memory allocation[END_REF], as well as all prior work on intra-array reuse [START_REF] De Greef | Memory size reduction through storage order optimization for embedded parallel multimedia applications[END_REF][START_REF] Quilleré | Optimizing memory usage in the polyhedral model[END_REF][START_REF] Lefebvre | Automatic storage management for parallel programs[END_REF], is based on the concept of conflicting (array) elements, i.e., the set of pairs of elements that should not be mapped to the same location. It is expressed as a binary relation , which is the counterpart, for array elements, of the well-known interference graph defined for register allocation.

In register allocation, vertices correspond to scalar variables, and edges indicate that two variables should not be mapped to the same register, so that graph coloring can be used to derive a valid register assignment. For intra-array reuse, this interference graph (the relation) is not expressed in extension, but in a symbolic way, e.g., with polyhedra specifying a set of conflicting indices for a given array A. This specification is then used to derive a valid mapping σ, i.e., a mapping such that σ( i) = σ( j) if i j (meaning that A( i) and A( j) should not be mapped to the same location).

For intra-array reuse, the construction of the conflicting indices requires some symbolic liveness analysis. Such array mappings can be defined in a post-scheduling phase, i.e., for a particular execution or schedule (scheduledependent mappings [START_REF] De Greef | Memory size reduction through storage order optimization for embedded parallel multimedia applications[END_REF][START_REF] Lefebvre | Automatic storage management for parallel programs[END_REF][START_REF] Quilleré | Optimizing memory usage in the polyhedral model[END_REF]), or before the final schedule is defined (scheduleindependent mappings [START_REF] Strout | Schedule-independent storage mapping for loops[END_REF][START_REF] Thies | A step towards unifying schedule and storage optimization[END_REF][START_REF] Yuki | Memory allocations for tiled uniform dependence programs[END_REF]), so that the mapping is valid for any further code transformation (or a subclass, such as loop tiling). The latter situation arises also when compiling programs, expressed in a parallel language such as OpenMP or X10, on top of a runtime system, in which case the exact schedule is not statically known.

Actually, all these cases are particular instances of the more general problem of defining liveness analysis and "simultaneously live" array elements in explicitlyparallel specifications. How to build such a relation is not our concern here: we assume it has been computed, possibly over-approximated, and that it expresses, for each array to be contracted, the set of pairs of conflicting indices (or their differences) in a compact symbolic form. For optimization purposes, we will mainly discuss the case where the conflicting differences are the integer points in a union of polyhedra K = P 1 ∪ • • • ∪ P r . But the theory can be used in more general situations, e.g., with formulas in Presburger arithmetic as available in the Integer Set Library [START_REF] Verdoolaege | isl: An integer set library for the polyhedral model[END_REF], or even polynomial expressions [START_REF] Feautrier | The power of polynomials[END_REF], as long as the corresponding optimizations can be carried out.

Modular Mappings

In lattice-based memory allocation, the mappings are restricted, both for optimization and code generation purposes, to modular mappings. A modular mapping (M, b), defined by a p × n integral matrix M and a positive integral vector b of dimension p, maps the index i of a n-dimensional array to σ( i) = M i mod b (the modulo is applied component-wise) in a p-dimensional array of shape b. The size of the mapping is the size of the resulting array, i.e., the product of the elements of b. Its dimension is the dimension of the resulting array, i.e., p. As modular mappings are affine, it is enough to work with the symmetric set of conflicting differences K = {( ij) | i j}. The initial theory of lattice-based memory allocation [START_REF] Darte | Lattice-based memory allocation[END_REF][START_REF] Darte | Lattice-based memory allocation[END_REF] focuses on the case where K is the set of integer points in a polyhedron, in which case lower and upper bounds on memory size can be given and more properties proved.

Although we focus on intra-array reuse in this paper, it is worth to point out that modular mappings can be used in other contexts, e.g., for bank allocation, to allow parallel accesses to memory [START_REF] Cilardo | Improving multibank memory access parallelism with lattice-based partitioning[END_REF][START_REF] Darte | Lattice-based memory allocation[END_REF], or more generally whenever renewable resources need to be shared.

Greedy Mapping and Lattice Constructions

In this section, we present our main contribution that extends lattice-based memory allocation. Our work is based on two dual approaches proposed by Darte et al. [START_REF] Darte | Lattice-based memory allocation[END_REF] that are equivalent when the set of conflicting differences K is a polyhedron. We show how we can extend them to handle the case where K is a union of polyhedra, and then to further optimize the selection of the basis (i.e., the matrix M in the modular mapping) to reduce the memory size. The resulting optimizations lead to two complementary greedy approaches, in which the rows of the matrix M are optimized in the opposite order. We then show in Section 4.3 how to combine them to try to get the best of the two worlds, the first one being well suited to handle parameters while the second one is more suitable to detect directions with constant (i.e., non-parametric) reuse.

Both approaches define mappings given a basis of Z n , i.e., how to choose the modulo vector b. The first approach directly works with the matrix M of the mapping while the second one works with its kernel, i.e., the set of vectors i such that σ( i) = 0. Such a set is a lattice of Z n , thus the name of the technique. Section 4.1 extends the former over the mapping space, Section 4.2 extends the latter over the lattice space. We combine the two in Section 4.3. Unimodular matrices (invertible in the integers) play an important role in this construction, both for the optimizations in the two approaches, and for their combination. In this case, the first approach builds the rows of M , while the second builds the column of M -1 .

Basis Selection in Mapping Space

The following mechanism [START_REF] Darte | Lattice-based memory allocation[END_REF] is the "successive modulo" principle [START_REF] Lefebvre | Automatic storage management for parallel programs[END_REF], generalized to any set of n linearly independent integral vectors. K is the 0-symmetric set of conflicting differences.

Heuristic 1.

• Choose n linearly independent integral vectors ( c 1 , . . . , c n ).

• Compute F * i ( c i ) = sup{ c i . z | z ∈ K, c j . z = 0, ∀j < i}, successively for all 1 ≤ i ≤ n.
• Define M the matrix with row vectors ( c i ) 1≤i≤n and b an integer vector such that b i > F * i ( c i ). For each i, F * i ( c i ) is the width along c i of the intersection of K and the orthogonal of the vector space defined by ( c 1 , . . . , c i-1 ). Our goal is to adapt this heuristic when K is described with a union of polyhedra and to design an optimization to choose a suitable basis.

Theorem 1. The modular mapping built by Heuristic 1 is a valid mapping for K, assuming that K is bounded, symmetric w.r.t. 0.

Proof. Let x ∈ K with M x mod b = 0. Since b 1 > F * 1 ( c 1 ) = sup{ c 1 . x | x ∈ K} = sup{| c 1 . x| | x ∈ K} (as K is 0-symmetric), c 1 . x = 0 mod b 1 implies c 1 . x = 0. Then, c 2 . x = 0 mod b 2 , but b 2 > F * 2 ( c 2 ) = sup{| c 2 . x| | x ∈ K, c 1 . x = 0}
, thus c 2 . x = 0. Continuing this process, we get c i . x = 0 for all i. Since the c i are n linearly independent vectors, this implies x = 0. This shows that the modular mapping σ = (M, b) is valid. Indeed, if i j, then ij ∈ K and thus, unless i = j, σ( i) = σ( j).

Theorem 1 shows that, although it was initially designed assuming that K is a polyhedron [START_REF] Darte | Lattice-based memory allocation[END_REF], Heuristic 1 is actually valid with weaker hypotheses. Also, the following important properties remain true:

• The different values F * i depend on the order in which the vectors ( c i ) 1≤i≤n are considered. Considering all n! orders can help reducing the size of the mapping. (In practice, n is small.)

• Increasing the values of b keeps the validity of the mapping. This can be used for example to restrict the values of b to power of 2. Also, if K is a union of different 0-symmetric pieces, one can compute a vector b i for each piece, each obtained with a particular order of the basis, and then take the maximum, component-wise, of the b i to get a valid b for the whole K. • The same proof as for Theorem 1 shows that the modular mapping

c 1 . x + b 1 ( c 2 . x + b 2 (• • • + b n-1 c n . x)) mod i b i is valid, of dimension 1,
and with same size. In other words, the dimension of a mapping is not related to its size, and some 1D mappings can even be found with a multi-dimensional approach.

Optimizing the Basis

We now show how the previous heuristic can be extended to define an optimized version, where the basis ( c i ) 1≤i≤n is built in a greedy fashion, so that F * i ( c i ) is minimized at each step:

c i = argmin{F * i ( c) | c ∈ Z n linearly independent with c j , j < i}
This is a min-max problem, which can be solved, when K is a union of polyhedra, thanks to either the duality theorem of linear programming or the affine form of Farkas lemma [START_REF] Schrijver | Theory of Linear and Integer Programming[END_REF], in the same way schedules with minimal latencies (∼ width) are built [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part II: Multi-dimensional time[END_REF][START_REF] Darte | Linear scheduling is nearly optimal[END_REF]. Let us detail the technique with Farkas lemma, which we recall here.

Lemma (Farkas, affine form). Let P be a non-empty polyhedron defined by affine inequalities

P = { x | Q x ≤ e}. Then c.
x ≤ δ for all x ∈ P iff there exists y ≥ 0 such that c = y.Q and y. e ≤ δ. Now, the width of P along a vector c can be computed as follows:

max{ c. x | x ∈ P } = min{δ | c. x ≤ δ, ∀ x s.t. Q x ≤ e} = min{δ | y ≥ 0, c = y.Q, y. e ≤ δ}
Similarly, the quantity F * i ( c) is equal to:

F * i ( c) = min{δ | y ≥ 0, c = y.Q + z.C i-1 , y. e ≤ δ} (1)
RR n°8840

where C i-1 is the matrix whose rows are c 1 , . . . , c i-1 , and the vector z does not need to be nonnegative. Finally, for a union of polyhedra K = ∪ 1≤j≤r P j where P j = { x | Q j x ≤ e j }, one just need to collect the different constraints expressed in Equation 1, with a y j and a z j for each P j , plus the additional constraints y j . e j ≤ δ for all 1 ≤ j ≤ r. It remains to find c such that F * i ( c) is minimized, which amounts to solve a linear program, with objective function δ, where all these different variables, c, δ, y j and z j , are unknowns.

Linear Independence and Unimodularity

With no additional constraint, the optimal solution is of course c = 0. But c should be restricted so that it is linearly independent with all c j with j < i, and nonzero if i = 1. One way to do it is to complete (e.g., by computing the Hermite normal form), at each step, the vectors ( c j ) j<i , into a n-dimensional basis, with additional vectors ( d j ) j≥i , and to write c = j<i λ j c j + j≥i λ j d j . Then, it is sufficient to solve one linear program for each j ∈ [i . . . n] with the additional constraint λ j ≥ 1, and to define c i as the best of these (n -i + 1) solutions. There is no need to check for λ j ≤ -1 since c andc lead to similar mappings. Another trick is to select a random integer vector r and to add a single constraint c. r ≥ 1 (or a similar reasoning with the λ j ). Unless bad luck, this would be enough to not miss the optimal with a single linear program. Alternate solutions are possible, e.g., by constructing a basis of the orthogonal of the vector space defined by ( c j ) j<i or a pseudo-inverse (see also the discussion on Pluto scheduler [START_REF] Bondhugula | Automatic transformations for communication-minimized parallelization and locality optimization in the polyhedral model[END_REF][START_REF] Acharya | PLUTO+: Near-complete modeling of affine transformations for parallelism and locality[END_REF] in Section 6.1).

Also note that, since

F * i ( c) = F * i ( c -j<i λ j c j )
, one can restrict the search to vectors c = j≥i λ j d j , i.e., such that λ j = 0 if j < i. The resulting basis ( c i ) 1≤i≤n (i.e., the matrix M ) will then be automatically unimodular (i.e., with determinant ±1). This property is not formally needed to define a mapping but it makes the construction easier. In particular, completing at each step the vectors into a basis is just one iteration of the Hermite form computation. Also, imposing c to have integer components is then equivalent to looking for integer values for the λ j . To see this, suppose that, before Step i of the heuristic, we have built a unimodular matrix U whose first i -1 rows (matrix C i-1 ) are the vectors ( c j ) j<i built so far, and the remaining rows (matrix D i ) are the vectors ( d j ) j≥i . At Step i, the chosen solution c i is such that λ j = 0 for j < i and the λ j , for j ≥ i, have no common divisor d (otherwise a better solution can be defined by dividing by d). Thus, there is a unimodular matrix V i of size (n -i + 1) such that (λ i , . . . , λ n ) is the first row of V i . Then:

I i-1 0 0 V i C i-1 D i = C i-1 V i D i
is a unimodular matrix whose first i rows are the ( c j ) j≤i . This technique can thus be used to enforce a final unimodular matrix but also to complete the basis at each step, on the fly.

Handling parameters

In practice, the set of conflicting pairs and the set of conflicting differences are parameterized by structure parameters from the program, such as loop or array bounds. Heuristic 1 has the nice property that it can easily handle parameters as long as they constrain K in an affine way. If

P j = { x | Q j x ≤ E j n + e j }.
The standard technique [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part II: Multi-dimensional time[END_REF] is to bound the width as an affine function of the parameters ∆ n + δ and to apply Farkas lemma considering that n as a variable too, possibly taking into account additional affine constraints on the parameters. We then get a set of constraints as before, this time with the variables δ and ∆, with the objective of minimizing ∆ n + δ. Classically, parameters are then ordered with some priority so as to optimize the width in a lexicographic manner w.r.t. this order. The width is then parametric (thus the corresponding modulo), but the vector c built at each step has constant components thus all mechanisms presented before for basis completion, linear independence, and unimodularity remain true.

Wrap Up

The optimized version of Heuristic 1 presented here as some similarities with multi-dimensional scheduling and tiling, as explained in Section 6.1. It can help finding better mappings than without basis optimization, especially when K is not a polyhedron. However, unlike multi-dimensional scheduling for maximal parallelism detection for which a greedy approach is asymptotically optimal [START_REF] Darte | Optimal parallelism detection in nested loops[END_REF], selecting the smallest width at each step can be sub-optimal, even in order of magnitude, because it may be better to select a direction with larger width if the next width is then smaller in the orthogonal. This situation arises for the example of Section 2, where the heuristic will simply select the canonical basis with both widths equal to N , i.e., with no array reuse. See more detailed explanations in Figure 2 which depicts the set K of differences of the live-out points in Figure 1. A more complex optimized variant of Heuristic 1 has been proposed [START_REF] Somashekaracharya | Automatic intra-array storage optimization[END_REF] to try to avoid this caveat, but as explained in Section 6.2, it can also fail to find directions of reuse with a small (constant) modulo. The lattice-based heuristic proposed in the next section is aimed to avoid this pitfall.

Basis Selection in Lattice Space

A dual approach to derive a modular mapping (M, b) is to build an integer strictly admissible lattice for K. Given n linearly independent vector ( a i ) 1≤i≤n , the lattice generated by ( a i ) 1≤i≤n is the set

Λ = { x | x = A u, u ∈ Z n }
where A is the matrix with column vectors ( a i ) 1≤i≤n . It is strictly admissible for K if λ ∩ K = { 0}. A modular mapping is valid if and only if its kernel is an integer strictly admissible lattice for K. Also, from such a lattice Λ, one can build a valid modular mapping whose kernel is Λ and whose size is equal to the determinant of Λ, i.e., det(A). This is the underlying idea of lattice-based memory allocation [START_REF] Darte | Lattice-based memory allocation[END_REF][START_REF] Darte | Lattice-based memory allocation[END_REF].

The following scaling mechanism [START_REF] Darte | Lattice-based memory allocation[END_REF] gives a way to build a strictly admissible lattice from a set of n linearly independent integer vector ( a i ) 1≤i≤n . Again, we want to check if it still works for any K, not just for a convex polyhedron K as initially formulated, then derive mechanisms to optimize the basis it works with.

Heuristic 2.

• Choose n linearly independent integral vectors ( a 1 , . . . , a n ). Figure 2: Optimized version of Heuristic 1 for the set K of differences of the live-out points in Figure 1. The selected basis is c 1 = (0, 1), then c 2 = (1, 0) (or the converse). The left figure shows the width F * 1 ( c 1 ) = N -1 computed for the first vector. In the next step of the heuristic, K is restricted to the orthogonal space of the first vector, and F * 2 ( c 2 ) = N -1 is computed as illustrated on the right. The resulting mapping is (x, y) → (y mod N, x mod N ), which is not a good allocation for this example. This is because the basis vectors are greedily selected to minimize the width at each dimension, missing the "hole" along (1, 1).

• Compute F i ( a i ) = inf{λ ≥ 0 | a i ∈ λK i }, successively for all 1 ≤ i ≤ n,
where K i = K + Vect( a 1 , . . . , a i-1 ).

• Define the lattice Λ generated by the vectors (ρ i a i ) 1≤i≤n where the ρ i are integers such that ρ i > 1/F i ( a i ).

For each i, µ a i with µ = 1/F i ( a i ) is the largest "multiple" of a i that belongs (if K is closed) to the extrusion K i of K along the vectors ( a j ) j<i , i.e., K i = K + Vect( a 1 , . . . , a i-1 ) = { x | x = y + j<i α j a j , y ∈ K}. In other words, µ a i -j<i α j a j ∈ K for some real numbers (α j ) j<i and, for any ρ > µ, ρ a i / ∈ K i .

Theorem 2. The lattice Λ built by Heuristic 2 is a strictly admissible integral lattice for K, if K is bounded, symmetric w.r.t. 0.

Proof. Let x ∈ Λ. If x = 0, one can write x = i j=1 u j ρ j a j for some integers (u j ) j≤i with u i = 0. Suppose x ∈ K. If K is 0-symmetric, one can assume u i ≥ 1 without loss of generality. Then a i = ( x-j<i a j )/(ρ i u i ), thus F i ( a i ) ≤ 1/(ρ i u i ) by definition of F i ( a i ), then F i ( a i ) ≤ 1/ρ i ,
which is impossible by definition of ρ i . Thus, Λ is a strictly admissible integral lattice for K.

Optimizing the Basis

We now assume that K is star-shaped, i.e., if x ∈ K, then λ x ∈ K for all 0 ≤ λ ≤ 1. In this case, ρ i is the smallest integer such that ρ i a i / ∈ K i . Theorem 2 shows that Heuristic 2, as for Heuristic 1, is still valid even if K is not a polyhedron. But how can we optimize the basis ( a i ) 1≤i≤n ? One could think that minimizing 1/F i ( a i ) is the right thing to do, i.e., that it is sufficient to pick RR n°8840 any vector outside the ith extrusion K i of K along the previously-built vectors (in which case ρ i = 1 can be chosen). This is of course not true: unlike for Heuristic 1, the size of the resulting mapping cannot be read directly from the ρ i , it depends also on the determinant of the a i . It seems thus more profitable to directly optimize the basis vectors l i = ρ i a i of Λ, e.g., by minimizing their norm, for some adequate norm to be defined. This is because the determinant, in absolute value, is bounded by the product of the Euclidian norms (and all norms have the same order of magnitude). Also, if the direction a i is chosen, the best solution is to choose l i to be the smallest integral vector out of the extrusion, in the direction of a i .

Another indication of why looking for small l i is more likely to be good is that each l i is a reuse vector (or occupancy vector in the UOV terminology [START_REF] Strout | Schedule-independent storage mapping for loops[END_REF]), i.e., A( x), A( x + l i ), . . . , A( x + k l i ) will reuse the same memory location. The mapping will exploit more reuse if a larger number of copies of l i can traverse the original array space thus, intuitively, a smaller vector will lead to more reuse. However, as for the UOV optimization [START_REF] Strout | Schedule-independent storage mapping for loops[END_REF], the best direction remains difficult to anticipate: it depends on the extent of the next extrusion K i+1 , i.e., of the other a j not yet defined. Nevertheless, despite this inaccuracy, minimizing the norm is interesting because it will select, in priority, non-parametric reuse vectors, thus favor the reduction of the mapping size in order of magnitude (if there is one large parameter N and p constant reuse vectors are found, the mapping size will be of order N n-p ). For all these reasons, at each step of the heuristic, we will look for an integral vector with minimal norm that is out of the extrusion. For computation reasons, we will choose a norm that can be minimized with linear programming, for example ||.|| ∞ (max of the absolute value of components) or ||.|| 1 (sum of the absolute value of components, i.e., Manhattan distance). To break ties, we can also use the norm proposed for the computation of AUOV [START_REF] Thies | A step towards unifying schedule and storage optimization[END_REF], which is a two-dimensional lexicographic optimization, first the Manhattan distance, then the sum of the absolute value of all differences of two components, which tends to lead to more "diagonal" vectors.

Linear Independence and Unimodularity

Unlike Heuristic 1 for which we perform an optimization based on an expression of K, here we want to find a short vector not in K i , or at least on the border of K i . We tried many different optimization schemes, using Farkas lemma or the duality theorem, but we found no better solution than working directly with the complement of K i , expressed as a union of polyhedra. Then, we just need to minimize the norm in each piece, with integer linear programming, and to pick the best solution l i . As any linear combination of ( a j ) j<i belongs to K i , linear independence is automatically satisfied. Now, let us see if we can enforce some unimodularity property and how it can be used to simplify the computations. Suppose the vectors ( a j ) j<i have been computed so that they can be completed into a unimodular matrix, thanks to additional vectors ( d j ) j≥i . We can look for the vector l i expressed in this basis: l i = j<i λ j a j + j≥i λ j d j . Then, we can either minimize the norm in the original basis (minimizing the components of l i ) or in this new basis (minimizing the λ j ). In both cases, since l i / ∈ K i then, by definition of K i , the same is true if we subtract to it any linear combination of ( a j ) j<i . This does not change F i ( a i ), the subsequent extrusions, and the determinant of the lattice. This amounts to assume that λ j = 0 for j < i. Then, if ρ i is the common divisor of the (λ j ) j≥i , we can select a i = ( j≥i λ j d j )/ρ i so that l i = ρ i a i . Finally, as we did for enforcing unimodularity in Heuristic 1, we can then complete ( a j ) j≤i into a unimodular matrix, with one simple computation of the Hermite form. The basis of the final lattice Λ is then given as a unimodular matrix (A) times a diagonal one (the ρ i ).

Things are simpler when unimodularity is enforced. A valid mapping can be easily defined by σ = (A -1 , ρ). Indeed, σ( x) = 0 iff A x = 0 mod ρ iff there exists y ∈ Z n such that A x = R y where R = diag(ρ 1 , . . . , ρ n ), i.e., x = A -1 R y, which means x ∈ Λ. When A is unimodular, it is also simpler to build the set of integral vectors in K i when K is, as in the isl library, expressed by a formula involving only integer variables and affine inequalities. Indeed, in this case,

K i = K ∩ { x | x = j<i α j a j , α j ∈ R} = K ∩ { x | x = j<i
α j a j , α j ∈ Z} is also defined with integer variables and affine inequalities, so expressible with isl. The set of integral vectors in its complement can then be computed by isl. It will also be described with integer variables and affine inequalities.

Handling Parameters

As we saw in Section 4.1, Heuristic 1 can be extended with no difficulties to the handle the case where K depends affinely on parameters. The final mapping may depend on these parameters, but not the matrix M , only the modulo vector b. This is more complicated for Heuristic 2. If all vectors ( l i ) 1≤i≤n are constant or equal to λ a i where only λ depends on the parameters, the previous construction can be performed and a valid mapping can be computed from the lattice. However, if this is not the case, i.e., if the direction of l i depends on the parameters, then computing K i+1 involves quadratic constraints (α i is multiplied by a parameter in the definition of K i+1 ), which is problematic. Similarly, defining a i and completing the vectors ( a j ) j≤i into a basis will need knowledge on the arithmetic properties (gcd) of the parametric components of l i .

An easy situation is when all vectors a i built during the optimization process are constant, except possibly the last one. Indeed, for i = n, there is no K i+1 to build and no basis completion to perform. This is the case for the example of Section 2, where we find successively l 1 = ρ 1 a 1 = (2, 2) T , then l 2 = ρ 2 a 2 = (2N -1, 0) T , which corresponds to the 2D mapping (x, y) → (y mod 2, xy mod (2N -1)). See detailed explanations in Figure 3. By construction, the 1D mapping σ (x, y) = xy + (2N -1)y mod 2(2N -1) is also valid [START_REF] Darte | Lattice-based memory allocation[END_REF].

Finally, let us point out that as long as, in the process, we seek a constant reuse vector l i , then l i belongs to K i for any value of the parameters, i.e., the intersection of all K i when the parameters vary, i.e., the complement of the union of all K i when the parameters vary. This means that we can first project out the parameters in K (equivalent to the union when the parameters vary), then compute the successive extrusions and complement sets. In the example of Section 2, projecting out the parameter N in K gives the non-parametric complement set K defined by {(x, y) | (x ≥ 2, y ≥ 2) or (x ≤ -2, y ≤ -2)}, from which we easily find l i = (2, 2). In general, we can thus first build a non-parametric set K, by projecting out the parameters, and work with it as long as we find a constant reuse vector. 1. The selected basis is a 1 = (1, 1) T , then a 2 = (1, 0) T , or equivalently (0, 1) T . The figure to the left shows the basis vector of the lattice l 1 = ρ 1 a 1 = (2, 2) T , which is the shortest vector (a multiple of a 1 ) that points outside of K. In the next step, the set K is extruded along a 1 and is now the infinite diagonal band shown in the right figure, leading to the second basis of the lattice: l 2 = ρ 2 a 2 = (2N -1, 0) T . The final mapping is given by the inverse of the matrix [ a 1 a 2 ] -1 = 0 1 1 -1 and the modulo factors defined by ρ 1 = 2 and ρ 2 = 2N -1, which is the mapping (x, y) → (y mod 2, xy mod (2N -1)).

Star Shaping

There remains one potential problem, if K is not naturally star-shaped. In this case, K i can have "holes" along a direction and the optimization of the norm in K i can lead to a small integral reuse vector l i ∈ K i such that λ l i ∈ K i for some positive integer λ, resulting in an invalid lattice. A possibility is to ignore this problem and to check, a posteriori, that this does not happen. This can be done with integer linear programming if l i is a constant vector. Another safer strategy is to first modify K into its star-shaped extension K * = { x | x = λ y, 0 ≤ λ ≤ 1, y ∈ K}. This however has to be done with care to avoid adding artificially integral points. Indeed, suppose that K = {(0, 0)} ∪ {(x, y) | 2 ≤ |y| ≤ 3}. The star-shaped extension of K is an open set {(x, y) | 0 < |y| ≤ 3} ∪ {(0, 0)}, while projecting naively λ in the expression of K * would add the line y = 0. But, as we are interested only in integral vectors, we want to work with K * = {(0, 0)} ∪ {(x, y) | 1 ≤ |y| ≤ 3}. This can be done as follows. First suppose that K is the integral points in a polyhedron P = { x | A x ≤ b}. We consider the set of integral points in K * = { x | ∃λ, 0 < λ ≤ 1, A x ≤ λ b} to which will be then added the vector 0. We then eliminate λ using the Fourier-Motzkin method, which leads to the following constraints:

K * = { 0} ∪              x (A x) i ≤ b i if b i > 0 (A x) i ≤ 0 if b i = 0 (A x) i < 0 if b i < 0 (A x) i b i ≤ (A x) j b j if b j < 0 < b i              RR n°8840
The inequality (A x) i < 0 is then modified in (A x) i ≤ -1 since A and x have integral components. As the star-shaped extension of a union of polyhedra is the union of the star-shaped extension of each polyhedron, the previous technique gives an algorithm to star-shape any union of polyhedra. Note however that, as for the convex hull, it does not work, in general, for a parametric set as the result is not always affinely parametric. But, as we use this procedure to find constant reuse vectors, we can first project out the parameters from K, then build the star-shaped extension K * of this non-parametric K as we just explained. To illustrate this principle, suppose that

K = {(0, 0)} ∪ {(x, y) ∈ Z 2 | 2 ≤ |y| ≤ 3, |x| ≤ N }. Eliminating N produces {(0, 0)} ∪ {(x, y) | 2 ≤ |y| ≤ 3}, i.e.
, the set we discussed previously. Finally, with Fourier-Motzkin, we get the star-shaped extension {(0, 0)} ∪ {(x, y) ∈ Z 2 | 0 < |y| ≤ 3}, and, restricting to integer points,

K * = {(0, 0)} ∪ {(x, y) ∈ Z 2 | 1 ≤ y ≤ 3} as expected.
We can then compute the set of integer points not in

K * as {(x, y) | |y| ≥ 4} ∪ {(x, y) | |x| ≥ 1, y = 0}.

Wrap Up

The optimized version of Heuristic 2 presented here can be viewed as a generalization of all approaches based on reuse/occupancy vectors: UOV [START_REF] Strout | Schedule-independent storage mapping for loops[END_REF], QUOV [START_REF] Yuki | Memory allocations for tiled uniform dependence programs[END_REF], AUOV [START_REF] Thies | A step towards unifying schedule and storage optimization[END_REF], and even (pseudo)-projection methods [START_REF] Quilleré | Optimizing memory usage in the polyhedral model[END_REF]. See more details in Section 6.3. It is a strict generalization in the sense that our technique is the first multi-dimensional reuse vector technique, thanks to the concept of extrusion.

Unlike Heuristic 1, it manipulates K, the complement of K, and not K itself. This duality was also identified in the context of liveness analysis. In some specialized contexts (as in all previous work on reuse vectors) K, or an over-approximation of it, can be built directly, and not as the complement of K.

In Figure 4, we describe how the different optimization criteria for the two heuristics favors one set of basis over the other. It is important to emphasize that the vectors for the two approaches have different meanings. One is the direction of the mapping, and the other is the reuse vector. Providing the equivalent mapping vectors that would be obtained from lattice-based approach, in the reverse order, i.e., c 1 = (1, -1) and c 2 = (0, 1), to the Heuristic 1 yields the same modulos, but computed as widths, thus the same mapping.

Currently, our technique can be fully implemented only when one can guarantee that the reuse vectors produced are constant, i.e., do not depend on parameters. In contrast, Heuristic 1 is naturally adapted to parametric optimization. This motivates the development of a combination of the two heuristics as exposed in Section 4.3: we will first look for constant reuse vectors with Heuristic 2, then complete the lattice/mapping by optimizing the dimensions in the opposite order thanks to Heuristic 1, constrained by the reuse vectors already found by Heuristic 2.

Note that when p constant vectors are found, one could re-optimize in this vector space, thanks to a suitable enumeration of lattices using the Hermite normal form [START_REF] Darte | Lattice-based memory allocation[END_REF], to find the best one in this space. However, such a search is potentially expensive if the determinant is large. Also, the problem of finding the best constant strictly admissible sub-lattice in general seems open because we do not know in which subspace the search needs to be done. favors, for the first vector, the canonical axis (1, 0) (or equivalently (0, 1)) over (1, -1) because the width along (1, -1) is wider. As shown in the figure to the left, the width is defined by the vector (N -1, 1 -N ) with a width equal to 2N -2, which is larger than N -1 obtained from the canonical basis. In the lattice space, i.e., for Heuristic 2, the canonical basis leads to a much longer reuse vector (it has to be multiplied by ρ = N ) as depicted in the right figure. Since a shorter (and constant) reuse vector is preferred in the basis selection, (2, 2) is selected.

Combining the Two Approaches

The two previous approaches can be combined by understanding the connection between a valid mapping and an admissible lattice, in particular when the lattice has a "scaled" unimodular basis, i.e., is given by a unimodular basis A multiplied by a diagonal matrix R = diag(ρ 1 , . . . , ρ n ). If ( a i ) 1≤i≤n are the columns of A, then the l i = ρ i a i form a basis of reuse vectors for the lattice. If C = A -1 is the inverse of A, with row vectors ( c i ) 1≤i≤n , then σ defined by σ( x) = C x mod b (following the notations of Heuristic 1) with b i = ρ i is a valid mapping. Furthermore, if K is the set of integral points in a union of polyhedra, then if the ρ i were computed, from i = 1 to n, following Heuristic 2 for the basis ( a i ) 1≤i≤n , they are also equal to the b i computed, from n to 1, following Heuristic 1 for the basis ( c i ) 1≤i≤n . This can be proved the same way as for a single polyhedron [START_REF] Darte | Lattice-based memory allocation[END_REF], thanks to duality in linear programming. Now, suppose that only the first i column vectors of A have been computed, for example because we were not able to find a constant vector l i+1 in the complement of K i+1 . How can we complete them with (n -i) parametric reuse vectors? These first i vectors do not fully constrain the final mapping, they only define a partial mapping with the first i vectors ( c j ) j≤i of the inverse of A. However, they fully determine the vector space they generate, as well as its orthogonal, i.e., the vector space where the missing vectors ( c j ) j>i should lie. The idea is then to use Heuristic 1, constrained to this orthogonal, to optimize them and complete the mapping. This gives rise to the following combined heuristic.

Heuristic 3.

• Project out the parameters in K to get a description K of the union of constraints valid for all parameters.

• Build the star-shaped extension K * of this non-parametric K as explained in Section 4.2.

• Use Heuristic 2 with K * to build optimized reuse vectors, until no reuse vector is found. If i is the number of vectors built, we get a unimodular matrix A with column vectors ( a j ) 1≤j≤n such that for all j ≤ i, l j = ρ j a j is the reuse vector that was built.

• Use Heuristic 1 with the initial set K, restricting the search to the orthogonal defined by ( a j ) j≤i , to get (n -i) optimized mapping vectors ( c j ) n≥j>i (numbered in this decreasing order). Let (w j ) n≥j>i be their corresponding successive widths.

• Define the matrix M with row vectors ( m j ) 1≤j≤n such that m j is the jth row of A -1 , if j ≤ i, and

m j = c j if j > i. Define b with b j = ρ j if j ≤ i and b j = w j if j > i.
Theorem 3. The mapping σ = (M, b) is a valid mapping for K.

Proof. Let x ∈ K such that M x = 0 mod b. With the same argument used in the proof of Theorem 1, and by definition of ( m j ) j>i in Heuristic 3 and of (b j ) j>i in Heuristic 1, we first get c n . x = 0, then successively c j . x = 0 for all j from n to (i + 1). Thus x belongs to the orthogonal of ( c j ) j>i , i.e., the vector space spanned by ( a j ) j≤i . As A is unimodular, we can write x = j≤i λ j a j with integers λ j . As the first i rows of M are those of the inverse of A, m j . x = λ j for j ≤ i. Thus λ j is a multiple of b j = ρ j as defined in Heuristic 2. With the same arguments used in the proof of Theorem 2, we conclude that x = 0, which means that the modular mapping defined by M and b is valid.

To force the search in the orthogonal of the ( a j ) j≤i , we can look for an integer linear combination of the last (n -i) rows of A -1 only. Then, using the same principles as for Heuristic 1, we can impose that the vectors ( c j ) j>i form a unimodular basis of this orthogonal, which, combined with the first i rows of A -1 , will form a n × n unimodular matrix M . There is thus in this case a complete correspondence between the b i found through Heuristic 1 for the row vectors of M , starting from the last one, and the ρ i found through Heuristic 2 for the column vectors of M -1 starting from the first one. But what is important here is that by enforcing the search of the mapping vectors ( c j ) j>i to the orthogonal of small reuse vectors ( a j ) j≤i , we find potentially larger widths for the first ( c j ) j>i . But this is on purpose, to make sure we do not miss a good directions of reuse later. This is what happens on the example of Section 2 if we apply the combined heuristic. We first find the reuse vector (2, 2), then we limit the search to the orthogonal direction. Here, in 2D, we directly get the mapping vector (1, -1). In general, this remains a heuristic as the orthogonal of the reuse vector may not be the right direction in terms of width, although it guarantees the right order of magnitude for the size of the mapping. This concludes the formal description of our combined optimization.

Evaluation

We validated our technique by running the different optimization steps, as a script, with the iscc calculator [START_REF] Verdoolaege | Counting affine calculator and applications[END_REF]. As iscc does not provide sufficient genericity, the script has to be tailored to each example but a generic implementation could be done using isl [START_REF] Verdoolaege | isl: An integer set library for the polyhedral model[END_REF].

Reverse-L Shaped Region and Optimizing Scripts

Our illustrating example can be solved using the following script: Kp represents K, with the exception that its parameters are variable (to circumvent iscc's impossibility of projecting parameters out). We then build K_1 by projecting those parameters. A star-shaping step should also take place (but again is not easy to do with iscc) but would give the same polyhedron as it is already star-shaped (and 0 is also in each basic set). We then compute nK_1 as the complement of K_1. The Norm map provides the ||.|| ∞ norm, which we want to minimize (the ||.|| 1 norm would give the same result here as the optimal solution is a vertex of nK_1). The inversion of the x and y coordinates is cosmetic and is so that we choose positive coefficients for equivalent vectors as a side effect of minimization. We then proceed to find the smallest vector, which, in our example, leads to l 1 = (2, 2) T as expected, providing our first lattice vector a 1 = (1, 1) T and its associated modulo 2. The search for the next vector will be unfructuous but is executed as follows:

Kp := {[N]->[x,
Extr_1 := {rat: [x,y]->[x',y']:
(exists e: x'=x+1*e and y'=y+1*e)};

Unimod_1 := range {[t]->[1*t,0*t]: t>0}; K_2 := coalesce Extr_1(K_1); nK_2 := {[x,y]}-K_2; a_2 := (Norm^-1)(lexmin (Norm(nK_2*Unimod_1)));
where Extr_1 is an extrusion operator along a 1 (notice the use of a rational map) and Unimod_1 is a constraint enforcing the unimodularity of the future basis. It is obtained by looking for vectors colinear to the vectors (here only one, (1, 0) T ) that complete a 1 into a unimodular basis. Here, K_2 is now the whole space, due to the parameter projection, so nK_2 is empty, which shows that we cannot find a second constant reuse vector to complete the first. If we did not project the parameters (and detected that K_2 was already star-shaped) then we would have found l 2 = (2N -1, 0), which is a valid parametric reuse vector (see again Figure 3). Now, let us assume that we did not find that last vector because we stopped Heuristic 2 as soon as it does not find a constant vector. We complete the vector a 1 into a unimodular basis as previously, for example with the vector RR n°8840

(1, 0), then we compute its inverse:

1 1 1 0

-1 = 0 1 1 -1
. This already gives the first vector c 1 = (0, 1) of the mapping (first row) and its associated modulo equal to 2. We then continue with Heuristic 1: The set K is the parameterized Kp and fK is the dual of K in the sense of Farkas lemma. Opt is used to order the coefficients in the correct order for the minimization: we want a small coefficient for N , then a small constant coefficient.

K := range ([n]->{[n]->[x,y]} * Kp); fK := (unwrap coefficients K)^-1; Opt := {[c_cst,c_n]->[c_n,c_cst]}; Ortho_1 := range {[t]->[-1*t,
Ortho_1 constraints the vector to be in the vector space that completes c 0 into a unimodular matrix (so here it is the second row of the inverse we computed). b combines those operators and is therefore a map that gives, for a given vector, the width of the symmetric conflict set in this direction. For example, lexmin (fK.Opt)([0,1]) would gives [1,0] which represents 1 * N + 0 in the Farkas space. We then compute the smallest width into b_min, and proceed to find a vector c 2 that attains it. The smallest modulo here is 2N -1 (as b_min = [2, -2] → 2N -2 and we need a modulo strictly larger).

Blur Filter

On the blur filter with an interleaved schedule [START_REF] Somashekaracharya | Automatic intra-array storage optimization[END_REF] we find the optimal allocation at the condition that we optimize both for the Maximum distance ||.|| ∞ and the Manhattan distance ||.|| 1 , as explained in Section 4.2.1. Indeed, as shown in Figure 5, (1, 2) is equivalent to (2, 2) for ||.|| ∞ and to (0, 3) for ||.|| 1 . With this optimization, we get the desired allocation of blurx[(y -2x) mod (2N + 1)].

x y Figure 5: Blur filter w. interleaved schedule and conflict differences

N →{ (x, y) | (x ≥ 1 ∧ x < N ∧ y ≤ 1 ∧ y ≥ 0) ∨ (x ≥ 0 ∧ x < N ∧ y ≤ -1 ∧ y ≥ -2 ∧ y > -N ) ∨ (x ≤ 0 ∧ x > -N ∧ y ≥ x ∧ y ≤ 2 ∧ y ≥ -1 ∧ y < N ) }

Diamond Tiling

Another interesting example is the diamond tiling example [START_REF] Somashekaracharya | Automatic intra-array storage optimization[END_REF] whose conflict set is shown in Figure 6. The technique of Bhaskaracharya et al., analyzed in Section 6.2, finds an allocation size 6B -5 with the mapping allocation A B [(t -3i) mod (6B -5)], where B is the tile size. Our algorithm finds the mapping A B [t mod 2, i mod (2B -1)], a mapping aligned with the canonical basis that Lefebvre and Feautrier [START_REF] Lefebvre | Automatic storage management for parallel programs[END_REF] thus also find (or equivalently Heuristic 1) if i is the first canonical dimension, and t the second. In our case, our combined heuristic does not fall in the trap: the first lattice vector found is (0, 2), i.e., a 1 = (0, 1) and a modulo of 2, then the second vector is a mapping vector (1, 0) with a modulo of 2B -1, so the canonical basis in the right order is found.

B →{ (i, t) | (i ≤ 0 ∧ t ≤ 1 -i ∧ t ≥ -1 + i ∧ t ≤ 2B + i ∧ t ≥ -2B -i) ∨ (i ≥ 0 ∧ t ≤ 1 + i ∧ t ≥ -1 -i ∧ t ≤ 2B -i ∧ t ≥ -2B + i) ∨ (i ≤ B ∧ i ≥ -B ∧ t ≤ 1 -i ∧ t ≥ -1 -i ∧ t ≤ B ∧ t ≥ -B) ∨ (i ≤ B ∧ i ≥ -B ∧ t ≥ -1 + i ∧ t ≤ 1 + i ∧ t ≤ B ∧ t ≥ -B) }

Other Examples

To compare our technique with the state of the art, we report in Table 1 the different mappings found on a set of examples previously introduced to evaluate the technique of Bhaskaracharya et al. [START_REF] Somashekaracharya | Automatic intra-array storage optimization[END_REF]. We provide the mappings found by the successive modulo (for the canonical basis) of Lefebvre-Feautrier [START_REF] Lefebvre | Automatic storage management for parallel programs[END_REF], by their technique, and by ours. The last column gives, in order of magnitude in terms of the parameters, the factor of improvement (the larger, the better), as compared to Lefebvre-Feautrier (thus it is 1 for the first row for each example). In all cases, we find the same or a better factor. Note that for heat-2d-tile, we took the tiling along the mapping (t, i, j) → (t, t + i, t + i + j) but this may not be the same.

Related Work

There are a number of existing techniques for memory allocation in the polyhedral model, i.e., for programs on which code analysis and optimizations based on manipulations of polyhedra and linear programming techniques can be applied [START_REF] Somashekaracharya | Automatic intra-array storage optimization[END_REF][START_REF] Darte | Lattice-based memory allocation[END_REF][START_REF] Lefebvre | Automatic storage management for parallel programs[END_REF][START_REF] Quilleré | Optimizing memory usage in the polyhedral model[END_REF]. This work extends the work by Darte et al. on lattice-based memory allocation [START_REF] Darte | Lattice-based memory allocation[END_REF][START_REF] Darte | Lattice-based memory allocation[END_REF]. As developed in Section 4, we extend this framework both to handle (non-convex) unions of polyhedra and to improve the choice of projections, i.e., mapping functions. The work by Lefebvre and Feautrier [START_REF] Lefebvre | Automatic storage management for parallel programs[END_REF] can be viewed as a special case of lattice-based allocation (more precisely in the form of Heuristic

LBM-D3Q19

LeFe

A[t mod 2, i mod N, j mod N, k mod N ] 1 BBA A[(i -2t) mod (N + 2), j mod N, k mod N ] 2 Lattice A[j mod N, (i -t) mod (N + 1), (k -t) mod (N + 1)] 2 LBM-D3Q27 LeFe A[t mod 2, i mod N, j mod N, k mod N ] 1 BBA A[(i -2t) mod (N + 2), j mod N, k mod N ] 2 Lattice A[(j -t) mod (N + 1), (i -t) mod (N + 1), (k -t) mod (N + 1)] 2 diamond-tile LeFe A[t mod B, i mod (2B -1)] 1 BBA A[(t -3i) mod (6B -5)] B/3 Lattice A[t mod 2, i mod (2B -1)] B/2 heat-2d-tiled LeFe A B [t mod B, i mod B, j mod B] 1 BBA A B [(i -2t) mod (4B -3), (t -2j) mod (5B -4)] B/20 Lattice A B [j mod (2B -1), i -t mod (3B -2)] B/6
i.e., in the space of mappings), where only a subset of the mapping space is used (and not optimized), i.e., it works for a fixed basis. The technique proposed by Quilleré and Rajopadhye [START_REF] Quilleré | Optimizing memory usage in the polyhedral model[END_REF], on the other hand, is more similar to Heuristic 2, in the sense that it explores the lattice space, characterizing legal projective allocations (possibly with modulo reuse), including projections along non-canonical directions. However, their primary objective is in minimizing the dimensionality of the resulting projection, and no algorithm is available to search among legal projections. Also, it can be used only with strong hypotheses, in particular for multi-dimensional schedules. Nevertheless, in this limited context, it does find constant reuse vectors when the conflict set is "flat", i.e., not fully-dimensional. Our optimized technique with successive extrusions generalizes this approach in a broader context.

Our technique has strong links with the search for multi-dimensional affine functions for scheduling, in particular for detecting tiling bands in nested loops (see Section 6.1). Our optimized version of Heuristic 1 uses similar linear programming techniques, as does the recent work of Bhaskaracharya et al. for intra-array reuse [START_REF] Somashekaracharya | Automatic intra-array storage optimization[END_REF]. This latter work has the same objectives as ours, but it uses a different approach to tweak Heuristic 1 and try to avoid a N 2 mapping in the main example of Section 2. We explain in Section 6.2 why this approach still has some weaknesses, at least for intra-array reuse. The way we tweak Heuristic 1, with the help of Heuristic 2, reveals new connections with early work on universal occupancy (reuse) vectors as we explain in Section 6.3.

Link with Multi-Dimensional Scheduling and Tiling

It is interesting to note the strong link with the algorithm used in the Pluto compiler [START_REF] Bondhugula | Automatic transformations for communication-minimized parallelization and locality optimization in the polyhedral model[END_REF][START_REF] Acharya | PLUTO+: Near-complete modeling of affine transformations for parallelism and locality[END_REF] to build code transformations enabling tiling. In Pluto, operations in nested loops, each captured by a textual statement S and a loop counter (or iteration) vector i, are reordered by defining affine mappings σ S ( i) such that σ T ( j) -σ S ( j) ≥ 0 whenever there is a dependence from (S, i) to (T, j). The objective function is to minimize the distance between these operations (this tends to improve data reuse), i.e., the maximum of all σ T ( j) -σ S ( i), as we do for the width computation in Heuristic 1.

However, there are three main differences compared to the scheduling algorithm. The first difference is that we assume a single σ, i.e., σ T = σ S , so that we can work with the difference ji. But we could also apply Heuristic 1 to map different arrays in the same memory space, each with a possibly different mappings, as explored in SMO [START_REF] Bhaskaracharya | SMO: An integrated approach to intra-array and inter-array storage optimization[END_REF]. The second difference is that, at each step, we remove all pairs such that σ( j) -σ( i) = 0 (as in the search for maximal parallelism [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part II: Multi-dimensional time[END_REF]) while, in Pluto, dependences need to be kept for defining other dimensions for tiling. The third difference is that there are no constraints such as σ( j) -σ( i) ≥ 0 as conflicting differences can have any sign. In other words, the technique is similar but we will get a smaller width at each step (because of fewer constraints).

We can import all tricks used in Pluto [START_REF] Bondhugula | Automatic transformations for communication-minimized parallelization and locality optimization in the polyhedral model[END_REF] and Pluto+ [START_REF] Acharya | PLUTO+: Near-complete modeling of affine transformations for parallelism and locality[END_REF], in particular for enforcing linear independence but, as our situation is simpler, the techniques we use are sufficiently cheap in our case. Nevertheless, this shows a strong link between multi-dimensional scheduling/tiling and array mapping: dependences are used for the first, pairs of conflicting elements for the latter, the main difference is that these pairs are not directed and can thus be "satisfied" by a mapping (i.e., σ( i) = σ( j)) either as a positive or a negative value.

Link with Bhaskaracharya et al. Intra-Array Reuse

Heuristic 1 is the natural extension of Lefebvre-Feautrier successive modulo technique [START_REF] Lefebvre | Automatic storage management for parallel programs[END_REF] with a greedy optimization of the basis. As we mention in the above, it has some similarities with multi-dimensional scheduling and tiling.

While optimizing the basis is usually not needed when K is convex (at least in order of magnitude) [START_REF] Darte | Lattice-based memory allocation[END_REF], this is not the case anymore when K is a union of polyhedra, and different bases can give different mapping sizes in order of magnitude. In particular, Heuristic 1 may suffer from the fact that, at each step, the optimization is equivalent to considering the convex hull of the intersection of K with the orthogonal of the previously-built vectors. Also, minimizing the width to get a small modulo may induce a bad choice for the next steps. This is our motivation, shared with Bhaskaracharya et al. [START_REF] Somashekaracharya | Automatic intra-array storage optimization[END_REF], to find a mechanism to force Heuristic 1 to not choose an hyperplane with smallest width, i.e., smallest modulo. Our result is Heuristic 2, and we now highlight the differences with respect to the work by Bhaskaracharya et al. [START_REF] Somashekaracharya | Automatic intra-array storage optimization[END_REF] in the following.

Their approach is formulated with a relation describing pairs of conflicting elements ( i, j) instead of conflicting differences ij as we do. For intraarray reuse, only one σ is searched, so this is equivalent. We thus rather explain their technique with conflicting differences to simplify the discussions and visualizations.

As we do, they assume that the conflict set K is given by the integer points in a union of polyhedra K = ∪ 1≤j≤r P j , but K is only half of the conflict set to exclude 0 (which can never be separated by any hyperplane) and to make it asymmetric1 . Then, instead of searching for a hyperplane with minimal width, they search for a hyperplane that "separates" as many polyhedra P j of the union as possible, thus fully on one side of the hyperplane, i.e., with either σ( x) ≥ 1 for all x ∈ P j or σ( x) ≤ -1 for all x ∈ P j . Among such solutions, the one with smallest width is chosen. If such a hyperplane is found for all j, it forms a valid one-dimensional allocation. When such a hyperplane does not exist, several hyperplanes are successively applied, focusing on the intersection of K with the orthogonal to the first hyperplane (as in Heuristic 1) until all the conflicts are resolved with the separating hyperplanes.

Although their algorithm may work well in many cases, careful look reveals certain situations where their approach misses good allocations, sometimes even using higher-dimensional array than necessary. The two main problems is that, as for Heuristic 1, minimizing the moduli in such a greedy fashion is not always good, and, more importantly, the resulting allocation is dependent on the way the conflict set is decomposed into a union of polyhedra, due to their primary objective defined at the granularity of the constituting polyhedra. This makes the approach quite unstable.

Figure 7 illustrates the problem with their heuristic, and how the decomposition may influence the allocation. In contrast, our approach finds better allocations for these examples by using a different objective function that does not rely on the way the conflict set is represented as a union of polyhedra. The key point is that, despite the link between multi-dimensional scheduling and affine mapping explained in Section 6.1, there is one fundamental difference. In scheduling, due to the constraint σ( i) -σ( j) ≥ 0 all dependences are made nonnegative (weakly separating hyperplane) and some are made positive σ( i) -σ( j) ≥ 1 (separating hyperplane). If a first schedule σ i separates P i and a second schedule σ j separates P j , then σ i + σ j is a also a schedule that separates both P i and P j . This is the reason why a greedy separating approach [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part II: Multi-dimensional time[END_REF] is optimal for detecting maximal parallelism. Here, this is not true because separation can be ≥ 1 or ≤ -1, which breaks the analogy.

Link with Universal Occupancy Vectors

An occupancy vector is a vector that points to another iteration that can safely reuse the memory location. Universal occupancy vectors (UOVs) are those that are legal for any legal schedule of a program [START_REF] Strout | Schedule-independent storage mapping for loops[END_REF]. For UOVs, the formulation is as follows. An iteration point v may overwrite a value produced at another iteration u if v (transitively) depends on all uses of u. Since all uses of u must be executed for v to execute without violating dependences, the allocation is valid for all possible schedules. Constraints for expressing UOVs are exactly the complement of constraints for expressing conflicts, in the context of scheduleindependent mappings. The theory of UOV is however rather limited, as only codes with uniform dependences and constant UOVs are captured. Otherwise, the set of legal UOV is hard to define. Thies et al. [START_REF] Thies | A step towards unifying schedule and storage optimization[END_REF] 7a, a hyperplane exists that does not intersect with the entire conflict polyhedra, satisfying their primary objective function. This leads to N extra storage as we illustrated in Section 2. Figure 7b is another possible decomposition where no hyperplane can partition the entire conflict set. In this case, their primary objective to maximize the number of polyhedra (in the disjoint union) fully on one side of the hyperplane, in combination with the secondary objective to minimize the width of a dimension in the final array, gives the vertical line. This produces N 2 storage in the end, which is worse by one full dimension that what we achieve. to handle affine dependences, but they had to restrict to the set of legal affine one-dimensional schedules. In this case, the set of all constant valid constant occupancy vectors (AUOV, for affine UOV) can be defined.

Although allocations legal for all possible schedules may seem too conservative, UOVs can still give efficient allocations for many programs. In fact, it is no coincidence that UOV-based allocations (and in particular QUOV-based allocations [START_REF] Yuki | Memory allocations for tiled uniform dependence programs[END_REF], designed for tiling) find good allocations for live-out values of tiled stencil programs [START_REF] Somashekaracharya | Automatic intra-array storage optimization[END_REF]. If we look at an allocation for a tile, the live-out values are those that must be preserved at the end of the execution of a tile. These values should never be overwritten since there are other uses outside of a tile. Thus, even if the tile itself as a specific schedule, live-out variables are captured well with schedule-independent mappings, e.g., with UOV-based allocations. These exploit possible reuse within the tile while keeping the live-out values.

Partially due to the restriction of UOVs that only use one projection (hence using d -1 dimensional array for d dimensional iteration space), the primary heuristic used in the search for a UOV is the length of the UOV [START_REF] Strout | Schedule-independent storage mapping for loops[END_REF]. The gcd of the UOV is directly connected to the memory consumption through modulo factors along the projection. Increases in the UOV length that do not influence the gcd often increase the memory usage by making the projection to be more steeply angled. For example, projecting a N × N domain along the vector (1, 0) T gives a line of length N , while a diagonal projection along (1, 1) T gives 2N -1. If such a projection is along some boundary of the domain, it may decrease the memory usage, but the length of the UOV is a good approximation otherwise. This shares the same idea with the heuristic used in our approach.

Conclusion and Future Work

In this paper, we have extended the lattice-based memory allocation in two important directions: unions of polyhedra, and better objective functions to find the basis vectors. The key insight is taken from the two dual approaches in the original lattice-based method, the first that works with hyperplanes on one side and the second with reuse vectors on the other side. We have shown that reuse vector based selection of the basis of the mappings finds more compact mappings with different examples.

Several research directions directly follow the work presented in this paper.

• One immediate future work is to explore how the Heuristic 2 can be made fully parametric. • Reasoning about the optimality of the resulting mapping. Can we guarantee that we obtain the mapping that is dimension-wise optimal? • While Heuristic 1 can be generalized to inter-array optimizations where different statements have different mappings, as done by Bhaskaracharya for their intra-array mapping [START_REF] Bhaskaracharya | SMO: An integrated approach to intra-array and inter-array storage optimization[END_REF], it is unclear how to extend our lattice-based method with reuse vectors. • Although we have focused on compactness of the allocation in this paper, it is not always the best for performance. All prior work on memory allocation in the polyhedral literature have never taken locality into account. An important direction to explore is how we can explore memory layout transformations for improved performance, e.g., through more cache-friendly layouts, and in particular cases where introducing redundant storage leads to better overall performance.

In this paper, we have separated the construction of conflict sets as an orthogonal component, since our method works for conflict sets that comes from many different inputs. In addition to the classical conflict sets computed based on multi-dimensional affine schedules, we may take those from other specifications, including explicitly-parallel specifications (such as OpenMP, X10, OpenStream, and so on), or when we consider kernel offloading with overlapped computation and I/O. Memory reuse analysis is important in many different scenarios, and now we have provided the theory that seamlessly generalizes across different cases.

Figure 1 :

 1 Figure1: An example to illustrate the key intuition. The reverse-L shape of width 2 may be viewed as the live-out values of a square tile of size N , when the data dependences (not depicted here) across tiles have a uniform length of 2 in each dimension. The memory allocation characterized by the reuse vectors in Figure1aleads to the optimal affine mapping (projection plus modulo 2) that we seek. Figure1billustrates another possible affine mapping (a projection) that may be found by other techniques, with N extra storage.

Figure 3 :

 3 Figure3: Optimized version of Heuristic 2 for the set K of differences of the live-out points in Figure1. The selected basis is a 1 = (1, 1) T , then a 2 = (1, 0) T , or equivalently (0, 1) T . The figure to the left shows the basis vector of the lattice l 1 = ρ 1 a 1 = (2, 2) T , which is the shortest vector (a multiple of a 1 ) that points outside of K. In the next step, the set K is extruded along a 1 and is now the infinite diagonal band shown in the right figure, leading to the second basis of the lattice: l 2 = ρ 2 a 2 = (2N -1, 0) T . The final mapping is given by the inverse of the matrix [ a 1 a 2 ] -1 = 0 1 1 -1 and the modulo factors defined by ρ 1 = 2 and ρ 2 = 2N -1, which is the mapping (x, y) → (y mod 2, xy mod (2N -1)).

Figure 4 :

 4 Figure4: How the basis selection techniques for the two heuristics favor one basis over the other. The basis selection for the mapping space (Heuristic 1) favors, for the first vector, the canonical axis (1, 0) (or equivalently (0, 1)) over (1, -1) because the width along (1, -1) is wider. As shown in the figure to the left, the width is defined by the vector (N -1, 1 -N ) with a width equal to 2N -2, which is larger than N -1 obtained from the canonical basis. In the lattice space, i.e., for Heuristic 2, the canonical basis leads to a much longer reuse vector (it has to be multiplied by ρ = N ) as depicted in the right figure. Since a shorter (and constant) reuse vector is preferred in the basis selection, (2, 2) is selected.

  y]: N>=2 and (-1<=-x,y<N or -1<=x,-y<N)}; K_1 := range Kp; nK_1 := {[x,y]} -K_1; Norm := { [x,y]->[m,-x,-y]: m >= x,y,-x,-y }; a_1 := (Norm^-1)(lexmin (Norm(nK_1)));

  1*t]: t>0}; b := (fK.Opt)*Ortho_1; b_min := lexmin range b; c_2 := sample (b^-1)(b_min);

Figure 6 :

 6 Figure 6: Diamond tiling, the conflict differences are

Figure 7 :

 7 Figure 7: How the technique by Bhaskaracharya et al. [2] encounters difficulties with the example from Figure 1. Figures 7a and 7b are two possible decompositions of the asymmetric version of the conflict set. In Figure7a, a hyperplane exists that does not intersect with the entire conflict polyhedra, satisfying their primary objective function. This leads to N extra storage as we illustrated in Section 2. Figure7bis another possible decomposition where no hyperplane can partition the entire conflict set. In this case, their primary objective to maximize the number of polyhedra (in the disjoint union) fully on one side of the hyperplane, in combination with the secondary objective to minimize the width of a dimension in the final array, gives the vertical line. This produces N 2 storage in the end, which is worse by one full dimension that what we achieve.

Table 1 :

 1 1, Comparison between the "successive modulo" technique (LeFe), Bhaskara. & al. algorithm (BBA) and our approach (Lattice).

	Example	Algo.

RR n°8840

How this is performed is not given. We assume that it is done by taking the intersection with the strict lexicographic order.RR n°8840

Acknowledgements

We thank Paul Feautrier for fruitful discussions, in particular for the Fourier-Motzkin trick to compute the star-shaped extension of K.