Hugo Thievenaz
email: hugo.thievenaz@inria.fr

Keiji Kimura

Christophe Alias
email: christophe.alias@inria.fr

Rephrasing Polyhedral Optimizations with Trace Analysis

Introduction

Compilers must restructure the application to use as best as possible computing and storage resources. In general, compiler optimizations are fragile and highly depend on the shape of the source code. The polyhedral model [START_REF] Feautrier | Some Efficient Solutions to the Affine Scheduling Problem, Part II: Multi-Dimensional Time[END_REF][START_REF] Feautrier | Polyhedron Model[END_REF][START_REF] Quinton | The mapping of linear recurrence equations on regular arrays[END_REF][START_REF] Sanjay | On synthesizing systolic arrays from Recurrence Equations with Linear Dependencies[END_REF] provides the mathematical foundations to develop domain-specific compiler optimizations focusing on computeintensive scientific kernels with affine loop kernels manipulating arrays. If polyhedral optimizations are precise and powerful, they often lack of scalability. For instance, maxfusion loop tiling [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF] freezes when applied directly to large loop nests, like those implementing deep learning algorithms [START_REF] Krizhevsky | Ima-geNet Classification with Deep Convolutional Neural Networks[END_REF]. Also, array contraction [START_REF] Alias | Bee+ Cl@ k: An implementation of lattice-based array contraction in the source-tosource translator Rose[END_REF] do not scale when applied to HLS, where hundreds of buffers needs to be allocated [START_REF] Alias | Data-aware process networks[END_REF]. An attentive analysis show that most of the time is spent in static polyhedral operations -projections, parametric ILP, subtractions or combinations of piecewise affine mappings.

In this paper, we defend the iconoclast idea that polyhedral optimizations might be computed without those operations, simply by applying a lightweight analysis on a few off-line execution traces. The main intuition being that, since polyhedral transformation are expressed as affine mappings, only a few points are required to infer the general mapping. Our hope is to compute those points from a few off-line execution traces thanks to a lightweight algorithms. We focus on array contraction, a well known technique to reallocate temporary arrays thanks to affine mappings

A[ ì i] → Â[σ A ( ì i, ì N
)] so the array size is reduced; and whose scalability needs to be improved. Our contributions to this problem are the following:

• We outline an off-line trace analysis to infer a contraction mapping of the form σ A ( ì i, ì N ) = ì i mod ì b( ì N ), which reproduces the results of the successive minima method [START_REF] Lefebvre | Automatic storage management for parallel programs[END_REF].

• In particular, we describe a liveness algorithm from an execution trace, and another to compute the maximum number of variables alive alongside a dimension, from which we get our scalar modular mappings. We show that a simple interpolation allow to infer the modulo mapping ì

N → ì b( ì N ).
This paper is structured as follows. Section 2 illustrates the array contraction problem. Section 3 outlines our trace analysis approach. Section 4 presents our preliminary results. Finally, Section 5 concludes this paper and draws future research directions.

Focus: array contraction

Array contraction consists in finding a mapping

A[ ì i] → Â[σ A ( ì i, ì N )]
reducing the size of A, where ì N is the vector of structure parameters. We illustrate array contraction on the Blur filter kernel:

f o r ( y = 0 ; y < 2 ; y + + ) f o r ( x = 0 ; x<N ; x + + ) blurx [ x , y ] = in [ x , y ] + in [ x + 1 , y ] + in [ x + 2 , y ] ; f o r ( y = 2 ; y<N ; y + + ) f o r ( x = 0 ; x<N ; x + + ) { blurx [ x , y ] = in [ x , y ] + in [ x + 1 , y ] + in [ x + 2 , y ] ; out [ x , y ] = blurx [ x , y -2] + blurx [ x , y -1] + blurx← [ x , y ] ; }
This kernel applies two consecutive elementary convolutions on the input signal in. The first convolution is bufferized to the array bluxr, in turn processed by the next convolution. With the canonical sequential schedule θ , blurx might be contracted with the mapping (x, y, N ) → (x mod N , y mod 3).

So far, that mapping was derived in two steps. First (step 1) an array liveness analysis is applied to compute the conflict relation ▷◁ θ , relating array cells whose liveness intervals intersect. Then (step 2), the conflict relation is analyzed to infer a minimum-range mapping such that conflicting arrays IMPACT'22, June 20, 2022, Budapest, Hungary Hugo Thievenaz, Keiji Kimura, and Christophe Alias cells are mapped to different places:

blurx( ì i) ▷◁ θ blurx( ì j) ∧ ì i ì j ⇒ σ blur x ( ì i, ì N ) σ blur x ( ì j, ì N )
Both steps require polyhedral manipulations.

Step 1 requires to check the emptiness of integer polyhedra. Depending on the contraction algorithm, step 2 may require rational projections ( [START_REF] Lefebvre | Automatic storage management for parallel programs[END_REF]), mixed-LP ( [START_REF] Alias | Bee+ Cl@ k: An implementation of lattice-based array contraction in the source-tosource translator Rose[END_REF]), or parametric lexico-minimization of integer polyhedra ( [START_REF] Somashekaracharya G Bhaskaracharya | Automatic storage optimization for arrays[END_REF]). Both steps are expensive. In this paper, we use an oracle which assumes the precomputation of the conflict relation. Thus, we focus on improving step 2.

Our approach

The following diagram outlines our approach.

Input program (P,θ )

Step 1: Log trace

Step 2: Get optimal mapping for trace

Step 3: Interpolation

Next parameter instance

Step 4: Check with Conflict Set Fail Array contraction mapping σ

Pass

We execute the input kernel on small parameter values (step 1). Then, we apply a liveness dataflow analysis on the trace to infer the conflict relation, summarized as a conflict difference set

D a = { ì i -ì j | a( ì i) ▷◁ θ a( ì j)} (step 2)
. Since D a is finite, the successive minima might be instanciated and computed with a small cost. Finally, we try to infer a parametrized modulo using an affine interpolation. If the mapping interpolated is correct (step 4), it is returned. Else we iterate on the next parameter value. Many open problems need to be investigated. How the parameters should be chosen? So far we considered programs with a single parameter N . In case several parameters are considered, one need to select linearly independant parameters. Other point, the contraction mapping is usually piecewise affine. With our method, we actually seek for the dominating value (on the biggest piece, not on the corners). How to guess parameters on that piece is still an open problem. Finally, we still depend on an oracle computed by means of polyhedral operations. We believe we can get rid of that oracle, providing we actually apply an instance of a polyhedral algorithm (here the successive minima).

Preliminary results

This section presents the preliminary results obtained with our approach. Our method have been implemented and evaluated on the following kernels, whose code is detailed in [START_REF] Thievenaz | Towards Trace-Based Array Contraction[END_REF]:

• pc-1D is a simple producer/consumer through a 1Darray. The producer and the consumer access the array cells in the same order with a phase shift of two steps. • pc-2D-single is a perfect loop nest operating over a 2D array with dependence vectors (1,0) and (0,1). • pc-2D-single is a perfect loop nest operating over a 2D array with dependence vectors (1,0) and (0,1). Only the last array column is live-out. • blur is the motivating example.

Our analyzer and the kernels have been compiled using gcc 9.3.0, the timings are obtained on an Intel Core i5-1135G7 CPU running at 2.40GHz. The results are depicted on the following table. For each kernel and each array we provide the original size, the mapping inferred from our algorithm. It. is the number of iterations before interpolating a correct mapping. For all kernel we tried N = 2, then N = 3, except for the blur kernel, which required N = 3 and 4, since N = 2 is a corner case. Must of the time is lost in the oracle. This must still be improved to provide competitive execution times.

Conclusion

In this paper, we have presented a first step toward a tracebased polyhedral model through array contraction. We have outlined a method able to infer a modulo mapping given a canonical basis from a few off-line executions traces. We have also shown that interpolation of the parameters from a subset of the possible traces is possible, and we have applied it to our examples. We have shown results that reproduce the original successive modulo method and answered positively to the question of generalization from a subset of traces.

In the future, we seek to apply the same process to more benchmarks of the polyhedral community, as our experimentations were limited by time. We also look forward to apply finer analysis on the modulo parametrization, like heuristics to choose parameter values that are relevant to both the conflict set checking and the interpolation, and more generally a finer approach to the choice of basis for any array allocation problem. The first parameter will also be investigated. Finally, we plan to get rid of the oracle and focus on pure trace analysis.


Hugo Thievenaz

Keiji Kimura

Christophe Alias

Lightweight Array Contraction by Trace-Based Polyhedral Analysis ⋆

Keywords: Compilation, Array Contraction, Polyhedral Model, Dynamic Analysis, Memory Allocation

Array contraction is a compilation optimization used to reduce memory consumption, by reducing the size of temporary arrays in a program while preserving its correctness. The usual approach to this problem is to perform a static analysis of the given program, creating overhead in the compilation cycle. In this work, we take a look at exploiting execution traces of programs of the polyhedral model, in order to infer reduced sizes for the temporary arrays used during calculations. We designed a four step process to reduce the storage requirements of a temporary array of a given scheduled program, in which we used an algorithm to deduce array access functions for which bounds are modulos of affine functions of parameters of the program. Our results show memory reductions of an order of magnitude on several benchmarks examples from PolyBench, a collection of programs from the polyhedral community. Execution time is compared to a baseline implementation of a static algorithm, and results show speed-up factors up to 20.

Introduction

The problem of temporary memory allocation is a challenge for programs meant to be running on platforms that have limited computing resources. Such temporary arrays manipulate results of intermediate computations that are disposed of at the end of the program. They are therefore sometimes oversized, when array cells are left unused and not overwritten by following computation despite their value no longer being used. Array contraction is a program transformation whose goal is to detect such unused array cells and replace a write to another cell to an unused one, in order to reduce the effective memory footprint of the array and shrink its maximum size, allocating less memory to the buffer while keeping program correctness intact.

Figure 1 introduces two direct applications of this method. Both the local memories used by the CPUs, and the buffer(s) used for communicating between local memory local memory CPU1 CPU2 bus buffer

Fig. 1: Two CPUs perform computations on local memories and communicate data through a buffer themselves, benefit from allocating smaller arrays. More precisely, many computing units require or may use on-the-fly memory sizing to exploit their structure (scratchpad memories, systems on chips). The class of programs studied for this optimization are affine loop kernels manipulating arrays (SCoPs, static control parts [START_REF] Bastoul | Putting polyhedral loop transformations to work[END_REF]), and this form of loop nest is the most common in many High-Performance Computing examples. The polyhedral model provides the necessary mathematical foundations to develop compiler optimizations such as array contraction, that focuses on compute-intensive scientific kernels containing such SCoPs. In this context, compilers rely on static analysis of the program to reduce the memory footprint of the program. Static analysis has been the basis of many works in the field [START_REF] Alias | Bee+ cl@ k: An implementation of lattice-based array contraction in the source-to-source translator rose[END_REF][START_REF] Bastoul | Putting polyhedral loop transformations to work[END_REF][START_REF] Bhaskaracharya | Automatic storage optimization for arrays[END_REF][START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF][START_REF] Darte | Lattice-based memory allocation[END_REF][START_REF] Simbürger | Polyjit: Polyhedral optimization just in time[END_REF]. However, dynamic analysis outclass static compiler ones when small execution traces can be efficiently produced and analysed. Static methods use polyhedral projections and Integer Linear Programming, which can be expensive depending on the shape of the code. It would seem that no approach to this problem, to our knowledge, has explored the option of using dynamic analysis of the program in order to infer compilation optimizations. In this work, we contradict this habit and study the problem of determining a buffer allocation function from analysis of several execution traces. The problem can be formulated as follows: given a program manipulating a temporary array A, we want to infer allocations functions σ A , of minimal image cardinal, such that any access A[i] can be safely replaced by Â[σ A (i)]. Our general approach is then to apply a lightweight analysis on a few offline execution traces, with the assumption that the input parameter instances chosen for those traces are small enough that the execution time is significantly smaller.

In this paper, we make the following contributions:

-We present a new method for storage optimization, a dynamic approach that uses offline execution trace analysis. In particular, we describe a liveness algorithm from such execution trace, and another to compute the maximum number of variables alive alongside a dimension, from which we get our scalar modular mappings.

-We show, through the use of interpolation, that we can identify parameters from said modulo and deduce a generalized mapping. -We implement this method on several benchmarks from the polyhedral community and show reductions both in implementation execution time and storage mappings deduced.

Our paper is structured as follows. Section 2 outlines the polyhedral model and the array contraction problem. Section 3 discusses related work. Section 4 describes our trace-based approach. Section 5 presents experimental results. Finally, Section 6 concludes this paper and draws research perspectives.

Background

We present the necessary background to the problem. We define the polyhedral model, and what is an usual polyhedral compilation flow. Then, we define the problem at hand, array contraction, and the related notions. The polyhedral model defined by [START_REF] Feautrier | Polyhedron model[END_REF] is an intermediate representation of a loop nest as a graph over points of Z n . The class of programs that can be represented in this model, and therefore are subject to polyhedral optimizations, is polyhedral programs. These are (sequences of, possibly nested) for loops where all loop bounds and conditions are affine functions of the surrounding loop iterators and program parameters. Each execution of a statement S, nested in a n-depth loop, namely an instance or operation, can be represented by ⟨S, i⟩ where i is a n-dimensional iteration vector of the surrounding loop indices. Its iteration domain D, the set of all possible iteration vectors for S, forms a graph over points of Z n .

Polyhedral model

f o r ( y =0; y <2; y++) f o r ( x =0; x<N ; x++) S : blurx [ x , y ] = in [ x , y ] + in [ x+1 ,y ] + in [ x+2 ,y ] ; f o r ( y =2; y<N ; y++) f o r ( x =0; x<N ; x++) T : blurx [ x , y ] = in [ x , y ] + in [ x+1 ,y ] + in [ x+2 ,y ] ; U : out [ x , y ] = blurx [ x , y -2] + blurx [ x , y -1] + blurx [ x , y ] ;
Running example We illustrate our algorithm on the 2D Blur filter illustrated with its iteration domains on Figure 2. This is a well-known example of the polyhedral community, that applies two consecutive elementary convolutions on the input signal in. We have represented the iteration domains of S, T and U as colored symbols. Polyhedral compilation flow Compilation flows usually produce an in-between form of the program at hand, named Intermediate Representations (IR), on which transformations can be applied more easily to optimize its execution. Polyhedral compilation is no different, and in the case of our method, we produce those IRs as a pre-analysis step, separated from the algorithm. Figure 3 describes our simplified polyhedral compilation flow. Source code gets transformed once to an intermediate representation through polyhedral raising, who is then subject to possibly multiple affine transformations in order to optimize execution of the target code. The focus of our algorithm, and its performance, is therefore focused entirely on the application of such polyhedral optimization, for which ours is Array Contraction.

Affine transformations. At the heart of a polyhedral compiler, code transformations are expressed by affine mappings specifying a new execution order: Definition 1 (Affine Scheduling). A schedule maps each execution instance ⟨S, i⟩ to an execution date θ S (i). In the polyhedral model, schedules are affine per statement, i.e. θ S (i) = A S i+b S , and dates are vectors of Z p ordered with the lexicographic ordering ≪. A schedule maps each iteration vector to its counterpart in the transformed, scheduled program.

A possible schedule for our motivating example is the canonical sequential schedule θ, which is the order specified by the original for statements of the program: θ S (y, x) = (0, y, x, 0), θ T (y, x) = (1, y, x, 0), θ U (y, x) = (1, y, x, 1).

Correctness. In the polyhedral model, data dependencies might be expressed between iterations. In the computation of the blur-interleaved example, for any instance where y ≥ 3, we have that an instance of ⟨U, x, y⟩ of the second convolution depends (flow) on the preceding instances ⟨T, x, y⟩, ⟨T, x, y -1⟩, ⟨T, x, y -2⟩ of the first convolution. Anti-and output-dependencies are expressed in the same way. The dependence relation is denoted by →. Of course, the schedule is constrained by data dependencies:

⟨S, i⟩ → ⟨T, j⟩ ⇒ θ S (i) ≪ θ T (j) (1)
This gives affine constraints which allow to compute affine schedules [START_REF] Feautrier | Polyhedron model[END_REF].

Array Contraction

The problem of array contraction, given a temporary array A, consists in finding a mapping

A[i] → Â[σ A (i)]
reducing or matching the unknown required size of A, minimal size for which the correctness of the program is intact. In our case, we seek memory mappings of the form σ A (i) = i mod b(N ), where b is an affine function of program's structure parameters N (e.g. array size).

Definition 2 (Conflict Relation).

A conflict relation ▷◁ θ is defined as the set of array cells whose lifetimes intersect during the execution of the program for the schedule θ.

The conflict relation induces a correctness condition on array contraction, as conflicting array cells might be mapped to different places:

a[i] ▷◁ θ a[j] ∧ i ̸ = j ⇒ σ a (i) ̸ = σ a (j) (2)
Running example (cont'd) With the original loop schedule, the temporary array blurx might be contracted with the mapping (y, x) → (y mod 3, x mod N ), when N ≥ 3. Indeed, blurx bufferizes the first convolution (S,T) before applying the next convolution (U) which only needs three rows y. This way, the footprint is reduced to 3 × N array elements. The successive minima technique [START_REF] Lefebvre | Automatic storage management for parallel programs[END_REF] is the state-of-the-art approach to compute such mappings. The method of this work by Lefebvre and Feautrier can be boiled down to the following process. The conflict relation is represented as a difference set ∆ a = {i-j | a[i] ▷◁ a[j]} for each array a; Then, for each array dimension k, the modulos are computed with b k

(N ) = 1+max{δ k | (δ 1 , . . . , δ n ) ∈ ∆ a }.
Finally, resolved conflicts are removed before iterating on the next array dimension :

∆ a := ∆ a ∩ {δ | δ k = 0}.
Our contribution, as we will show later, consists in a lightweight instantiation of this algorithm on several small execution traces, followed by an interpolation to obtain a general mapping. We show experimentally that our results are obtained way faster than with the Lefebvre-Feautrier method, the latter dealing with costly parametric Integer Linear Programming (ILP). More fundamentally, this work is a proof-of-concept that costly polyhedral analysis might be rephrased as lightweight trace analysis. This opens new perspectives to scale polyhedral compilers.

Related work

We quickly go over the multiple works related to our subject. We first present defining works on the subject. We then go over closely related work on the subject of array contraction. Finally, we present loosely related work on trace manipulation and analysis, but no work on dynamic array contraction by trace analysis has crossed our eyes.

Affine array contraction As described by [START_REF] Bhaskaracharya | Automatic storage optimization for arrays[END_REF][START_REF] Lefebvre | Automatic storage management for parallel programs[END_REF], and again in this work, the successive modulo technique seeks to reduce the memory storage requirements of an already scheduled program, by performing static analysis in order to construct a conflict set. The array dimensions are reduced by finding contraction moduli along the array's axes. While recalling that the method of Lefebvre and Feautrier [START_REF] Lefebvre | Automatic storage management for parallel programs[END_REF] obtains on the example blur-interleaved a storage mapping (y, x) → (y mod 3, x mod N ), the more advanced work of Bhaskaracharya et al. [START_REF] Bhaskaracharya | Automatic storage optimization for arrays[END_REF] infers a more refined mapping (y, x) → 2x -y mod (2N + 1), because their approach consider the change to a better basis for the contraction vectors.

Inter-array and intra-array contraction The type of optimization we are looking for in this paper is an intra-array optimization as designed by [START_REF] Bhaskaracharya | Automatic storage optimization for arrays[END_REF], and references such as [START_REF] Alias | Bee+ cl@ k: An implementation of lattice-based array contraction in the source-to-source translator rose[END_REF][START_REF] Feautrier | Polyhedron model[END_REF][START_REF] Quilleré | Optimizing memory usage in the polyhedral model[END_REF] focus on this intra-array analysis. This means that the analysis performed is done on a per-array basis. [START_REF] Bhaskaracharya | Automatic storage optimization for arrays[END_REF][START_REF] Bhaskaracharya | Smo: An integrated approach to intra-array and inter-array storage optimization[END_REF] build a technique for intra-array as well as inter-array optimization, a technique that consider the reduction of multiple temporary arrays, allowing them to find even more reduced mapping by changing (often reducing) the dimensionality of the array(s) considered for the analysis. [START_REF] Clauss | Symbolic polynomial maximization over convex sets and its application to memory requirement estimation[END_REF] calculate the memory requirements of a program by approaching them as a polynomials in the parameters of the program, but their method has to relax the solutions as rational instead of integer, and only give an upper bound of the memory consumption.

Trace analysis In terms of trace analysis, some work has already covered similar topics such as loop recognition and trace prediction [START_REF] Ketterlin | Prediction and trace compression of data access addresses through nested loop recognition[END_REF] and trace-based affine loop reconstruction [START_REF] Rodríguez | Trace-based affine reconstruction of codes[END_REF]. The former compresses traces (as sequences of scalars) and constructs a loop nest producing such sequence of numbers. The latter focuses on reconstructing loops based on their predictable affine behavior, from the addresses of the memory accesses, and presents a terminating algorithm to reconstruct the loop function entirely. These works, therefore, focus only on rebuilding incomplete traces, and not on the usage of traces in a compilation process.

Our approach

This section presents the contributions made to the problem of array contraction, and detail our method of offline trace execution and analysis to infer a general mapping. 4 depicts our approach. We start with an input Intermediate Representation, comprised of the program, its schedule θ, and its Dependence Graph. We also input the code generated from θ (Template code), which will allow to produce traces.

Overview of the approach

First, we compute the input parameter instances on which running the program to obtain interesting traces. We also instrument the template code to prepare the trace generation. Then, for each input parameter instance N , we generate the trace (Log trace) and we apply a lightweight instance of the successive minima method (Get trace mapping). We end up with a collection of trace mappings. Finally, we infer the general mapping (working for any input parameter) from an interpolation between each input parameter instance and its corresponding output modulo mapping (Interpolation). All these steps are detailed thereafter.

Generating input parameter instances

Our trace analysis operates on execution traces of programs, meaning we have to instantiate our kernel program with scalar values for its parameters N . Since we would like to interpolate modulos as an affine forms of parameters N → b k (N ), we need |N | + 1 parameter instances. Also, the parameter instances must be independent to enforce a unique interpolation. We first explain how the first parameter instance is computed. Then, we explain how we get the remaining parameter instances.

First parameter instance, O We derive the set of parameters covering all the dependencies by projecting the dependency constraints of the program on each of the parameters. Then, we compute a minimum value for each of the parameters. Usually the constraints are of the form N ≥ ℓ with ℓ some constant lower bound. Hence, we may infer a lower bound for each parameter with a simple syntactic heuristic without using expensive linear programming techniques. This gives the first parameter instance, denoted as O (for "origin").

This heuristic assumes that in the main program execution, all the dependencies are reached. However, this is not always the case, and if we deal with imperfect loop nests, then the result is an upper bound, giving potential additional overhead by operating on traces with bigger parameter instances than needed. How to extend this heuristic to the general case is left to future work. In our running example, the intersection of all the dependence constraints boils down to the polyhedron {N |N ≥ 3}. Hence O = (3), denoting the parameter instance N = 3.

Remaining parameter instances One set of parameter values is not sufficient to establish an interpolation. For each parameter N i , 1 ≤ i ≤ p, we create a new parameter instance, linearly independent from the rest. A straightforward way to compute such new instances is to build the set of increments by each canonical vector e i : I = (O, O + e 1 , . . . , O + e p ), where O + e i = (N 1 , ..., N i + 1, ..., N p ). These parameter instances will lead to a unique affine interpolation, as we will see later. On our example, we would obtain [START_REF] Bhaskaracharya | Automatic storage optimization for arrays[END_REF][START_REF] Bhaskaracharya | Smo: An integrated approach to intra-array and inter-array storage optimization[END_REF], denoting the parameter instances N = 3 and N = 4.

Inferring a mapping on a trace

The following algorithm 1 describes our lightweight instance of the successive minima method to compute a mapping from an execution trace.

We apply a direct liveness analysis on the trace to compute the difference set ∆ a we defined earlier, then each modulo is computed as the maximum difference measured alongside each array dimension i, following the lines of the successive minima algorithm. Since we deal with finite (and small) integer sets, no ILP is required.

This method is an instance of the Lefebvre-Feautrier algorithm [START_REF] Lefebvre | Automatic storage management for parallel programs[END_REF] since we produce the same mapping while operating on a trace, i.e. an instantiated program. Since the mapping is affine, we can directly apply a linear interpolation to deduce a generalized parametrized mapping, as we will describe in the next section.

Algorithm 1: Find the mapping for the array a on the trace T

Result: mapping i → i mod m function GetMapping(T ,a) (In, Out) ← liveness(T ) CS ← p {(a[i], a[j]) | a[i], a[j] ∈ In(p)} ∆a ← {i -j | (a[i], a[j]) ∈ CS}
for each array dimension i, starting from 0, in increasing order do mi ← 1 + max{δi | (0, ..., 0, δi, ...) ∈ ∆a} end Running example (cont'd) For our blur-interleaved example, such analysis would show, on the trace for N = 3, for the blurx array, that the biggest width alongside the y axis is 2, so there are a maximum of 2+1 array cells in conflict at any given control point p. The observation on x is similarly done, and again we measure a width of 2, and so a number of conflicts of 3. Hence, we obtain (y, x) → (y mod 3, x mod 3). This is repeated for the trace with N = 4 where we obtain (y, x) → (y mod 3, x mod 4).

Interpolation

From the mapping instances deduced, we show how to interpolate a generalized mapping that depends on program parameters. We retrieve mappings of the form i → i mod b(N ) by a direct affine interpolation from the pairs of inputs (parameter instances) and outputs (modulo scalars found). We realise this by solving the following systems of equations. Let p be the number of program parameters (N = (N 1 , . . . , N p )), and k the number of array indices (i = (i 1 , . . . , i k )). Then:

σ a (i) =    i 1 mod f 1 (N 1 , ..., N p ) . . . i k mod f k (N 1 , ..., N p )    (3)
This system of equation ( 3) defines the f ℓ functions that we are determining. Expecting to deal with affine functions, each f ℓ can be written in the homogeneous form:

f ℓ (N ) = τ ℓ • N 1 (4)
where τ ℓ is a vector of size p + 1, τ i being the coefficient of N i , for 1 ≤ i ≤ p, and τ p+1 being the constant coefficient.

We have collected sample values from these affine functions, for each array index ℓ, represented by:

     f ℓ (O) = m 0 . . . f ℓ (O + e p ) = m p (5)
Which can be written as Aτ ℓ = m:

     O 1 O + e 1 1 . . . . . . O + e p 1      τ ℓ =    m 0 . . . m p    (6)
We now show that this system has always a unique solution in Z p+1 :

Theorem 1. A is unimodular.

Proof. We apply the Gaussian elimination method to express the determinant of A. We may substract from each of the first p columns that we label each a i , the last column a p+1 multiplied by N i without changing the determinant. The resulting matrix is as such:

det A = 0 • • • • • • 0 1 1 0 • • • 0 1 0 1 • • • 0 1 . . . . . . . . . . . . . . . 0 • • • 0 1 1
It immediately follows that the determinant of this matrix is (-1) p+1 × det I p , the permutation of the p + 1-th and the first column leading to the (-1) p+1 factor, and det I p the determinant of the lower-left matrix which is the identity. Therefore, det A = ±1 and so A is unimodular.

⊓ ⊔

Because A is unimodular, the linear equation system always has integer solutions. Therefore, for any given program with p parameters, p + 1 traces are both necessary and sufficient to produce an interpolation.

Running example (cont'd) On traces, we obtained the trace mappings (y, x) → (y mod 3, x mod 3) for N = 3, and (y, x) → (y mod 3, x mod 4) for N = 4. Denoting (y, x) → (y mod b 1 (N ), x mod b 2 (N )) the general mapping, we have b 1 (3) = 3 and b 1 (4) = 3. Hence we solve:

Aτ 1 = 3 1 4 1 τ 1 = 3 3 from which we deduce τ 1 = 0 3 . Hence b 1 (N ) = τ 1 • N 1 = 3.
Also, b 2 (3) = 3 and b 2 (4) = 4. Hence we solve:

Aτ 2 = 3 1 4 1 τ 2 = 3 4 from which we deduce τ 2 = 1 0 . Hence b 2 (N ) = τ 2 • N 1 = N .
Hence, we get the parametrized program mapping (y, x) → (y mod 3, x mod N ).

Experimental Results

This section presents our implementation and the results obtained with our approach, and make a comparison with the successive minima approach.

Experimental setup

We have implemented our method as an automatic code generator in C++ named PoLi. Our tool takes as input an intermediate representation of a kernel and first outputs a C program where its statements have been swapped out with calls to trace generation methods. Then, the lifetime analysis is performed on several execution traces, from a remote compilation and execution of the modified kernel. Finally, the deduced mapping is directly applied by modifying the access functions of the temporary arrays to the ones deduced.

The baseline implementation, used to compare our analysis time and storage requirement measurements with, is an implementation of the successive modulo technique [START_REF] Lefebvre | Automatic storage management for parallel programs[END_REF]. The C kernels have been compiled using gcc 9.3.0 with flags "-fPIC -O3", while the implementation itself has been compiled using g++ 9.3.0 using flags "-O3 -ldl -lstdc++fs". Every compilation and execution of the kernels, and so their time measurements have been done on an Intel Core i5-1135G7 CPU running at 2.40GHz. No HPC computer is required, as we deal with compilation, not execution. The LF method is compiled with the same directives. We list the examples present in our benchmarks, which are part of the PolyBench test suite3 :

-fibonacci computes the n-th term of the fibonacci sequence. It showcases very low runtime because of a very simple single loop nest. -pc-2d and pc-2d-line, two examples of a producer-consumer mechanic in two dimensions, respectively without and with the last 2 rows of A explicitly being output dependencies. Those show the relevance of the method to only temporary memory. -blur-interleaved and blur-tiled, two examples of the 2D blur filter, respectively with producer-consumer statements interleaved (motivating example), and tiled scheduling. Together, they highlight the versatility of the method, matching the Lefebvre-Feautrier approach for the interleaved case, but outpaces it when the loop nest gets more complex with more loop counters added for the tiling. -2mm, example of two successive matrix multiplication and assignment. This example shows that the Lefebvre-Feautrier method also suffers from the multiplicity of arrays in the program, which skyrockets its runtime compared to our approach.

Results

Table 1 depicts the kernels and their targeted temporary array, alongside its original size, and the mapping found is the reduced size inferred from our algorithm. Parameter instance represents the starting parameter values chosen for the analysis. The execution times shown are not the ones of the modified kernels' executions, as the mappings found are the same for both methods. Rather, the first average runtime describes, for our method PoLi, the sum of the measured 1: Mappings and runtimes obtained using our approach (PoLi) compared to the baseline successive modulo method (LF) [START_REF] Lefebvre | Automatic storage management for parallel programs[END_REF] time spent on the generation of the parameter instance, the time spent instantiating the trace and the time spent interpolating the resulting mappings. This is compared to the baseline runtime which represents the time spent applying the instance of the Lefebvre-Feautrier approach we have implemented, and we show the speed-up factor between the two methods ran successively. We can observe that the fibonacci example has a speedup of more than 20, explained by the small trace parameters chosen, as the runtime of PoLi on this example is noticeably the lowest out of all. blur-tiled and 2mm have respectively bigger iteration dimension and a greater overall number of arrays, meaning the Lefebvre-Feautrier approach irremediably takes more time projecting over those several dimensions, whereas our method takes advantage of the smallness of the parameter instances selected and suffers way less from more arrays and array dimensions. The complexity of the LF method is directly tied to the iteration dimension in an exponential fashion, while our approach is less sensitive to it. blur-interleaved and pc-2d-line both present smaller speed-up factors, as our parameter instance generation gives an upper bound too big, while the dependencies can still be respected with lower parameter values. Therefore, more carefulness is required in the selection of the starting parameter instance, meaning that a better method to infer parameter instances is also of the essence. On these two examples, our approach still manages to match or outperform the Lefebvre-Feautrier method while having unnecessarily large starting parameter instances.

Conclusion

In this paper, we have presented a novel lightweight method for array contraction in the polyhedral model. This work is the very first step towards a new paradigm of trace-based analysis to scale polyhedral compilation and demonstrate a promising proof of concept on the array contraction problem. We design and implement an automatic array contraction tool, that takes as input the source code of the kernel and outputs optimized target code in regards to storage space consumption. We present the algorithms and methodology used in our tool. Execution times are compared to those of the Lefebvre and Feautrier method and shows promising speed-up factors. Results answers positively to the question of the possibility of generalization from a subset of execution traces.

In the future, we seek to apply another methodology to the starting parameter instance deduction, in order to choose minimal parameters regardless of the form of the loop. We also look forward to deduce more complex mappings, of the form i → Ai mod b(N ), similarly to [START_REF] Alias | Bee+ cl@ k: An implementation of lattice-based array contraction in the source-to-source translator rose[END_REF][START_REF] Bhaskaracharya | Automatic storage optimization for arrays[END_REF]. More generally, we seek to apply the paradigm of trace analysis to other problems of the polyhedral compilation, to further study the potential yield of trace analysis in the compilation process. Finally, we plan to address the interaction of array contraction and other optimizations passes by integrating our implementation in an automatic parallelizer.

Fig. 2 :

 2 Fig. 2: Motivating example: 2D Blur filter

Fig. 3 :

 3 Fig. 3: Simplified polyhedral compilation flow of our method

Fig. 4 :

 4 Fig. 4: Our approach

  Table

	Kernel	Mapping found	Parameters	PoLi time (ms)	LF time (ms)	Speed-up
	fibonacci	i mod 2	N = 2, 3	0.00103	0.024221	23.5
	pc-2d	i mod N	N = 2, 3	0.00284	0.045513	16.0
		j mod N
	pc-2d-line	i mod 2	N = 3, 4	0.01022	0.064114	6.3
		j mod N
	blur-2d	y mod 3	N = 5, 6	0.15636	0.187037	1.2
	blur-tiled	x mod N y mod 3	N = 5, 6	0.166067	4.041242	24.3
	2mm	x mod 4 i mod N	N = 2, 3	0.096936	2.228872	23.0
		j mod N

available at https://web.cse.ohio-state.edu/ pouchet.2/software/polybench/


Ilham Lasfar

Christophe Alias

Matthieu Moy

Rémy Neveu

Alexis Carré

Affine Multibanking for High-Level Synthesis

Keywords: High-Level Synthesis, Multibanking ACM Reference Ilham Lasfar, Christophe

come   L'archive ouverte pluridisciplinaire

Introduction

Since the end of Dennard scaling, the energy efficiency (flop/J) of computers has become a major challenge as soon as the energy budget is limitated. The best solution is to rely on specialized circuits, which ultimately trade energy efficiency for programmability. In the last decade, FPGAs appeared as a credible alternative for big data and high-performance computing applications. However, programming an FPGA is tedious: given a function to implement, the circuit configuration must be built from scratch by the developer. Hence the emergence of high-level circuit compilers (high-level synthesis, HLS) [START_REF]Nios II C2H Compiler User Guide[END_REF][START_REF] Catapultc | Mentor CatapultC High-Level Synthesis[END_REF][START_REF] Gaut | Gaut: High-Level Synthesis tool From C to RTL[END_REF][START_REF]Impulse-C, Accelerate Software using FPGAs as Coprocessors[END_REF][START_REF] Ugh | Ugh: User-Guided High-Level Synthesis[END_REF], able to translate a C program to an FPGA circuit configuration. Unlike software parallelisation, there is no parallel runtime to place the computation and the data among processing elements: all the parallelization decisions must be taken at compile-time.

In this short paper, we outline a source-to-source transformation to address the multibanking problem: data are mapped to distinct memories (banks) so parallel access will target distinct banks to avoid data transfer serialization. Given a program and a schedule prescribing parallelism, we are able to infer a complete reorganization of the data into banks and to generate the transformed program accordingly. This problem has been investigated for simple kernels with perfect loop IMPACT'22, June 20, 2022, Budapest, Hungary 2022.

nests [START_REF] Cong | Automatic memory partitioning and scheduling for throughput and power optimization[END_REF][START_REF] Wang | Memory partitioning for multidimensional arrays in high-level synthesis[END_REF][START_REF] Wang | An integrated and automated memory optimization flow for FPGA behavioral synthesis[END_REF], usually on convolution-like kernels with differents patterns. In a different context, approaches for false sharing removal [START_REF] Zinenko | Modeling the conflicting demands of parallelism and temporal/spatial locality in affine scheduling[END_REF] expose interesting ideas -though not directly transposable. We outline a general HLS algorithm which subsumes these approaches and makes the following contributions:

• We propose a novel and general formalization of the multibanking problem, which subsumes the previous approaches. • We propose a complete algorithm to compute our multibanking transformation using the polyhedral model. • Our approach reduces the number of banks and the maximal bank size, without hindering parallel accesses.

This paper is structured as follows. Section 2 illustrates the problem of multibanking on a motivating example. Section 3 outline the main ideas of our algorithm for multibanking. For a complete description, the reader is refered to [START_REF] Lasfar | Affine Multibanking for High-Level Synthesis[END_REF]. Section 4 outlines our preliminary results. Finally, Section 5 concludes this paper and outlines future work.

Multibanking

We illustrate the problem on the 2D convolution product, depicted in Figure 1.(a). The loop is assumed to be executed in sequence. For each iteration, array accesses are done in parallel. Because of memory limitation, parallel references must be mapped to different banks. A solution depicted on Figure 1.(b) is to put each reference in(i, j) into memory bank bank in (i, j) = 3i + j mod 9 at offset offset in (i, j) = i mod N . In general, we seek affine mappings bank a :

ì i → ϕ a ( ì i) mod σ ( ì N ) and offset a : ì i → ψ a ( ì i) mod τ ( ì N )
for each array a where ϕ a ,ψ a , σ and τ are affine functions. These functions define an affine multibanking. Note that those reallocation functions map different arrays to a common array memory organized with outer bank dimensions and inner offset dimensions. The size of that common memory is

Π d σ d ( ì N ) × Π d τ d ( ì N ), where σ d ( ì N ) denotes the d-th dimen- sion of vector σ ( ì N ).

Overview of our multibanking algorithm

This section outlines our multibanking algorithm, we present the intuitions behind the derivation of the bank mapping and the offset mapping. For a deeper level of detail, reader is refered to [START_REF] Lasfar | Affine Multibanking for High-Level Synthesis[END_REF]. A bank mapping bank a ( ì i) = ϕ a ( ì i) mod σ ( ì N ) is correct w.r.t. the parallel execution prescribed by a schedule θ iff two different memory cells a( ì i) and a( ì j) accessed at the same time belong to different banks:

for ( i =1; i <N -1; i ++) for ( j =1; j <N -1; j ++) out [i , j ] = in [i -1 ,j -1]+ in [i -1 , j ]+ in [i -1 , j +1]+ in [i ,j -1] + in [i , j ] + in [i , j +1] + in [ i +1 ,j -1]+ in [ i +1 , j ]+ in [ i +1 , j +1
a( ì i) ∥ θ b( ì j) ∧ ì i ì j ⇒ bank a ( ì i) bank b ( ì j)
We relax the formulation to postpone the computation of the modulo:

a( ì i) ∥ θ b( ì j) ∧ ì i ≪ ì j ⇒ ϕ a ( ì i) ≪ ϕ b ( ì j) ϕ a ( ì i) ≪ ϕ b ( ì j)
means that there exists a dimension d such that ϕ d a ( ì i) < ϕ d b ( ì j) and both vectors are identical above:

ϕ ℓ a ( ì i) = ϕ ℓ b ( ì j) for all ℓ < d.
Hence, the dimensions of ϕ might be computed incrementally across dimensions d, as an affine schedule would be. Once a dimension ϕ d is found, we focus on the unresolved parallel conflicts (still in the same bank) with

ϕ d a ( ì i) = ϕ d b ( ì j).
And we iterate on the next dimension until all conflicts are resolved.

The number of banks might be estimated with maximum difference ϕ d b ( ì j) -ϕ d a ( ì i) for conflicting cells a( ì i) and b( ì j). It gives the modulo value σ d , along dimension d:

a( ì i) ∥ θ b( ì j) ∧ ì i ≪ ì j ⇒ ϕ d b ( ì j) -ϕ d a ( ì i) ≤ σ d ( ì N )
The coefficients of the affine form σ d might be minimized lexicographically, under the constraints of correctness exposed above. Iterating the process for each dimension yield a general, correct and efficient bank mapping. These constraints are turned to existentially guarded affine constraints thanks to the affine form of Farkas lemma, following the lines of [START_REF] Feautrier | Some Efficient Solutions to the Affine Scheduling Problem, Part II: Multi-Dimensional Time[END_REF], formalized as a domain-specific language in [START_REF] Alias | fkcc: the Farkas Calculator[END_REF].

Offset mapping Finding the offset in a bank might be achieved by the same algorithm, on different constraints.

An offset mapping offset a ( ì i) = ψ a ( ì i) mod τ ( ì N ) is correct w.r.t. a schedule θ iff two array cells a( ì i) and a( ì j) mapped to the same bank and whose liveness conflict (a( ì i) ▷◁ θ a( ì j)) are mapped to different offsets:

bank a ( ì i) = bank b ( ì j) ∧ a( ì i) ▷◁ θ b( ì j) ∧ ì i ì j ⇒ offset a ( ì i) offset b ( ì j)
The dimensions of ψ a and τ ( ì N ) are computed incrementally in the same way as the bank mapping, with similar modulo minimization constraints. All the details are given in [START_REF] Lasfar | Affine Multibanking for High-Level Synthesis[END_REF].

Preliminary results

This section presents the preliminary results obtained with our multibanking approach.

We have applied our algorithm using the fkcc scripting tool [START_REF] Alias | fkcc: the Farkas Calculator[END_REF] on the motivating kernel, assuming a sequential execution and the arrays references to be accessed in parallel for each iteration. The synthesis results were obtained using VivadoHLS 2019.1, targeting a Kintex 7 FPGA (xc6k70t-fbv676-1). We transformed each array reference A[u( ì

i)] as Â[bank A (u( ì i)), offset A (u( ì i))]
. Then, we used the VivadoHLS array partitioning pragma on the bank dimension.

Our results are depicted on Table 1. The tool indicates the memory contention for the base kernel (base) are resolved on the transformed kernel (opt). The additional resources are due to the multibanking circuitry (steering logic). The overall latency is reduced from 79380 cycles to 31763 cycles. We suspect the speedup to be mitigated by the cost of the multibanking circuitry.

Conclusion

In this paper, we have outlined a unified, general HLS algorithm for multibanking, using the polyhedral model. Our approach reduces the overall size of memory banks, without hindering parallel memory accesses. Preliminary results encourage to pursue with this approach.

In the future, we plan to extend the field of experimental validation to more general linear algebra kernels under pipelining constraints. Also, we plan to explore how to minimize bank size separately as well as the trade-off surface/gain through a unified parametric formulation.

Figure 1 .

 1 Figure 1. Motivating example

Table 1 .

 1 Synthesis results

	Kernel	Latency BRAM DSP FF LUT
	conv2D-simple original	79380	0	0	471 981
	conv2D-simple opt	31763	0	0	1013 1692

Acknowledgments

This work has been partially funded by the fond recherche ENS-Lyon through the HLSIMU project.


Alec Sadler

Christophe Alias

Hugo Thievenaz

A Polyhedral Approach for Scalar Promotion

Keywords: compiler optimization, polyhedral model, array scalarization Compas

Memory accesses are a well known bottleneck whose impact might be mitigated by using properly the memory hierarchy until registers. In this paper, we address array scalarization, a technique to turn temporary arrays into a collection of scalar variables to be allocated to registers. We revisit array scalarization in the light of the recent advances of the polyhedral model, a general framework to design optimizing program transformations. We propose a general algorithm for array scalarization, ready to be plugged in a polyhedral compiler, among other passes. Our scalarization algorithm operates on the polyhedral intermediate representation.

In particular, our scalarization algorithm is parametrized by the program schedule possibly computed by a previous compilation pass. We rely on schedule-directed array contraction and we propose a loop tiling algorithm able to reduce the footprint down to the available amount of registers on the target architecture. Experimental results confirm the effectiveness and the efficiency of our approach.

Introduction

Using properly memory hierarchy until registers is of prime importance to improve the performances of a program, especially with the increasing gap between the peak rate of processing arithmetic units and the memory bandwidth. This trend in computer architecture, called the memory wall, boils down to the invention of memory hierarchy, and its counterpart in automatic code optimization. Array scalarization, or scalar promotion, [START_REF] Callahan | Improving register allocation for subscripted variables[END_REF][START_REF] Jiménez | Register tiling in nonrectangular iteration spaces[END_REF][START_REF] Renganarayana | Compact multi-dimensional kernel extraction for register tiling[END_REF][START_REF] Domagała | Register allocation and promotion through combined instruction scheduling and loop unrolling[END_REF] consists in transforming an array into a group of scalar variables, to be allocated to registers. In addition to reduce the memory traffic, hence the overall performances, it generally improves the precision of compiler optimizations, as dependences resolved through a register might be finely analyzed. In particular, register tiling [START_REF] Jiménez | Register tiling in nonrectangular iteration spaces[END_REF][START_REF] Domagała | Register allocation and promotion through combined instruction scheduling and loop unrolling[END_REF] splits a computation into blocks where register pressure make possible scalar promotion. Most of these approaches are monolithic, they are designed as endto-end optimizations without taking account of the scheduling constraints induced by previous compilation passes. In this paper, we focus on the polyhedral model [START_REF] Rajopadhye | On synthesizing systolic arrays from recurrence equations with linear dependencies[END_REF][START_REF] Quinton | The mapping of linear recurrence equations on regular arrays[END_REF][START_REF] Feautrier | Dataflow analysis of array and scalar references[END_REF][START_REF] Feautrier | Some efficient solutions to the affine scheduling problem. Part I. onedimensional time[END_REF][START_REF] Feautrier | Some efficient solutions to the affine scheduling problem. Part II. Multidimensional time[END_REF][START_REF] Feautrier | Polyhedron model[END_REF], a general framework to design loop transformations and data remapping for code optimization. Polyhedral compilers makes possible to reason about programs and their transformations thanks to a powerful geometric abstraction. We propose to rephrase array scalarization as a generic polyhedral compilation pass, parametrized by an input schedule -the result of a previous polyhedral compilation pass. We exploit array contraction [START_REF] Lefebvre | Automatic storage management parallel programs[END_REF][START_REF] Bhaskaracharya | Automatic storage optimization for arrays[END_REF] to expose array-level data reuse, and we propose an additional loop tiling algorithm to reduce the memory footprint of temporary arrays to a tunable constant size. At the end, we expose a minimum amount of scalar variables ready to be assigned a register. Specifically, we make the following contributions :

-We propose a general algorithm for array scalarization, ready to be plugged in a polyhedral compiler. In particular, our algorithm is parametrized by the program schedule which might be the result of a previous polyhedral pass. -Our transformation reduces as much as possible the code size for array scalarization and exposes directly the scalar variables to be put in distinct registers. This way, the work of the register allocator is dramatically reduced compared to seminal approaches for scalarization. -We propose a loop tiling algorithm able to reduce to footprint of some temporary arrays to a constant value. This algorithm is used on demand, when required. -We present a complete experimental validation showing the effectiveness and the efficiency of our approach. The remainder of this paper is structured as follows. Section 2 presents the required notions in polyhedral compilation. Section 3 describes our scalarization algorithm Section 4 presents our exprimental validation. Finally, Section 5 concludes this paper and draws research perpectives.

Preliminaries

This section outlines the concepts of polyhedral compilation used in this paper. In particular, we define the polyhedral intermediate representation of a program.

Polyhedral model

The polyhedral model [START_REF] Rajopadhye | On synthesizing systolic arrays from recurrence equations with linear dependencies[END_REF][START_REF] Quinton | The mapping of linear recurrence equations on regular arrays[END_REF][START_REF] Feautrier | Dataflow analysis of array and scalar references[END_REF][START_REF] Feautrier | Some efficient solutions to the affine scheduling problem. Part I. onedimensional time[END_REF][START_REF] Feautrier | Some efficient solutions to the affine scheduling problem. Part II. Multidimensional time[END_REF][START_REF] Feautrier | Polyhedron model[END_REF] is a general framework to design loop transformations, historically geared towards source-level automatic parallelization [START_REF] Feautrier | Polyhedron model[END_REF] and data locality improvement [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF]. It abstracts loop iterations as a union of convex polyhedra -hence the name -and data accesses as affine functions. This way, precise -iteration-level -compiler algorithms may be designed (dependence analysis [START_REF] Feautrier | Dataflow analysis of array and scalar references[END_REF], scheduling [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem. Part II. Multidimensional time[END_REF] or loop tiling [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF] to quote a few) . The polyhedral model manipulates program fragments consisting of nested for loops and conditionals manipulating arrays and scalar variables, such that loop bounds, conditions, and array access functions are affine expressions of surrounding loops counters and structure parameters (input sizes, e.g., N)). Thus, the control is static and may be analysed at compile-time. With polyhedral programs, each iteration of a loop nest is uniquely represented by the vector of enclosing loop counters ⃗ i. The execution of a program statement S at iteration ⃗ i is denoted by ⟨S, ⃗ i⟩ and is called an operation or an execution instance. The set D S of iteration vectors is called the iteration domain of S.

S = {(y, x) | 0 ≤ y < 2 ∧ 0 ≤ x < N}, D T = D U = {(y, x) | 2 ≤ y < N ∧ 0 ≤ x < N}
for the 2D blur filter presented later.

Polyhedral intermediate representation (IR)

In polyhedral compilers, the intermediate representation (IR) usually consists of a program P summarized as a set of statements S and their iteration domains D S , a schedule θ (typically the original sequential order), an optional tiling ϕ (a reindexing transformation which groups iteration into tiles to be executed atomically) and an optional array contraction function σ (as arrays might be remapped with an allocation function a[ ⃗ i] → a opt [σ a ( ⃗ i)], usually with a smaller footprint) .

Compas'2022 : Parallélisme / Architecture/ Système MIS/UPJV -Amiens France, 5-8 juillet 2022

f o r ( y = 0 ; y < 2 ; y + + ) { f o r ( x = 0 ; x<N; x ++) { S : b l u r x [ x ] [ y ] = i n [ x ] [ y ] + i n [ x + 1 ] [ y ] + i n [ x + 2 ] [ y ] ; } } f o r ( y = 2 ; y<N; y++) f o r ( x = 0 ; x<N; x ++) { T : b l u r x [ x ] [ y ] = i n [ x ] [ y ] + i n [ x + 1 ] [ y ] + i n [ x + 2 ] [ y ] ; U: out [ x ] [ y ] = b l u r x [ x ] [ y -2] + b l u r x [ x ] [ y -1] + b l u r x [ x ] [ y ] ; } }

Our Approach

This section outlines our approach on a running example.

Running example

We illustrate our scalarization approach on the 2D blur filter kernel depicted in Figure 1. The computation is divided into two steps. First, an horizontal filter (statements S and T ) is applied to the input picture in and stores the result into the array blurx. Then, a vertical filter (statement U) is applied to blurx and stores the final result to the array out. The whole might be seen as a producer/consumer through the temporary array blurx. Since blurx is a temporary array, it might be contracted and then scalarized, provided array contraction leads to a constant (nonparametrized by N) size. We point out that the array in cannot be scalarized directly in statement S, since it is not a temporary array. Nonetheless, a temporary version of in produced by a loop at the beginning of the program could perfectly be contracted and then scalarized, with a register pressure depending on the time shift between the producer and S. This preprocessing is used on some of our experimental results. Our scalarization algorithm is intended to be used in a polyhedral compilation chain. Hence a schedule might be imposed by the previous compilation steps. In the following, we consider two scheduling scenarios : the original execution order and a loop permutation.

Scenario 1. Original execution order

With the original schedule θ S (y, x) = (0, y, x), θ T (y, x) = (1, y, x, 0), θ U (y, x) = (1, y, x, 1), 3 iterations of x must be completed before the execution of U. Indeed, the second filter applied by U required three vertical cells of blurx, in particular the three first, for each x. Hence the allocation σ blurx (x, y) = (x mod N, y mod 3), with the non-constant (parametrized) footprint 3N. In that case, blurx cannot be directly scalarized. We propose to tile the iteration domain to limit the conflicting cells in the x direction. With that tiling, illustrated in Figure 1.(b), the footprint becomes 3h with h the tile size in the x direction. On a x86-64 machine with 14 general registers, we would set the tile size to h = [14/3] = 4.

Scenario 2. Loop permutation

We now assume that the outcome of the previous polyhedral compilation steps is a permutation of the loops x and y. This is described with the schedule θ S (y, x) = (0, x, y), θ T (y, x) = (1, x, y, 0), θ U (y, x) = (1, x, y, 1). In that case, we directly have the allocation σ blurx (x, y) = (x mod 1, y mod 3) with a constant footprint 3. Hence scalarization might be applied directly, without the need to apply further loop tiling.

Our algorithm

We now present our main algorithm (scalarization) and all its subroutines (tiling, unroll_factors and code_generation). They can be found in the appendix. Also, theorical proof of the algorithm might be found in the companion research report [START_REF] Sadler | A Polyhedral Approach for Scalar Promotion[END_REF].

Our main algorithm, Algorithm 1 inputs the intermediate representation of the program (P, θ) and an optional loop tiling ϕ. It outputs the polyhedral IR of the scalarized program (P out , θ out ), which might feed the next polyhedral compilation pass until the final code generation.

Array Contraction. First, we contract temporary arrays with the original schedule and tiling (step 2). Input and output arrays are ignored, since they cannot be contracted. As mentioned in section 3.1, input adn output arrays might be scalarized at the price of adding copy-in and copy-out code to temporary buffers. This might be addressed by a preprocessing polyhedral pass and will not be discussed further in this paper.

Adjusting Register Pressure. When the contraction fails to produce only temporary arrays with constant size (step 3) and no loop tiling is imposed, we try to tile the program in such a way the footprint is reduced to a constant, non-parametrized, size (step 5, Algorithm 2). Then, the arrays are recontracted (step 6). At this point, the tile size is adjusted so the product of σ modulos fits the available amount of registers. This is simply done by iterating step 6 on tile size ⃗ S from size (1, . . . , 1), incrementing each tile size component at each iteration, until the temporary arrays with constant contracted size all have a footprint (modulo product) tightly less than the available amount of registers. Arrays which still have a parametric size are skipped (step 8). When no array remains, meaning that the tiling failed to restrict at least one array to a constant size, our algorithm stops and returns the original program.

Code Emission. Finally, we scalarize the arrays with constant size. First, we compute the unrolling factors for the loops formally described by θ (step 13, Algorithm 3). These are the loops produced after the final polyhedral code generation for P under the scheduling constraint θ. Of course, we do not have syntactically these loops at this point of the polyhedral compilation, and we have to reason directly on θ. Then, we produce the polyhedral IR for the final scalarized program (step 14). We apply the unrolling (and our tiling ϕ when step 5 was required) with respect to θ and we generate the program statements with scalar variables to be allocated to registers.

Tiling Algorithm

We now describe our tiling procedure depicted in Algorithm 2. Our goal is to tile the program to bound the parametric terms of the array allocation σ. From now, we consider the running example, scenario 1. Recall that we obtained σ blurx (x, y) = (x mod N, y mod 3), hence the need to tile the iteration domain on the x direction to restrict the number of conflicting array

Compas'2022 : Parallélisme / Architecture/ Système MIS/UPJV -Amiens France, 5-8 juillet 2022 cells to a constant value. Actually, there is two notions of direction : a parametric direction in the array index domain, clearly identified : x, from which we deduce a parametric direction in the iteration domain, which happens to be the same, here. More precisely, given a statement S and an array reference a[u( ⃗ i)], we want to infer a variation ⃗ ∆ k in the iteration domain D S of S which incurs a variation in the direction ⃗ δ k (vector with 1 at position k, 0 elsewhere) in the direction k of the array index domain (here k = 1, δ 1 = (1, 0)). If σ a (⃗ c) = A⃗ c mod s( ⃗ N), this amounts to solve :

A • u( ⃗ ∆ k ) = ⃗ δ k
This affine equation is classically solved thanks to standard linear algebra techniques (lines 8 to 15). Note that Q -g denotes the generalized inverse of Q. The outcome is the set P S of directions ⃗ ∆ k of the iteration domain D S of statement S for which at least one reference a[u(i)] makes a step in a parametric direction ⃗ δ k according to σ a . Then, a tiling is computed (line 19) using the pluto algorithm [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF]. Finally, we keep only the hyperplanes going into a parametric direction. We point out that our algorithm will lead to a contraction of temporary arrays to a constant size if hyperplanes do not cross dependences hold by those arrays. Otherwise, a copy of sources should be kept along complete slices of the iteration domain. Note that the pluto algorithm tends to avoid that pitfall by pushing the resolution of dependences to innermost hyperplanes.

Experimental Results

This section presents our experimental results on several polyhedral programs.

Experimental setup

We have implemented our scalarization algorithm. The final code was generated using the iscc polyhedral code generator [START_REF] Verdoolaege | Counting affine calculator and applications[END_REF]. We have applied our algorithm to the following kernels :

-2D-blur-filter. Our running example, applying a 2D blur filter to an input.

fibonacci. This kernel generates the N first fibonacci numbers, and returns the last one.

pc-2d-interleaved. Producer/consumer throught a 2D array, where the consumer executes 2 iterations after the producer. pc-1d. Same as before for an array of one dimension pc-2d. This kernel applies a stencil pattern on a 2D array, with dependence vectors (1, 0) and (0, 1). cnn. Simple CNN with a convolutive layer followed a ReLU layer. -2mm. Multiplication of three matrices together (A × B × C).

gemm. BLAS kernel computing C := αA × B + βC. On the experiments, A and B where chosen as N × N matrices. poly. Multiplication of monovariate polynomials P and Q of degree N, represented by their array of coefficients. Kernels cnn, 2mm, gemm and poly were preprocessed to enable the contraction of input/output arrays, along the lines described in Section 3.1. Benchmarks were done by executing both the default and scalar program with different array sizes. Executions were made on a single-x86_64 intel CPU, with 14 registers. The CPU features 4 cores, with 64KB of cache L1, 512 KB of cache L2 and 4MB of cache L3. Compilation was done with GCC11 -O0 to measure exactly the impact of our optimization.

Results

Figure 2 depicts our results. Every graph shows runtimes for both default and scalar version, as well as the speed-up, for multiple array size. For every example, similar behaviours can be observed, such as cache effects when the memory footprint gets large enough. Cache memory becomes saturated, and another phase of the curve starts.

For almost example, speed quite a lot the program. On 2D blur filter, it is interesting to note that the scalarized version show a bigger growing rate compared to the default version, which with the data size, unlike fibonacci, pc-2D-interleaved, cnn, 2mm, and gemm, which exhibits a constant speed-up. On pc-2d, we observe instabilities on both curves with the ratio slighty going down. On gemm and the poor performances are explained by the number of conditional branches in the target program to handle corner-cases, that we to cause many branch misprediction. This is the main weakness of direct polyhedral code generation.

Conclusion

In this paper, we have proposed a complete algorithm for array scalarization as a composable pass in a polyhedral compiler. Our algorithms features a loop tiling to reschedule the input kernel so the footprint of the temporary arrays may be tuned to fit into the registers of the target architecture. We have also provided a complete correctness proof of our approach, completed with an experimental validation on a set of representative polyhedral kernels used in linear algebra and signal processing applications.

In the future, we would like to investigate how to improve the polyhedral code generation to reduce the conditional branches, which bound unexpectedly our speed-ups on some kernels. MIS/UPJV -Amiens France, 5-8 juillet 2022

Figure 1 .

 1 (b) provides the iteration domains D

FIGURE 1 -

 1 FIGURE 1 -Running example : 2D Blur filter

Compas' 2022 :

 2022 FIGURE 2 -Experimental results

Algorithm 2 : 2 foreach 3 4 P S ←-∅ 5 foreach 6 7 begin 8 if u is non-singular then 9 12 Write 13 ⃗ 14 ⟨⃗ e 1 ,

 2234567891213141 TILING Data: Program (P, θ), allocation σ Result: Scalarization-aware tiling ϕ 1 begin reference S : . . . a[u(⃗ i)] . . . do Write σ a (⃗ c) = A⃗ c mod s( ⃗ N) k s.t. s( ⃗ N)[k] is parametrized do Add a basis of ⃗ ∆ k s.t. A • u( ⃗ ∆ k ) = ⃗ δ k : Add ⃗ ∆ k = u -1 • A -1 ( ⃗ δ k ) to P S A • u( ⃗ ∆ k ) = Q ⃗ ∆ k + ⃗ r / * get a solution * / ∆ 0 ←-Q -g ( ⃗ δ k -⃗ r) / *add a solution basis * / . . . , ⃗ e p ⟩ ←ker Q 15 Add each ⃗ e i + ⃗ ∆ 0 to P S ( ⃗ i) = T ⃗ i + ⃗ u 23 foreach line vector ⃗ ℓ j of T do 24 if ⃗ ℓ j • ⃗ ∆ ̸ = 0 for some ∆ ∈ P S

Compas'2022 : Parallélisme / Architecture/ Système MIS/UPJV -Amiens France, 5-8 juillet 2022 Appendices 2D blur filter, scenario 1, code generated void b l u r _ k e r n e l ( double * in , double * out , i n t N) { r e g i s t e r double blurx_00 , blurx_01 , blurx_02 , blurx_10 , blurx_11 , blurx_12 , blurx_20 , blurx_21 , blurx_22 , blurx_30 , blurx_31 , b l u r x _ 3 2 ;

f o r ( i n t c0 = 0 ; c0 <= f l o o r d (N -1 , 4 ) ; c0 += 1 ) { b l u r x _ 0 0