Konrad Trifunovic

Dr Albert Cohen

M John Cavazos

M Philippe Clauss

M Albert Cohen

M Ayal Zaks

In Partial Fulfillment of the Requirements for the Degree of DOCTOR OF PHILOSOPHY Discipline: Computing science

Introduction

The responsibility of a compiler is to translate the programs written in a language suitable for use by human programmer into the machine language of the target machine. Naturally, the compiler is required to produce the correct translation to the machine code.

In 1944 John von Neumann has laid the foundation to almost all widely used architectures today -that of decoupling the program memory from data memory. This also implied a sequential execution model, in which the machine instructions are executed sequentially, changing the global state of the memory through stores and loads.

Besides producing the correct code, compilers were enhanced with optimizations, so that the translated machine code performance could match the performance of the hand-crafted one. Contributions in scalar optimizations, interprocedural analysis, instruction scheduling and register allocation mainly achieved this goal.

The major driving forces for the sustained computing performance increase were the advancements in instruction level parallelism (ILP) and ever increasing clock rates of uniprocessors -according to the Moore's law. This obviated the need for advanced compiler optimizations, since the performance boost could have been obtained by running non-modified sequential codes on the new machines. Advanced compiler optimizations were considered as a domain of specialized, scientific computing community, specializing in the expensive parallel and vector machines programmed mainly in FORTRAN.

End of uniprocessor era

The trends have changed in the early 2000s, since it was noticed that increasing the ILP brought diminishing returns, altogether with physical limitations and power dissipation that have put an end to the increase in the clock-speeds of uniprocessors. The capacity to increase the number of transistors has since then been used in another direction: providing multiple cores on the single chip.

The era of multi-core

Multi-core processors are now in widespread use in almost all areas of the computing: in desktop or laptop machines, accelerators like the Cell Broadband Engine, GPUs, GPGPUs and in mobile embedded devices implemented as MP-SoCs (Multiprocessor System-on-Chip).

In the scientific computing field, where specialized multiprocessors have been used since 1960's, the commodity multicore and GPGPUs are entering very rapidly, as main building blocks of the computational clusters.

INTRODUCTION

Optimizing compilers are again in demand

To harness the power of multiple cores and complex memory hierarchies, the need for powerful compiler optimizations is now again in high demand. Compilers are now required to efficiently map the semantics of the, mostly imperative, high-level programming languages into the efficient use of target architecture resources [START_REF] Hall | Compiler research: the next 50 years[END_REF].

Increasing the number of cores on a chip does not improve the performance of the legacy, singlethreaded applications. What is more, the sequential programming paradigm, due to the simplicity of the von Neumann computational model, is still used to program the new applications. Those applications do not benefit from the performance improvements of the new architectures.

In order to utilize the peak performance of the current multi-core architectures, the optimizing compiler is responsible for finding the intrinsic parallelism within the computations of the source program. Also, the memory hierarchy has to be utilized carefully, in order to narrow down the gap between latencies of the main memory and the computing core -the memory wall problem.

Semantic gap

There is a notorious semantic gap between the semantics of the currently used imperative programming languages and the target platform execution model. Nevertheless, while more than 200 parallel programming paradigms and environments have been proposed, few of them are in a widespread usethe majority of the applications are still written in imperative sequential languages such as C, C++ or FORTRAN.

Most current applications are written using a single thread of execution. There is a good reason for this -writing single-threaded code is intuitive and practical: one has to manage only a single state space, it is easily analysable, the proof of correctness could be obtained in an inductive way, the execution trace is predictable and reproducible.

On the other hand, writing the parallel programs introduces the problems such as communication, synchronization and data transfers. Writing hand-crafted parallel code is costly, error-prone and, mostly importantly, impossible to debug, since the race conditions in parallel program are generally non-reproducible. Hand-crafted parallel programs are non-portable, since the synchronization mechanisms, communication costs, scheduling policies and other issues are highly platform specific. The non-portability manifests itself in two ways: platform non-portability and performance non-portability.

Automatic program transformations

Very promising solution to the mentioned semantic gap problem is to enable the automatic transformation of sequential code into the form that best fits the target execution architecture.

The task of the compiler is now to understand the semantics of the input sequential source program and to perform the transformation, while preserving the semantics, into the form that best fits the target platform execution model, e.g., parallel threads, vectorized instructions, memory accesses clustered into blocks. This is a challenging task, since all the responsibility of translating the semantics into performance now relies on the compiler.

Loop Transformations

Focusing the program transformations on loop nests is particularly interesting, since the majority of compute-intensive applications spend most of their execution time in loops. Loop transformations have a long history and they have been used in restructuring compilers for enabling automatic parallelization and vectorization of the scientific FORTRAN codes.

Loop transformations potentially have a very powerful performance impact, since they might enable parallelism, vectorization or memory locality. The major challenge for efficiently employing the loop transformations is to find the best composition of loop transformations and to predict the performance impact of those transformations.

The Polyhedral Model

The polyhedral model is a powerful algebraic abstraction for reasoning about loop transformations expressed as schedules. The polyhedral model might be viewed as an abstract interpretation framework for reasoning about dynamic program execution statically -at compile time. The foundations of the polyhedral model have their roots in operational research and systolic array design.

An execution trace of the sequential program is modelled as an ordered sequence of statement instances. It is assumed that each statement of the input program might be executed multiple times in a loop, giving rise to statement instances -one instance for each execution of the loop.

Single-Shot Optimization Approach

Traditionally, the program transformation using the polyhedral model is done in a sequence of three steps: [START_REF] Agakov | Using machine learning to focus iterative optimization[END_REF] getting the program into the polyhedral model, (2) computing the scheduling transformation expressed in the polyhedral model, and (3) generating the transformed code from the polyhedral model.

The crucial step is the scheduling step that provides the actual program transformation. The automatic scheduling algorithms [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part I: one dimensional time[END_REF][START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF] based on integer linear programming exist, and they are based on the best effort heuristics.

The advantage of the best effort scheduling heuristic is that the transformation is computed within a single step. Hopefully, the best possible transformation within the space of legal transformations is chosen.

The space of legal (loop) transformations is huge, and selecting the best transformation within this space is still an open problem. The best effort heuristics rely on a limited performance predicting costfunctions, abstracting away the architectural details of the target platform. The optimal transformation might easily be mispredicted.

Iterative Optimization Approach

Due to the intrinsic difficulty of selecting the best program transformation within a single step, a viable alternative approach is the iterative feedback-directed compilation. Iterative optimization resorts to testing the different program transformations. For each transformation, the code is generated. The code is then executed and the runtime impact of the transformation is evaluated.

A carefully crafted heuristic for generating the multiple program transformations is the key component of this approach. Experimental evidence confirms that iterative approach can be much better at adapting the semantics of the program to the execution platform than the single-shot, model based heuristics. Obviously -by trying more times, there is a better chance of getting the optimal transformation.

While the iterative approach might give better results in terms of the output program performance, it has one fundamental drawback: the whole process might be very time consuming. It is not uncommon for the iterative approach to take hours or days to converge to the best solution.

The iterative approach breaks the traditional view of the compilation: the compiler is now a component of the optimizing feedback loop, and not a self-standing tool that transforms the input source code to the binary. While the iterative compilation is a good approach to specialized program optimization, the two characteristics of this approach -that of having an unpredictable running time and reliance on 1. INTRODUCTION the actual generated program execution -preclude it from being incorporated into the general purpose compiler.

Problem Statement and Proposed Solution

Our view is that the single-shot, best effort heuristics for loop optimizationsare not precise enough to model the complex performance characteristics of the modern architectures. On the other hand, the iterative compilation approach, due to its unpredictable runtime and reliance on experimental program evaluation is not directly employable in the general purpose compiler. The problem statement could be summarized as follows:

We want to keep the effectiveness of the model-based transformation approach, while bridging the gap between transformation prediction precision between single-shot and iterative approaches.

In order to achieve this, we propose a solution in the form of a new, search-based transformation approach, that could be summarized as follows:

A new transformation search approach is based on the controlled enumeration of the transformation candidates, each transformation candidate being evaluated at the compilation time, according to a precise cost function based on the target machine model.

We postulate that the new proposed solution brings us closer to solving the stated problem. The new approach takes the benefits of both model-based and iterative approaches. The effectiveness is achieved by having an expressive but limited search space, whose complexity could be controlled. The fact that the cost function is evaluated at compile-time, and not through experimental evaluation, brings the predictability of the running time of the transformation search, allowing it to be incorporated into the general purpose compiler.

The precision is achieved by the fact that we use a machine dependent cost function that directly evaluates the execution cost on a specified machine. This function is much more precise than the simplistic linear cost functions, that are used in the current model-based approaches. Though, this precision comes with a price: non-linear cost functions cannot be optimized using linear programming machinery, as it is done in the current model-based approaches -our solution uses a search based strategy to evaluate those non-linear cost functions.

Contributions

The main contributions of this thesis are to provide algorithms, methodologies and the compilation framework that enable our new search-based strategy. The relevant contributions are detailed in the respective chapters.

Our contributions are built on top of the polyhedral model, which serves as the theoretical foundation of the work. The polyhedral model is used to express the data dependences, the cost-model functions and schedules.

Contrary to the current source-to-source polyhedral frameworks, we have investigated the direct polyhedral compilation on the low-level three address code. While being a challenging task, this low-level of abstraction is the foundation for our precise performance cost-modelling that is out of reach of the current source-to-source polyhedral compilers.

We take a fresh look at the problem of finding the program transformations. Contrary to the current linear cost-function based scheduling algorithms, we propose a search strategy that evaluates the complex cost-model function on the discrete set of legal transformations represented as decision diagrams. By doing so, we take the benefits of the iterative optimization, while not requiring the generation of the transformed programs and measuring their actual runtime.

Thesis overview

The thesis is organized in 3 parts and a conclusion. Contributions and related work are provided in each chapter.

The first part of the thesis gives the background material, starting from the basic mathematical definitions. The polyhedral model is defined, together with the detailed discussion of this abstract representation. Later, we discuss the topics of representing the program semantics, polyhedral transformations and legality of those transformations. The mathematical background is included so that the manuscript is self-contained -a reader familiar with the topic of linear optimization might directly skip to the sections devoted to the polyhedral program model.

The second part of the thesis is a detailed description of our three-address-code polyhedral compilation framework -GRAPHITE. The crucial design decisions were explained, together with our contributions in not so well investigated topics. A novel approach for efficiently handling the restrictions imposed by memory-based dependences is explained as well.

The third part is the core subject of the thesis. It is dedicated to the description of the search-based transformation methodology. A detailed explanation of our precise, machine-level cost-model is provided. The final part of the thesis is our proposal for the new search strategy, starting with a review of the most important related work in this field. We show how our search strategy is built on an enumeration of the discrete sets of legal transformations, represented as decision diagrams.

Our techniques are implemented in GRAPHITE polyhedral framework, that is a part of the widely used GCC compiler.

Part I Basic Concepts

Chapter 2

Background

The polyhedral model representation is based on linear algebra, linear programming and convex optimization theory.

In this chapter we present the polyhedral model and the mathematical notations we use in the rest of the dissertation. The basic linear algebra and convex optimization background is given as well.

The Section 2.1 gives the necessary mathematical background and it is expected to be self-contained. A complete coverage of the material could be found in [START_REF] Schrijver | Theory of linear and integer programming[END_REF]. The Section 2.2 gives an introduction to the polyhedral model and places it in the context of the compilers. The Section 2.3 goes into the details of the polyhedral representation and covers the basic components of this representation. An extensive coverage of the literature and an historical overview is given in Section 2.4.

Mathematical background

We denote by R, Z and Q a field of real, integer and rational numbers respectively. A column vector w is denoted by w = (w 1 , . . . , w m ) T .

Given two vectors a ∈ R n and b ∈ R m we will denote a concatenation of two vectors: c = (a, b) T ∈ R n+m . A concatenated vector c if formed by taking n components from a vector a and m components from a vector b: c = (c 1 , . . . , c n , c n+1 , . . . , c n+m ) T = (a 1 , . . . , a n , b 1 , . . . , b m ) T

The vector concatenation notation will be used very often to distinguish several parts of a vector. In a similar way, given a matrix A ∈ R k×n and a matrix B ∈ R k×k , we can compose two matrices to get: C = (A|B). A composed matrix C is obtained by taking n columns from the matrix A and m columns form the matrix B:

C = (C •,1 . . .C •,m+n )
where C •,i denotes i-th column of the matrix C. In the same way C k,• denotes k-th rows of the matrix C. Definition 2.1.5 (Convex polyhedron). The set P ⊆ R n is convex polyhedron if it can be represented as an intersection of a finite number of affine half-spaces of R n . Each half-space is called a face of the polyhedron. The set of all affine inequalities representing faces of the polyhedron can be compactly represented by an matrix A and a vector b:

P = {x ∈ R n |Ax + b ≥ 0}
where A ∈ R m×n and b ∈ R m , m being the number of affine inequalities.

Definition 2.1.6 (Parametric polyhedron). The set P (g) ⊆ R n is parametric polyhedron, parametrized by a vector g ∈ R k , if it is a convex polyhedron that can be represented as:

P (g) = {x ∈ R n |Ax + Bg + b ≥ 0}
where A ∈ R m×n , B ∈ R m×k and b ∈ R m . The number of parameters in a parameter vector is k and m is the number of affine inequalities.

Definition 2.1.7 (Polytope, bounded polyhedron). A polyhedron P ⊆ R n is called a bounded polyhedron if there exists a λ ∈ R + such that the polyhedron P is contained in a bounding box:

P ⊆ {x ∈ R n | -λ ≤ x j ≤ λ, 1 ≤ j ≤ n}
For the program analysis problems that we are interested in, we will consider only integer solutions, so we restrict all vectors to contain only integer components: x ∈ Z n . Also, A ∈ Z n×m and b ∈ Z n . For more details please refer to [START_REF] Schrijver | Theory of linear and integer programming[END_REF]. An efficient implementation of the presented concepts is given by [START_REF] Bagnara | The parma polyhedra library: Toward a complete set of numerical abstractions for the analysis and verification of hardware and software systems[END_REF].

The polyhedral model

Classical compiler internal representations are based on a syntax and semantics of imperative programs: abstract syntax trees, control-flow graph, def-reach chains, SSA form, three-address-code, basicblocks [3,[START_REF] Muchnick | Advanced Compiler Design and Implementation[END_REF]. Those representations capture a syntax, control and data dependences in a program.

In such reduced representations of the dynamic execution trace, each statement of a high-level program occurs only once, even if it is executed many times (e.g., when enclosed within a loop). Representing a program this way is not convenient for aggressive loop optimizations which operate at the level of dynamic statement instances.

Compile-time constraints and lack of adequate algebraic representation of loop nest semantics prevent traditional compilers from adapting the schedule of statement instances of a program to best exploit architecture resources. For example, compilers typically cannot apply loop transformations if data dependences are non-uniform [START_REF] Wolfe | High performance compilers for parallel computing[END_REF] or simply because profitability is too unpredictable.

A well known alternative approach, facilitating complex loop transformations, represents programs in the polyhedral model. This model is a flexible and expressive representation for loop nests with statically predictable control flow. Such loop nests, amenable to algebraic representation, are called static control parts (SCoP) [START_REF] Feautrier | Array expansion[END_REF][START_REF] Girbal | Semiautomatic composition of loop transformations for deep parallelism and memory hierarchies[END_REF].

The polyhedral model captures an exact execution semantics of a statically predictable control flow program. We can think of the polyhedral model as an abstract interpretation framework [START_REF] Cousot | Automatic discovery of linear restraints among variables of a program[END_REF] in which we can precisely reason about a dynamic execution of a program at a compilation time.

Static Control Parts

Providing a static analysis of a dynamic execution trace of an arbitrary program is undecidable problem [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part I: one dimensional time[END_REF]. In order to apply the polyhedral model to real computing problems, we have to restrict the scope of the analysis to program parts with statically determinable control flow. Those program parts are contained in affine loops. Definition 2.2.1 (Affine loops). Affine loop nests are nests of perfectly or imperfectly 1 nested loops with constant strides, having their loop bounds and memory accesses expressed as affine functions of outer loop induction variables and global parameters.

The scope of the polyhedral program analysis and manipulation are Static Control Parts containing sequences of affine loop nests and conditionals with boolean expressions of affine inequalities. This class of programs includes non-perfectly nested loops, non-rectangular loops and conditionals with boolean expressions of affine inequalities [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part I: one dimensional time[END_REF]. Global parameters are symbolic variables that are invariant inside a loop.

We will use an abbreviation SCoP (Static Control Part) to denote the parts of a program amenable for polyhedral analysis [START_REF] Girbal | Semiautomatic composition of loop transformations for deep parallelism and memory hierarchies[END_REF]. In Chapter 4 we will show an algorithmic approach that is used in our tool to extract SCoP parts from an arbitrary imperative program. An example of the program part that fits in SCoP category is shown in Figure 2.1. The global parameters for this affine loop nest are N and K.

Polyhedral representation

The polyhedral program representation is based on the concept of dynamic statement instances, which is the basic atomic unit of a program execution. The set of those instances is modelled in a linear algebraic framework, which we are going to present in the subsequent parts of this section.

We will define a concept of the polyhedral statement and the three components of the polyhedral representation:

-Iteration domains -capturing the set of dynamic statement instances -Data access functions -capturing the memory access patterns expressed as affine functions of loop induction variables -Scheduling functions -capturing the relative execution order of the statement instances Building on those primitives, we will show in Section 3.2 how to statically extract the semantics of the SCoP programs. Also, we will show how to express the transformations in the polyhedral representation in Section 3.3.

Polyhedral statements

The polyhedral model is an abstraction that expresses the program execution trace as an execution of the dynamic polyhedral statements.

This abstraction is built on top of some internal representation used in a compiler. Depending on the level of abstraction, one polyhedral statement could correspond to:

1. syntactic statement -a source-level statement in a source program.

basic block of a low-level internal compiler representation 3. three-address code instruction

The polyhedral model is mostly used in the context of source-to-source compilers [START_REF] Girbal | Semiautomatic composition of loop transformations for deep parallelism and memory hierarchies[END_REF][START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF][START_REF] Kelly | Optimization within a unified transformation framework[END_REF][START_REF] Meister | Productivity via automatic code generation for pgas platforms with the r-stream compiler[END_REF][START_REF] Irigoin | Semantical interprocedural parallelization: An overview of the pips project[END_REF]. Thus, the first case, where a polyhedral statement directly corresponds to a source-level statement in a source program, is the most common. It is also the most natural correspondence for presentation purposes.

In this dissertation, we investigate the application of the polyhedral model to the low-level, three address code compilation. Thus, in Chapter 4, we will also use an abstraction level where a polyhedral statement corresponds to a basic block of low-level internal compiler representation.

In the rest of this dissertation we will implicitly use a term statement to refer to an abstract polyhedral statement of the polyhedral model. Depending on the context, the term statement might refer to a polyhedral statement corresponding to the source-level statement2 , or to the basic block in a low-level internal representation. We will try to make this distinction clear in the following chapters.

f o r ( v = 0 ; v < N; v ++) f o r ( h = 0 ; h < N; h ++) { S1 : s = 0 ; f o r ( i = 0 ; i < K ; i ++) f o r ( j = 0 ; j < K; j ++) S2 : s += image [ v+ i ] [ h+ j ] * f i l t e r [ i ] [ j ] ; S3 : o u t [ v ] [ h ] = s >> f a c t o r ; } Figure 2.1 -Main loop kernel in Convolve

Iteration domains

Once we have defined an abstract polyhedral statement, we are interested in analyzing dynamic statement instances. An execution of an abstract sequential program might be seen as a totally-ordered interleaving of statement instances. Each statement might be executed multiple times in a loop. Definition 2.3.1 (Iteration vector). The iteration vector i of a statement S consists of values of the induction variables of all loops surrounding the statement.

Each statement S is executed once for each possible value of a surrounding loop induction variables We denote a statement by S and a statement instance by (S, i), where i is an iteration vector of a statement.

If a statement S belongs to a SCoP then the set of all iteration vectors i relevant for S can be represented by a parametric polytope:

D S (g) = i ∈ Z n | A • i + B • g + b ≥ 0
which is called the iteration domain of S, where g is the vector of global parameters whose dimension is n g . Global parameters are invariants inside the SCoP, but their values are not known at compile time (parameters representing loop bounds for example).

We will combine matrices A, B and a column vector b in one matrix D = (A | B | b). We will also concatenate vectors i, g to form an homogeneous column vector (i, g, 1) T , where last component is a scalar constant 1. We can then provide a more compact representation of a polyhedron:

D S (g) = i ∈ Z n | D • (i, g, 1) T ≥ 0
Let us consider a computational kernel given in Figure 2.1. We will name the components of the iteration vector according to loop induction variables and global parameters given in a source code. The resulting iteration domain polyhedron for statement S 1 is then represented as:

    1 0 00 0 -1 0 10 -1 0 1 00 0 0 -1 10 -1           v h N K 1       ≥ 0 v ≥ 0 v ≤ N -1 h ≥ 0 h ≤ N -1
We can represent the same iteration domain polyhedron in a more concise way:

D S 1 = (v, h) T ∈ Z 2 | 0 ≤ v, h ≤ N -1
The iteration domains for statements S 2 and S 3 are represented as follows:

D S 2 = (v, h, i, j) T ∈ Z 4 | 0 ≤ v, h ≤ N -1 ∧ 0 ≤ i, j ≤ K -1 D S 3 = (v, h) T ∈ Z 2 | 0 ≤ v, h ≤ N -1

Data access functions

Given an affine program, all accesses to memory locations must be expressed through array references. We will treat accesses to scalar values as accesses to single-element arrays A[0].

A data reference to an array A ∈ A is denoted by R = A, f , where f is a subscript function. Together, they describe an array access of the form A[ f (i)]. A subscript function is a function that maps an iteration vector i to an array subscripts:

f (i) = F • (i, g, 1) T
where F ∈ Z m×(n+n g ) is a matrix containing m rows -m being the number of array subscripts, n and n g being the iteration vector and global parameters vector sizes respectively.

In affine programs all subscript functions must be expressed as affine functions. A given statement S can contain multiple data references. A set of write and read data references of a statement S is denoted as W S and R S respectively.

As an example, let us take the statement S 2 from the example in Figure 2.1:

S2: s += image[v+i][h+j] * filter[i][j];
The set of write and read references is the following:

W S 2 = { s, f (i) = 0 } R S 2 = { s, f (i) = 0 , image, f 1 , filter, f 2 }
The access functions that map the iteration vectors to array subscripts have the following form:

f 1 ((v, h, i, j) T ) = v + i h + j f 2 ((v, h, i, j) T ) = i j

Schedules

An execution trace of a sequential program could be described as a total order relation on a set of statement instances {(S, i) : i ∈ D S }. In order to impose a total execution order on a set of statement instances, each statement instance (S, i) is assigned a timestamp denoted by a vector t.

The t vector represent a multidimensional time. In order to impose a temporal ordering of time vectors, an ordering relation between two time vectors has to be defined. We use a lexicographic order between two vectors to impose their temporal order. Given two vectors a = (a 1 , . . . , a n ) T and b = (b 1 , . . . , b n ) T , we define a lexicographic order relation a ≪ b as follows:

(a 1 , . . . , a n ) T ≪ (b 1 , . . . , b n ) T ⇐⇒ ∃i : 1 ≤ i < n, a i = b i ∧ a i+1 < b i+1
In other words, there exist an integer i such that all vector components up to i are equal, while the vector component at the i + 1 position of the vector b is greater than respective component of the a vector. The idea of using multidimensional timestamp vectors was proposed by Feautrier [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part II: multidimensional time[END_REF] and later by Kelly and Pugh [START_REF] Kelly | A unifying framework for iteration reordering transformations[END_REF].

For each polyhedral statement S we define an affine scheduling function as a mapping between an iteration vector i and a timestamp vector t:

t = θ S (i) = Θ S • (i, g, 1) T
where Θ S ∈ Z p×(n+n g +1) is a scheduling matrix having p rows, where p is a dimensionality of a timestamp vector t.

Once we have defined the scheduling function and the lexicographic ordering of the timestamp vectors, we can formally define an execution order of statement instances: Definition 2.3.2 (Execution order). A given statement instance (S, i) is executed before a statement instance (S ′ , i ′ ) if and only if:

θ S (i) ≪ θ S ′ (i ′ )

Related work and an historical overview

The polyhedral model is well-established theoretical foundation for reasoning about program semantics and transformations of static-control programs.

The seminal work of Karp, Miller and Winograd [START_REF] Karp | The organization of computations for uniform recurrence equations[END_REF] on scheduling the computations of uniform recurrence equations can be considered as an origin of the polyhedral model theory. It introduced many concepts such as dependences, dependence graph, scheduling, iteration vectors.

Later work on automatic systolic array synthesis by Rajopadhye introduced the terms of timing function, allocation function, affine mappings and polyhedral domains. The PolyLib library for performing operations on rational polyhedra represented as a system of affine inequalities was conceived by Le Verge and Wilde [START_REF] Wilde | A library for doing polyhedral operations[END_REF][START_REF] Loechner | Parameterized polyhedra and their vertices[END_REF].

Different authors have used different notations to refer to the same concepts relating to polyhedral compilation. The most recent attempt to unify the notation and introduce the canonical form was that of Girbal [START_REF] Girbal | Semiautomatic composition of loop transformations for deep parallelism and memory hierarchies[END_REF]. When referring to the polyhedral model we are referring to the notation introduced in this chapter, which is mainly based on the notation coming from the work of Girbal [START_REF] Girbal | Semiautomatic composition of loop transformations for deep parallelism and memory hierarchies[END_REF].

A group at INRIA [START_REF] Girbal | Semiautomatic composition of loop transformations for deep parallelism and memory hierarchies[END_REF] has demonstrated that one can consider the polyhedral model as a fullycapable internal representation of the compiler. WRaP-IT [START_REF]Optimisation d'applications -Composition de transformations de programme: modèle et outils[END_REF] project has shown an internal representation based on the polyhedral model incorporated into production quality compiler Open64 [START_REF] Coarfa | Open-source compiler technology for source-to-source optimization[END_REF].

As the last note, we remind that the polyhedral analysis is restricted to affine programs -programs that have static control and data access patterns described by the system of affine constraints. There are works that try to widen the scope of applicability of the polyhedral model to programs going beyond static-control [START_REF] Benabderrahmane | The polyhedral model is more widely applicable than you think[END_REF][START_REF] Barthou | Fuzzy array dataflow analysis[END_REF]. They mainly rely on the conservative approximations of the program semantics.

Summary

In this chapter we have given a self-contained mathematical background that forms a foundation of the polyhedral model. We have formally defined the affine functions, affine hyperplanes, polyhedra and polytopes. We have shown the place of the polyhedral model in the general compilation context. We have defined the class of the affine programs for which the polyhedral model representation could be obtained.

In Section 2.3 the basic building block of the model are defined : 1) iteration domains, 2) data access functions and 3) schedules.

Given the basic components of the polyhedral model, we define the way to express the multidimensional program transformations and to assess their legality by introducing the data dependence concepts in the next chapter.

Chapter 3

Program Transformations in the Polyhedral Model

Canonical Form of the Scheduling Functions

We have given the definition of the scheduling function in Section 2.3.4, defining it as as a function that expresses a mapping from statement iteration vectors to timestamp vectors. While the representation of the iteration domains and access functions (shown in Sections 2.3.2 and 2.3.3) follows straightforwardly from their definition, the question of representing the scheduling functions leaves many degrees of freedom.

One might restrict the scheduling functions to one-dimensional timestamps. But as Feautrier [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part II: multidimensional time[END_REF] has shown, one-dimensional timestamps cannot express the schedule of some class of the programs. Feautrier [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part II: multidimensional time[END_REF] has also shown that multidimensional schedules are expressive enough to encompass the full class of affine programs.

Within the class of multidimensional schedules, the following question arises: how many dimensions should a timestamp vector contain? Feautrier [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part II: multidimensional time[END_REF] provides a scheduling algorithm that minimizes the necessary number of dimensions of the timestamp vector. Bondhugula [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF] provides a scheduling approach that requires at least d scheduling dimensions, d being the maximal loop depth of the statements within SCoP.

For the compiler construction purposes, it is desirable to have a canonical form of the scheduling functions. Girbal et al. [START_REF] Cohen | A polyhedral approach to ease the composition of program transformations[END_REF][START_REF] Bastoul | Putting polyhedral transformations to work[END_REF][START_REF] Cohen | Facilitating the search for compositions of program transformations[END_REF] have proposed the canonical form of the scheduling functions that unifies the format of the scheduling matrices representing those functions.

We will discuss the motivation for providing the canonical form of the scheduling matrices, then we will define the form and discuss some of its properties.

Motivation

Given the scheduling function θ S (i) represented as θ S (i) = Θ S • (i, g, 1) T , the goal is to provide a canonical form of the scheduling matrix Θ S .

Even if a matrix Θ S ∈ Z d×(n+n g +1) has a full rank (requiring d ≥ n), a direct mapping from an iteration vector i to a timestamp t is not enough to distinguish instances of two different statements. Let us take an example in Figure 3.1. Given the two statements S 1 and S 2 , if we simply take a full-rank identity matrix as a scheduling matrix for both statements, we would get:

θ S 1 ((i 1 , i 2 ) T ) = (i 1 , i 2 ) T θ S 2 ((i ′ 1 , i ′ 2 ) T ) = (i ′ 1 , i ′ 2 ) T f o r ( i 1 = 0 ; i 1 <N ; i 1 ++) { f o r ( i 2 = 1 ; i 2 < N ; i 2 ++) { S1 : X[ i 1 ] [ i 2 ] = X[ i 1 ] [ i 2 ] -X[ i 1 ] [ i 2 -1] * A[ i 1 ] [ i 2 ] / B[ i 1 ] [ i 2 -1 ] ; S2 : B[ i 1 ] [ i 2 ] = B[ i 1 ] [ i 2 ] -A[ i 1 ] [ i 2 ] * A[ i 1 ] [ i 2 ] / B[ i 1 ] [ i 2 -1 ] ; } } Figure 3.

-An excerpt from ADI numerical kernel

There is no way to express that a given instance of statement (S 1 , (i 1 , i 2 ) T ) happens before an instance

(S 2 , (i ′ 1 , i ′ 2 ) T ) for (i 1 , i 2 ) T = (i ′ 1 , i ′ 2 )
T , since both statement instances are mapped to the same timestamp vector.

An easy solution to this problem is to extend the timestamp vectors with constant entries that disambiguate different static statements in an original program. Those constant entries corresponds to relative execution order of statement instances inside a common loop nest. Now we have:

θ S 1 ((i 1 , i 2 ) T ) = (i 1 , i 2 , 0) T θ S 2 ((i ′ 1 , i ′ 2 ) T ) = (i ′ 1 , i ′ 2 , 1) T and it is always true that (i 1 , i 2 , 0) T ≪ (i ′ 1 , i ′ 2 , 1) T for (i 1 , i 2 ) T = (i ′ 1 , i ′ 2 )
T , since the lexicographic order is disambiguated at the last component of the timestamp.

Since the polyhedral model allows non-perfectly nested loops, a similar problem might occur when different statements, say S 1 and S 2 have different loop depths.

Aforementioned problems were limiting factors in early works on polyhedral program transformations [START_REF] Ancourt | Scanning polyhedra with do loops[END_REF][START_REF] Li | A singular loop transformation framework based on non-singular matrices[END_REF][START_REF] Ramanujam | Beyond unimodular transformations[END_REF][START_REF] Xue | Automating non-unimodular loop transformations for massive parallelism[END_REF]. This was mainly due to the limitations of code generation algorithms that required the scheduling matrices to be unimodular [START_REF] Ancourt | Scanning polyhedra with do loops[END_REF] or at least invertible [START_REF] Li | A singular loop transformation framework based on non-singular matrices[END_REF][START_REF] Ramanujam | Beyond unimodular transformations[END_REF][START_REF] Xue | Automating non-unimodular loop transformations for massive parallelism[END_REF]. The recent developments [START_REF] Kelly | Code generation for multiple mappings[END_REF][START_REF] Griebl | Code generation in the polytope model[END_REF][START_REF] Bastoul | Code generation in the polyhedral model is easier than you think[END_REF][START_REF] Vasilache | Polyhedral code generation in the real world[END_REF] in the polyhedral code generation have obviated those restrictions.

Canonical form of the scheduling matrix

Feautrier [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part II: multidimensional time[END_REF] and later Kelly and Pugh [START_REF] Kelly | A unifying framework for iteration reordering transformations[END_REF] have proposed a timestamp encoding for characterizing an execution order of statement instances within a non-perfectly nested loops.

Girbal et al. [START_REF] Cohen | A polyhedral approach to ease the composition of program transformations[END_REF][START_REF] Bastoul | Putting polyhedral transformations to work[END_REF][START_REF] Cohen | Facilitating the search for compositions of program transformations[END_REF] have defined the canonical form of the scheduling matrix Θ S . This encoding is generalized to handle arbitrary compositions of affine transformations. The canonical form encodes a mapping from iteration vector i to timestamp vector t, as well as a static statement order within the common loop. The scheduling matrix Θ S in the canonical form has the following format:

Θ S =             0 • • • 0 0 • • • 0 β S 0 A S 1,1 • • • A S 1,n Γ S 1,1 • • • Γ S 1,n g Γ S 1,n g +1 0 • • • 0 0 • • • 0 β S 1 A S 2,1 • • • A S 2,n Γ S 2,1 • • • Γ S 2,n g Γ S 2,n g +1 . . . . . . . . . . . . . . . . . . . . . A S n,1 • • • A S n,n Γ S n,1 • • • Γ S n,n g Γ S n,n g +1 0 • • • 0 0 • • • 0 β S n             (3.1)
The scheduling matrix Θ S is composed of three components: -Component A is an invertible (full rank) matrix capturing the relative ordering of iteration vectors. Changing coefficients of this component corresponds to loop interchange, skewing and other unimodular transformations.

-Column β reschedules statements statically, at all nesting levels. It expresses code motion, loop fusion and fission. -Component Γ captures loop shifting (pipelining) effects. We call the A and Γ matrix components a dynamic schedule components -since those components are used to reschedule dynamic statement instances, while β vector component is called a static schedule component -since it corresponds to a static statement order inside a program.

The matrix Θ S ∈ Z (2n+1)×(n+n g +1) has 2n + 1 rows, n being the loop depth of the statement S. Odd rows correspond to static execution order of statements that are enclosed in a common loop nest. Each such a row contains a single integer constant. All the odd rows could be summarized by β vector having n + 1 components. We call this vector a static scheduling vector, since it represents the constant part of the timestamp.

Rows at even positions represent an affine mapping from iteration vectors i and global parameters g to their respective timestamp component. Affine mappings enable an instance-wise multidimensional affine scheduling [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part I: one dimensional time[END_REF][START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part II: multidimensional time[END_REF][START_REF] Lim | Communication-free parallelization via affine transformations[END_REF][START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF][START_REF] Pouchet | Loop transformations: convexity, pruning and optimization[END_REF]. Changing the coefficients of A and Γ matrix components enables the representation of arbitrary affine schedules which encompass classical loop transformations like loop interchange, skewing, shifting or unimodular transformations [5,[START_REF] Allen | Optimizing Compilers for Modern Architectures[END_REF][START_REF] Wolfe | High performance compilers for parallel computing[END_REF][START_REF] Banerjee | Unimodular transformations of double loops[END_REF][START_REF] Darte | Loop shifting for loop compaction[END_REF] and much more.

An example

The scheduling matrix representing the original program schedule has the A matrix initialized to the identity matrix, the Γ matrix is all zeros and the β vector encodes relative statement order inside a program source-code (or other program intermediate representation).

As an example, let us take a non-perfectly nested loop shown in Figure 2.1. The scheduling matrix that corresponds to an original execution order of statement S 1 is given as follows:

Θ S 1 =       0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0      
this canonical scheduling matrix is decomposed into following components:

A S 1 = 1 0 0 1 Γ S 1 = 0 0 0 0 β S 1 =   0 0 0  
The matrix Θ S 1 represents the scheduling function:

θ S 1 ((v, h) T ) = (0, v, 0, h, 0) T .
In the same way, the scheduling matrix Θ S 2 is decomposed into the following components (the statement S 2 is nested within 4 nested loops):

A S 2 =     1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1     Γ S 2 =     0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0     β S 2 =       0 0 1 0 0      
that represent the scheduling function: θ S 2 (v, h, i, j) T = (0, v, 0, h, 1, i, 0, j, 0) T . The scheduling function for statement S 3 is:

θ S 3 ((v, h) T ) = (0, v, 0, h, 2) T .
Note that the timestamp components at odd positions are integer constants, whereas the components at even positions correspond to affine mapping of loop induction variables to their respective timestamp components. The fifth components of all time vectors are integer constants 0, 1, 2 -they correspond to textual order of statements S 1 , S 2 and S 3 respectively. Multidimensional time vectors are compared lexicographically. If we take statements S 1 and S 2 with their respective scheduling functions, we can say that:

θ S 1 ((v, h) T ) ≪ θ S 2 ((v ′ , h ′ , i ′ , j ′ ) T ) i.f.f. v < v ′ ∨ (v = v ′ ∧ h < h ′ ) ∨ (v = v ′ ∧ h = h ′ ∧ 0 < 1)
In other words, an instance of statement S 1 occurs before an instance of statement S 2 when one of the following conditions is satisfied:

v < v ′ , (v = v ′ ∧ h < h ′ ) or (v = v ′ ∧ h = h ′ ∧ 0 < 1)
. The first two conditions are obvious: if a statement instance is scheduled at an earlier iteration of the outermost loop then it happens before. Likewise, if a statement instance is scheduled at the same outermost loop iteration, but at an earlier innermost loop iteration, then it happens before. The last condition states: if two statement instances are scheduled at the same iteration of their common outer loops, a statement instance (S 1 , (v, h) T ) is going to be executed earlier than a statement instance (S 1 , (v ′ , h ′ , i ′ , j ′ ) T ), because it is statically scheduled to happen earlier at the fifth timestamp component (0 < 1).

Capturing semantical constraints

We have shown a powerful way to express the program transformations -by changing the schedule of the executed statement instances. Obviously, if we want to perform a legal program transformation, the semantics of the transformed program must match the semantics of the original program.

There are numerous ways to define the program semantics. But for the purpose of the polyhedral transformations, it is enough to define a semantical equivalence of the transformed program with respect to an original program. To determine this equivalence it is not necessary to define the exact semantics of each computational step.

The classical polyhedral model does not model the actual computations that are performed inside the polyhedral statements S ∈ S. Polyhedral statement is treated as a black box that reads from and/or stores to a memory -thus changing the visible state. A given sequential program can be seen as a sequence of reads from and writes to a memory, irrelevant of the computations that are performed within the statements [START_REF] Allen | Optimizing Compilers for Modern Architectures[END_REF].

The transformation is legal if the transformed program computes the same values as the original program. This definition allows the reordering of the statement instances, as long as the transformed program would provide the equivalent results as the original one.

We use a well defined concept of data dependence as a semantical constraint that has to be preserved in order for a transformation to be defined as a legal. [START_REF] Bernstein | Analysis of programs for parallel processing[END_REF][START_REF] Allen | Optimizing Compilers for Modern Architectures[END_REF].

A good survey of data dependence abstractions was given in [START_REF] Yang | Minimal data dependence abstractions for loop transformations[END_REF][START_REF] Darte | On the optimality of allen and kennedy's algorithm for parallel extraction in nested loops[END_REF][START_REF] Vivien | On the optimality of feautrier's scheduling algorithm[END_REF][START_REF] Banerjee | Dependence Analysis for Supercomputing[END_REF]. Naturally, the more precise the data dependence abstraction is, the more legal transformations are exploitable. The most precise data dependence abstraction is that of Feautrier [START_REF] Feautrier | Dataflow analysis of scalar and array references[END_REF] expressed as a dependence polyhedra. The similar approach was proposed by Pugh [START_REF] Pugh | The omega test: a fast and practical integer programming algorithm for dependence analysis[END_REF][START_REF] Pugh | Eliminating false data dependences using the omega test[END_REF][START_REF] Pugh | An exact method for analysis of value-based array data dependences[END_REF].

A practical and efficient implementation of instance-wise dependence analysis was provided by Vasilache [START_REF] Vasilache | Violated dependence analysis[END_REF]. Barthou [18] has extended the scope of the application of the polyhedral model to non-affine programs by providing the fuzzy array dataflow analysis -an approach that handles the non-affine program parts in a conservative way.

We will proceed by presenting the data dependence representation and analysis. Later we will introduce a refinement of dependence analysis that captures only the dataflow dependences.

Data dependences

Two statement instances (S i , i S i ) and (S j , i S j ) are in a data dependence relation if they access the same memory location and at least one of those accesses is a write memory access. We designate a statement instance that happens earlier as the source instance, and a statement instance that happens later as the target instance. If a source instance is a write and a target instance is a read then the data dependence is a true dependence. If a source instance is a read and a target instance is a write then the dependence is a antidependence. Lastly, if a source instance is a write and a target instance is a write then the data dependence is an output dependence. All the possible types of dependences are summarized in Table 3.2.

It is not practical to compute data dependences between all statement instances in a program. What is more, if a loop iteration count depends on a parameter that is unknown at the compilation time, it is not even possible to enumerate all the data dependences in a program.

A compact way to summarize all the data dependences between all statement instances in a program is to use a structure called the Reduced Data Dependence Graph (RDDG) [START_REF] Darte | Scheduling and Automatic Parallelization[END_REF].

A RDDG is a directed multi-graph, G = (V, E), whose vertex set V is a set of all program statements, i.e., V = S. There is an edge e ∈ E from a vertex S i to a vertex S j if there is a data dependence relation between an access in some instance of S i and an access in some instance of S j . Each edge e ∈ E is labelled by a dependence polyhedron P e . The dependence polyhedron P e that labels an edge from S i to S j represents the set of iteration vector pairs denoting statement instances that are in data dependence relation:

P e = {(i S i , i S j ) | (S i , i S i
) and (S j , i S j ) are data dependent }

In order to compute such an instancewise dependence polyhedron, one has to construct an intersection of the following affine constraints [START_REF] Vasilache | Violated dependence analysis[END_REF]:

Conflict condition both statement instances access the same memory location. This is equivalent to stating that their subscript functions have all their subscripts equal:

f S i (i S i ) = f S j (i S j )
Causality condition the instance (S i , i S i ) happens before the instance (S j , i S j ) in the original program execution order:

θ S i (i S i ) ≪ θ S j (i S j )
Execution condition both statement instances are actually executed, i.e., both instances belong to statement iteration domains:

i S i ∈ D S i ∧ i S j ∈ D S j
Since all the conditions could be expressed as a system of affine equalities or inequalities, the resulting affine constraint system represents an affine relation that can be easily manipulated by an integer linear programming and linear algebra tools. The causality condition is a lexicographical comparison of vectors and could be represented as a disjunction of N S i S j constraint components, where N S i S j is a common loop depth of statements S i and S j .

An example

As an example, let us take a matrix-vector multiplication kernel shown in Figure 3.3. There is a true data dependence between an instance of the statement S 1 and an instance of the statement S 2 , denoted As the source of the dependence we consider the statement S 1 containing a write to an array b[i] whose corresponding subscript function is f S 1 (i) T = (i). As the target of the dependence we consider the statement S 2 containing a read from an array b[i] whose corresponding subscript function is f S 2 (i, j) T = (i). We summarize the necessary affine constraints in Table 3.1. Taking an intersection of those affine constraints, we end up with the following dependence polyhedron 1 :

f o r ( i = 0 ; i < N; i ++) { S1 : b [ i ] = 0 ; f o r ( j = 0 ; j < N; j ++) S2 : b [ i ] = b [ i ] + A[ i ] [ j ] * x [ j ] ; }
(i) T ∈ D S 1 0 ≤ i < N execution condition (i ′ , j ′ ) T ∈ D S 2 0 ≤ i ′ < N ∧ 0 ≤ j ′ < N conflict condition f S 1 ((i) T ) = f S 2 ((i ′ , j ′ ) T ) i = i ′ causality condition θ S 1 ((i) T ) ≪ θ S 2 ((i ′ , j ′ ) T ) i ≤ i ′
P e :         1 -1 0 0 0 -1 1 0 0 0 0 1 0 0 0 0 -1 0 1 -1 0 0 1 0 0 0 0 -1 1 -1               i i ′ j ′ N 1       ≥ 0 i ≥ i ′ i ≤ i ′ i ′ ≥ 0 i ′ ≤ N -1 j ′ ≥ 0 j ′ ≤ N -1
A dependence polyhedron P e is a relation that describes all pairs of iteration vectors (i) and (i ′ , j ′ ) of statement instances that are in a dependence relation. The pairs of iteration vectors are concatenated in a single vector (i, i ′ , j ′ ), so the relation is a subset of Cartesian product space:

P e = {(i, i ′ , j ′ ) ∈ Z dim(S 1 )+dim(S 2 ) |(S 2 , (i ′ , j ′ )) is dependent on (S 1 , (i))}
thus a dependence polyhedron could be concisely represented as:

P e 1 = {(i, i ′ , j ′ )|i = i ′ ∧ 0 ≤ i ′ ≤ N -1 ∧ 0 ≤ j ′ ≤ N -1}
In our example, the pair is decomposed in the first component (i) T that corresponds to the source iteration, while the second component (i ′ , j ′ ) T corresponds to the target iteration. In other words, if the value was stored in b[i] in the i-th iteration of the statement S 1 , it is read by all the iterations of the statement S 2 for which i ′ = i. There is no constraint on the j ′ , since for i = i ′ , the statement instance S 2 is always executed after an instance of the statement S 1 .

A true data dependence between S 1 and S 2 is not the only data dependence that exists in a program given in Figure 3.3. For example, there exists a write-after-write (output) dependence between statements S 1 and S 2 , described by the same polyhedron as a true dependence. The reason for the output data dependence is the fact that the write to an array b[i] at S 1 is subsequently overwritten by a write to an array

Memory based dependences

It is a well known fact [START_REF] Lefebvre | Automatic storage management for parallel programs[END_REF] that anti and output data dependences can always be removed from the dependence graph by a technique called an array expansion [START_REF] Maydan | Array-data flow analysis and its use in array privatization[END_REF][START_REF] Calland | On the removal of anti-and output-dependences[END_REF]. The anti and output dependences are called memory based dependences [START_REF] Pugh | An exact method for analysis of value-based array data dependences[END_REF] since they are induced by the fact that the same memory cell is reused for storing values. One can remove those dependences by attributing a separate memory cell to each statement instance, so that the same memory cell will not be reused -an array expansion.

On the contrary, the true dependences have to be preserved, since they capture the relation between the producer and a consumer of the data. The subset of the true dependences, named dataflow dependences has to be always preserved -those dependences cannot be removed from the dependence graph.

In the next subsection we will discuss the technique for computing the set of dataflow dependencesan array dataflow analysis. The whole Chapter 5 will be devoted to a technique for efficiently handling memory based dependences.

Array data-flow analysis

Standard data dependence relation gives an information on whether there exists an aliasing of memory locations. If a statement instance (S i , i S i ) writes to a memory location A[ f S i (i S i )], and a statement instance (S j , i S j ) reads from the same memory location, they are considered to be in a true data dependence relation, even if there was an intervening write coming from some other statement instance (S k , i S k ).

This notion captures the fact of the potential memory conflict in the case the statement instances are rescheduled, but it does not provide an information on the flow of values: a read in the statement instance (S j , i S j ) will only see a value that was stored in statement instance that was scheduled as the latest before statement instance (S j , i S j ).

In order to capture the flow of values for each array element, an instancewise array data-flow analysis based on the polyhedral model was provided by Feautrier [START_REF] Feautrier | Dataflow analysis of scalar and array references[END_REF] and Pugh [START_REF] Pugh | An exact method for analysis of value-based array data dependences[END_REF]. It has also been studied by Maydan and Lam [START_REF] Maydan | Array-data flow analysis and its use in array privatization[END_REF].

The purpose of array data-flow analysis is to give an exact characterization of the source statement instance that defines value stored in a given array element that is seen by the given read statement instance. array data-flow analysis is also known under the name of value-based array data dependence analysis [START_REF] Pugh | An exact method for analysis of value-based array data dependences[END_REF], since it captures the data dependences induced by the flow of values, and not merely by the aliasing of memory locations.

We will define a true data-flow dependence between two statement instances as a true data dependence satisfying an additional constraint stating that an element written in a source statement instance must not be overwritten by some other statement that is executed between source and target statement instances. The additional constraint is called a liveness condition. If we put together all the affine constraints together, we will get the following set of affine constraints:

Conflict condition both statement instances access the same memory location. This is equivalent to stating that they access the same memory array A and their subscript functions have all their subscripts equal:

f S i (i S i ) = f S j (i S j )
Causality condition the instance (S i , i S i ) happens before the instance (S j , i S j ) in the original program execution order:

θ S i (i S i ) ≪ θ S j (i S j )
Execution condition both statement instances are actually executed, i.e., both instances belong to statement iteration domains:

i S i ∈ D S i ∧ i S j ∈ D S j
Liveness condition there is no overwriting statement instance that happens between source and target statement:

∀S k ∈ S kill : ¬∃i S k : θ S i (i S i ) ≪ θ S k (i S k ) ≪ θ S j (i S j ) ∧ f S i (i S i ) = f S k (i S k )
S kill is the set of statements containing array references that might write to the array A. Note that both S i and S j are included in S kill .

While this definition is a denotational description of how to compute value-based (data-flow) dependences, it does not express those constraints directly as a set of affine constraints. Indeed, an existential quantifier and a negation operator might induce a solution that contains a set of non-convex polyhedra. Pugh [START_REF] Pugh | An exact method for analysis of value-based array data dependences[END_REF] shows different heuristics that control the complexity of the solution.

An ILP based formulation

Featurier [START_REF] Feautrier | Dataflow analysis of scalar and array references[END_REF] has shown that the previously mentioned liveness condition might be equivalently expressed as the following problem:

Given a read statement instance (S j , i S j ) containing a read data reference A[ f S j (i S j )], compute an unique statement S i and its iteration vector i S i that form a source statement instance

(S i , i S i ) being the source of the value read in A[ f S j (i S j )].
This provides a notion of source function H : A[ f S j (i S j )], (S j , i S j ) → (S i , i S i ): an argument is a data reference A[ f S j (i S j )] and its associated statement instance (S j , i S j ) while the result of the function is source statement instance (S i , i S i ) that produced a value in A[ f S j (i S j )].

Among the set of possible source statement instances only the one that is executed as the latest is the true source of the given value, since all values written by other sources are overwritten by the latest executed instance.

In order to compute the unique value of the source function for each possible data reference, Feautrier [START_REF] Feautrier | Dataflow analysis of scalar and array references[END_REF] proposes to use a PIP [START_REF] Feautrier | Parametric integer programming[END_REF] (Parametric Integer Programming) algorithm to compute the lexicographically largest element of the parametric set Q.

The parametric set Q A[ f S j (i S j )], (S j , i S j ) = {(S k , i S k )}, contains the set of all executed statement instances (S k , i S k ) that write to the given array element A[ f S j (i S j )] and that happen before a statement

(S 1 , (i) T ) b[i] = 0; (S 2 , (i, 0) T ) b[i] = b[i] + A[i][0] * x[0]; (S 2 , (i, 1) T ) b[i] = b[i] + A[i][1] * x[1]; . . . . . . (S 2 , (i, N -1) T ) b[i] = b[i] + A[i][N-1] * x[N-1];
Figure 3.5 -A slice of the execution trace of matvect kernel instance (S j , i S j ) (instances that satisfy conflict, execution and causality condition). An unique solution is obtained by solving the following parametric integer programming problem:

max ≪ θ(S k , i S k ), (S k , i S k ) ∈ Q A[ f S j (i S j )], (S j , i S j )
In our approach, we will use an implementation of the array dataflow analysis as proposed by Feautrier and implemented in an ISL [START_REF] Verdoolaege | isl: An integer set library for the polyhedral model[END_REF] library. An extension of array data-flow analysis to non-affine programs proposed by Barthou [18] follows the same principle, but it introduces an additional predicate variables that model the non-affine access functions and non-affine control flow in the program.

An Example

As an illustration, let us consider a kernel in Figure 3.3. Figure 3.5 illustrates a slice of the sequential execution trace for some fixed outer loop induction variable i. For a given outer loop induction variable i one instance of statement S 1 is executed: (S 1 , (i) T ). Subsequently N -1 instances of S 2 are executed: first an instance (S 2 , (i, 0) T ) is executed, next an instance (S 2 , (i, 1) T ) and so on.

By definition, the alias based dependence analysis considers all the pairs of dependences between statements S 1 and S 2 : (S 1 , (i) T ) → (S 2 , (i, 0) T ), (S 1 , (i) T ) → (S 2 , (i, 1) T ), . . . , (S 1 , (i) T ) → (S 2 , (i, N -1) T ). Those dependences are induced by a write b[i] = 0 in statement S 1 and subsequent reads ... = b[i] in statement instances of statement S 2 .

The array data-flow analysis will only consider the true value-based dependence (S 1 , (i) T ) → (S 2 , (i, 0) T ). There is no value-based dependence (S 1 , (i) T ) → (S 2 , (i, 1) T ), since the value written in statement instance (S 1 , (i) T ) is overwritten in statement instance (S 2 , (i, 0) T ) and thus not visible in instances (S 2 , (i, 1) T ), . . . , (S 2 , (i, N -1) T ).

Transformations and legality

A program transformation in the polyhedral model is expressed as an reordering transformation on statement instances. One expresses a transformation as a set of scheduling matrices. For each statement S i ∈ S one provides a scheduling matrix θ S i that completely describes a multidimensional schedule of statement instances. Definition 3.3.1 (Dependence satisfaction). A RDDG edge e ∈ E from S i to S j , labelled by dependence polyhedron P e , is satisfied by schedules θ S i and θ S i iff:

∀(i S i , i S j ) T ∈ P e : θ S i (i S i ) ≪ θ S j (i S j ) (3.2)
A dependence satisfaction condition for a dependence edge e ∈ E simply states that for each pair of dependent statement instances, the source instance (S i , i S i ) is scheduled before the target instance (S j , i S j ). This is expressed as a lexicographical comparison of their timestamps which are transformed through their respective scheduling functions θ S i (i S i ) and θ S j (i S j ).

A program transformation is expressed as the set of scheduling matrices:

T = {θ S 1 , . . . , θ S p }
where p is the number of statements inside the program statement set S. We say that a transformation T is legal for a given RDDG G=(V, E) if all the edges e ∈ E are satisfied by their respective statement schedules.

A priori versus a posteriori legality

Depending on the usage scenario, there are two major approaches to assuring the legality of transformations in the polyhedral model: a posteriori checking the legality of the given transformation. A transformation might be given by an user or obtained semi-automatically [START_REF] Chen | CHiLL: A framework for composing high-level loop transformations[END_REF][START_REF] Girbal | Semiautomatic composition of loop transformations for deep parallelism and memory hierarchies[END_REF]. A violation analysis [START_REF] Vasilache | Violated dependence analysis[END_REF] is performed to check the legality of the proposed transformation.

a priori computing the automatic transformations that are guaranteed to be legal. All the affine scheduling approaches fall into this category [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part I: one dimensional time[END_REF][START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part II: multidimensional time[END_REF][START_REF] Griebl | Automatic parallelization of loop programs for distributed memory architectures[END_REF][START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF][START_REF] Lim | Maximizing parallelism and minimizing synchronization with affine transforms[END_REF][START_REF] Lim | An affine partitioning algorithm to maximize parallelism and minimize communication[END_REF].

A computational procedure for a posteriori legality check is based on a legality condition given in Equation 3.2. But in order to express this condition in a way that could be handled by the polyhedral techniques, it has to be rephrased into an equivalent condition:

∃(i S i , i S j ) T ∈ P e : θ S j (i S j ) ≪ θ S i (i S i )
This set represents those pairs of statement instances that are in dependence relation, but whose order is inversed (please note that lexicographical comparison is an inversion of the one shown in Equation 3.2). If this set is empty, there are no illegal inversed statement instances, and the given dependence relation is preserver. The same check is performed for each dependence edge e ∈ E.

An emptiness check is done by performing a Fourier-Motzkin elimination [START_REF] Schrijver | Theory of linear and integer programming[END_REF] restricted to integer solutions only [START_REF] Pugh | The omega test: a fast and practical integer programming algorithm for dependence analysis[END_REF]. An efficient implementation of the violation check was proposed by Vasilache [START_REF] Vasilache | Violated dependence analysis[END_REF].

In the case there is some pair of dependent statement instances that is inverted it is not a dead-end: one might resort to correct the schedule by applying a corrective shifting. This technique was proposed by Vasilache [START_REF] Vasilache | Automatic correction of loop transformations[END_REF].

An another approach is to automatically build the legal schedules a priori. This is the approach taken by all the affine scheduling algorithms [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part I: one dimensional time[END_REF][START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part II: multidimensional time[END_REF][START_REF] Griebl | Automatic parallelization of loop programs for distributed memory architectures[END_REF][START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF][START_REF] Lim | Maximizing parallelism and minimizing synchronization with affine transforms[END_REF][START_REF] Lim | An affine partitioning algorithm to maximize parallelism and minimize communication[END_REF].

One can build the space of legal transformations by putting together the dependence satisfaction constraints for all the dependence edges e ∈ E in RDDG. As a solution, one gets the space of legal scheduling matrices, expressed as the set of possible matrix coefficients that give the legal transformation.

Building and solving such a constraint system is based on affine form of Farkas lemma [START_REF] Schrijver | Theory of linear and integer programming[END_REF] (to be defined in Chapter 7) and Fourier-Motzkin elimination [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part I: one dimensional time[END_REF]. The related problem that arises is that of selecting the best schedule among the set of possible legal schedules. That problem is the core of the Part III of this dissertation.

An Example

As an example, consider a matvect kernel in Figure 3.3. The original scheduling matrices Θ S 1 and Θ S 2 are shown in Figure 3.6. By definition, the original scheduling is always legal -indeed a RDDG is computed from the original scheduling matrices. Let us consider a true data dependence edge from S 1 to S 2 , whose dependence instances are summarized by the following dependence polyhedron:

Θ S 1 =   0 0 0 1 0 0 0 0 0     i N 1   =   0 i 0   Θ S 2 =       0 0 0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0           i j N 1     =       0 i 1 j 0      
P e = {(i, i ′ , j ′ ) T | 0 ≤ i < N ∧ 0 ≤ i ′ < N ∧ 0 ≤ j ′ < N ∧ i = i ′ }.
An illustration of the original execution order and data dependence instances for a dependence edge S 1 → S 2 is shown in Figure 3.7.

As an instance of a possible reordering transformation, let us change the scheduling matrix of statement S 2 by interchanging the original scheduling matrix columns. A program transformation is represented by a set of new scheduling matrices: T = {Θ ′S 1 , Θ ′S 2 }. The scheduling matrix of statement S 1 is left intact. The scheduling is illustrated in Figure 3.8.

The transformed schedule does not satisfy the dependence edge e from S 1 to S 2 summarized by the polyhedron P e . This is illustrated in Figure 3.11 showing a new (rescheduled) execution order, together with an instance of dependence pair that is not satisfied. Indeed, a true data dependence from S 1 to S 2 requires that an instance (S 1 , (2) T ) is executed before an instance (S 2 , (2, 0) T ), but in a transformed schedule an instance (S 2 , (2, 0) T ) is executed first -a data dependence is inversed.

On the other hand, Figure 3.9 shows a scheduling transformation that is legal. As in the previous example, the scheduling matrix of the statement S 1 is left intact. The scheduling matrix of the statement S 2 is formed by putting an integer constant in the first row and leaving the rest of the rows intact. In this way, we have changed the β (see Section 3.1) component of the scheduling matrix. .10 illustrates the new execution order -first, all statement instances of the statement S 1 are executed, and then all instances of S 2 are executed. This execution order preserves all dependence pairs for a true data dependence from S 1 to S 2 .

Θ ′S 1 =   0 0 0 1 0 0 0 0 0     i N 1   =   0 i 0   Θ ′S 2 =       0 0 0 0 0 1 0 0 0 0 0 1 1 0 0 0 0 0 0 0           i j N 1     =       0 j 1 i 0       Figure 3.8 -Illegal transformation matrices Θ ′S 1 =   0 0 0 1 0 0 0 0 0     i N 1   =   0 i 0   Θ ′S 2 =       0 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0           i j N 1     =       1 i 0 j 0      
Please note that in order to prove legality of a transformation, one has to prove that all dependence pairs are satisfied. On the other hand, in order to prove that a transformation is illegal it is enough to demonstrate that there exists at least one dependence pair that is inversed.

Related work

A loop transformation theory based on data dependences was developed by researchers working on optimizing compilers for high performance computing [START_REF] Banerjee | Data dependence in ordinary programs[END_REF][START_REF] Kuck | Dependence graphs and compiler optimizations[END_REF][START_REF] Allen | Automatic translation of fortran programs to vector form[END_REF]. The data dependence graph was used as a central concept for ensuring the legality of transformations. The loop transformations were performed on a syntax representation.

The first attempts to model the imperative program loop transformations through linear algebraic representation were those of unimodular transformations [START_REF] Lamport | The parallel execution of do loops[END_REF][START_REF] Banerjee | Dependence Analysis for Supercomputing[END_REF][START_REF] Wolf | A loop transformation theory and an algorithm to maximize parallelism[END_REF]. The loop transformations ware represented by unimodular matrices. A sequence of several loop transformations is represented in a single unimodular matrix and transformed source code is generated automatically. The drawback of unimodular transformations is the fact that they are amenable to perfectly nested loops only. Some works were trying to circumvent the limitation of unimodular transformations so that they could be applied to non-perfectly nested loops [START_REF] Li | A singular loop transformation framework based on non-singular matrices[END_REF][START_REF] Xue | Unimodular transformations of non-perfectly nested loops[END_REF]. But they have limitations in code generation phase -requiring the transformation matrices to be invertible.

A unifying theory that could handle non-perfectly nested loops was developed [START_REF] Pugh | Uniform techniques for loop optimization[END_REF][START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part I: one dimensional time[END_REF][START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part II: multidimensional time[END_REF][START_REF] Lim | Communication-free parallelization via affine transformations[END_REF]. It relies on multidimensional schedules of statement instances. They also show how to automatically search for a schedule that optimizes program parallelization or data locality.

Feautrier [START_REF] Feautrier | Array expansion[END_REF] has demonstrated the feasibility of handling static control program analysis by using integer linear programming, and he proposed a Parametrized Integer Programming (PIP) algorithm. Later, Feautrier has shown how to handle the problem of instancewise array dataflow analysis [START_REF] Feautrier | Dataflow analysis of scalar and array references[END_REF] and minimal latency scheduling problem [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part I: one dimensional time[END_REF] -two seminal works in the history of theoretical polyhedral compilation.

Another technique developed by Pugh [START_REF] Pugh | The omega test: a fast and practical integer programming algorithm for dependence analysis[END_REF] is based on an extension of Fourier-Motzkin variable elimination which led to the development of library for manipulating integer sets and relations represented as systems of affine constraints called Omega [START_REF] Kelly | The omega library interface guide[END_REF]. Pugh has also solved a problem of array dataflow analysis under the name of value based array dependence analysis [START_REF] Pugh | An exact method for analysis of value-based array data dependences[END_REF].

More recent advances in the scalable and general code generation techniques [START_REF] Bastoul | Efficient code generation for automatic parallelization and optimization[END_REF][START_REF] Bastoul | Code generation in the polyhedral model is easier than you think[END_REF][START_REF] Vasilache | Polyhedral code generation in the real world[END_REF] enabled the polyhedral model to be integrated into the real production compiler.

Summary

In this chapter we have summarized the state of the art in program transformations expressed in the polyhedral model. The concepts of data dependences and array dataflow analysis were introduced. The notion of transformation legality was briefly discussed. Some implementational issues of array dataflow analysis and transformation legality check were presented.

We will extend the concepts that were presented in this chapter with our new contributions. In Chapter 5 we will show that we can relax the conditions for determining the legality of the transformation. In Chapter 6 we will show that we can construct a precise cost models based on the polyhedral representation, provided that we do not treat the polyhedral statement as a black box. Chapter 7 provides our contribution to the schedule construction problem.

Part II

The framework Chapter 4

Graphite framework

Despite several decades of research into the polyhedral model, there is still no general-purpose production compiler using the polyhedral model internally. The situation is changing with the demonstration of the scalability of polyhedral algorithms and with the widespread dissemination of multicore processors and hardware accelerators.

This chapter describes GRAPHITE framework that incorporates the polyhedral analyses and transformations into GCC1 compiler. GCC compiler is one of the most widely used open-source compilers for a variety of imperative programming languages (C, C++, ADA, Fortran) and platforms (x86, ARM, PowerPC, MIPS).

The original motivation for GRAPHITE project is the study of loop optimizations on a low-level three-address code [3] intermediate representation. Polyhedral representation is extracted from the GCC three-address code intermediate representation that is in the static single assignment (SSA) [START_REF] Cytron | An efficient method of computing static single assignment form[END_REF] form. This is a major difference with traditional source-to-source polyhedral compilers that operate on high-level abstract syntax level. Operating directly on the three-address code brings in new challenges but also new opportunities: we can leverage existing scalar analyses in the compiler and interact directly with a wealth of optimizations for extracting coarse and fine-grained parallelism and improving memory locality.

We will first discuss the related work on different approaches to polyhedral compilation. A general overview of the compilation flow in GRAPHITE will be presented. Later, a discussion of the relevant and interesting design issues follows, together with the conclusion.

Related Work

There have been many efforts in designing an advanced loop-nest transformation infrastructure. Most loop restructuring compilers introduced syntax-based models and intermediate representations. ParaScope [START_REF] Cooper | The ParaScope parallel programming environment[END_REF] and Polaris [START_REF] Blume | Parallel programming with Polaris[END_REF] are dependence-based, source-to-source parallelizers for Fortran. KAP [START_REF] Kap C/ | OpenMP for Tru64 UNIX and KAP DEC Fortran for Digital UNIX[END_REF] is closely related to these academic tools.

SUIF [START_REF] Hall | Maximizing multiprocessor performance with the SUIF compiler[END_REF] is a platform for implementing advanced compiler prototypes. PIPS [START_REF] Irigoin | Semantical interprocedural parallelization: An overview of the pips project[END_REF] is one of the most complete loop restructuring compilers, implementing polyhedral analyses and transformations (including affine scheduling) and interprocedural analyses (array regions, alias). Both of them use a syntax tree extended with polyhedral annotations, but not a unified polyhedral representation.

The MARS compiler [START_REF] O'boyle | MARS: a distributed memory approach to shared memory compilation[END_REF] unifies classical dependence-based loop transformations with data storage optimizations. However, the MARS intermediate representation only captures part of the loop in-formation (domains and access functions): it lacks the characterization of iteration orderings through multidimensional affine schedules.

The first thorough application of the polyhedral representation was the Petit tool [START_REF] Kelly | Optimization within a unified transformation framework[END_REF], based on the Omega library [START_REF] Kelly | Code generation for multiple mappings[END_REF]. It provides space-time mappings for iteration reordering, and it shares our emphasis on per-statement transformations, but it is intended as a research tool for small kernels only. We also use a code generation technique that is significantly more robust than the code generation in Omega [START_REF] Bastoul | Code generation in the polyhedral model is easier than you think[END_REF].

Semi-automatic polyhedral frameworks have been designed as building blocks for compiler construction or (auto-tuned) library generation systems [START_REF] Kelly | Optimization within a unified transformation framework[END_REF][START_REF] Cohen | Facilitating the search for compositions of program transformations[END_REF][START_REF] Vasilache | Automatic correction of loop transformations[END_REF][START_REF] Chen | CHiLL: A framework for composing high-level loop transformations[END_REF][START_REF] Tiwari | A scalable autotuning framework for computer optimization[END_REF]. They do not define automatic methods or integrate a model-based heuristic to construct profitable optimization strategies.

The GRAPHITE project was first announced by Pop et al. in 2006, [125]. The design of GRAPHITE is largely borrowed from the WRaP-IT polyhedral interface to Open64 and its URUK loop nest optimizer [START_REF] Girbal | Semiautomatic composition of loop transformations for deep parallelism and memory hierarchies[END_REF]. The CHiLL project from Chen et al. revisited the URUK approach, focusing on source-to-source transformation scripting [START_REF] Chen | CHiLL: A framework for composing high-level loop transformations[END_REF][START_REF] Tiwari | A scalable autotuning framework for computer optimization[END_REF].

Unlike URUK and CHiLL, GRAPHITE aims at complete automation, resorting to iterative search and cost modeling of the profitability of program transformations. In addition, unexpected design and algorithmic issues have been discovered, partly due to the design of GCC itself, but mostly due to the integration of the polyhedral representation in a three-address code in SSA form.

To the best of our knowledge, two proprietary polyhedral compilers based on a low-level internal representation are currently in development: the R-Stream compiler from Reservoir Labs [START_REF] Meister | Productivity via automatic code generation for pgas platforms with the r-stream compiler[END_REF], and IBM's polyhedral extension of its XL compiler suite [START_REF] Renganarayana | Compact multi-dimensional kernel extraction for register tiling[END_REF]. Little work has been published on compilation issues of three-address code based polyhedral model abstraction.

An overview of the compilation flow

Traditionally the polyhedral model is used in the source-to-source [START_REF]PLUTO: A polyhedral automatic parallelizer and locality optimizer for multicores[END_REF]124] translators and may be viewed in terms of the three step process: (1) extraction of the polyhedral representation of the static control program part, (2) transformation on the polyhedral abstraction and (3) generation of the equivalent source code fragment from the abstraction.

On the other hand, GRAPHITE is the polyhedral transformation engine that is only a small part of the complete compiler toolchain -it is implemented as an optimization pass of GCC compiler. As such, GRAPHITE operates on a three-address code internal representation. A place of GRAPHITE inside GCC compiler is illustrated in Figure 4.1.

GCC compiler is decoupled into three principal components: front-end, middle-end and back-end, as shown in Figure 4.1. Front-end is language dependent. It is responsible for parsing the source language 2 into GENERIC [START_REF] Fsf | GNU Compiler Collection (GCC) Internals Manual[END_REF] high-level abstract syntax tree. GENERIC representation is lowered into GIMPLE [START_REF] Fsf | GNU Compiler Collection (GCC) Internals Manual[END_REF] intermediate three-address code code used in compiler middle-end.

Middle-end is language independent optimization stage. Inside this stage multiple optimization passes are executed 3 . Intraprocedural analyses and optimizations are performed. GIMPLE three-address code is transformed into SSA form. Various optimizations based on SSA form are performed, including DCE (Dead Code Elimination), forward propagation and copy propagation. Loop based optimizations such as loop invariant motion, loop unswitching (moving conditional statements out of the loop), PRE(Partial Redundancy Elimination) are also performed.

GRAPHITE is scheduled as a loop optimization pass operating on a transformed GIMPLE code in SSA form. GRAPHITE itself is composed of several steps that will be explained in the following section. After middle-end optimization passes, the transformed GIMPLE code is in turn lowered into RTL (Register Transfer Language) [START_REF] Fsf | GNU Compiler Collection (GCC) Internals Manual[END_REF] which is used in the GCC back-end. Back-end is target machine dependent and can produce machine code for different architectures (x86, x86_64, PowerPC, MIPS, ARM).

Further in this work we will not consider neither front-end neither back-end, since we are operating on a language and architecture neutral middle-end stage. The only architecture dependent part of the GRAPHITE loop transformation search algorithm is the dependence on the target architecture cost model as shown in Figure 4.1.

GIMPLE internal representation

GRAPHITE is operating on a three-address code internal representation in SSA (Single Static Assignment) form. An internal representation called GIMPLE is composed of basic blocks containing three-address code statements. Each statement has no more than 3 operands. Necessary temporaries are introduced to store intermediate values needed to compute complex expressions. Each statement could contain at most one data access to a memory location, either read or write. All the high level control structures are expressed as conditional jumps and the lexical scopes are removed [START_REF] Fsf | GNU Compiler Collection (GCC) Internals Manual[END_REF].

GCC keeps track of control flow between basic blocks in CFG (Control Flow Graph) graph, whose vertices are the basic blocks and whose directed edges show the control-flow. Loops are strongly connected components of the CFG.

SSA [START_REF] Cytron | An efficient method of computing static single assignment form[END_REF] form ensures uniqueness of a variable definition: each variable is assigned only once. When transforming into SSA form, each original variable is split into versions if there are multiple definitions of the same variable. So, for example, if there are two distinct assignments to the variable x, it would be split into two variables x 1 and x 2 .

GIMPLE in SSA form contains φ-functions [START_REF] Cytron | An efficient method of computing static single assignment form[END_REF] of the form φ(x 1 , x 2 , . . . , x n ), where x 1 , x 2 , . . . , x n are variable versions and the number of arguments of the φ-function is the number of distinct control flow paths that reach the φ-function. φ-functions are inserted at the places where multiple distinct variable versions x 1 , x 2 , . . . , x n might reach the use, depending on the control flow.

An excerpt of GIMPLE internal representation and the corresponding source code is shown in Figure 4.7.

Inside GRAPHITE

GRAPHITEis the polyhedral compilation framework that operates on GIMPLE code in SSA form. As such, it has to first extract the polyhedral representation out of three-address code. Also, since it is a GCC pass, it has to produce GIMPLE three-address code.

This two tasks are more complicated when operating on three-address code than when operating on high-level abstract syntax trees. In this section we will discuss an approach of mapping three-address GIMPLE code into the polyhedral model, and mapping the polyhedral model back into the three-address GIMPLE code. Actual polyhedral transformations are the topic of subsequent chapters. GRAPHITE is itself composed of several sequential stages.

GRAPHITE pass

symbol table

Figure 4.2 -Inside Graphite

-Data dependence analysis. Data-dependence analysis is performed to construct RDDG. Special treatment for scalar dependences and reductions is applied. -Transformation search. A search for a legal loop transformation sequence based on a cost model is performed. This is the topic of chapters 6 and 7. -GIMPLE code generation. Three-address GIMPLE code corresponding to the transformation expressed in the polyhedral model is generated. We will discuss the technical details of GRAPHITE steps in the subsequent sections. We will not be covering the actual loop transformations, since the search for efficient loop transformations is extensively covered in Chapters 6, 5 and 7. We will start by briefly discussing an induction variable recognition, which is an essential building block for the operation of GRAPHITE.

Induction variable recognition

Induction variable analysis is essential part of any loop optimizing compiler [START_REF] Engelen | Efficient symbolic analysis for optimizing compilers[END_REF]. GCC implements induction variable analysis by constructing TREC(Tree of RECurrences) [START_REF] Pop | Induction variable analysis with delayed abstractions[END_REF] -a closed form expressions that capture the evolution of induction variable as a function of iteration counts:

χ(i) = χ(i 1 , i 2 , . . . , i n )
TREC expressions are either constants χ = c or they are defined recursively as χ = {χ a , +, χ b } k , where χ a and χ b are trees of recurrences, c is an integer constant or variable name, and subscript k is the loop dimension along which we are evaluating an induction variable [START_REF] Pop | Induction variable analysis with delayed abstractions[END_REF]. If we consider an evolution of induction variable in the loop at level k, χ a is an initial value and χ b is an increment at each loop iteration.

TREC expressions can represent linear, affine, polynomial and exponential functions as well. We will not provide details on construction and evaluation of those expressions, for the reference please see [START_REF] Pop | Induction variable analysis with delayed abstractions[END_REF].

The use of TREC expressions in GRAPHITE is threefold: (1) detection of static control part regions, (2) construction of the polyhedral information for the given CFG subgraph and (3) code generation.

SCoP outlining

The first step in GRAPHITE is the extraction of the maximal subgraphs of CFG that have static control described by the affine loop bounds, affine conditionals and affine data access patterns. Each such a subgraph is called SCoP(Static Control Part) in GRAPHITE terminology.

Since GRAPHITE is operating on a low level three-address code scattered in a CFG, high level syntactical information is lost: loop structure, loop induction variables, loop bounds, conditionals, data accesses and reductions. Extraction of maximal SCoP regions proceeds in two steps:

1. Construct a region tree that represents the nesting of reducible [START_REF] Carter | Folklore confirmed: reducible flow graphs are exponentially larger[END_REF] regions in CFG.

2. Traverse the tree, starting from the topmost root region. For each region (in the traversal order):

(a) Check if all the basic-blocks belonging to the region satisfy the affine program conditions. If yes:

i. Mark the region as a SCoP else i. Mark the region as nonSCoP. Repeat the procedure for child regions.

In order to check whether a basic block satisfies affine program conditions several checks are performed:

1. if a basic block is a loop header block: If some basic-block does not satisfy affine program conditions, the region that contains this basic block is marked as a non-SCoP. In such a case, the sub-regions (child nodes in a region tree) are checked recursively.

The most favorable case is when the outermost region is proved to be a SCoP. This is obvious, especially if the outer SCoP contains the outermost loop in a function. In a less favorable case, some child nodes in the region tree are marked as a SCoP and some as a non-SCoP. It is less favorable, since inner SCoPs contain inner-loops and also because inner SCoPs are fragmented. The least favorable case happens when a SCoP region is composed of a single basic block only -such SCoPs are called degenerate and they do not contain any loop. Degenerate SCoPs are not optimized.

As an example, a CFG with marked reducible regions is shown in Figure 4.3. Regions themselves are either nested within each other or they are ordered in a sequence. A relative order and nesting of regions is represented as a region-tree and it is shown in Figure 4.4. Regions containing whole loops are marked (L1 and L2 respectively). Regions that contain basic blocks numbered 2 and 8 are degenerate regions, since they are not part of any loop. The region marked as L1 includes a subregion L2 and three subregions executed in sequence. If a region L1 is a SCoP then it would contain two nested loops, which is the best On the other hand, let us consider an example in Figure 4.6. The corresponding region tree is shown in Figure 4.5. An original CFG represents a single nest loop with one basic block scheduled before the loop, and another one scheduled after the loop. The loop header basic block (number 3) contains nonaffine loop bounds, so it is marked as not satisfying affine conditions. Thus, the region that contains the loop is not marked as a SCoP. The remaining regions form three separate SCoPs: SCoP1, SCoP2 and SCoP3. Those SCoPs are degenerate and are not amenable to loop optimizations.

Mapping three-address intermediate representation into polyhedra

After outlining SCoP regions, each SCoP is processed so to map the three-address code contained inside region into the polyhedral representation. Only SCoPs that are containing at least one loop are processed.

As shown in Figure 4.7, the strongly connected components of the CFG subgraph that belongs to the region form the loops. Given any two loops, one of them is either completely nested in the other or they are disjoint. Each basic block is enclosed in at least one loop region. The nesting of the loop regions determines the loop depth of the basic block.

Three components (Section 2.3) of the polyhedral model representation are extracted from threeaddress code IR: The dimensionality of the iteration space is equal to the number of loops belonging to a SCoP and enclosing the basic block.

2. Schedules. The initial scheduling function θ B for each basic block B belonging to a SCoP is built. The initial scheduling function encodes the execution order of the basic block in an original program -before a transformation is applied.

3. Data references. Each basic block might contain memory data accesses in its statements. In a region complying with affine program constraints the only memory accesses are through array accesses whose access functions could be expressed as affine functions of loop induction variables.

For each basic block B the sets of write and read data references, W B and R B respectively, are constructed.

We will discuss the technical details of the polyhedral model representation extraction in the rest of this subsection. A discussion is based on an example CFG shown in Figure 4.7, obtained during the compilation of a matrix vector multiplication kernel shown in the same figure.

Reconstructing loop domains

SCoP extraction detects the static control region containing two nested loop regions: L1 and L2. Basic block B 4 is the loop header of the loop L2 in which the variable j is used to controls the exit condition of the loop. Since the code is in SSA form, variable j has two versions: j 1 and j 2 , corresponding to different assignments to the original j variable.

In order to detect the evolution of the variable j in the consecutive loop iterations of the loop L2, the TREC construction algorithm is used. TREC expression construction is done by following def-use chains of the induction variable. In SSA form, def-use chains are explicit, since there is an unique definition for each use.

The variable version j 1 that is used in a loop exit condition is defined in a statement S 7 : j 1 = j 2 + 1. This statement uses a variable version j 2 that is in turn defined in a φ-function in the statement S 1 : j 2 = φ( j 1 , 0). The statement S 1 uses the variable version j 1 that is defined in the statement S 7 . This defuse chain closes a cycle and a TREC expression of the form {0, +, 1} 2 is deduced. This is illustrated by dashed arrows inside the basic block B 4 in Figure 4.7.

A TREC expression of the form {0, +, 1} 2 means that the induction variable initial value is the integer constant 0 and that the induction variable is incremented by the step of 1 at each iteration of the loop L2.

The variable version j 1 is used in the conditional statement S 8 : j 1 < N. Combining this (affine) condition with a TREC an affine constraint on loop induction variable could be deduced: 0 ≤ j ≤ N -1. The same induction variable analysis is performed for the loop L1.

Affine expressions defining loop domains for all the basic blocks of the SCoP are shown in Figure 4.7.

Reconstructing data references

The property of the three-address code used in GCC is that each statement might contain at most one read or write data reference to a memory location. Furthermore the affine program constraint mandates that all the memory accesses could only be performed through arrays whose subscript functions are affine functions.

For each basic block that belongs to the SCoP all the statements are scanned. For each memory access statement the subscript functions are reconstructed by evaluating the evolution of induction variables that are used as array indices.

For the reconstruction of access functions the TREC expressions of the array indices are evaluated. For the basic block B 3 there is one memory access statement S 2 : b[i 2 ] = 0.0. The basic block B 3 belongs to the loop L1 and the evolution of the array index variable i 1 is expressed as a TREC {0, +, 1} 1 . This expression is interpreted as an affine access function f (i) = i. The complete set of write references of the basic block B 3 is:

f o r ( i = 0 ; i < N; i ++) { S1 : b [ i ] = 0 ; f o r ( j = 0 ; j < N; j ++) S2 : b [ i ] = b [ i ] + A[ i ] [ j ] * x [
j_1 = j_2 + 1; j_2 = PHI <j_1, 0> S1 S2 S3 S4 S5 S6 S7 S8 S1 S2 S3 S2 S1 B3 B4 B5 B6 goto <bb 4> loop L2 i_2 = PHI <i_1, 0> b[i_2] = 0.0; D0 = A[i_2][j_2]; D1 = x[j_2]; D2 = D1 * D0; prephi = PHI <D3, 0.0> D3 = D2 + prephi; return; B8 loop L1 if (i_1 < N) goto <bb 7>; else goto <bb 8>; i_1 = i_2 + 1; b[i_2] =D3; goto <bb 3>; B7 D B4 = {(i, j) | 0 ≤ i ≤ N -1 ∧ 0 ≤ j ≤ N -1} θ B4 (i, j) T = (0, i, 1, j, 0) T W B4 = ∅ R B4 = { A, f ((i, j) T ) = (i, j) T , x, f ((i, j) T ) = j } R B3 = ∅ W B3 = { b, f (i) = i } θ B3 (i) T = (0, i, 0) T D B3 = {(i) | 0 ≤ i ≤ N -1} D B6 = {(i) | 0 ≤ i ≤ N -1} θ B6 (i) T = (0, i, 2) T W B6 = { b, f (i) = i } R B6 = ∅
W B 3 = { b, f (i) = i }.
Symbol b is a symbolic reference to the base address of the b array.

The basic block B 4 belongs to the loop L2 (with iteration vector is (i, j) T ). Its statement S 3 :

D0 = A[i 2 ][ j 2 ]
represents a read memory array access. The TREC of the first subscript is {0, +, 1} 1 , while the TREC of the second subscript is {0, +, 1} 2 . Combining those two expression we get a multidimensional affine access function : f ((i, j) T ) = (i, j) T . The complete set of read references of B 4 is:

R B 4 = { A, f ((i, j) T ) = (i, j) T , x, f ((i, j) T ) = j }.
Writes/reads to/from scalar variables are not represented. Writes and reads to temporary scalar variables are treated as an internal state of the basic block. Since they are not producing visible side-effects they do not have to be captured. We will discuss this further in the section dealing with data dependence analysis. In Figure 4.7 the represented data references are marked (the statements that access the memory are emphasized).

Construction of the original schedules

The scheduling functions θ B encoding the original program execution order are constructed for each basic block belonging to a SCoP region of a CFG.

As it was defined in section 2.3.4, the scheduling matrix Θ B is decomposed into dynamic (A and Γ matrices) and static (β vector) scheduling components. For the scheduling matrices Θ B corresponding to the original execution order, their dynamic scheduling components, A and Γ matrices, are an identity matrix and all-zero matrices respectively. The static scheduling β vector for each basic block B is deduced from the region tree. The child nodes of each region are numbered. By following a path from the SCoP root region down to the given basic block, we construct a β vector that consists of the numbers labelling the child positions. This idea is described in [START_REF] Bastoul | Code generation in the polyhedral model is easier than you think[END_REF]. For example, the β vector for the basic block B 3 is (0, 0). The β vector for the basic block B 4 is (0, 1, 0). The construction of static scheduling components (β vectors) based on a region tree guarantees an uniqueness of the scheduling time-stamps for each basic block instance.

As a side note: the basic blocks that do not contain any computational statements are not considered for the polyhedral representation. Empty basic blocks are simply discarded from the region tree. In our example, the basic blocks4 B 5 and B 7 are not considered as a part of the polyhedral representation. Since the polyhedral model captures all the necessary information to reconstruct the control-flow according to the new schedule, the basic blocks containing only goto statements are ignored, since they are going to be recreated in the code generation stage of the GRAPHITE.

Mapping IR symbols into the polyhedra

The GIMPLE intermediate representation statements contain symbolic variable names. Variables are either scalars or arrays. On the other hand, the polyhedral representation is a mathematical abstraction representing sets of integer vectors 5 . The mapping between symbolic names of the three-address code and the variables of the polyhedra has to be maintained. For that purpose a mapping, similar to compiler symbol table, is used. An example of this correspondence is seen in Figure 4.7.

Differences from the source-to-source compilers

The source-to-source compilers [124,[START_REF]PLUTO: A polyhedral automatic parallelizer and locality optimizer for multicores[END_REF] directly map the source code level statements into the polyhedral model. For example, the polyhedral representation for the source code in Figure 4.7 would contain two domains D S 1 and D S 2 corresponding to the two syntactical statements S 1 and S 2 . Also, the scheduling functions are provided per source level statement:

θ S 1 , θ S 2 .
After the source code is translated into GIMPLE three-address code, the high-level statements are broken down into several three-address code statements -temporary variables are inserted to hold the intermediate results. Inside GRAPHITE we are storing the polyhedral representation per each basic block and not per each GIMPLE statement. Thus, we have three domains, one for each basic block: D B 1 , D B 2 and D B 3 . The same holds true for scheduling functions. We use a superscript B i to emphasize this fact.

The fact that we store the scheduling functions per basic block and not per three-address code statement contained inside basic block, means that all the statements inside one basic block share the same scheduling function. Statements inside the basic block are not rescheduled. In section 4.5.1 we will discuss the implications of this decision.

In the source code of a typical imperative language (C, fortran), one assignment statement contains one write (left hand side) and several reads (right hand side). After transforming the source level statement into three-address code the original write and read data references are scattered into several statements. The original data references might be scattered among several basic blocks, or even optimized out completely.

For example, the original source code shown in Figure 4.7 contains one write data reference (b[i]) in textual statement S 1 . The corresponding data reference in the three-address code (shown in Figure 4.7) is contained inside basic block B 3 in the statement S 2 . The two read data references contained in the statement S 2 (A[i][j] and x[j]) have their corresponding references inside the basic block B 4 (in statements S 3 and S 4 ). On the other hand, the write data reference (b[i]) from source-level statement S 2 appears in the basic block B 6 as statement S 1 . The original read from array b inside source-level statement S 2 does not appear in three-address code at all -it has been optimized out and replaced by a scalar access 6 .

D B 3 = (i) | 0 ≤ i ≤ N -1 θ B 3 (i) T = (0, i, 0) T W B 3 = { b, f (i) = i } R B 3 = / 0 D B 4 = (i, j) | 0 ≤ i ≤ N -1 ∧ 0 ≤ j ≤ N -1 W B 4 = / 0 R B 4 = { A, f ((i, j) T ) = (i, j) T , x, f ((i, j) T ) = j } θ B 4 (i, j) T = (0, i, 1, j, 0) T D B 6 = (i) | 0 ≤ i ≤ N -1 W B 6 = { b, f (i) = i } R B 6 = / 0 θ B 6 (i) T = (0, i, 2) T
D S 1 = (i) | 0 ≤ i ≤ N -1 θ S 1 (i) T = (0, i, 0) T W S 1 = { b, f (i) = i } R S 1 = / 0 D S 2 = (i, j) | 0 ≤ i ≤ N -1 ∧ 0 ≤ j ≤ N -1 θ S 2 (i, j) T = (0, i, 1, j, 0) T W S 2 = { b, f (i) = i } R S 2 = { A, f ((i, j) T ) = (i, j) T , x, f ((i, j) T ) = j , b, f (i) = i }

Putting it all together

The complete polyhedral representation of the three-address code fragment contained inside the SCoP from Figure 4.7 is shown in Table 4.9. Figure 4.7 shows the same information attached to the respective basic blocks.

As a comparison, Table 4.10 shows the (hypothetical) polyhedral representation of the respective source-level code in Figure 4.7.

The data dependence analysis of three-address code

The polyhedral model used in the source-to-source compilers [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF], [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part I: one dimensional time[END_REF] only allows memory references through array accesses. An access to a scalar value is treated as an access to an one-element array: a = ... is rewritten as a[0] = .... While this approach is practical for the source-to-source compilers, transforming the numerical kernels where array accesses constitute the majority of data references, it has prohibitive cost if used in the three-address code compilation context.

The lowering process from the source code down to the three-address code introduces (many) temporary scalar variables. Treating all the scalar temporary variables as single-element arrays and taking them into an account for the dependence analysis might lead to an explosion in the size of data dependence graph. This is of particular concern, if we take into an account that the scheduling algorithms [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF][START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part I: one dimensional time[END_REF], based on the polyhedral model, have a complexity that is dependent on the number of edges of a data dependence graph G = (V, E).

In order to efficiently manage the data dependences between scalar variables in the three-address code, the smallest necessary subset of the scalar variables is taken into a consideration for the dependence analysis.

Internally in GRAPHITE, those scalars whose data dependences must be exposed are marked as if they were single-element arrays. As a preprocessing step, GRAPHITE is marking those scalars that must to be processed by the dependence analysis as single-element arrays. Scalars that are not marked are not taken into the consideration for the dependence analysis.

In order to mark the necessary scalars for the dependence analysis, GRAPHITE is classifying the scalar data dependences into four categories:

1. Intra basic block scalar data dependences 2. Cross basic block scalar data dependences 3. Inductions

Reductions

Only the cross basic block and the reduction scalar dependences have to be explicitly represented. This is done by introducing the shadow single-element arrays for those scalars that are involved in dependences. Scalar data dependences that are induced by the computation of the induction variables and scalar data dependences that are contained within a single basic block are not represented explicitly in the polyhedral data dependence graph.

Given the three-address code in Figure 4.11, the effect of the preprocessing step that introduces the shadow arrays is shown.We will provide the brief discussion of the details related to handling each kind of scalar dependences in GRAPHITE. The illustration shown in Figure 4.11 will serve as common base example for all the discussions.

Intra basic block data dependences

As mentioned in Subsection 4.3.2, the basic block is the basic unit of scheduling. The statements inside a basic block are not rescheduled. The transformation process from the source-code to the lower level three-address code introduces temporary scalar variables to store the intermediate computation results.

Definition 4.4.1 (Intra basic block dependence).

A scalar data dependence between the scalar value X defined within a basic block B at iteration i and its subsequent uses within the same iteration i of the same basic block B are called intra basic block dependences Since the relative order of the statements inside a basic block is not changed by a transformation, intra basic block dependences are always preserved, regardless of the basic block schedules θ B . As a consequence, intra basic block dependences are not explicitly represented in the polyhedral dependence graph.

For example, consider scalar variables D0, D1, D2 defined in B 4 -they are defined and used within the same basic block and within the same iteration. The relative ordering of the statements S 3 , S 4 , S 5 and S 6 will not be changed in the polyhedral transformation. The flow dependence from statements S 3 and S 4 to statement S 5 would always be preserved. Please note that this approach is similar to the explicit privatization used for dependence removal [START_REF] Allen | Optimizing Compilers for Modern Architectures[END_REF].

L1 loop: B6 i_2 = PHI <i_1, 0> b[i_2] = 0.0; phi_out_of_ssa[0] = b[i_2]; j_1 = j_2 + 1; if (j_1 < N) goto <bb 5>; pre.3 = phi_out_of_ssa[0]; j_2 = PHI <j_1, 0> if (i_1 < N) goto <bb 7>; i_1 = i_2 + 1; D0 = A[i_2][j_2]; D1 = x[j_2]; D2 = D1 * D0; D3 = D2 + pre.3; phi_out_of_ssa[0] = D3; cross_bb_dep[0]=D3 D4 = cross_bb_dep[0]; b[i_2] =D4;

Cross basic block data dependences

Since the different basic block instances might be rescheduled independently, the data dependences that occur between scalar values belonging to two different basic blocks must be captured in order to restrict the possible rescheduling of the basic block instances.

If some scalar variable X is assigned in one instance of the basic block and if it is subsequently used in another basic block or subsequent iteration of the same basic block, that scalar value will get the shadow array that exposes the data dependence to the dependence analyser. Those scalars are treated as single element arrays, so that dependence analysis could capture them in the dependence test.

An example is given in

Inductions

Scalar variables representing the loop counters or auxiliary induction variables that could be expressed as TREC expressions are ignored for dependence analysis. The data dependences induced by computation of those variables could be ignored, since they could be represented by scheduling invariant expressions of the form:

χ(i) = χ(i 1 , i 2 , . . . , i n ).
Those expressions are reconstructed in the GIMPLE during the code generation step, after the transformation has been applied. Since the expression could be correctly reconstructed, regardless of the applied scheduling, there is no need to explicitly keep those dependences in the dependence graph.

Please note that in the source-to-source compilers, the induction variables are implicit and their computation is not taken into an account for the dependence analysis. We base our reasoning on the same principle, but in our case we have to explicitly discover those inductions inside a three-address code.

Reductions

Reductions [START_REF] Redon | Scheduling reductions[END_REF] form a special computational pattern used very often in the scientific codes. Indeed, an effective parallelization of loop kernels relies heavily on an effective detection and parallelization of reductions. Mostly often it is used to perform a summation over the values of an array.

Data dependence theory puts a constraint on the relative scheduling of statement instances. But some operations, if proved to be associative and commutative algebraic operations, might be scheduled in an arbitrary order.

If the operator along which the reduction happens can be proved to be commutative and associative, then the dependences that are induced by variables involved in a reduction are marked as belonging to reduction.

As an illustration, the data dependences that form a reduction cycle are shown in Figure 4.11. Indeed, the scalar D3 accumulates the values computed in the scalar D2 in each iteration of the loop. An addition operator '+' is commutative and associative, if performed on integer data type, so the order in which the summation is performed does not influence the correctness of the final results.

Putting it all together

Taking into an account all the dependences that have to be explicitly represented, we come up with the data dependence graph shown in Figure 4.12. Please note that the set of vertices is composed from the three-address code statements, but those statements are contained within basic blocks, which is rep- As a comparison, the dependence graph that is built by a source-to-source compiler is shown in Figure 4.13. The dependence graph of the three-address code contains, besides true data dependences, many anti and output data dependences -those are jointly classified as memory based dependences. These dependences might significantly constrain the space of legal transformations. Later in Chapter 5 we will show how to efficiently deal with those kind of dependences, such that the space of legal polyhedral transformations on the three-address code is equivalent to the source-level polyhedral compilation.

B3:S2 B6:S2 B4:S8 B6:S1 B3:S3 B4:S2 -true -output -anti B4:S7 δ o δ o δ o δ δ -1 δ -1 δ δ -1 δ -1 δ δ δ o δ δ o δ -1

Transformations

The loop transformations in the GRAPHITE framework are performed on the polyhedral representation. In order to maintain the composability of loop transformations, we perform all the loop transformations through the scheduling matrix of the polyhedral model. This is in contrast to the traditional [START_REF] Cohen | Facilitating the search for compositions of program transformations[END_REF][START_REF] Vasilache | Violated dependence analysis[END_REF][START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF] approach to transformations in the polyhedral model, where some transformations are expressed through scheduling matrices (affine transformations), where the loop tiling and strip-mining transformations require the modification of the iteration domains. We will later show how we do overcome this non-homogeneity in expressing the loop transformations.

Once the analysed program is brought into the polyhedral model representation and after the desired transformations are selected, the output code is generated from the transformed schedules. The automatic search for the transformation is conducted according to the cost-model driven search approach shown in chapters 7, 5 and 6.

The state of the art polyhedral code generator CLooG [START_REF] Bastoul | Code generation in the polyhedral model is easier than you think[END_REF][START_REF] Bastoul | Cloog: The chunky loop generator[END_REF] generates the loops that scan the integer vectors contained in polyhedra. The integer vectors are scanned in the lexicographic order specified by the scheduling function. The code generator is not responsible for maintaining the legality of the transformations -iterations are scanned according to the transformed scheduling matrices. It is the duty of the violated dependence analysis, discussed in Chapter 5, to guarantee the legality of the transformations.

The CLooG code generator generates an AST (abstract syntax tree) representing the newly generated loop nests, conditionals and affine loop bounds and conditions. This structure is traditionally used in the source-to-source compilers [START_REF]PLUTO: A polyhedral automatic parallelizer and locality optimizer for multicores[END_REF]124] to generate the target source code. In the case of GRAPHITE the AST is used to generate the output CFG that corresponds to new loops. The symbol table (as shown in Figure 4.2) that maps the original three-address code into polyhedral model variables is now used to regenerate the target three-address code.

All the loop transformations, including loop tiling, are expressed seamlessly through scheduling functions and it is the duty of CLooG to generate the appropriate loop structures. GRAPHITE is interacting with parallelizer and vectorizer (as shown in Figure 4.1) passes of GCC compiler: if the loop is detected to be a parallel DOALL loop, it is marked as a such in the code generation phase. After GRAPHITE, the marked loop will be processed by the parallelizer to generate the multithreaded code. The similar interaction happens with vectorizer pass.

We will now discuss several design issues related to the transformations and code generation in our framework.

Granularity of scheduling

As we have mentioned in Chapter 2 (subsection2.3.1), the basic building block of the polyhedral representation could be chosen according to the abstraction level on which the polyhedral model is to be built.

In the source-to-source compilers, the basic, and the most natural, atomic unit of scheduling is a source-level statement. In a three-address-code based polyhedral compiler, such as GRAPHITE, it is not practical to build a polyhedral representation for each three-address code statement. The most reasonable choice is to have a basic block of the control flow graph as the basic scheduling unit.

Nevertheless, the choice of having a basic block as the basic unit of scheduling is not always optimal and may be way too restrictive. Let us consider an illustration in Figure 4.14. The source-level code for the ADI 7 kernel is shown. It contains two syntactic statements: S 1 and S 2 . The corresponding threeaddress code is shown as well. One can observe that all the three-address statements are clustered in one basic block, the basic block B 4 .

Since the design choice is to have a scheduling function per basic block, all the three-address code statements contained within the same basic block would get the same schedule. Obviously, this is not desirable, especially in the case when we want to apply per statement affine scheduling [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF]. But having one scheduling function per each three-address code statement is non-practical, for the scheduling scalability reasons.

A solution that we have applied is to split such a basic block into clusters, such that each cluster would contain one write statement (highlighted in Figure 4.14) and all the necessary temporary variables that are produced within the same basic block. This could easily be achieved by following use-def chains in the SSA form. After this preprocessing, each cluster becomes a separate basic block that could be scheduled independently. [START_REF] Allen | Optimizing Compilers for Modern Architectures[END_REF]. taken from Polybench benchmark suite

f o r ( i 1 = 0 ; i 1 <N ; i 1 ++) { f o r ( i 2 = 1 ; i 2 < N ; i 2 ++) { S1 : X[ i 1 ] [ i 2 ] = X[ i 1 ] [ i 2 ] -X[ i 1 ] [ i 2 -1] * A[ i 1 ] [ i 2 ] / B[ i 1 ] [ i 2 -1 ] ; S2 : B[ i 1 ] [ i 2 ] = B[ i 1 ] [ i 2 ] -A[ i 1 ] [ i 2 ] * A[ i 1 ] [ i 2 ] / B[ i 1 ] [ i 2 -1 ] ; } } goto <bb 4> return; loop:l1 l2 loop: # i2_38 = PHI <i2_22(3), 1(7)> D.5548_5 = X[i1_37][i2_38]; D.5549_6 = i2_38 + -1; D.5550_7 = X[i1_37][D.5549_6]; D.5551_8 = A[i1_37][i2_38]; D.5552_9 = D.5550_7 * D.5551_8; D.5553_11 = B[i1_37][D.
# i1_37 = PHI <0(2), i1_23(6)> i1_23 = i1_37 + 1; if (i1_23 != 500) goto <bb 6>; else goto <bb 8>; goto <bb 7>; S1: X[i1][i2] = X[i1][i2] -X[i1][i2-1] * A[i1][i2] / B[i1][i2-1]; B[i1][i2] = B[i1][i2] -A[i1][i2] * A[i1][i2] / B[i1][i2-1]; S2: B7 B4 B3 B5 B8 B6

Functions represented as relations

The classical definition (Chapter 2, subsection 2.3.3) of the polyhedral model components states that the data accesses and schedules are expressed as data access functions and scheduling functions:

f : i ∈ D → s ∈ D θ : i ∈ D → t ∈ D
The data access function is the one-to-one mapping from an iteration space (i ∈ D) to a data space of array subscripts (s ∈ D). The scheduling function is the one-to-one mapping from an iteration space (i ∈ D) to a time domain of multidimensional timestamps (t ∈ D).

In GRAPHITE we extend this definition with a notion of access relation and scheduling relation:

F : (i, s) ∈ D × D θ : (i, t) ∈ D × D
The mapping from the iteration vector to the subscript function is not anymore restricted to one-toone mappings. This enables us to represent the full access regions -when the data reference information could not be expressed as a direct affine function. This is the case if only the approximation of the memory access is available, coming from an interprocedural analysis for example. The notion of access regions is a known term used in interprocedural compilation [START_REF] Triolet | Direct parallelization of call statements[END_REF]. The access relation is now expressed by an access polyhedron:

F = (i, a, s) | F × (i, a, s, g, 1) T ≥ 0
An additional component, the alias set number 'a', captures points-to information (pointer aliasing); it allows us to extend the notion of array memory access to pointers accesses, provided that the precise aliasing information [START_REF] Pearce | Efficient field-sensitive pointer analysis of C[END_REF][START_REF] Andersen | Program analysis and specialization for the c programming language[END_REF] could be obtained and that the original access functions (or regions) could be reconstructed, as shown in Section 4.3.2.

Similarly, the scheduling is represented as a scheduling relation. 8 represented by an scheduling polyhedron:

θ = (t, i) | Θ × (t, i, g, 1) T ≥ 0 .
Having the schedules expressed as relations, we can express the tiling transformation with a scheduling relation alone, without modifying the statement iteration domain D S . This is of particular importance for the transparent integration of analytical cost-modelling approach discussed in Chapter 6 and lazy memory expansion scheme shown in Chapter 5.

Adding new scheduling dimensions -loop tiling

Loop tiling or loop blocking is an extremely important loop transformation [START_REF] Wolfe | Iteration space tiling for memory hierarchies[END_REF][START_REF] Irigoin | Supernode partitioning[END_REF][START_REF] Ahmed | Tiling imperfectly-nested loop nests[END_REF] used both for coarsening the grain of parallelism and data locality improvement [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF]. An example of the tiling applied to a 2d iteration domain is shown in Figure 4.15. The new iterators ii and jj are introduced to express the tile-space, while the intra-tile iterations are constrained to the boundaries of the containing tile.

In the classical polyhedral framework, this transformation is expressed both as an iteration domain transformation (increasing the loop nest depth) and scheduling transformation(permuting the scheduling dimensions) [START_REF] Cohen | Facilitating the search for compositions of program transformations[END_REF][START_REF] Vasilache | Violated dependence analysis[END_REF][START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF]. But this formulation breaks the principle of transformation composabilitythe ability to express an arbitrary sequence of loop transformation through scheduling function only.

We postulate that it is possible to express the loop tiling transformations by modifying the scheduling relation θ alone. Given a timestamp variable t that corresponds to some scheduling dimension of the dynamic scheduling component of the scheduling matrix (Section 3.1.2), we introduce a new time variable tt, called supernode [START_REF] Irigoin | Supernode partitioning[END_REF] time variable, by extending the scheduling relation with a new dimension and introducing an inequality: 8. also called scattering relations following CLooG [START_REF] Bastoul | Code generation in the polyhedral model is easier than you think[END_REF] code generator terminology

B • tt ≤ t ≤ B • tt + (B -1)
where B is the compile-time 9 known loop blocking factor.

The new supernode time variable tt denotes the B-sized blocks of the original iteration domain, while the original time variable t is now constrained to scan the iterations within the block denoted by tt.

As an example, we will consider the representation of the strip-mining, which is the basic building block for implementing tiling -in order to get a loop tiling, the strip-mining has to be performed along each dimension considered for tiling. Consider a code snippet:

f o r ( i = 0 ; i < N ; i ++) S1 ;
The original scheduling matrix in the canonical form could be decomposed into the following components:

A S 1 = 1 Γ S 1 = 0 β S 1 = 0
The concise representation of the scheduling relation is:

θ S 1 = (t 1 ,t 2 ,t 3 , i) | t 1 = 0 ∧ t 2 = i ∧ t 3 = 0
The strip-mining of the t 2 time dimension with a factor B = 64 is expressed as follows:

θ S 1 stripmined = (t 1 ,tt 2 ,t 2 ,t 3 , i) | t 1 = 0 ∧ t 2 = i ∧ 64tt 2 ≤ t 2 ≤ 64tt 2 + 63 ∧ t 3 = 0
The corresponding code generated by CLooG [START_REF] Bastoul | Code generation in the polyhedral model is easier than you think[END_REF] code generator is shown:

9. Thus, we do not support parametric tiling as [START_REF] Hartono | Parametric multi-level tiling of imperfectly nested loops[END_REF] The floord stands for the floor function of integer division. Please note that we show the syntactic code, but only for the presentation purposes. Internally in GRAPHITE, the loops are generated as basic blocks within a CFG. Please note that the original iteration domain D S 1 of the statement S 1 is not modified in any way, even though the iteration domain of the generated code is changed, but this change was expressed solely through scheduling relation.

Conclusions and future work

In this chapter we have shown a detailed design of the GRAPHITE polyhedral compilation framework which is currently a standard component of the general purpose GCC [START_REF]Gnu compiler collection (gcc[END_REF] compiler suite. The author is a direct contributor to this project.

Two proprietary polyhedral compilers based on a low-level internal representation are currently in development: the R-Stream compiler from Reservoir Labs [START_REF] Meister | Productivity via automatic code generation for pgas platforms with the r-stream compiler[END_REF], and IBM's polyhedral extension of its XL compiler suite [START_REF] Renganarayana | Compact multi-dimensional kernel extraction for register tiling[END_REF]. But little has been published on the specific problems and trade-offs that arise when introducing the polyhedral model based compilation flow into the realm of low-level three-address code based compiler.

Our work is the first widely published contribution [START_REF] Trifunovic | Graphite two years after: First lessons learned from real-world polyhedral compilation[END_REF][START_REF] Trifunovic | Enabling more optimizations in graphite: ignoring memory based dependences[END_REF][START_REF] Trifunovic | Polyhedral-model guided loop-nest auto-vectorization[END_REF] on the specifics of the direct polyhedral compilation on a three-address code representation.

Benefits of direct manipulation of the three-address code

The traditional, source-to-source based approach to the polyhedral compilation regards the syntactical statement as the basic atomic unit of the computation, abstracting away the internal state of the computation and giving relevance to the loads/stores of the statement only.

On the other hand, the intrinsic nature of the three-address code is the explicit visibility of the statethe intermediate results of all the computations are visible through intermediate temporary variables.

This intrinsic property of the three-address code is the source of many difficulties in applying the polyhedral model concept, but it also gives many degrees of freedom that could be exploited. Those degrees of freedom go along two axes: the granularity of the scheduling (statement versus basic-block) and the granularity of visible state (intra versus inter basic data block dependences).

The low-level nature of the three-address code enables us to consider a precise modelling of the actual instruction costs -a property that is out of reach of the source-to-source compiler [124,[START_REF]PLUTO: A polyhedral automatic parallelizer and locality optimizer for multicores[END_REF]. We will heavily exploit this property in Chapter 6, where we present the first-of-a-kind approach to precise cost modelling of SIMD vectorization that is low-level, machine specific in nature.

We will summarize the benefits of providing the direct polyhedral compilation of the three-address code:

Additional degrees of freedom. Since the original high-level syntactic computations are broken down into the atomic three-address code instructions, there are additional degrees of freedom with respect to the scheduling granularity (which statements should be clustered together) and exposing data dependences (whether to keep the data dependence internal to the basic block, or to externalize it).

Tight interaction with a compiler. This benefit is more of a technical nature: since we are operating within the pipeline of the compiler optimization passes, all the analysis information like aliasing, pointer analysis, data alignment, induction variable analysis, interprocedural analysis is available for the use within the polyhedral transformation framework.

Precise modelling of instruction costs. The three-address code is much closer to the machine code of the target architecture than an abstract syntax tree of a source program. Also, the recognition of the intrinsics (for vectorization) and idioms [START_REF] Pottenger | Idiom recognition in the polaris parallelizing compiler[END_REF] is readily available at the level of three-address code. The actual instruction costs are available from the backend of the compiler, since they are anyway used in the instruction scheduling phase. We will heavily rely on this property in Chapter 6.

Semantical transparency. The three-address code is expressed in SSA [START_REF] Cytron | An efficient method of computing static single assignment form[END_REF] form which captures the essential flow of the scalar data in the program, and could be regarded as the functional representation of the computation [START_REF] Appel | Ssa is functional programming[END_REF]. Also, the essential scalar optimizations based on the SSA form are performed before the GRAPHITE stage is reached.

This leads to a degree of semantical transparency with respect to the form of the input program : two semantically equivalent kernels would be represented in the same way in SSA form, when they reach the polyhedral framework. This is in contrast with source-to-source compilers, which are fragile with respect to syntactical details -if the programmer forgot to bring the program in the form expected by the syntactical parser, the program would simply be ignored.

Future work

The aforementioned benefits of having the polyhedral transformations within the three-address code compiler bring many new research opportunities.

We have used GRAPHITE framework to conduct a research on the iterative search and cost modeling of the profitability of the loop transformations, as explained in subsequent chapters of this dissertation.

But many new possible and interesting research opportunities emerged. Some problems were left unsolved, since their in-depth investigation was not on the main line of this thesis. We will summarize several of those research problems:

Interplay with other compiler optimizations. The GRAPHITE polyhedral framework is used to perform the loop transformations expressed as scheduling relations. But the loop transformations are just a part of the bigger picture. GRAPHITE is a part of the complete compilation pipeline consisting of more than 200 optimization passes -scalar optimizations, instruction selection and scheduling, register allocation, loop unrolling and much more. An in-depth study of the possible interactions between loop optimizations and other optimization seems very interesting, though intractable in its full generality .

Some optimizations might conflict with the polyhedral transformations. A notorious example is loop unrolling, that might be performed before GRAPHITE. Obviously, unrolling a loop before going into the polyhedral representation can have a dramatic effect on the scalability of the polyhedral optimizations, and it is better avoided at all. Data layout transformations. The expressive power of the polyhedral model is used in the context of data-layout transformations [START_REF] Lu | Data layout transformation for enhancing data locality on nuca chip multiprocessors[END_REF] as well. At the same time, those transformations are studied independently, as the components of other optimization stages of the compiler [START_REF] Golovanevsky | Trace-based data layout optimizations for multi-core processors[END_REF]. Combining both might be beneficial, as the polyhedral model precisely captures the execution order, while the internal representation of the compiler has a low-level details about the data-layout. Also, the controlled memory expansion schemes like [START_REF] Lefebvre | Automatic storage management for parallel programs[END_REF] are a key missing part for fully enabling the potential of automatic parallelization with GRAPHITE.

Extending the scope of the analyzable programs. GRAPHITE restricts its scope of analyzable and transformable programs to the static control programs, much in the same way as the traditional polyhedral compilation approaches. Recent works [START_REF] Benabderrahmane | The polyhedral model is more widely applicable than you think[END_REF] on extending the applicability of the polyhedral model to non-static control programs (loop bounds and conditionals that depend on the input data) seem a promising research direction for the GRAPHITE. This is particularly true for the three-address code Chapter 5

Lazy memory expansion scheme

A program transformation needs to be safe -the semantics of the original imperative program cannot be changed. In order to preserve the legality, data-dependences (Section 3.2.1) need to be analyzed. There are essentially two types of data-dependences: data-flow dependences and memory-based dependencesas discussed in Section 3.2.2. While preserving the data-flow dependences is always mandatory, spurious memory-based dependences can be removed, or simply ignored in certain conditions. So called true or dataflow [START_REF] Lefebvre | Automatic storage management for parallel programs[END_REF] dependences are imposing ordering constraints between write and read operations -they should always be preserved, since this preserves the right producer-consumer ordering, which in turn guarantees the correctness of the computations.

Memory-based dependences 1 are induced by the reuse of the same memory location to store multiple, temporary values. Spurious scalar dependences not only increase the total number of dependences in the RDDG graph, but, most importantly, they reduce the degrees of freedom available to express effective loop transformations and parallelization.

Memory-based dependences could be removed by introducing new memory locations, i.e. expansion of the data structures [START_REF] Cohen | Parallelization via constrained storage mapping optimization[END_REF]. While the expansion approaches can remove spurious memory-based dependences, they have to be avoided whenever possible due to their detrimental impact on cache locality and memory footprint -expanding a scalar value into an array requires extra memory storage, proportional to the size of the iteration space.

Designing a polyhedral compilation framework on a three-address code exacerbates the problem of spurious memory dependences even further, since the gimplification 2 process introduces many temporary variables in the internal representation.

In this chapter, we show a technique that detects those memory-based dependences that could be ignored when checking for a legality of the transformation. If a memory-based dependence could not be ignored -it is violated in another words -it can be removed proposes an expansion or privatization that would enable a transformation.

Using a proposed technique, we can get rid of many memory-based dependences, either those induced by the lowering of a source program into three-address-code or those introduced by a programmer. Only if a dependence cannot be ignored, it is removed through scalar/array expansion.

1. anti and output dependences [START_REF] Allen | Optimizing Compilers for Modern Architectures[END_REF] 2. gimplification is a GCC jargon term denoting the lowering of a high-level AST into a low-level GIMPLE internal GCC representation

State of the art

The problem of eliminating the memory-based (anti and output) dependences has been studied since a long time in the automatic parallelization community [START_REF] Kuck | Dependence graphs and compiler optimizations[END_REF][START_REF] Padua | Advanced compiler optimizations for supercomputers[END_REF]. It is a well known fact that memorybased dependences hinder the opportunities for automatic parallelization and other loop transformations. Darte [START_REF] Calland | On the removal of anti-and output-dependences[END_REF] gives a good overview on the removal of the anti and output data dependences.

The simple and most natural solution is to perform a memory expansion -assigning a separate memory location for each executed iteration [START_REF] Feautrier | Array expansion[END_REF]. A similar approach is that of privatization -assigning a private copy to each thread executing in parallel [START_REF] Tu | Automatic array privatization[END_REF].

While the approach of memory expansion or privatization helps the parallelism, it has an inherent cost -an additional memory footprint or additional copy in/copy out operations. The cost of those additional overheads might be prohibitive. What is more -there is no guarantee that removing the memory-based dependences will help extract more parallelism.

In order to control the cost of the array expansion footprint, the two general approaches are:

-Perform a maximal expansion [START_REF] Barthou | Maximal static expansion[END_REF], apply a transformation, and then do an array contraction [START_REF] Darte | New complexity results on array contraction and related problems[END_REF]4,[START_REF] Darte | Lattice-based memory allocation[END_REF][START_REF] Gao | Collective loop fusion for array contraction[END_REF] which minimizes the memory footprint. Approaches like [START_REF] Maydan | Array-data flow analysis and its use in array privatization[END_REF][START_REF] Lefebvre | Automatic storage management for parallel programs[END_REF][START_REF] Cohen | Storage mapping optimization for parallel programs[END_REF] fall into this category. This approach gives the maximal degree of freedom for parallelization or loop transformation, but an array contraction phase is not always capable of optimizing the memory footprint. -Control the memory expansion phase by imposing constraints on the scheduling. Approaches like [START_REF] Cohen | Parallelization via constrained storage mapping optimization[END_REF], [START_REF] Thies | A unified framework for schedule and storage optimization[END_REF] fall into this category. This category of approaches tries to optimize the memory footprint, but it might restrict schedules, thus loosing optimization opportunities.

Our contribution

There is a trade-off between parallelization and memory usage: if we expand maximally, we will get the maximal degree of freedom for parallelization and loop transformations, but with a possibly huge memory footprint. If we choose not to expand at all, we will save memory, but our parallelization or loop transformation possibilities would be limited.

Our approach takes the lazy-expansion strategy: we do not expand memory before transformation. When checking for transformation legality, we simply ignore all memory based dependences. Only after applying a given transformation, we perform a violation analysis to check which memory based dependences might have been violated, and we propose to expand memory or to change a schedule.

By taking our approach, we are combining the best from the two mentioned state of the art approaches: we do not perform a full expansion before a transformation and we do not restrict the transformation too early.

The lazy-expansion strategy is not used to compute transformations automatically, as it is done in techniques using linear programming approach [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part II: multidimensional time[END_REF][START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF], instead, it is used in an iterative enumeration of possible schedules as it is done in [START_REF] Trifunovic | Polyhedral-model guided loop-nest auto-vectorization[END_REF].

Motivating example

Consider a matrix multiplication numerical kernel given in a Figure 5.1. If we analyze this kernel at the source-code level, we will get the dependence graph shown in Figure 5.7. It contains both dataflow 3 , and memory-based 4 dependences. Those data dependences do not prevent the rescheduling transformation that permutes the 'i' and 'j' loops.

3. true, a.k.a. read-after-write 4. write-after-write and write-after-read f o r ( i = 0 ; i < N ; i ++) f o r ( j = 0 ; j < N ; j ++) { S1 : illegal execution order

A[ i ] [ j ] = 0 ; f o r ( k = 0 ; k < N ; k ++) S2 : A[ i ] [ j ] += B[ i ] [ k ] * C[ k ] [ j ] ; } Figure 5.1 -matrix multiplication kernel f o r ( i = 0 ; i < N ; i ++) f o r ( j = 0 ; j < N ; j ++) { S1 : t = 0 ; f o r ( k = 0 ; k < N ; k ++) { S2 : t += B[ i ] [ k ] * C[ k ] [ j ] ; } S3 : A[ i ] [ j ] = t ; }
(2,2) (2,1) (1,2) (1,3)
(2,3)

Figure 5.4 -Illegal execution order

If the same kernel was written in an optimized way, shown in Figure 5.2, the dependence graph would now contain new output dependences due to the reuse of the scalar variable t for storing intermediate results.

The similar thing would happen if a compiler, such as [START_REF]Gnu compiler collection (gcc[END_REF], performs a lowering of the high-level source code into a three-address code internal representation. After the source code is transformed into low-level form, the compiler performs many optimization passes. One of those passes is PRE [START_REF] Muchnick | Advanced Compiler Design and Implementation[END_REF] (Partial Redundancy Elimination) which does the following scalar optimization: instead of accumulating a values into an array, it initializes a scalar value, accumulates values into that scalar and then stores the scalar into an array element. Conceptually, the idea is the same as if the user has rewritten the code as shown in Figure 5.2. An equivalent three-address code seen in GCC compiler is shown in Figure 5.5.

A data dependence graph corresponding to code in Figure 5.5 is shown in Figure 5.6. After introducing a scalar into the loop, a new write-after-write dependence on statement S 1 has been introduced: δ WAW S 1 →S 1 . This dependence stems from the fact that the same temporary scalar value is overwritten in each iteration of the containing loop.

This output dependence has to be respected, which forces the sequential execution in the original scheduling order. Figure 5.3 shows that if we execute the code in a sequential manner, according to the original loop nesting (loop i as outermost, loop j as innermost), then the dependences would be preserved. If we try to interchange loops i and j, we would invert a dependence constraint, thus violating the write-after-write dependence on the scalar t. This is shown in Figure 5. 4.

But an intuition tells us that it is legal to interchange loops i and j and still have a correct output code. An essential observation is that some memory based dependences (write-after-write and write-after-read) could be ignored when performing some transformations. But how do we determine when it is safe to On the other hand, let us consider an example where an user has written a code in Figure 5.2 manually and provided it to the compiler. As an example let us consider the legality of loop distribution transformation shown in Figure 5.8. This transformation would not be legal, and an output dependence on the scalar t could not be ignored. Though, an expansion of the scalar t could be performed, as shown in Figure 5.9. Expansion removes an output dependence on the scalar t and allows the distribution transformation to be legal. But an expansion introduces an additional cost -that of transforming the scalar t into an array t[N][N]). This is a situation that should be rather avoided.

δ W AR S2→S1 δ RAW S2→S6 δ W AR S6→S2 δ W AW S1→S1 δ W AW S5→S5 δ RAW S1→S2 δ RAW S5→S7 δ W AR S7→S5 δ W AW S8→S0 δ W AW S6→S6 Figure 5.6 -Data Dependence Graph S1 S2 δ W AW S 2 →S 2 δ RAW S 1 →S 2 δ RAW S 2 →S 2 δ W AR S 2 →S 2 δ W AW S 1 →S 2
In the following section we will show how to formally prove which memory-based dependences, given a desired transformation, could be ignored (as in the example of loop interchange ) and which must be removed by a memory expansion. f o r ( i = 0 ; i < N; i ++) f o r ( j = 0 ; j < N; j ++) t = 0 ; f o r ( i = 0 ; i < N; i ++) f o r ( j = 0 ; j < N; j ++)

{ f o r ( k = 0 ; k < N; k ++) { t += B[ i ] [ k ] * C [ k ] [ j ] ; } A[ i ] [ j ] = t ; } Figure 5.8 -illegal -loop distribution without an ex- pansion f o r ( i = 0 ; i < N ; i ++) f o r ( j = 0 ; j < N; j ++) t [ i ] [ j ] = 0 ; f o r ( i = 0 ; i < N ; i ++) f o r ( j = 0 ; j < N; j ++) { f o r ( k = 0 ; k < N; k ++) { t [ i ] [ j ] += B [ i ] [ k ] * C[ k ] [ j ] ; } A[ i ] [ j ] = t [ i ] [ j ] ; }

Framework

We will first define the concept of live range and we will show the way of computing the closed form expression that summarizes the set of live ranges. The core of the approach will be shown in subsection 5.3.2, where we discuss a computational procedure for checking the potential violation of value live range instances: if a value live range instance of some memory location M is violated, then the memory based dependence induced by the reuse of the memory location cannot be ignored. A detailed example is shown in subsection 5.3.3.

Live ranges

An execution trace of a sequential program can be seen as an interleaving of read and write instructions. During an execution of a loop, values are produced and stored into memory locations. Those values are then read by subsequent instructions. The value stored in a memory location is live until its last read, before it is destroyed by a subsequent write to the same memory location.

We are interested in formally analyzing and modelling value live ranges inside loops, using the polyhedral model. Given a code in Figure 5.2 we can represent the execution trace and instances of live ranges as shown graphically in Figure 5.10. A similar concept, named utility span of a value, was introduced in [START_REF] Lefebvre | Automatic storage management for parallel programs[END_REF].

Consider a value V written into a memory cell M[b] by an instruction instance w = (S W , i W ). We can compute a set of instruction instances R = {(S R , i R )} such that there is a direct data-flow of value V from instance w to instances in set R.

Each instance of a value live range is a tuple, describing a write and read instruction instances. Let us consider a set of value live range tuples:

L = {< (S W , i W ), (S R , i R ) >}
We want to have a closed form expression that summarizes all instances in this set. We can decompose the set L into a finite number of convex polyhedra, where each polyhedron describes live range instances for a pair of statements:

λ S W →S R = {(i W , i R ) | Λ × (i W , i R , g, 1) T ≥ 0}
Each convex polyhedron λ S W →S R represents instances of statements that form a definition/use pairs. This polyhedron is constructed by enforcing the following conditions: conflict condition -write and read statement instances refer to the same memory location:

F S W (i W ) = F S R (i R )
causality condition -a read instruction is scheduled after a write instruction:

θ S w (i W ) ≺ θ S r (i))
liveness condition -there is no intervening write w k = (S KW , i KW ) happening between (S LW , i LW ) and (S R , i LR ) instruction instances.

Described polyhedra satisfies exactly the same conditions that an array dataflow analysis (Section 3.2.3) requires. It is computed by the array dataflow algorithm of Feautrier, detailed in [START_REF] Feautrier | Dataflow analysis of scalar and array references[END_REF]. Given a memory location M, all the read accesses and their source functions (Section 3.2.3) are combined into the set L. A practical implementation of this computation is given in [START_REF] Trifunovic | Enabling more optimizations in graphite: ignoring memory based dependences[END_REF].

Live range violation analysis

After applying a transformation, the relative execution order of statement instances might change. This change might cause live ranges described in the set L, mapped to the same memory location, to interfere.

The classical dependence analysis theory states that a source instruction instance must execute before a sink instruction instance:

θ ′ S W (i W ) ≺ θ ′ S R (i R ).
Once dataflow dependences are computed, we use a simple dependence violation analysis [START_REF] Vasilache | Violated dependence analysis[END_REF] to check this condition on all dataflow dependences.

Since we propose to ignore the memory based dependences, how do we assure that a transformation would not reschedule some statement instance, say (S KW , i KW ), to overwrite a value V , produced at instance (S W , i W ), before it is read by a statement instance < S R , i LR >, as it was expected in the original program?

In other words: we have to check that such a case may not happen after a transformation. We form the following system of equations, modelled using a polyhedral model:

Vio S W →S KW →S R = (i W , i R ) ∈ λ S W →S R : θ ′S W (i W ) ≺ θ ′S KW (i KW ) ≺ θ ′S R (i R )∧ F S W (i W ) = F S R (i R )
New schedules θ ′ represent a polyhedral transformation that we check for legality. If a violation set Vio S W →S KW →S R is empty, there is no violation after a transformation. We have to perform this check for each dataflow dependence polyhedron δ S j →S i and for each possible write statement S KW .

Legality of loop parallelization:

Previously mentioned check takes into an account the fact that that schedule functions θ ′ are describing a sequential execution order. Sequential execution order has a property that no two different statement instances could be scheduled at the same time: so either θ

′ S i (i) ≺ θ ′ S j (j) or θ ′ S j (j) ≺ θ ′ S i (i).
If we consider loop parallelization transformation, we need to take into an account that some statement instances might be executed simultaneously θ ′ S i (i) = θ ′ S j (j). Thus, value live range legality violation check used for sequential transformations could not be used.

Given a loop at level k that we consider for parallelization 5 , we proceed as follows: we check that there are no loop carried [START_REF] Allen | Optimizing Compilers for Modern Architectures[END_REF] dataflow dependences at level k. Then, we check that no two distinct statement instances θ ′S i (i) and θ ′S j (j) sharing the same prefix up to depth k -1 write to the same memory location (instances associated with the same iteration of the parallel loop but different iterations of some inner loops can still write to the same memory location).

Supporting array/scalar expansion:

Our approach is compatible with well known array/scalar expansion approaches. If a transformation produces at least one violated live range (the set Vio is not empty) then we can choose to expand the variable M whose live ranges are violated. A precise characterization of violated live range instances in a set Vio could be used to drive the needed degree of expansion as in [START_REF] Lefebvre | Automatic storage management for parallel programs[END_REF]. If we do not want to perform any expansion, we can iteratively choose a new schedule that does not violate any live ranges.

Supporting privatization: Privatization is a common concept in the loop parallelization community. We can use our framework to automatically detect which scalars/arrays need to be privatized to enable loop parallelization transformation. This is actually how GRAPHITE is driving autopar -a GCC loop parallelizer.

An example

Let us take the GIMPLE code from Figure 5.5 and consider a memory location phi_out_of_ssa. Figure 5.10 shows an interleaving of writes and reads to this memory location. A slice of the execution trace, for a finite number of iterations, is shown. Value live ranges are shown as well. Some value live range instances contained in a set L phi_out_of_ssa : < (S 1 , (0, 0)), (S 2 , (0, 0, 0)) > < (S 6 , (0, 0, 0)), (S 2 , (0, 0, 1)) > < (S 6 , (0, 0, 1)), (S 2 , (0, 0, 2)) > . . . < (S 6 , (0, 0, N -2)), (S 2 , (0, 0, N -1)) > < (S 1 , (0, 1)), (S 2 , (0, 1, 0)) > . . . After array dataflow analysis, we come up with two closed form expressions: λ S 1 →S 2 and λ S 6 →S 2 . These two polyhedra summarize value live range instances between statements S 1 and S 2 , and between S 6 and S 2 respectively. They have the following form:

λ S 1 →S 2 = {< (i, j), (i ′ , j ′ , k ′ ) >: i ′ = i ∧ j ′ = j ∧ k ′ = 0 ∧ 0 ≤ i < N ∧ 0 ≤ j < N} λ S 6 →S 2 = < (i, j, k), (i ′ , j ′ , k ′ ) >: i ′ = i ∧ j ′ = j ∧ k ′ = k + 1∧ 0 ≤ i < N ∧ 0 ≤ j < N ∧ 0 ≤ k < N -1
For the purpose of this example, we would like to check whether interchanging loops i and j is a transformation that preserves non-conflicting condition on all live range instances. Referring again to Figure 5.5, we see that we are interested in the schedule of statements S 1 , S 2 and S 6 . Their scheduling functions in an original program are as follows:

θ S 1 (i, j) T = (0, i, 0, j, 0) T θ S 2 (i, j, k) T = (0, i, 0, j, 1, k, 1) T θ S 6 (i, j, k) T = (0, i, 0, j, 1, k, 8) T
If we perform a loop interchange transformation, we will get the following transformed scheduling functions:

θ ′S 1 (i, j) T = (0, j, 0, i, 0) T θ ′S 2 (i, j, k) T = (0, j, 0, i, 1, k, 1) T θ ′S 6 (i, j, k) T = (0, j, 0, i, 1, k, 8) T
By applying the sequential version of the violation check, we find that the violation sets are empty:

Vio S 1 →S 1 →S 2 = / 0 Vio S 1 →S 6 →S 2 = / 0 Vio S 6 →S 1 →S 2 = / 0 Vio S 6 →S 6 →S 2 = / 0
proving that there are no violations of live range intervals, after an interchange transformation. The conclusion is that the output dependences induced by the scalar variable phi_out_of_ssa.6[0] (shown as δ WAW S 1 →S 1 and δ WAW S 6 →S 6 in Figure 5.6), could be safely ignored when checking the legality of the transformation.

This check has to be performed for other memory accesses as well. In addition, we perform a scheduling violation analysis (Section 3.3) on true dataflow (read-after-write) dependences only: δ S 1 →S 2 and δ S 6 →S 2 , to check that relative writes/reads are executed in a correct order.

A classical dependence violation analysis [START_REF] Vasilache | Violated dependence analysis[END_REF], that takes into an account the full dependence graph, would consider this transformation to be illegal, thus precluding (obviously) legal interchange.

Performance benefits

The original motivation for coming up with the lazy expansion scheme was the fact that a threeaddress-code internal representation of the compiler that we were developing (Chapter 4) inherently contains many scalar variables.

It is a known fact that memory based write-after-write dependences on scalar variables hamper any effort on automatic parallelization. Indeed, all the codes that contain a reduction operation, like the matrix matrix multiplication shown in Figure 5.1, are internally transformed into the form where a scalar value is used to accumulate the values. Also, if an user writes such a code, such as the version of the matmult shown in Figure 5.2 the parallelization is not possible, since there is a loop carried output dependence [START_REF] Allen | Optimizing Compilers for Modern Architectures[END_REF]. Our scheme enables the transparent handling of all the non-harmful output dependences, also in the case of loop parallelization. We ran measurements on a suite of computationally intensive numerical kernels [START_REF]Polybench, the polyhedral benchmark suite[END_REF]. We compared the parallelizability of the kernels using both dependence models: the classical, all data-dependences analysis, and our approach of combined array data-flow analysis with live range analysis.

Classical

After integrating the new approach into the GRAPHITE polyhedral framework of the GCC compiler, we found that the new approach allowed for more loops -and mostly importantly, more outer loops -to get parallelized, as depicted in Table 5.1. This translates naturally into the speedups obtained by running those benchmarks on chip multiprocessor 5.11.

The baseline is the w/o dataflow column, standing for a sequential run when the current transformation violation analysis is used. The dataflow column represents a sequential run when the new, lazy live range violation analysis is used. As expected, changing the analysis doesn't influence the runtime of sequential runs, since we perform no sequential code transformation, which is why these columns show the same runtimes for all kernels.

Enabling parallelization on top of the current dependence analysis produces only a mild speedup on jacobi-2d-imper, while it does not impact or sometimes even slows down the others. This is mainly due to the low speedup potential of inner loop parallelization.

When autoparallelization is enabled with the new data-flow and live range analysis, we achieve substantial speedups of up to 3.8x. We experience one degradation in atax kernel, which we relate to the overhead incurred by autopar, unrelated to which data dependence model we use, since it happens for both approaches.

These are only the preliminary results, showing the potential of lazy memory expansion scheme combined with privatization and automatic parallelization. The integration of the GRAPHITE framework with automatic parallelization phase of the GCC compiler was discussed thoroughly in Chapter 4.

Summary

We have shown how a value live range analysis can be used to enable a lazy memory expansion scheme -memory is expanded only if absolutely necessary for the correctness of the transformation.

Polyhedral compilation traditionally takes as an input a dependence graph with all dependences, constructs a legal transformation and generates code. The other class of so called violation analysis based compilation flow [START_REF] Vasilache | Automatic correction of loop transformations[END_REF] takes as an input a dependence graph with all dependences as well, it constructs a transformation (without legality check), and only then it checks whether the transformation is legal. If it is, then it proceeds with code generation. If it is not, then it reiterates and proposes a next transformation until it finds a legal one.

We take the violation analysis approach a step further: we do not take into the account all the dependences. We split the dependence graph into those dependences that are true data-flow dependences and those that are memory based. We do not check memory-based dependences for scheduling violation violation.

By not taking the memory-based dependences into an account for scheduling violation analysis we cannot guarantee that original values written into specific memory location are not destroyed before they could be read by an appropriate instruction. In order to maintain this property, we first define live ranges for all values generated in an original program. Then we check whether a transformation would respect those value live ranges by performing the live range violation analysis. If we prove that all value live ranges are still respected, then the transformation is legal.

If a transformation destroys value live range set for a memory location, it is still possible to correct this: we could choose to replicate (memory expand) those memory locations. This leads to a controllable memory expansion scheme -we only expand when there is a need for expansion. A complete scheme of our approach is shown in Figure 5.12.

Conclusions and perspectives

We presented an algorithmic framework to maximize the effectiveness of loop transformations while maintaining per-hardware-thread memory usage invariant. Unlike traditional array expansion and contraction techniques, dependence removal can be implemented with guaranteed memory footprint constraints. This framework is a promising trade-off between array expansion and the the degrees of freedom for affine transformation.

In the context of a low-level three-address-code polyhedral compilation solving this problem is even more urgent, as evidenced in the experimental evaluation section 5.4. The general compilation flow that includes the lazy memory expansion scheme, presented in 5.5, fits well into the compilation flow of GRAPHITE framework (Chapter 4) and a transformation search method shown in Chapter 7.

Part III

Towards a search strategy

Chapter 6

Analytical cost model

The critical aspect of modern architectures is the need to efficiently coordinate the exploitation of multiple forms of parallelism provided by platforms, while carefully coordinating the use of the memory hierarchy efficiently.

Systematic solutions to harness the interplay of multi-level parallelism and locality are emerging, by advances in automatic parallelization and loop nest optimization [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF][START_REF] Pouchet | Iterative optimization in the polyhedral model: part ii, multidimensional time[END_REF]. These rely on the polyhedral model of compilation to facilitate efficient exploration and application of very complex transformation sequences. However, the analytical cost models employed by those solutions are not able to model the low-level, machine-dependent aspects of the target architectures.

One such an aspect is exploiting short SIMD parallelism by vectorization [START_REF] Nuzman | Outer-loop vectorization -revisited for short SIMD architectures[END_REF]. Naturally, such a kind of vectorization has not been represented within the polyhedral model due to its low-level, machinedependent nature. As a result, there remains a gap in providing a combined framework for exploring complex loop transformation sequences together with vectorization.

The problem of providing the analytical cost-models for predicting the performance is deemed difficult [START_REF] Pouchet | Iterative optimization in the polyhedral model: Part i, one-dimensional time[END_REF], since it seems infeasible to provide that static analytical model that could predict the performance of the architectures that are ever-growing in their complexity.

An alternative solution, based on an iterative compilation, has demonstrated the superior performance over the static compilation [START_REF] Agakov | Using machine learning to focus iterative optimization[END_REF][START_REF] Almagor | Finding effective compilation sequences[END_REF]. Nevertheless, the major drawback of the iterative compilation schemes is their running time, precluding them from being incorporated into the general-purpose, single-pass compilers such as Intel ICC, IBM XL or GCC [START_REF]Gnu compiler collection (gcc[END_REF].

Our goal is to provide a hybrid approach: we provide a precise and parameterizable analytical performance cost-predictor and we use it within a transformation search loop within the polyhedral model framework.

We present the cost-model function that is sensitive to the the actual loop transformation expressed as a schedule (Chapter 3). By exploring a finite and restricted space of schedules (as explained in Section 6.3.1), the cost-model function is evaluated for different loop transformations. In this chapter we present a simple scheme for the search space construction, while the much more expressive search space is presented in Chapter 7.

This search based cost-model function evaluation scheme still incurs some overhead in terms of the compilation time, since the search space has to be constructed and traversed. But as we show, this space is reduced to a reasonable size, and an evaluation of the cost-model function for each point in the search space is effective.

We achieve this effectiveness by enumeratively evaluating the non-linear cost function for different loop transformations, without any code generation nor running the result code, contrary to a feedback directed iterative search [START_REF] Pouchet | Iterative optimization in the polyhedral model: Part i, one-dimensional time[END_REF]. Only after the best schedule is determined, the code is generated only once and the final, optimized executable is obtained. This scheme is illustrated in Figure 6.3.

The most recent work of Park et al. [START_REF] Park | Predictive Modeling in a Polyhedral Optimization Space[END_REF] is aiming at solving the similar problem: constructing a predictive cost-model that drives the selection of the compiler transformations without resorting to the iterative search during the compilation time. But the fundamental difference, compared to our approach, is that they use the machine learning technique to (still, iteratively) learn the optimization sequences for the selected suite of benchmarks, based on the premise that the learned sequences are good predictors for the rest of programs. Our model is not based on learning, but analyzing the behavior of the program in order to predict its performance.

We are not trying to build a cost-model function that models all the possible aspects of the target architecture. We focus ourselves on the aspect of SIMD vectorization, which so far has not been considered in the polyhedral model. Yet this is a low-level and highly target specific optimization that has to be modelled precisely in order to yield the best results. We will give an overview of the SIMD vectorization in Subsection 6.1.1.

In a summary: in this chapter we present a precise analytical, performance predicting cost-model that captures the low-level details of the target architecture in the polyhedral model. In particular, we help bridge the aforementioned gap between analytical and feedback directed methods by incorporating the low-level (vectorization) considerations into a polyhedral model.

A note to the reader: the term cost-model is used to refer to the analytical function used to predict the total execution cost of a given program with applied transformation. The term polyhedral model is used to refer to the complete program analysis framework, as explained in Part I. The actual cost-model function is implemented within the polyhedral model framework.

Motivation

Fine-grain data level parallelism is one of the most effective ways to achieve scalable performance of numerical computations. This effectiveness translates into area and power savings for processor architectures, and complexity and performance benefits for applications. These are well known advantages of SIMD over MIMD architectures and programming models, and are best illustrated by the computational density of graphical processing units (GPUs). In this chapter we concentrate on improving the exploitation of short-SIMD parallelism available in modern instruction sets and vector-processing extensions (including Altivec, Cell SPU, and SSE). We will hereafter refer to such fine-grain data parallelism as subword parallelism.

SIMD Vectorization

Automatic vectorization for modern short-SIMD instruction sets, such as Altivec, Cell SPU and SSE, has been a popular topic, with successful impact on production compilers [START_REF] Wu | An integrated Simdization framework using virtual vectors[END_REF][START_REF] Bik | The Software Vectorization Handbook. Applying Multimedia Extensions for Maximum Performance[END_REF][START_REF] Bik | Automatic intra-register vectorization for the Intel architecture[END_REF][START_REF] Nuzman | Autovectorization in GCC -two years later[END_REF]. Exploiting subword parallelism in modern SIMD architectures, unlike automatic vectorization for traditional vector computers [START_REF] Allen | Optimizing Compilers for Modern Architectures[END_REF], however suffers from several limitations and overheads [START_REF] Nuzman | Outer-loop vectorization -revisited for short SIMD architectures[END_REF][START_REF] Nuzman | Auto-vectorization of interleaved data for simd[END_REF][START_REF] Naishlos | Autovectorization in gcc[END_REF] (involving alignment, redundant loads and stores, support for reductions and more) which complicate the optimization dramatically.

Automatic vectorization was also extended to handle more sophisticated control-flow restructuring including if-conversion [START_REF] Shin | Superword-level parallelism in the presence of control flow[END_REF] and outer-loop vectorization [START_REF] Nuzman | Outer-loop vectorization -revisited for short SIMD architectures[END_REF]. Classical techniques of loop distribution and loop interchange [START_REF] Allen | Optimizing Compilers for Modern Architectures[END_REF] can dramatically impact the profitability of vectorization. To be successful, it is vital to avoid inapt strategies that incur severe overheads, which are quite common in complex memory access patterns.

Motivating Example

Loop vectorization can adversely affect performance, mainly due to data management overheads. The first and foremost goal of a vectorization cost-model is to avoid performance degradations while not missing out on improvement opportunities. In addition, a cost model should also drive the selection of a vectorization strategy, assuming there exist a profitable one. The problem is particularly acute on modern subword SIMD architectures and compilers due to the interplay among compilation passes, the interactions among micro-architecture components, and the dependence of run-time information (array alignment, loop trip count).

Given a loop-nest, a compiler needs to choose which loop to vectorize, and at which position, employing one of several strategies (innermost-or outer-loop vectorization, in-place or based on innermosting, as explained below). This in-turn brings to play other loop-transformations, most notably loop-interchange but also loop-peeling and others. Searching among all these alternatives becomes a non-trivial problem. This problem is especially apparent in computations featuring loop nests that are (1) deep, (2) can be vectorized in several ways, (3) are amenable to other loop optimizations, and (4) are sensitive to underlying architectural features (e.g. alignment handling mechanisms, sensitivity to data locality). Figure 6.1 introduces the Convolve kernel -a simple example of a loop nest exhibiting the above features. Convolve performs a 2D-convolution on a NxN 16-bit pixel image. Such convolutions are vastly used in DSP programs for various applications, from edge-detection using 3x3 filters (as in the UTDSP suite [START_REF] Lee | Utdsp benchmarks[END_REF]), to 17x17 polyphase filters. 1There are d possible vectorization alternatives for a loop-nest of depth d (in our case d = 4), without involving any other loop-transformation: we can vectorize any of the j, i, h or v loops "in-place" -i.e. in their original position (loop-level). For instance, we can vectorize the j-loop in its innermost position (as shown in Figure 6.2a), which is the common practice of vectorizing compilers. Employing loopinterchange to permute one loop (inwards or outwards) into a specific position within the loop nest and vectorizing it there (keeping the other loops intact), increases the search-space to d • d possibilities. Figures 6.2b and 6.2c show examples of vectorizing the h-loop after permuting it to the innermost and next-to-innermost positions, respectively. Using loop-permutation more aggressively to reorder the loops of a nest according to a specific permutation, and then vectorizing one of them, results in a total of d(d!) combinations. If we also employ loop peeling to align memory accesses, the search space grows to d(d!)V F where V F is the Vectorization Factor (number of elements operated upon in parallel in a vector). In our case this amounts to 768 alternatives.

The search space becomes quite large even for modest depths and only few transformations (interchange and peeling), as shown, and can easily reach much higher volumes if deeper loop nests and/or more loop-transformations are considered. Also note that programs often contain many loop-nests, where

f o r ( v = 0 ; v < N ; v ++) f o r ( h = 0 ; h < N ; h ++) { s = 0 ; f o r ( i = 0 ; i < K ; i ++) { v s [ 0 : 7 ] = { 0 , 0 , . . . , 0 } ; f o r ( v j = 0 ; v j < K ; v j +=8) { v s [ 0 : 7 ] += image [ v+ i ] [ h+ v j : h+ v j + 7 ] * f i l t e r [ i ] [ v j : v j + 7 ] ; } s += sum ( v s [ 0 : 7 ] ) ; } o u t [ v ] [ h ] = s >> f a c t o r ; } (a) j-loop vectorized at level 4 f o r ( v = 0 ; v < N ; v ++) f o r ( h = 0 ; h < N ; h ++) o u t [ v ] [ h ] = 0 ; f o r ( v = 0 ; v < N ; v ++) f o r ( i = 0 ; i < K ; i ++) { f o r ( j = 0 ; j < K ; j ++) { c = f i l t e r [ i ] [ j ] ; v f i l t e r [ 0 : 7 ] = { c , c , . . . , c } ; f o r ( vh = 0 ; vh < N ; vh +=8) { o u t [ v ] [ vh : vh + 7 ] += v f i l t e r [ 0 : 7 ] * image [ v+ i ] [ vh+ j : vh +7+ j ] ; } } } f o r ( v = 0 ; v < N ; v ++) f o r ( h = 0 ; h < N ; h ++) o u t [ v ] [ h ] = o u t [ v ] [ h ] >> f a c t o r ; (b) h-loop vectorized at level 4 f o r ( v = 0 ; v < N ; v ++) f o r ( h = 0 ; h < N ; h ++) o u t [ v ] [ h ] = 0 ; f o r ( v = 0 ; v < N ; v ++) f o r ( i = 0 ; i < K ; i ++) { f o
] ; } } f o r ( v = 0 ; v < N ; v ++) f o r ( h = 0 ; h < N ; h ++) o u t [ v ] [ h ] = o u t [ v ] [ h ] >> f a c t o r ;
(c) h-loop vectorized at level 3 Approaches that generate each alternative and rely on its (possibly simulated) execution or on performance evaluation at a later low-level compilation stage, are competitive in terms of accuracy but are significantly inferior in terms of scalability to analytical approaches that reason about costs and benefits without actually carrying out the different loop transformations beforehand. Operating on the polyhedral representation itself, rather than relying on code generation, is therefore a key ingredient. Hybrid approaches can provide a more practical solution by combining the feedback-based approach with classical analytical models to narrow the search space. The modest compile-time requirements of our purely analytical approach (about 0.01s to build the model and search for the optimal vectorization strategy for Convolve) facilitates its integration in a production compiler.

Having potentially very large search-spaces is only one aspect of the vectorization decision-making problem. A complementary challenge to dealing with the very large search-spaces, is how to evaluate the costs and benefits associated with each alternative efficiently and accurately. Some trade-offs are clearly visible in Figure 6.2. For example, variants (b,c) use loop-permutation, which in this case incurs an overhead of extra memory traffic to/from the out array. On the other hand variant (a) incurs a reduction epilogue overhead (see sum operation) in each iteration of the i-loop. Outer-loop vectorization (vectorizing a loop other that innermost-loop) is used in (c), implying that more code is vectorized. The innermost j-loop in this case continues to advance sequentially, operating simultaneously on values from V F = 8 consecutive h-loop iterations. On the other hand (b) has better temporal locality ( f ilter[i][ j] is invariant in the innermost loop) and the misalignment is fixed (this is explained in more detail later). Overall the speedup factors obtained by transformations a, b, c on PPC970 (relative to the original sequential version shown in Figure 6.1) are 2.99, 3.94, 3.08 respectively. On the Cell SPU the respective speedups are 2.59, 1.44, 3.62.

The following sections describe our approach and demonstrate how our cost model computes its predictions within the analytical polyhedral-based model, considering different loop transformations and metrics. Final cost-model predictions for Convolve and analysis of the speedups are given in Section 6.3.3, where we show that the cost model is able to correctly predict the best vectorization option for both PPC and SPU.

Polyhedral Modelling of Vectorization Metrics

Several key costs impact the expected performance of vectorized code, including: strides of accesses to memory, memory access alignment, loop trip counts, reduction operations across loop iterations and more. These factors depend on the modified scheduling θ S ′ and on the modified iteration domain D ′S of each statement.

The underlying assumption of vectorization is that the kernel of a loop usually executes faster if vectorized than if not, but that associated overheads may hinder the vectorized version, diminishing its speedup compared to the original scalar version, and more so for loops that iterate a small number of times. Indeed, if the number N of iterations of a loop is smaller than its vectorization factor V F, there is no potential for speeding it up using vectorization; on the contrary, vectorizing such a loop may only slow down its execution due to additional preparatory actions and checks. Furthermore, even if N is larger than V F, the number of iterations of the vectorized loop, although positive, may not suffice to out-weigh the overheads incurred by vectorization.

Modelling the Access Patterns

Recall from Section 2.3.3 that in the classical polyhedral framework memory access functions for array references are represented as a vector of affine expressions:

f (i) = F × (i, g, 1)
T but there is no notion of the data layout of an array.

One may combine this access function with the data layout of the array. For each array reference, one may form a linearized memory access function ℓ, capturing the stream of memory access addresses as a function of the iteration vector:

ℓ(i) = b + (L i |L g |ω) × (i, g, 1) = b + L i i + L g g + ω (6.1)
where b is the base address2 of the array and (L i |L g |ω) is the row vector of coefficients that encodes the layout information (assuming row-major data layout). This vector is composed of three parts: L i is scheduling-dependent, L g depends on global parameters, and ω is the constant offset part.

Assuming that matrix M is defined as

M[r 1 ][r 2 ] . . . [r m ],
we can construct the vector R encoding the strides along each subscript:

R = m-1 ∏ i=1 r i , m-1 ∏ i=2 r i , . . . , m-1 ∏ i=m-1 r i , r m-1 , 1
Then the following equation holds:

(L i |L g |ω) = R × F (6.2)
where the matrix F defines the access function f . For example, taking the kernel in Figure 6.1 and assuming array image is defined as image [START_REF] Shin | Compiler-controlled caching in superword register files for multimedia extension architectures[END_REF][144], the linearized access for the array image in the statement S 2 of Figure 6.1 can be represented as:

(144, 1) × 1 0 1 0 0 0 0 0 1 0 1 0 0 0 = (144, 1, 144, 1, 0, 0, 0) meaning that linearized access function is:

ℓ(v, h, i, j) = b + 144v + h + 144i + j.
The linearized access function is crucial for computing the cost of vectorized data load/store instructions, which constitutes the majority of the vectorizer overhead.

Access Pattern Sensitivity to Scheduling

Based on the canonical scheduling matrix representation defined in Section 3.1 the rescheduled timestamp vector t is expressed as follows (for modelling purposes we can ignore β):

t = (A|Γ) × (i, g) = Ai + Γg (6.3)
thus the original iteration vector is i = A -1 (t -Γg) which together with Equation (6.1) gives us the new, transformed linearized access function:

ℓ ′ (t) = b + L i A -1 t + (L g -L i A -1 Γ)g + ω. (6.4)
with a new vector of coefficients:

L ′ = (L i A -1 |L g -L i A -1 Γ|ω).
Taking as example the kernel in Figure 6.1, the linearized access functions for arrays image and filter are as follows:

ℓ image (i) = b + [ 144 1 144 1 0 0 ω ] × (i, g, 1) T ℓ filter (i) = b + [ 0 0 144 1 0 0 ω ] × (i, g, 1) T
After applying loop interchange, by swapping the columns 3 and 4, we obtain the following transformation matrix (only A part is shown):

A ′S 2 =     1 0 0 0 0 1 0 0 0 0 0 1 0 0 1 0    
the new access functions become as follows:

ℓ image (t = (v, h, j, i)) = b + 144v + h + i + 144 j ℓ ′ image (t) = b + ( 144 1 1 144 0 0 ω ) × (i, g, 1) T ℓ ′ filter (t) = b + ( 0 0 1 144 0 0 ω ) × (i, g, 1)
T Notice that the memory access strides with respect to the new scheduling dimensions have changed. This has a dramatic impact on the performance of the vectorized code. Indeed, if we chose to vectorize the innermost level, the vectorized code will suffer from a very costly memory access operations: the stride is 144, so the elements are not consecutive and a vector cannot be loaded within one vector-load instruction.

This shows that linearized access functions, on which the total vectorization cost depends, is transformed automatically with the scheduling transformations. Thus, we do not need to generate code in order to compute the vectorization cost after applying a set of loop transformations -the vectorization cost is the direct function of scheduling matrix. For the rest of presentation we focus on L i part of the linearized access function coefficient vector.

Cost Model Function

Our cost model is based on modelling the total execution time of all statement instances, given the modified iteration domain D ′S and the modified schedule θ S ′ of each statement S. We compute the cost function for statement S as follows:

c(D ′S , θ ′S ) = |D ′S | VF ( ∑ c vect_instr ) + ∑ m∈(W S ) (c a + |D ′S | VF (c vect_store + f m )) + ∑ m∈(R S ) (c a + |D ′S | VF (c vect_load + c s + f m ))
where |D ′S | denotes the integer cardinality of the iteration space (total number of dynamic instances of S) and VF is the vectorization factor. We currently support the simplest case, when loop bounds form a rectangular polyhedral iteration space:

|D S | = ∏ dim(S) i=1 (UB i -LB i ).
In a more general case, algorithms to computing the number of integer points inside polyhedra could be used [START_REF] Barvinok | A polynomial time algorithm for counting integral points in polyhedra when the dimension is fixed[END_REF][START_REF] Clauss | Counting solutions to linear and nonlinear constraints through Ehrhart polynomials: applications to analyze and transform scientific programs[END_REF][START_REF] Pugh | Counting solutions to presburger formulas: how and why[END_REF].

Given the desired scheduling dimension d, we can compute a finite difference ∆ d of the linear memory address with respect to the time dimension:

∆ d = ℓ(i 1 , . . . , i d + 1, . . . , i d S ) -ℓ(i 1 , . . . , i d , . . . , i d S ) = L i d
∆ d is a memory access stride w.r.t. a schedule dimension; we will simply call it the stride. It can be determined directly from linearized access vector L i by looking at its d-th component.

Factor c s considers the penalty of load 3 instructions accessing memory addresses with a stride across the loop being vectorized. Accesses to non-unit strided addresses require additional data unpack or pack operations, following and/or preceding vector load or store instructions, respectively [START_REF] Nuzman | Auto-vectorization of interleaved data for simd[END_REF].

For example, VF scalar accesses to memory addresses with stride ∆ d v across the loop being vectorized may require ∆ d v vector loads (each with cost c 1 ), followed by ∆ d v -1 vector extract odd or extract even instructions (each with cost c 2 ), to produce one vector holding the desired VF elements. On the other hand, if several accesses to the same address are vectorized together (i.e. ∆ d v = 0), a vector "splat" instruction is often required to propagate the loaded value across all elements of a vector (with cost c 0 ). Factor f s is computed as a function of the stride ∆ d v :

c s =    ∆ d v = 0 : c 0 ∆ d v = 1 : 0 ∆ d v > 1 : ∆ d v • c 1 + (∆ d v -1) • c 2    (6.5)
Factor f a considers the alignment of loads and stores. Typically, accesses to memory addresses that are aligned on VF-element-boundaries are supported very efficiently, whereas other accesses may require loading two aligned vectors from which the desired unaligned VF elements are extracted (for loading) or inserted (for storing). This alignment overhead may be reduced considerably if the stride ∆ of memory addresses accessed across loop levels d v + 1..d S is a multiple of VF, since the misalignment remains constant inside the vectorized loop. If this is the case there is the opportunity to reuse loaded vectors and use invariant extraction masks. By having the transformed linearized access function:

ℓ(i) T = b + L i 1 i 1 + . . . + L i d v i d v + • • • + L i dim(S)
i dim(S) + L g g + ω it is easy to check if misalignment inside the vectorized loop remains constant -the coefficients from L i d v +1 to L i dim(S) (corresponding to strides of all inner loops of the vectorized loop) have to be a multiple of VF.

If the misalignment is constant inside the vectorized loop we also check if the base address which is accessed on each first iteration of the vectorized loop (d v ) is known to be aligned on VF-elementboundary; if so then there is no need for re-aligning any data: f a = 0. This is done by considering strides across outer-loops (enclosing the vectorized loop, if exist), and initial alignment properties such as array alignment. In order to check the alignment in outer loops, we need to check if coefficients from L i 1 to L i d v -1 are multiple of VF.

By putting together all considerations for alignment, cost can be modelled as:

f a =        aligned : 0 var. misalign. : |D S |(c 1 + c 3 + c 4 ) fixed misalign. : |D S 1..d v -1 |(c 1 + c 3 )+ |D S |(c 1 + c 4 )        (6.6)
3. Vector store operations with strided access is not yet implemented in GCC.

where c 3 represents the cost of building a mask based on the misalignment amount, c 4 is representing the cost of extraction or insertion and c 1 is the vector load cost. |D S 1..d v -1 | denotes the number of iterations around the vectorizer loop level.

The vectorization factor VF of a loop is determined according to the size of the underlying vector registers and the smallest data-type size operated on inside the loop. Each individual vector register will thus hold VF values of this small size. However, if there are variables in the loop of larger size, storing VF copies of them will require multiple vector registers, which in turn implies that the associated instructions need to be replicated. Factor f m records the extra overhead that is associated with this replication. Additional factors that depend on the specific machine resources available may also impact the performance of vectorization, such as the size of register files, available ILP, and complex vector instructions.

By applying different loop interchange transformations and choosing different loops to vectorize, the performance of the resulting vectorization varies considerably. Our model was able to predict the best possible combination of loop interchange and outer/inner vectorization strategy.

Vectorization Profitability Metrics

To summarize, the above vectorization profitability metrics can be classified into the following three classes:

Scheduling invariant metrics are not affected by changing the execution order of the statement instances. Thus, they are invariant with respect to a scheduling function θ S . Vector to scalar reduction cost and multiple type support costs fall into this category.

Scheduling sensitive metrics are affected by changing the execution order of the statement instances. Those metrics have the greatest impact on the resulting cost, since they vary with the change of the scheduling order and thus are affected by loop transformations such as interchange. Changing of the scheduling function θ S directly affects the cost of strided memory accesses, non-aligned vector loads and spatial locality.

Code generation dependent metrics depend on the actual code-generation strategy implemented in a compiler 4 . These metrics go against the design of our cost model, as they need to regenerate the compiler's internal structures from the modified polyhedral representation. Yet no significant performance factor falls into this category. There is a practical difficulty however: idiom recognition (e.g., saturated arithmetic, absolute differences) are typically a source for scheduling-sensitive metrics, as the identification of algorithmic idioms can be performed on the array data flow [START_REF] Redon | Scheduling reductions[END_REF]. Yet vectorizers currently rely on syntactically fragile pattern-matching techniques and do depend on the code generation.

Evaluation

We present a systematic method to evaluate the cost function, within a selected transformation search space. Later, we will show the experimental evaluation of the predicted speedups, compared to the actual execution on two different architectures for a selected benchmark, altogether with a detailed discussion.

Driving the search process

To select an optimal vectorization strategy one needs to construct and traverse the relevant search space. We want to do so without generating different syntactic versions of the code and then vectorizing each of them, which is inefficient and sometimes infeasible.

φ(x) = n ∑ i=1 c(D S i , θ S i ). (6.7)
The parameters for the cost function for the single Statement S i are its iteration domain D ′S i (number of dynamic instances of a Statement depends on its iteration domain) and its scheduling θ ′S i (cost of accessing memory by a Statement instance depends on execution order of other instances). Section 6.2.3 describes this cost function in detail. The optimization goal is to search for a vector of transformations x min that minimizes the cost function φ(x):

x min = min x∈X φ(x) (6.8)
Vector x min represents an optimal program version in the polyhedral model. After extracting the SCoPs and building the polyhedral representation of all statements by the GRAPHITE framework, we perform the optimization of each SCoP according to Algorithm 1: first we compute the base cost for the unmodified (input program) representation, by computing the cost of executing all dynamic instances of all statements S i in the original scheduling order. The current optimal cost is stored in cost min and is updated incrementally by applying different transformations (skipping over infeasible ones) on the polyhedral model (stored in vector x) and computing the new costs using cost function φ(x). Besides the schedule transformation, performed by permuting (PERMUTE) the columns of the compo-nent A of the schedule, each possible level v is strip-mined, which is the way to model the vectorization at level v. At the end of this process the optimal scheduling representation is available in x min .

Note that Algorithm 1 shows only one possible way of constructing a search space. We chose to consider all combinations of loop interchanges due to their impact on vectorization. This small (yet expressive) search space makes it compatible with the constraints of a production compiler. Much more expressive search space will be explained in Chapter 7.

Algorithm 1 Driver for search space exploration d ← level of a deepest Statement S in a SCoP n ← number of Statements in a SCoP

x min ← [θ S 1 , . . . , θ S n , D S 1 , . . . , D S n ]
⊲ Start with the original schedules and iteration domains cost min ← φ(x min ) for all σ ∈ (set of d-element permutations) do for i = 1 to n do We evaluate our approach by introducing our model into the polyhedral framework of GCC production compiler and comparing its performance estimates for different loop interchanges and vectorization alternatives against actual execution runs of a set of benchmarks. Table 6.1 summarizes the main relevant features of the kernels used in our experiments: a rate 2 interpolation (interp), block finite impulse response filter (bkfir), an 8 × 8 discrete cosine transform (dct [START_REF] Lee | Utdsp benchmarks[END_REF]), 2D-convolution (convolve), a kernel from H.264 (H264), video image dissolve (dissolve), weight-update for neural-nets training (alvinn) and a 16 × 16 matrix-matrix multiply (MMM) (including a transposed version MMM T ). The first four columns of Table 6.1 show the number of iterations N i of loops nested within the main loop-nest of each benchmark, starting with N 1 for the outermost loop and moving inwards. Loop nests with less that 4 nested loops have empty entries (e.g., N 2 refers to the innermost loop in the doubly-nested bkfir). For example, convolve has a 4-nest loop whereas bkfir has only a doubly-nested loop and thus N 2 refers to its inner-most loop. Multiple values in an entry represent multiple distinct loop nests.

θ ′S i ← PERMUTE(σ, θ S i ) d S i ← level of loop nesting for Statement S i for v = 1 to d S i do D ′S i ← STRIPMINE(v, D S i ) x ← [θ ′S 1 , . . . , θ ′S n , D ′S 1 , . . . , D ′S n ] if φ(x) < cost
Similarly, the next four columns of Table 6.1 show the strides ∆ i of the memory references across each of the nested loops, with multiple values in an entry representing the strides of different memory references. For example, strides of 8, 512 and 16 are found in the innermost loops of dct, alvinn and MMM respectively, where columns of 2D arrays are scanned resulting in strides at the length of the rows. Lastly, zero strides imply that duplication of a single value across a vector is required.

We first evaluate the cost-model qualitatively, demonstrating that the scores it computes are consistent using one detailed example (subsection 6.3.3). The following subsection (subsection 6.3.4) evaluates the model relative to actual experiments on a set of kernels, analyzing the mispredictions and showing that overall the relative performance ordering of the variants is largely preserved.

Qualitative Evaluation

We use the convolve kernel qualitatively (see Figure 6.1, Section 6.1.2). We study only a small subset of the search space described in Chapter 7, restricting our attention to the d × d combinations of moving each loop inwards/outwards and vectorizing it there, plus the option to vectorize each of the d loops without any interchange. Note however that our technique opens up a much larger transformation space, as shown in the search methodology described in Chapter The results of running our model against the d × d = 16 combinations, estimating the performance of each combination for a Cell/SPU and a PowerPC system are shown in Table 6.2 and Table 6.3 respectively. The loops are numbered in the tables from 1 (outer-most) to 4 (inner-most). Entry (i, j) shows the estimated speedup over the sequential version, obtained by moving loop j to position i followed by vectorizing loop j at new position i. Thus entries along the diagonal refer to vectorization with no interchange. Entries (4,4), (4,2), (3,2) (in bold) correspond to the vectorization alternatives shown in Figures 6.2a, 6.2b, 6.2c respectively.

While the estimated speedups for these versions are a little too optimistic compared to actual measured speedups (see Section 6.1.2), the relative ordering of the speedups for both platforms is accurate 5and the cost model is able to identify the best choice among the three.

The convolve entry in Table 6.1 reveals the key factor for the performance degradations predicted for loops v, i (columns 1 and 3) -there are very large strides along these loops (∆ 1 , ∆ 3 = 128). The overhead involved in vectorizing these loops and strides is described in Section 6.2.3. The remaining candidate loops for vectorization are therefore loops 2 and 4 (h and j). The best speedup is predicted for entry (3,4) which corresponds to using outer-loop vectorization to vectorize the j-loop after moving it to level 3. The original i-nest is a perfect nest (there are no operations outside the innermost loop within that nest) and so there are no overheads incurred by this interchange (as opposed to interchanging an imperfect-nest like the h-loop, e.g. as in cases (4,2),(3,2)/Figures 6.2b,6.2c, which involve scalar expansion and loop-distribution costs). In addition, outer-loop vectorization avoids reduction-epilogue costs and also increases the portion of the code that is being vectorized compared to vectorizing the j-loop in its original innermost location. Note that this choice is different from the traditional approach: compilers usually either apply inner-most loop vectorization (entry (4, 4) in the tables) or apply innermosting (entries (4, * )).

Partial experimental evaluation of convolve confirms these predictions. In Figure 6.6 we show the obtained speedups relatively to the cost model estimations (denoted exp,model respectively) for PPC970 and Cell SPU for entries (3, 2), (4, 2), (3,4) and (4, 4) in the tables. When a given vectorization technique cannot be applied due to limitations of our current implementation of vectorization in the GCC compiler, the scalar performance is reported. This happens in some cases of strided accesses that are not yet fully supported (and that would certainly degrade performance).

We evaluate the relative speedup of four different vectorization alternatives: innermost-loop vectorization (inner), interchange followed by innermost-loop vectorization (innermosting), and in-place outer-loop vectorization, with and without optimized realignment using unrolling (outer and outer-opt).

The experiments were generated automatically using an enhanced version of GCC. Speedup are measured over the sequential version of the benchmark, compiled with the same optimization flags. Interchange, when used, was applied manually. Time is measured using the getrusage routine on pow-erpc970, and the decrementer utility on the SPU. Experiments were performed on the IBM PowerPC PPC970 processor with Altivec, and an SPU of the Cell Broadband Engine. Both architectures have 128 bit wide vector registers, and similar explicit alignment constraints.

The first set of kernels (interp, bkfir, dct and MMM) is expected to gain most from in-place outerloop vectorization with realignment optimization, as they consist of imperfect loop-nests (and therefore get penalized for interchange), and exhibit high data-reuse opportunities across the (vectorized) innerloop that can be exploited by the unrolling optimization. They also have inner-loop reductions (which are generally done more efficiently using outer-loop vectorization), and two of the benchmarks in this set (dct and MMM) also have large strides in the innermost loop (as the access is column-wise). Alvinn has a perfect nest and no reuse opportunities, and therefore in-place outer-loop vectorization should not gain over traditional interchange, but innermost loop vectorization should be avoided due to the large stride. The last group of benchmarks (MMM T , dissolve and H264) have consecutive access in the innermost loop, but strided access in the outer-loop, and so for these we expect inner-loop vectorization to be the best technique. This behavior can be clearly observed in the SPU speedups in Figure 6.4, where overall the exp and model graphs are largely consistent, with the preservation of the relative performance ordering of the variants. Exceptions are due to low-level target-specific factors that our model does not take into account. Most notable is the misprediction in the first set of benchmarks, where bkfir and dct are the only benchmarks for which outer-loop vectorization is inferior to innermost loop vectorization due to an SPU specific issue (unhinted branch).

Target-specific issues come to play also on the PPC970 (Figure 6.5). The most significant one appears in the fixed-point bkfir and interp where inner-loop vectorization employs a specialized Altivec instruction to compute a dot-product pattern. We have not yet incorporated idioms into the cost-model and so it does not anticipate this behavior. The model also does not try to estimate register pressure, and therefore does not predict the degradation in performance incurred by the unrolling optimization on interp due to register spilling (this problem does not occur for SPU having 128 vector registers, compared to the 32 Altivec registers of PowerPC970). Inaccurate modelling of spatial locality in our current implementation is the reason why the cost model misses the improved data cache locality when interchanging loops in alvinn (this problem does not occur on the Cell SPU as it doesn't have a cache). Lastly, in some cases interchange can be done with smarter scalar-expansion (hoisting), whereas the model estimates the associated overhead of a naive scheme. This sometimes pessimizes the predicted speedup of interchanged versions both on PPC and the SPU.

Related Work

Vectorization Cost-Model Related Work. Leading optimizing compilers recognize the importance of devising a cost model for vectorization, but have so far provided only partial solutions. Wu et al. conclude [START_REF] Wu | An integrated Simdization framework using virtual vectors[END_REF] regarding the XL compiler that "Many further issues need to be investigated before we can enjoy the performance benefit of simdization ... The more important features among them are ... the ability to decide when simdization is profitable. Equally important is a better understanding of the interaction between simdization and other optimizations in a compiler framework". Likewise, Bik stresses the importance of user hints in the ICC vectorizer's profitability estimation [START_REF] Bik | The Software Vectorization Handbook. Applying Multimedia Extensions for Maximum Performance[END_REF], to avoid vectorization slowdowns due to "the performance penalties of data rearrangement instructions, misaligned memory references, failure of store-to-load forwarding, or additional overhead of run-time optimizations to enable vectorization."; on the other hand opportunities may be missed due to overly conservative heuristics.

These state-of-the-art vectorizing compilers incorporate a cost model to decide whether vectorization is expected to be profitable. These models however typically apply to a single loop or basic-block, and do not consider alternatives combined with other transformations at the loop-nest level. This work is the first to incorporate a polyhedral model to consider the overall cost of different vectorization alternatives in a loop-nest, as well as the interplay with other loop transformations.

Loop-nest auto-vectorization in conjunction with loop-interchange has been addressed in prior art [START_REF] Allen | Automatic translation of fortran programs to vector form[END_REF][START_REF] Allen | Optimizing Compilers for Modern Architectures[END_REF][START_REF] Wolfe | High performance compilers for parallel computing[END_REF]. This however was typically in the context of traditional vector machines (such as Cray), and interchange was employed as a preprocessing enabling transformation. Overheads related to short-SIMD architectures (such as alignment and fine-grained reuse) were not considered.

Costs of specific aspects of short-SIMD vectorization were addressed in more recent works. Realignment and data-reuse were considered together with loop-unrolling [START_REF] Shin | Compiler-controlled caching in superword register files for multimedia extension architectures[END_REF], but in the context of straightline code vectorization, and not for the purpose of driving loop vectorization. A cost model for vectorization of strided-accesses was proposed in [START_REF] Nuzman | Auto-vectorization of interleaved data for simd[END_REF], but it does not consider other overheads or loop transformations.

Polyhedral-Model Related Work. Bondhugula et al. [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF] integrate inner-loop vectorization as a post-pass of their tiling heuristic, and leverage the interchangeability of inner loops to select one that is vectorizable. Their method does not take into consideration the respective vectorization overheads, nor does it model reductions. Nevertheless, their tiling hyperplane and fusion algorithm can serve as a complementary first pass for our technique, favoring the extraction of interchangeable inner loops.

Pouchet et al. demonstrate how one can systematically study the interplay of loop transformations with backend optimizations (including vectorization) and complex microarchitectures by constructing huge search spaces of unique, valid transformation sequences [START_REF] Pouchet | Iterative optimization in the polyhedral model: part ii, multidimensional time[END_REF]. These search spaces are tractable using carefully crafted heuristics that exploit the structure of affine schedules. An analytical performance model capable of characterizing the effect of such complex transformations (beyond loop interchange, and accommodating for large-scale locality effects) does not currently exist. Known analytical cache models for loop transformations are quite mature in some domains, loop tiling in particular [START_REF] Boulet | Pen)-ultimate tiling?[END_REF], yet remain sensitive to syntactic patterns and miss key semantical features such as loop fusion effects [START_REF] Cascaval | Compile-time based performance prediction[END_REF][START_REF] Fraguela | Probabilistic miss equations: Evaluating memory hierarchy performance[END_REF].

Conclusions and future work

This contribution narrows the gap between analytical cost-model based static compilation techniques [START_REF] Cascaval | Compile-time based performance prediction[END_REF] and fully iterative search based compilation flow [START_REF] Agakov | Using machine learning to focus iterative optimization[END_REF][START_REF] Pouchet | Iterative optimization in the polyhedral model: Part i, one-dimensional time[END_REF][START_REF] Pouchet | Iterative optimization in the polyhedral model: part ii, multidimensional time[END_REF], while still maintaining the running time of the compilation acceptable. We have focused on the single aspect of the machine specific cost modelling -that of automatic vectorization on SIMD architectures.

We have presented the first-of-a-kind approach to low-level machine-specific analytical cost model, used for prediction of the performance of vectorization within the polyhedral model framework. Our contributions of this work are the following:

-Polyhedral modelling of subword parallelism. We demonstrate how to leverage the polyhedral compilation framework naturally and efficiently to assess opportunities for subword parallelism in combination with complex loop transformation sequences. Our integrated approach relies on abstract matrix manipulations instead of code generation, thereby shortening the evaluation time of the cost model dramatically, compared to iterative optimization approaches [START_REF] Pouchet | Iterative optimization in the polyhedral model: Part i, one-dimensional time[END_REF] -Evaluation in a production compiler. Our model is fully automated and implemented based on GCC [START_REF]Gnu compiler collection (gcc[END_REF] 4.5 compiler. -Studying the interplay between loop transformations. We provide a thorough empirical investigation of the interplay between loop interchange with array expansion and loop nest vectorization of both inner and outer loops on modern short-SIMD architectures.

Future work

The experimental results are promising, and the running time of the search space exploration is acceptable for incorporation into the general purpose compiler like GCC. Nevertheless, the cost model is somewhat simplistic and does not incorporate the modelling of cache locality issues [START_REF] Fraguela | Probabilistic miss equations: Evaluating memory hierarchy performance[END_REF], interplay with thread-level parallelism nor the incorporation of memory traffic optimizing loop transformations [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF] like loop interchange. We see the future directions that could fully exploit the potential of this approach:

Interplay with other transformations. The method could easily be extended and employed in the future, to consider the effects of additional transformations within the polyhedral framework. The search space that is proposed in Chapter 7, which is in turn based on [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF], could be used to combine the search for thread-level parallelism, memory locality optimization and optimal vectorization strategy in one, combined manner.

Precise data locality model. Our cost model includes a very rough estimation of the data localitythat of computing the data access strides. While it can capture the spatial locality of the given memory access pattern, it cannot precisely model the temporal locality. The model could be extended with cache miss equations, an approach proposed by Fraguela [START_REF] Fraguela | Probabilistic miss equations: Evaluating memory hierarchy performance[END_REF]. The preliminary integration efforts are ongoing, at the time of writing this dissertation.

Machine learning of target specific factors. The analytical cost function shown in subsection 6.2.3 relies on the machine specific atomic instruction costs. While we have obtained those costs for the architectures that we were interested in, mainly through microbenchmarking and knowledge of the instruction latencies, obtaining those coefficients for the new architectures might be. In order to facilitate the porting of the cost-model to the new architectures, we would like to investigate the possibility of employing machine learning, like in [START_REF] Cavazos | Automatic performance model construction for the fast software exploration of new hardware designs[END_REF], to obtain the coefficients automatically for the new architectures. Please note that this does not change our approach: a learning phase is done only once, to obtain the coefficients. Later, the ready cost-model function is used in a single-pass compiler directly, without any extra penalty.

Chapter 7

Transformation search strategy

In the previous chapters we have discussed the basic components needed in the polyhedral model based optimization compiler such as static analysis, data dependence analysis, cost model, and code generation. But the core part of each polyhedral model based compilation flow is an algorithm for obtaining the actual transformation. The aim of this chapter is to provide our contribution to this key part of the compilation flow.

While numerous works on program transformations -and polyhedral transformations in particularhave been published, there is no general consensus on whether they are practical and efficient for use in a general purpose optimizing compilers. In this chapter we propose a transformation search methodology, instead of proposing a single specialized transformation algorithm.

The chapter starts with introducing the related work in Section 7.1. A necessary notation, based on the state of the art approaches, is introduced in Section 7.5. The core of the method is explained in Sections 7.6 through 7.10.

Related work

Optimizing for various architectural features requires complex program transformations to exploit thread-level parallelism, vectorization opportunities and memory hierarchy in a combined manner.

Pioneering works [START_REF] Allen | Automatic translation of fortran programs to vector form[END_REF][START_REF] Lamport | The parallel execution of do loops[END_REF][START_REF] Banerjee | Unimodular transformations of double loops[END_REF][START_REF] Banerjee | Loop Transformations for Restructuring Compilers[END_REF][START_REF] Wolf | A loop transformation theory and an algorithm to maximize parallelism[END_REF] on loop transformations look for the best transformation according to some ad-hoc criteria. Furthermore, they are restricted to a very limited subset of program loops.

Loop optimizations expressed in the polyhedral model subsume all previous works on loop transformations. Optimizations in the polyhedral model are expressed as scheduling problems. There are as many scheduling algorithms as there are different optimization criteria, but they all share the common property: the computed schedule has to provide a legal transformation. The state of the art polyhedral transformation frameworks generally fall in two categories:

-best effort, cost model based scheduling approaches -iterative search based approaches The works of Feautrier [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part I: one dimensional time[END_REF][START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part II: multidimensional time[END_REF], Lim and Lam [START_REF] Lim | Maximizing parallelism and minimizing synchronization with affine transforms[END_REF][START_REF] Lim | An affine partitioning algorithm to maximize parallelism and minimize communication[END_REF], Griebl [START_REF] Griebl | Automatic parallelization of loop programs for distributed memory architectures[END_REF] and Bondhugula [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF] fall into the first category, while the works of Pouchet [START_REF] Pouchet | Iterative optimization in the polyhedral model: Part i, one-dimensional time[END_REF][START_REF] Pouchet | Iterative optimization in the polyhedral model: part ii, multidimensional time[END_REF] and Vasilache [START_REF] Vasilache | Automatic correction of loop transformations[END_REF] fall into the second category.

The basic principle underlying the best effort scheduling approaches is a well defined cost model function. The problem is cast as a linear programming problem where the data dependences provide the set of constraints, and an objective is to find a solution that minimizes the given cost model function, while preserving the legality of the problem.

Obviously, different cost functions provide different solutions. An objective function in the approach of Feautrier is the latency. The goal of the scheduling is to provide minimum latency schedules. Lim and Lam provide an objective function that minimizes the order of required synchronizations and maximize the number of parallel loop dimensions.

The state of the art cost model based approach of Bondhugula [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF] provides a cost function that minimizes the dependence distance. It turns out that this simple cost function is very powerful, since it maximizes the parallelism in the outermost loop levels and minimizes the volume of memory communication.

While the best effort scheduling approaches provide only one solution -the best solution according to the cost function -the other class of the scheduling approaches is based on the iterative search.

Iterative search based scheduling approaches do not provide a single solution. Indeed, they provide a space of the legal and distinct schedules. Each distinct schedule corresponds to one program transformation. The transformations within a space are assessed according to their 'quality' -mostly often the speedup provided.

Iterative approach in the polyhedral model was pioneered by Pouchet [START_REF] Pouchet | Iterative optimization in the polyhedral model: Part i, one-dimensional time[END_REF][START_REF] Pouchet | Iterative optimization in the polyhedral model: part ii, multidimensional time[END_REF]. Pouchet shows the systematic way to generate the space of legal and distinct schedules expressed in the polyhedral model and the way to traverse those schedules. For each distinct schedule a code is generated and the transformed program is run on the target hardware. The schedule corresponding to program version with the best runtime is selected.

Since the search space of all legal schedules might be huge, Pouchet proposes to use machine learning techniques to speedup the search [START_REF] Pouchet | Iterative optimization in the polyhedral model: part ii, multidimensional time[END_REF].

Limitations of the current approaches

The best effort cost model based techniques are limited in their ability to find the best transformation adapted to a particular architecture.

The approach of Feautrier [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part I: one dimensional time[END_REF][START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part II: multidimensional time[END_REF] finds the minimal latency schedule. But this cost model is based on an imaginary machine model where synchronization costs and memory locality issues are ignored. The algorithm is optimal on an idealized PRAM machine, but it is unlikely to be optimal on any real-world multiprocessor machine. This fact was empirically proven by Bondhugula [START_REF] Bondhugula | Effective automatic parallelization and locality optimization using the polyhedral model[END_REF].

The approach of Lim and Lam [START_REF] Lim | Maximizing parallelism and minimizing synchronization with affine transforms[END_REF][START_REF] Lim | An affine partitioning algorithm to maximize parallelism and minimize communication[END_REF] has a similar limitation. While it takes into an account the cost of synchronization, it does not model the memory traffic costs. The optimization process is based on basic linear algebra rather than on linear programming. This has an implications on the quality of obtained solution, since the method can reduce the degree of synchronization and maximize the degree of parallelism, but it can miss an exact optimal solution.

The state of the art approach of Bondhugula [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF] has a simple, but very powerful cost model that matches the key performance factors relevant to modern architectures. It takes into an account both parallelism and memory traffic optimization.

While being well adapted to extracting parallelism and optimizing for memory locality through loop tiling, the approach Bondhugula [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF] does not take into account some lower level architectural features, such as SIMD vectorization and possible cache interferences.

As evidenced in [START_REF] Pouchet | Combined iterative and model-driven optimization in an automatic parallelization framework[END_REF], the mentioned problem of Bondhugula's approach could be alleviated if combined with iterative search strategy. The method shown in [START_REF] Pouchet | Combined iterative and model-driven optimization in an automatic parallelization framework[END_REF] does not really into any of the categories discussed so far, since it is a combination of both.

The iterative search methods try to overcome the lack of the precision of the cost model based approaches but with an extra cost -that of generating multiple versions of program transformations.

The iterative search method of Pouchet [START_REF] Pouchet | Iterative optimization in the polyhedral model: Part i, one-dimensional time[END_REF][START_REF] Pouchet | Iterative optimization in the polyhedral model: part ii, multidimensional time[END_REF] builds a search space of all legal and distinct program transformations expressed as one-dimensional [START_REF] Pouchet | Iterative optimization in the polyhedral model: Part i, one-dimensional time[END_REF] or multi-dimensional [START_REF] Pouchet | Iterative optimization in the polyhedral model: Part i, one-dimensional time[END_REF][START_REF] Pouchet | Iterative optimization in the polyhedral model: part ii, multidimensional time[END_REF] affine schedules.

The main problem of iterative search methods is the problem of bounding the search space of legal program versions. In their full generality, the space of affine schedules is unbounded. One has to bound the possible values of the scheduling coefficients to some restricted range in order to obtain the finite-size search space.

Pouchet [START_REF] Pouchet | Iterative optimization in the polyhedral model: Part i, one-dimensional time[END_REF] proposes some practical solutions to the problem of bounding the search spaces to make the iterative search methods practical. We might add a note that the cost model based approaches do not face the problem of bounding the coefficient values, since those values are computed directly from the analytical cost model.

The feedback directed iterative search methods based on empirical program evaluation are obviously not suitable for integration into the general purpose optimizing compiler because of their unpredictable running time -even for simple program kernels, a millions of legal transformations might be evaluated [START_REF] Pouchet | Iterative optimization in the polyhedral model: Part i, one-dimensional time[END_REF].

Since the search space explored by the iterative search method has to be bounded upfront, some possible solutions might be missed, even if an exhaustive search is performed.

Obviously, neither class of the methods is better than the other. While the best effort, cost model based approaches provide the desired solution in the single step, they might not precisely model the target architecture and they might provide suboptimal results. On the other hand, the feedback directed, iterative search based methods might provide the optimal transformations -missed by cost model based approaches -but at the cost of an exhaustive search of the transformation space.

The common limitation of all the methods discussed so far is the lack of control over the complexity of the generated code. To the best of our knowledge, there are no works addressing the problem of finding the good program schedule that also generates the least complex code, though Darte [START_REF] Darte | Loop shifting for loop compaction[END_REF] has pointed out this problem. This problem is very relevant for enabling the integration of the scheduling algorithm in a general purpose compiler. We will elaborate on this problem Section 7.3.

All the mentioned methods rely on the representation of the data dependence graph and satisfaction of all the data dependences. As we have shown in Chapter 5, some of the data dependence constraints could be relaxed and some dependence constraints could be completely removed with an extra cost of memory expansions. Removing the dependence constraints enables more optimization choices.

While the techniques for removing data dependence constrains have been extensively studied [START_REF] Feautrier | Array expansion[END_REF][START_REF] Tu | Array privatization for shared and distributed memory machines (extended abstract)[END_REF][START_REF] Calland | On the removal of anti-and output-dependences[END_REF][START_REF] Barthou | Maximal static expansion[END_REF][START_REF] Knobe | Array ssa form and its use in parallelization[END_REF], they were treated separately from the scheduling problem itself. The only work that addresses those two problems in a combined manner is [START_REF] Cohen | Parallelization via constrained storage mapping optimization[END_REF]. The related problems were discussed extensively in Chapter 5. We summarize the limitations of the current best effort, cost model based approaches:

-Provide only a rough, linear analytical model based on an imaginary target execution platform -Focus on a particular performance issue, while not taking into an account the holistic effect of multiple program transformations -The optimization process takes a greedy optimization approach that can stick to local minimum and miss a global optimal solution -No control over the complexity of the generated code -No way to integrate the relaxation of the scheduling constraints through dependence removal The limitations of the iterative search based approaches are summarized:

-Exhaustively searching the space of all transformations is expensive in terms of time -The complete search space is unbounded and a bounding solution has to be provided -Even if an exhaustive search is performed, the optimal solution might be missed -Does not help in understanding the program transformations -No control over the complexity of the generated code

Motivation

Having on mind the limitations of the current scheduling approaches, our goal is to propose the transformation search methodology that could provide the best trade-off between cost model based scheduling methods and search based techniques.

The core idea is to combine the enumeration search method with an evaluation of the specialized analytical cost model function. In this way, we combine the best of both worlds: the precision of iterative enumeration and effectiveness of the cost model based approaches.

The specialized cost function is tailored to the specific architecture and models the selected performance phenomena. Such a cost model function was shown in Chapter 6.

Previous cost model based approaches could not handle the complex non-linear objective functions. The reason is simple: they use (integer) linear programming techniques that only can operate on linear functions.

We cannot use the techniques proposed by Feautrier and Bondhugula to evaluate and choose the best transformation according to the non-linear cost function. Instead, we resort to enumerative evaluation of the non-linear cost function, in a similar way as the iterative approach.

Feedback directed iterative optimization [START_REF] Pouchet | Iterative optimization in the polyhedral model: Part i, one-dimensional time[END_REF][START_REF] Pouchet | Iterative optimization in the polyhedral model: part ii, multidimensional time[END_REF] relies on the feedback, for each generated program version, obtained from measuring the performance of the transformed program. But instead of relying on the empirical feedback, we plug in the non-linear objective cost function to provide the assessment of each transformation version.

Provided that we have a precise cost function tailored to a target architecture and that we can construct a bounded but expressive search space, we can perform an enumerative search of the cost function evaluations and pick the best point in the search space.

The success of our approach relies on providing the following key properties:

1. construction of an expressive search space of affine transformations 2. ability to control the size of the search space 3. ability to efficiently evaluate the cost function for each point in the search space Based on our new search methodology, we want to address the problems that were not addressed in an integrated manner in the previous approaches. Mostly importantly, we aim at providing the transformation engine that could be integrated as an optimization pass of a general purpose compiler.

Two properties of our method enable its easy integration into the general purpose compiler: 1) embedding of the precise cost model tailored to a given architecture 2) not relying on the empirical feedback as the current iterative optimization approaches [START_REF] Pouchet | Iterative optimization in the polyhedral model: Part i, one-dimensional time[END_REF][START_REF] Pouchet | Iterative optimization in the polyhedral model: part ii, multidimensional time[END_REF] An enumeration based nature of our approach enables us to control the complexity of the generated code by ignoring some transformations (loop skewing for example) that might generate the transformed code that is too costly [START_REF] Darte | Loop shifting for loop compaction[END_REF].

Also, an enumerative approach for the schedule construction goes in tandem with a lazy memory expansion scheme discussed in Chapter 5. There is a trade-off between the cost of the memory expansion and the expressiveness of the transformation. To the best of our knowledge, there is no work that shows how to capture this trade-off in a purely best effort based methods [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part II: multidimensional time[END_REF][START_REF] Lim | An affine partitioning algorithm to maximize parallelism and minimize communication[END_REF][START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF]. An enumerative based approach can capture this trade-off in its search space.

Problem statement

The goal of the transformation search strategy is the enumeration of the search space of a class of legal multidimensional schedules. The space of the legal solutions L is obtained. Each legal solution T ∈ L have the form:

T = {θ S 1 , . . . , θ S p }
Each legal solution describes the complete scheduling matrices for all statements S = {S 1 , . . . , S p }.

For each scheduling dimension of the legal solution a set of scheduling properties (parallel, sequential, permutable) is determined.

Notation and definitions

Here we give the necessary definitions and notations used throughout the rest of the chapter. Given a scheduling matrix Θ S for a statement S, we will denote the k-th row of the scheduling matrix Θ S as Θ S k . This row represents an affine scheduling hyperplane θ S k . The similar notation is used in works relating to program transformations and scheduling in the polyhedral model, bearing different names or meanings -schedules, mappings or partitions. [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF][START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part I: one dimensional time[END_REF][START_REF] Kelly | Minimizing communication while preserving parallelism[END_REF][START_REF] Lim | Communication-free parallelization via affine transformations[END_REF].

Recalling the defined canonical form of the scheduling matrix in section 2.3.4, the scheduling matrix Θ S could be split in three components: A, Γ and β. Thus, a k-th row of the matrix that represents an affine function has the following form:

θ S i k = a S i k,1 i 1 + a S i k,2 i 2 + . . . + a S i k,n i n + γ S i k,1 g 1 + . . . + γ S i k,n g g n g + β S i k
For the purpose of this chapter, we will ignore the Γ part and consider it to be set to all zeros. 1 Taking this into an account, we can represent a single row of the matrix as:

θ S i k = a S i k,1 i 1 + a S i k,2 i 2 + . . . + a S i k,n i n + β S i k ( 7.1)
Definition 7.5.1 (Affine form of Farkas Lemma). Given a nonempty polyhedron P = {x ∈ R n |Ax + b ≥ 0} then an affine form f : R n → R is non-negative at each point in a polyhedron P iff it is a positive affine combination of the polyhedron faces:

f (x) = λ 0 + m ∑ k=1 λ k (A k,• x + b k ), λ 0 , λ 1 , . . . , λ m ≥ 0 (7.2)

Dependence satisfaction

We will briefly recall the definition given in Section 3.3. Given a dependence graph G = (V, E), a valid schedule has to satisfy all the dependence edges e ∈ E. Each dependence edge e ∈ E is labelled by a dependence polyhedron P e . A given dependence edge e from S i to S j is satisfied iff:

∀(i, j) T ∈ P e : θ S i (i) ≪ θ S j (j) (7.3)
where ≪ is the lexicographical comparison operator. Since we are assuming multidimensional schedules, a given dependence edge e has to be satisfied at some level. We formally define this condition: Definition 7.5.2 (Dependence satisfaction at level). A dependence edge e from statement S i to statement S j , labelled by a dependence polyhedron P e , is satisfied at level l if it is the first level to satisfy the following:

∀(i, j) ∈ P e : ∀ 1≤k<l : θ S j k (j) -θ S i k (i) = 0 ∧ θ S j l (j) -θ S i l (i) ≥ 1
We say that the given dependence is strongly satisfied at level l, if the condition θ S j l (j) -θ S i l (i) ≥ 1 holds for all statement instance pairs that are in dependence. If only the weaker condition θ S j l (j) -θ S i l (i) ≥ 0 is satisfied, we say that a given dependence is weakly satisfied at a given level. Both definitions assume that the dependence is at least weakly satisfied at all previous levels 1 ≤ k < l. Definition 7.5.3 (Outer parallelism enabling hyperplanes). Given the set of statements S = {S 1 , . . . , S p } we say that the set of hyperplanes {θ S 1 1 , . . . , θ S p 1 } enable parallelism at the outermost level (l = 1) if the following condition is satisfied:

∀(i, j) ∈ P e : θ S j 1 (j) -θ S i 1 (i) = 0 , ∀e ∈ E
In other words, for all dependence edges ∀e ∈ E and for all pairs of statement iterations in a dependence relation P e , the source and the sink of the dependence are mapped to the same scheduling hyperplane. If we treat the scheduling hyperplane as a partition, as in [START_REF] Lim | Maximizing parallelism and minimizing synchronization with affine transforms[END_REF], it means that all the dependent iterations are executed within the same partition, while the partitions themselves could be executed in parallel. This kind of parallel execution is known as the DOALL or synchronization free parallelism.

Definition 7.5.4 (Permutable scheduling band). The set of scheduling hyperplanes {θ S 1 , . . . , θ S p } at levels l, l + 1, . . . , l + s -1 form a permutable scheduling band iff:

∀e ∈ E l , ∀k(l ≤ k ≤ l + s -1) : ∀(i, j) ∈ P e , θ S j k (j) -θ S i k (i) ≥ 0 where E l is the set of dependence edges not satisfied up to level l -1. The loops that correspond to scheduling hyperplanes within the permutable scheduling band could be permuted freely, without violating any dependence. The similar definition is used in several works [START_REF] Lim | An affine partitioning algorithm to maximize parallelism and minimize communication[END_REF][START_REF] Darte | Combining retiming and scheduling techniques for loop parallelization and loop tiling[END_REF][START_REF] Irigoin | Supernode partitioning[END_REF], where the goal is to find the maximal bands of fully permutable loops.

We have given the basic definitions that form the building blocks of our search methodology. Please note that some definitions are local -they apply for one dependence edge e at the time. In our transformation search method we will combine those constraints, step by step, into the global solution that satisfies the full dependence graph G = (V, E).

The space of legal affine transformations

As stated before, the goal of any program transformation based on the polyhedral model is to find the coefficients of the scheduling matrices: T = {θ S 1 , . . . , θ S p }.

Taking the full dependence graph G = (V, E), each dependence edge e ∈ E is described by a dependence polyhedron P e . A dependence edge e from S i to S j mandates that a condition shown in Equation 7.3 must be strongly satisfied. If we now consider the set of one-dimensional schedules, this condition becomes:

∀(i, j) T ∈ P e : θ S j (j) -θ S i (i) > 1 (7.4)
But the scheduling matrices that encode the scheduling functions θ S i and θ S j have unknown coefficients. If we use the representation of the single row of the matrix as an affine function, as shown in Equation 7.1, we will get the following form of the constraint:

∀(i, j) T ∈ P e : a S j 1,1 j 1 + . . . + a S j 1,n j n + ω S j 1 -(a S i 1,1 i 1 + . . . + a S i 1,n i n + ω S i 1
) > 1 But this constraint is nonlinear, since both loop induction variables and coefficients are unknown. This obstacle can be avoided by linearizing the constraints by using the Farkas lemma. All the affine scheduling algorithms 2 rely on the Farkas lemma [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part I: one dimensional time[END_REF][START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part II: multidimensional time[END_REF][START_REF] Pouchet | Iterative optimization in the polyhedral model: part ii, multidimensional time[END_REF][START_REF] Lim | Improving Parallelism and Data Locality with Affine Partitioning[END_REF][START_REF] Pouchet | Loop transformations: convexity, pruning and optimization[END_REF][START_REF] Kelly | Finding legal reordering transformations using mappings[END_REF] to perform the linearization.

2. An alternative approach is to use an equivalent vertex method as in [START_REF] Gupta | Scheduling in the z-polyhedral model[END_REF][START_REF] Vivien | On the optimality of feautrier's scheduling algorithm[END_REF] To perform the linearization, one has to represent the dependence polyhedron P e as an intersection of affine inequalities: P = {x ∈ R n |Ax + b ≥ 0}. By applying Equation 7.2 and substituting the dependence satisfaction constraint on the left side of the equation and the affine combination of the faces of the dependence polyhedron P e on the right side, one will get the Farkas multipliers -the unknowns λ 0 , λ 1 , . . . , λ m ≥ 0. Afterwards, one can project out those multipliers and obtain the system of affine inequalities describing the set of transformation coefficients that provide a legal schedule w.r.t. the dependence polyhedron P e -this set itself is a polyhedron L P e .

As an example, let us take the code shown in Figure 7.4. In this case we have only one statement and one dependence edge e : S 1 → S 1 . The dependence is the true (read-after-write) dependence described by the polyhedron:

P e = {(i, j, i ′ , j ′ )|i ′ = i + 1 ∧ j ′ = j + 1 ∧ 0 ≤ i ≤ N -2 ∧ 0 ≤ j ≤ N -2}
Our goal is to find the unknown coefficients a S 1 1,1 , a S 1 1,2 , ω S 1 1 of the one-dimensional affine scheduling function θ S 1 1 satisfying the dependence constraints:

∀(i, j, i ′ , j ′ ) ∈ P e : a S 1 1,1 i ′ + a S 1 1,2 j ′ + ω S 1 1 -(a S 1 1,1 i + a S 1 1,2 j + ω S 1 1 ) > 1 (7.5)
In the case of uniform dependences, application of the Farkas lemma could be simplified. 3 The distance between source and sink iteration is constant -i ′ = i + 1 and j ′ = j + 1. If those equalities are put into the previous constraint system, the loop induction variables are canceled out. After simplifying we get the constraint:

a S 1 1,1 + a S 1 1,2 ≥ 1.
This inequality defines an unbounded polyhedron L P e containing the possible values for scheduling coefficients. This set is an unbounded set -there is no upper nor lower bound on the coefficient values. Please note that coefficient ω S 1 1 is completely unconstrained. The legal solution polyhedron is described as follows:

L P e = {(a S 1 1,1 , a S 1 1,2 , ω S 1 1 )|a S 1 1,1 + a S 1 1,2 ≥ 1} (7.6) Figure 7.1 shows an iteration domain, and uniform dependences for the program in Figure 7.4. One of the possible choices for the coefficients is a S 1 1,1 = 1, a S 1 1,2 = 1, ω S 1 1 = 0, which gives the skewed scheduling hyperplane θ S 1 1 = i + j shown in Figure 7.1. The skewed scheduling hyperplanes are executed sequentially, while the iterations that belong to the same hyperplane might be executed in parallel. But the simpler solution is to take a S 1 1,1 = 1, a S 1 1,2 = 0, ω S 1 1 = 0, which gives an orthogonal scheduling hyperplane -the one that corresponds to the original loop iterator i.

The goal of all the affine scheduling algorithms [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part I: one dimensional time[END_REF][START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part II: multidimensional time[END_REF][START_REF] Pouchet | Iterative optimization in the polyhedral model: part ii, multidimensional time[END_REF][START_REF] Lim | Improving Parallelism and Data Locality with Affine Partitioning[END_REF][START_REF] Pouchet | Loop transformations: convexity, pruning and optimization[END_REF][START_REF] Kelly | Finding legal reordering transformations using mappings[END_REF] is to choose the best coefficient values among those contained inside the legal solution polyhedron L P e . The best coefficients are chosen according to some well defined criteria -maximal parallelism [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part I: one dimensional time[END_REF][START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part II: multidimensional time[END_REF], minimal synchronization [START_REF] Lim | Improving Parallelism and Data Locality with Affine Partitioning[END_REF] or minimal memory traffic [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF] for an instance. Another approach is to exhaustively search for all the possible schedules, given some bounds on the coefficient values. That is the approach that iterative compilation takes [START_REF] Pouchet | Iterative optimization in the polyhedral model: Part I, one-dimensional time[END_REF][START_REF] Pouchet | Iterative optimization in the polyhedral model: part ii, multidimensional time[END_REF]. Pouchet shows that it is feasible to exhaustively search the space of all the legal and distinct schedules, provided that the scheduling coefficients are bounded. As shown empirically by Pouchet, bounding the coefficients to integer values in the interval [-1, 1] gives the schedules that are usually expressive enough [START_REF] Pouchet | Iterative optimization in the polyhedral model: Part I, one-dimensional time[END_REF].

In the example we have shown, there is only one dependence edge. One dependence edge gives one set of the constraints on the coefficients expressed as a polyhedron L P e . But in the general case, the dependence graph contains several dependence edges. To obtain the global solution, one has to take the 3. For a full example of applying a Farkas lemma, a reader is could refer to [START_REF] Pouchet | Iterative optimization in the polyhedral model: Part I, one-dimensional time[END_REF].

Multidimensional schedules

One-dimensional schedules are not always obtainable. In practice, there are many codes for which the set of legal coefficients for the one-dimensional schedules is empty. As shown by Feautrier [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part II: multidimensional time[END_REF], one needs to extend the notion of one-dimensional schedules with the multidimensional schedules. It has been proven [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part II: multidimensional time[END_REF] that there exists a multidimensional schedule for each dependence graph induced by the loop nests coming from imperative programs. Indeed, the multidimensional schedule corresponding to the original loop nest is always legal, and it is one of the possible solutions.

Having the one-dimensional schedule, it is always mandatory to strongly satisfy all the dependences at the first (and only) dimension. Having the multidimensional schedule, there is a combinatorial choice of the depth of the dimension that strongly satisfies the dependence.

Feautrier decouples the problem of selecting the dependence satisfaction depth from forming the set of feasible schedule coefficient values. By introducing the 0 -1 variable x e k ∈ {0, 1}, one can model whether the given dependence edge e ∈ E is satisfied weakly (x e k = 0) or strongly (x e k = 1):

∀(i, j) ∈ P e : θ S j k (j) -θ S i k (i) ≥ x e k
Each dependence edge e ∈ E has to be satisfied strongly at some depth l. Let us consider the set of scheduling hyperplanes θ S 1 k , . . . , θ S p k at different dimensions, k ∈ 1..m. If the dependence edge e is satisfied strongly (x e l = 1) at depth l, this edge could be removed from the system of constraints when considering deeper levels k > l. Putting it more formally, for a given edge e : S i → S j , the following is a sufficient condition for the legality of the multidimensional schedule w.r.t. the dependence edge e:

∃1 ≤ l ≤ m, x e l = 1 ∧ ∀k < l, x e k = 0 ∧ ∀k ≤ l, ∀(i, j) ∈ P e , θ S j k (j) -θ S i k (i) ≥ x e k (7.8)
This reasoning comes from the simple fact on lexicopositivity of the difference of the timestamp vectors [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part II: multidimensional time[END_REF]. 4 The fact that there is a combinatorial choice of the weak/strong satisfaction for each edge e ∈ E and for each scheduling level 1 ≤ k ≤ m has led to the heuristics that try to structuralize this decision problem. Feautrier [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part II: multidimensional time[END_REF] has proposed to maximize the number of strongly satisfied dependences at each level: . This is the greedy approach that tries to satisfy as much dependences as possible in the outer scheduling dimensions, and to expose the parallelism in the inner loops. Also, it aims at multidimensional schedules that have the minimal dimensionality.

Bondhugula [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF][START_REF] Bondhugula | Automatic transformations for communication-minimized parallelization and locality optimization in the polyhedral model[END_REF] aims at different objective: that of exposing the outer loop parallelism. Thus, the strong dependence satisfaction is moved to the higher dimensions (inner loops), and the parallelism and tiling opportunities are exposed in the outermost loop levels.

As an example of the impact of the choice of dependence satisfaction (strong/weak), let us again take an example in Figure 7.4. If we consider the first dimension of the multidimensional schedule, after applying the Farkas lemma and simplifying, the constraint is:

a S 1 1,1 + a S 1 1,2 ≥ x e 1 , x e 1 ∈ {0, 1}
There is a choice between satisfying the dependence strongly x e 1 = 1 or weakly x e 1 = 0. According to Feautrier's heuristic, one maximizes the number of the strongly satisfied dependences at each level, thus x e 1 = 1. This leads to the constraint a S 1 1,1 + a S 1 1,2 ≥ 1, which resolves the schedule within the first dimension. This enables the innermost loop parallelism, since the second scheduling dimensions is completely unconstrained (after strong satisfaction at the first dimension, the dependence edge is removed from consideration at deeper levels).

Figure 7.2 shows the minimum latency orthogonal hyperplanes θ S 1 1 = i providing the first scheduling dimension. The iterations belonging to the same scheduling hyperplane could be executed in any orderin parallel in particular. With Bondhugula's approach [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF] the goal is to minimize the dependence distance starting from the outermost scheduling dimensions. This translates to the heuristic which prefers to weakly satisfy the dependence at the outer levels and to strongly satisfy the dependence at the inner levels. This leads to the constraint a S 1 1,1 + a S 1 1,2 ≥ 0. In particular, it is preferable to obtain the outer parallel enabling scheduling hyperplane: a S 1 1,1 + a S 1 1,2 = 0. Since the first scheduling dimension is not strongly satisfying the edge, the satisfaction has to occur at deeper dimension. This forces the choice of the strong satisfaction at the second scheduling level a S 1 2,1 + a S 1 2,2 ≥ 1. By assigning the values to the scheduling coefficients for the first level, one gets a possible parallel enabling scheduling hyperplane θ S 1 1 = ij (a S 1 1,1 = 1, a S 1 1,2 = -1). For the second dimension, the scheduling hyperplane is θ S 1 2 = i (a S 1 2,1 = 1, a S 1 2,2 = 0). This multidimensional schedule enables the outermost parallelism -the outermost loop is a DOALL loop. This is depicted in Figure 7. 3. Depending on the underlying architecture, either the inner parallelism (Figure 7.2) or outer parallelism (Figure 7.3) might be more beneficial.

All the mentioned approaches stick to the one scheme of selecting the dimensions that are to be strongly satisfied. A new approach of Pouchet et al. [START_REF] Pouchet | Loop transformations: convexity, pruning and optimization[END_REF] tries to express the convex space of all the possible choices for the x e k variables and to enumerate those choices. While this seems a promising approach, the amount of the possible solutions is prohibitively huge.

In the later sections of this chapter we will show the approach of enumerating the possible decisions in a specially designed binary decision tree. This enables us to enumerate the different decisions, while still maintaining the reasonable problem sizes.

Discrete sets of one-dimensional legal solutions

While the previously mentioned affine scheduling algorithms consider a very expressive search space of all legal affine transformations, it is not always necessary nor desirable to consider the full space of affine transformations. What is more, Feautrier's affine scheduler [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part I: one dimensional time[END_REF][START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part II: multidimensional time[END_REF] or Bondhugula's [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF] Transformation name Matrices involved Interchange, skewing (unimodular), reversal A Strip-mining, tiling D, θ Pipelining (multidimensional) β

Parametric pipelining (multidimensional) Γ Reversal (unimodular) A, Γ Motion, fusion, fission (distribution) β Table 7.1 -Classical loop transformations and their polyhedral matrix representation combined parallelism and memory locality heuristic are best effort heuristics -they provide one and only one solution chosen by some cost function.

As evidenced by Pouchet [START_REF] Pouchet | Combined iterative and model-driven optimization in an automatic parallelization framework[END_REF], the best effort heuristics are not well adapted to different architectural features -the optimal solution according to heuristics might be optimal for one architecture and suboptimal for other architecture. The other problem is

State of the art

We take the idea of Darte [START_REF] Darte | Loop shifting for loop compaction[END_REF] stating that arbitrary affine transformations are not always desirable. Narrowing down the scope of affine transformations might result in the better control over the generated code.

Our goal is to enable the exhaustive search of the possible affine transformations while still narrowing down the search space to the reasonable size. We have already shown in Section 7.6 that the space of all the legal coefficients for affine transformation might be unbounded.

Bounding the coefficients to the [-1, 1] interval was proposed by Pouchet [START_REF] Pouchet | Iterative optimization in the polyhedral model: Part i, one-dimensional time[END_REF][START_REF] Pouchet | Iterative optimization in the polyhedral model: part ii, multidimensional time[END_REF] in his work on iterative compilation within the polyhedral model.

Pouchet also shows [START_REF] Pouchet | Iterative optimization in the polyhedral model: Part i, one-dimensional time[END_REF] that some coefficients of the affine transformation are less critical than others. This fact was investigated empirically by assessing the impact of different transformation coefficients on the performance of the transformed output program. That fact has led to the decoupled heuristic [START_REF] Pouchet | Iterative optimization in the polyhedral model: Part i, one-dimensional time[END_REF] -first the coefficients corresponding to the loop induction variables are found, and later the rest of the coefficients are completed.

Girbal has shown [START_REF] Girbal | Semiautomatic composition of loop transformations for deep parallelism and memory hierarchies[END_REF] that the different parts of the canonical scheduling matrix format correspond to different classical loop optimizations [START_REF] Allen | Optimizing Compilers for Modern Architectures[END_REF]. This is summarized in Table 7.1.

Narrowing the search space

If one wants to enumerate all the possible transformations within an affine transformation space, one has to bound the values of the transformation coefficients.

Instead of bounding the values of the coefficients of affine transformation, as done by Pouchet [START_REF] Pouchet | Iterative optimization in the polyhedral model: Part i, one-dimensional time[END_REF], we restrict the possible space of transformations right away at the beginning of our search. We restrict our search space to the following transformations:

-loop permutation -non-parametric shift -motion, fusion, distribution Strip-mining and tiling transformations are treated specially and they are applied afterwards. Restricting the search space to loop permutations only, we do not look for full unimodular transformation space (skewing, loop reversal). This enables us to restrict the search space considerably, but the resulting search space is non-convex. f o r ( i = 1 ; i < N; i ++) f o r ( j = 0 ; j < N; j ++) S1 : The motivation would be clear after following the subsequent sections. Loop permutation defines only the orthogonal hyperplanes. While this decision might restrict some possible transformations (like loop skewing for extracting fine grained parallelism) it is powerful enough for the purpose of finding outer and inner loop parallelism in many cases. The problem is similar to Kennedy's [5] loop selection problem but it is expressed in the much more general framework of the polyhedral model: we support the non-perfectly loop nests, non-uniform dependences, shifting for transformation correction and integrate this into the loop distribution/fusion decisions.

A[ i ] [ j ] = B [ i -1][ j ] ; S2 : B[ i ] [ j ] = A[ i -1][ j ] ;

Building non-convex sets of one-dimensional legal schedules

We have to take into an account the restrictions on the possible schedules that we have imposed upfront when considering the affine scheduling functions. Considering the A part of the scheduling matrix for each statement S i , for each row k we constrain the scheduling coefficients to the following form:

m ∑ s=1 a S i k,s = 1, a S i k,s ∈ {0, 1} (7.9)
In other words, each row is composed of all zeros except one position -say s -that contains an integer constant 1. This corresponds to selecting the loop s [5] at the scheduling level k. Let us consider an example in Figure 7.4. The general form of the affine scheduling function for the statement S 1 is:

θ S 1 1 = a S 1 1,1 i + a S 1 1,2 j + ω S 1 1
A respective affine scheduling legality constraint (shown in Section 7.6 ) is a S 1 1,1 + a S 1 1,2 ≥ 1. Taking into an account the constraints we have imposed in Equation 7.9, we get the following table of the possible solutions: solution no. a S 1 1,1

a S 1 1,2 ω S 1 1 1. 0 1 unconstrained 2.
1 0 unconstrained The solution is trivial: there are only two possible loop interchanges and all of them are legal. The shifting factor ω S 1 1 is unconstrained 5 . The solution set is non-convex. The actual solution set is the enumeration of legal loop permutations alongside with legal shifting factors for each selection of loop permutations.

Let us take a more complex example, containing two statements, S 1 and S 2 , and two dependence edges e 1 : S 1 → S 2 and e 2 : S 2 → S 1 . The dependence polyhedra describing those two dependence edges are the following:

P e 1 = {(i, j, i ′ , j ′ )|i ′ = i + 1 ∧ j ′ = j ∧ 1 ≤ i ≤ N -2 ∧ 0 ≤ j ≤ N -1} P e 2 = {(i, j, i ′ , j ′ )|i ′ = i + 1 ∧ j ′ = j ∧ 1 ≤ i ≤ N -2 ∧ 0 ≤ j ≤ N -1} 5.
The shifting transformation expresses the relative shift of two statements. It makes sense only if we have more than one statement. In the case of a single statement we will assume ω S 1 1 = 0.

a S 1 1,1 a S 1 1,2 a S 2 1,1 a S 2 1,2 convex part 1 0 1 0 ω S 2 1 -ω S 1 1 ≥ 0 0 1 0 1 ω S 2 1 -ω S 1 1 ≥ 1 Table 7.2 -Solution set L e 1 a S 1 1,1 a S 1 1,2 a S 2 1,1 a S 2 1,2 convex part 1 0 1 0 ω S 1 1 -ω S 2 1 ≥ 0 0 1 0 1 ω S 1 1 -ω S 2 1 ≥ 1 Table 7.3 -Solution set L e 2
The prototype schedules from both statements are:

θ S 1 1 = a S 1 1,1 i + a S 1 1,2 j + ω S 1 1 (7.10) θ S 2 1 = a S 2 1,1 i + a S 2 1,2 j + ω S 2 1 (7.11)
The legality condition is expressed as: θ S 2 1 (i ′ , j ′ ) -θ S 1 1 (i, j) ≥ 1 Since we have the uniform dependences, we can substitute i ′ = i + 1 and j ′ = j + 1 in the legality condition expression. After simplification we get:

(a S 2 1,1 -a S 1 1,1 )i + (a S 2 1,2 -a S 1 1,2 ) j + ω S 2 1 -ω S 1 1 + a S 2
1,1 ≥ 0 After applying Farkas lemma (the details could be found in Appendix A), we get the following system of constraints on the scheduling coefficients:

     a S 2 1,1 -a S 1 1,1 ≥ 0 a S 2 1,2 -a S 1 1,2 ≥ 0 ω S 2 1 -ω S 1 1 + 2a S 2 1,1 -a S 1
1,1 -1 ≥ 0 By enumerating the possible values of the coefficients that in addition satisfy the Equation 7.9, we get the set of solutions shown in Figure 7.2.

The solution set is now split into the non-convex part and convex part. The non-convex part corresponds to the enumeration of the legal loop selections -encoded within A part of the scheduling matrix. For each legal assignment of the a S i 1,k coefficients, the remaining constraints on the shifting part ω S i 1 are expressed in the convex polyhedron form.

The similar solution set is obtained for the dependence edge e 2 : S 2 → S 1 . The solution set is shown in Table 7.3.

While the solution sets represent the space of legal coefficients for single dependence edges, the aim is to obtain the global solution, as in Equation 7.7. For that purpose, the intersection of the solution sets has to be provided. After combining the solution sets L e 1 and L e 2 , we obtain the global solution L shown in Table 7.4.

a S 1 1,1 a S 1 1,2 a S 2 1,1 a S 2 1,2 convex part 1 0 1 0 ω S 2 1 -ω S 1 1 + 1 ≥ 0 ∧ ω S 1 1 -ω S 2 1 + 1 ≥ 0 0 1 0 1 ω S 2 1 -ω S 1 1 ≥ 0 ∧ ω S 1 1 -ω S 2 1 ≥ 0 Table 7.4 -Solution set L = L e 1 L e 2 A S
% % J J J J J J J J J J J J A S 2 1,• 1 2 % % K K K K K K K K K K K K A S 2 1,• 2
% % J J J J J J J J J J J J A S 2 1,• 1 2 % % K K K K K K K K K K K K A S 2 1,• 2 1
y y t t t t t t t t t t t t t At each depth k of the tree, we consider the nodes of the form A S k 1,• that represent the decision variables that encode the possible rows of the A part of the scheduling matrix for each statement S k . A given node A S k 1,• might have up to dim(S k ) legal child nodes. The legal child node is the node that does not lead to the ⊥ terminating node.

ω S 1 1 -ω S 2 1 ≥ 0 ⊥ ω S 1 1 -ω S 2 1 ≥ 1
Initially, we represent the sets for the pairs of statements only, so each tree would have a depth of two. When combining the trees into the global solution, we will obtain a tree that encodes the global schedule with a depth p, where p is the number of the statements whose schedules we are constructing. We will now show how to perform the basic set operations on n-ary decision treesthat we construct.

A S 0 1,• 1 2 A S 1 1,• 1 y y t t t t t t t t t t t t 2 % % J J J J J J J J J J J J A S 2 1,• 1 2 % % K K K K K K K K K K K K A S 2 1,• 2 1 y y t t t t t t t t t t t t t ω S 2 1 -ω S 1 1 ≥ 0 ⊥ ω S 2 1 -ω S 1 1 ≥ 1
Figure 7.9 -The result of the operator EMBED(L e 1 , {A S 0 1,• }) on two sets represented as n-ary decision trees it is mandatory that they have the same sets of decision variables.

Given a decision tree with a set of variables {X1, . . . , Xn}, an embed operation EMBED extends this set with one decision variable, say X n+1 , while not changing the elements contained in the original set.

As an example, let us take the original set in Figure 7.7. Figure 7.9 shows the result of the embed operation on the set L e 1 with a (hypothetical) decision variable A S 0 1,• . Now the decision tree is extended with one more decision variable -A S 0 1,• . The set of the solutions remains the same -the value of the A S 0 1,• variable does not impact the set of valid results represented in a set. The ordering of the decision variables has an impact on the size of the decision tree. Embedding the variable at the deeper levels increases the number of the nodes, since it requires the replication of that decision node.

Intersection operator

The intersection operator performs the set intersection. The result is the intersection of two solution sets, say L ′ and L ′′ , so that the resulting set L = L ′ ∩ L ′′ contains the solutions that are both legal in L ′ and in L ′′ .

The operator requires that both arguments contain the same set of decision variables and that their ordering is the same. If this is not the case, then the EMBED is used to introduce the common decision variables in both sets.

The intersection operation L ′ ∩ L ′′ = INTERSECT(L ′ , L ′′ ) proceeds by recursively following the decision variable nodes from top to bottom for both arguments and finding the set of common terminating nodes. For those terminating nodes that have legal solutions, the intersection of convex sets is found.

An example is shown in Figure 7.12 which is the result of applying the intersection operator on trees shown in Figures 7.10 and 7.11. A common set of legal terminating nodes for sets L ′ and L ′′ is the terminating node for X 1 = 2, X 2 = 2. Subsequently, the common terminating node is labelled by the intersection of convex sets P 1 and P 3 .

Projection operator

The projection operator PROJECT(L, X) takes the decision tree L containing the decision variable X and removes that decision variable from the representation while keeping the set of legal solutions

X 1 1 }} } } } } } 2 A A A A A A A X 2 1 2 A A A A A A A A X 2 2 1 }} } } } } } } P 1 ⊥ P 2 Figure 7.10 -set L ′ X 1 1 }} } } } } } 2 A A A A A A A X 2 1 2 @ @ @ @ @ @ @ @ X 2 2 1 ~~~~~~~⊥ ⊥ P 3 Figure 7.11 -set L ′′ X 1 1 }} } } } } } 2 # # G G G G G G G G G X 2 1 2 @ @ @ @ @ @ @ @ X 2 2 1
{ { w w w w w w w w w w

⊥ ⊥ P 3 ∩ P 2 Figure 7.12 -INTERSECT(L ′ , L ′′ ) X 1 1 }} } } } } } 2 A A A A A A A X 2 1 2 A A A A A A A A X 2 2 1 }} } } } } } } P 1 ⊥ P 2 Figure 7.13 -set L X 2 1 }} } } } } } 2 @ @ @ @ @ @ @ @ X 2 2 A A A A A A A 1 ~~~~~~~P 1 ⊥ ⊥ P 2
Figure 7.14 -Split tree after removing X 1 node selected by remaining decision variables 6 .

While removing the decision variable nodes corresponding to the variable X that is projected out, the decision tree is split into several subsets. To obtain the final solution, the union of those subsets is taken and represented with the remaining decision variables.

As an example, let us take the set in Figure 7.13. The goal is to project out the X 1 variable. After removing the node corresponding to decision variable X 1 , the tree is split in two parts shown in Figure 7.14. Those two parts are combined into the final solution by taking their union.

6. This operation is the 'non-convex' counterpart of the projection operation performed on convex sets.

A note on complexity

By splitting the solution set representation into the non-convex part represented by n-ary decision treesand convex part, operated on by using the regular convex optimization theory, we have a better control over the complexity of operations on those sets.

Initially we build |E| decision trees representing the sets L e , ∀e ∈ E. Each such a decision tree represents the legal solutions w.r.t to one dependence edge. Those trees have the depth of two -one level for each statement. Given an edge e : S i → S j , the tree that represents the legal solutions for this dependence might have at most dim(S i ) • dim(S j ) nodes. Each node might be labelled by one convex polyhedron. The worst case is rarely attainable, even in the case of simple uniform dependences.

Having the full dependence graph, one needs to build |E| decision trees. Let d be the maximal depth of the statement in a program. The worst case complexity of building the constraints for all the dependence edges is O(|E| • d 2 ).

In order to build the global solution, one needs to combine all those pairwise constraints into the global system. This is achieved by performing the sequence of intersect operations, which in turn might require an embedding operation. This simple scheme is shown in the form of pseudo-code in Figure 7. [START_REF] Banerjee | Loop Transformations for Restructuring Compilers[END_REF].

The complexity of the operations is sensitive to the order in which the intersection operations are performed. The most desirable is to intersect two sets, say L ′ , L ′′ , that have exactly the same set of decision variables -no embedding operation is needed in this case. The complexity of such an intersection is O(n 1 n 2 ), where n 1 and n 2 are the number of tree nodes.

Sets of multi-dimensional legal solutions

So far we have discussed the representation of one-dimensional solution sets. While the one-dimensional solution sets are the starting point for the exploration of transformations, the ultimate goal is to provide the set of legal multidimensional solutions, since one-dimensional schedules are not possible for some class of the programs.

Representing weak/strong satisfaction

When we were considering the sets of one-dimensional schedules only, the strong satisfaction is required in the first dimension -all the dependences have to be satisfied within the first and only dimension. If we extend the set of solutions to multidimensional schedules, the combinatorial problem becomes more complex. There is a combinatorial choice of whether to strongly or weakly satisfy the dependence edges e ∈ E at each scheduling level.

A S 1 1,• 1 y y t t t t t t t t t t t t 2 % % J J J J J J J J J J J J A S 1 1,• 1 y y s s s s s s s s s s s s 2 % % J J J J J J J J J J J J A S 2 1,• 1 2 & & L L L L L L L L L L L L L A S 2
P ω S 1 1 -ω S 2 1 ≥ -1 ω S 1 1 -ω S 2 1 > -1 ω S 1 1 -ω S 2 1 ≥ 0 ω S 1 1 -ω S 2 1 > 0 Figure 7.
% % L L L L L L L L L L L L L A S 2 1,• 1 2 ' ' O O O O O O O O O O O O O O O A S 2 1,• 2 1
x x q q q q q q q q q q q q q x e 1

' ' O O O O O O O O O O O O O x e 2 1 0 1 > > > > > > > > -1 ≤ ω S 1 1 -ω S 2 1 ≤ 1 -1 < ω S 1 1 -ω S 2 1 ≤ 1 -1 ≤ ω S 1 1 -ω S 2 1 < 1 -1 < ω S 1 1 -ω S 2 1 < 1 ω S 1 1 -ω S 2 1 = 0 ⊥
= 0 = 0 A S 1 1,• = (01) A S 2 1,• = (01) ω S 1 1 -ω S 2 1 = 0 Both edges strictly weakly satisfied. Enables outermost DOALL loop. > 0 > 0 A S 1 1,• = (10) A S 2 1,• = (10) -1 < ω S 2 1 -ω S 1 1 < 1
= 0 > 0 A S 1 1,• = (10) A S 2 1,• = (10) ω S 1 1 -ω S 2 1 = 1
Edge e 1 strictly weakly satisfied, edge e 2 strongly satisfied. Schedule is not complete at level 1, requires completion at second level.

=> = 0 A S 1 1,• = (10) A S 2 1,• = (10) ω S 2 1 -ω S 1 1 = 1
Edge e 2 strictly weakly satisfied, edge e 1 strongly satisfied. Schedule is not complete at level 1, requires completion at second level. In order to get the complete solution, one needs to complete the schedules at scheduling dimensions 1 < k ≤ dim(S i ). The complete schedule has the full-ranked A part of the scheduling matrix 7 . In order to get the complete, full-ranked scheduling for each statement S i , each statement needs to be provided with a full ranked A S i matrix, with the dimensionality of dim(S i ) × dim(S i ).

A S 1 = 0 1 1 0 A S 2 = 0 1 1 0 ω S 1 1 -ω S 2 1 = 0 -1 < ω S 2 2 -ω S 1 2 < 1
Bondhugula [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF] is using the concept of orthogonal sub-spaces so that the solutions found at the following rows of the scheduling matrix are linearly independent w.r.t. rows that have already been found. This guarantees the full-ranked matrix.

In our search space, the problem of finding the independent solutions is easier. We restrict the rows of A matrix to satisfy Equation 7.9. If we have selected the position s of the integer constant 1 at the first row of the scheduling matrix, then this position could not be used in subsequent rows.

In general, if we have selected independent positions s 1 , . . . , s k for levels 1 ≤ l ≤ k, those positions cannot be used when completing the subsequent rows k + 1 ≤ l ≤ dim(S i ). This is equivalent with stating that the A represents a loop permutation matrix. An example of the full ranked permutation matrix is shown in Figure 7.20. Obviously, the solution shown in Figure 7.21 is not a full-ranked multidimensional schedule.

Given the complete solution set for the first scheduling dimension represented as a decision tree, each valid leaf node could be completed with the next scheduling level that gives the independent solution. While it is possible to provide a completion for each leaf node of the decision tree, this procedure might lead to combinatorial explosion, since each subsequent level is a also a decision tree. In Section 7.10 we will discuss the systematic search methodology that enables us to control the complexity of completing the schedule at deeper scheduling levels.

If we want to complete the second scheduling dimension, we have to find a solution set that satisfies some subset of the dependence edges E ′ ⊂ E, while providing a row A S i 2,• , for each statement S i , that is linearly independent from a row A S i 1,• . We will denote this operation as ORTHOGONAL(E ′ , L, k), where E ′ is the subset of the edges that have to be satisfied at a new scheduling level k and L is the solution set of the solutions found at levels 1 ≤ l ≤ k -1.

Let us take a decision tree shown in Figure 7.19 that represents the solution set L 1 . We want to complete the scheduling matrices by providing the second scheduling dimension. Let us take one of the solutions in the set L 1 and provide an orthogonal completion for it. Figure 7.23 shows one chosen solution from the first dimension (emphasized with double arrows) and the orthogonal completion in the second dimension. After orthogonal completion in the second dimension, we have one complete solution:

If we now assign some concrete values for shifting factors:

ω S 1 1 = 0, ω S 2 1 = 0, ω S 2 2 = 0, ω S 1 2 = 0, the complete solution is: θ S 1 1 = j, θ S 1 2 = i
7. This is also necessary for correct code generation from the polyhedral model

Finding global solutions

In previous sections we have discussed the way to build the pairwise (per dependence edge e) solutions sets, the way to combine then in the global solution for one scheduling level, and a technique to complete the scheduling with independent solutions at further dimensions. Also, we have shown how to encode the decision on whether to strictly/weakly satisfy a given dependence edge.

On top of those operators and data structures, we build a tunable search methodology that enumerates the possible multidimensional schedules. It might be exhaustive -it will enumerate all the possible schedules and edge satisfaction conditions, or it might be goal driven -it will filter out only those schedules that have a certain property (to be defined later), giving an empty solution set if there is no such a schedule.

In addition to enumerating the multidimensional schedules, the method reports the different properties for each multidimensional schedule. For our purposes, the properties are:

1. A given dimension (loop) level could be executed in parallel 2. A given dimension (loop) must be executed sequentially

A given band of dimensions forms a permutable loop band

Some notation

Some extra notation is introduced for the sake of the search algorithm:

The solution set of legal transformations is denoted by L. The solution set satisfying dependence edge e 1 is denoted by L e 1 . The solution set satisfying the set of dependence edges E is denoted by L E . Each single solution T ∈ L of the solution set L consists of part that encodes the A part of the rows A k and satisfaction decision variables x e∈E k . We denote those the A part of the solution as T A . If a particular T A solution is selected in a solution tree, the remaining part T X encodes the strong/weak satisfaction of the edges e ∈ E. Each assignment of the decision variable x e∈E k that leads to a non-empty solution selects some subset E S ∈ E of the edges that could be strongly satisfied. Indeed, the full decision tree of the decision variables x e∈E k determines the family X of the possible strong satisfy subsets E S ∈ X . As an example of the solution set, let us take a solution set L E shown in Figure 7.19. Here an edge set is E = {e 1 , e 2 }. One particular solution is T

= {A S 1 1,• = 1, A S 2 1,• = 1, x e 1 1 = 0, x e 2 1 = 0}. The A part of the solution is T A = {A S 1 1,• = 1, A S 2 1,• = 1}. If the A part T A of
the tree is selected, one gets the the sub-tree of decision variables. In this example, if one selects T A = {A S 1 1,• = 1, A S 2 1,• = 1} one will obtain the tree rooted at x e 1 1 that represents the family X of possible edge satisfaction sets E S . In this example, the family is equivalent to the powerset 2 E of all edges.

Recursive procedure

Here we show a complete recursive procedure for computing the multidimensional schedules. For simplicity, we only show the procedure for obtaining one solution, but it could be easily adapted to explore the space of more complete schedules -it is a matter of changing the termination condition of the main loop and implementing a way of controlling the desired number of candidate solutions at each scheduling depth.

The overall recursive procedure for enumerating the multidimensional schedules is shown in Figure 7.24. An initial set of solutions for the first scheduling dimension is obtained:

L ← COMBINEINITIALSOLUTIONS(E)
The procedure for computing the initial solution set was shown in Figure 7.16. The recursive procedure is called with the following arguments:

COMPLETESCHEDULE (L, E, d = 1)
The initial solution set is L, the full set of the dependence edges E is to be satisfied, and the building process is starting with the first scheduling dimension d = 1.

In the case the procedure is invoked with d = dim max (S ), dim max (S ) being the maximum domain dimensionality of the statements, 8 one complete multidimensional schedule is obtained. This is the terminating condition for the recursion. At step (1.1) the schedule properties are provided and the complete schedule is saved at step (1.2). By returning Trueat the step (1.3), it is signalled to the calling procedure that the complete schedule is obtained.

A loop starting at step (3) is executed until at least one solution is obtained. As mentioned earlier, this loop could be changed so that more than one solution is obtained if desired.

At each iteration of the loop the next solution is selected at step (3.1) and its A part is stored as T A . At step (3.2) an orthogonal subspace for the transformation T A is obtained.

A given transformation T A determines an orthogonal subspace at the deeper levels. The core part of the search algorithm is to determine the solution set L ′ for the next scheduling level such that H ∩ L ′ = / 0.

If no such a set could be obtained, then the given transformation T A at the level d could not be completed with independent solutions at deeper levels.

There is a combinatorial number of ways to provide the set L ′ . Indeed the set L ′ = LE -E S . There is a combinatorial number of ways to select the strong satisfaction edge set E S . In a special case when E S = / 0, the L ′ = L E = L. This is the case when we can have a permutable loop band.

In order to avoid useless backtracking steps, a quick test is performed to determine whether a given T A could have a non-empty completion at deeper levels. First, a set of edges E must that must be removed at the current level d is determined at step (3.3). We will show the way to perform this check in Subsection 7.10.2. Later, at steps (3.4) and (3.5), it is determined whether this set could be strongly satisfied by a selected transformation T A . If it cannot be strongly satisfied, then the set T A is not taken into the consideration anymore.

At step (3.6) the family of strong satisfaction sets is restricted to family X ′ that contains only the satisfaction sets E S that are supersets of those edges that must be satisfied: E S ⊃ E must . This steps narrows down the search space of strong edge satisfactions that is traversed at further steps.

The loop at step (3.9) performs the traversal of the strong satisfaction sets E S . A naive approach is to traverse all the E S ∈ X ′ , in any order. The order in which those sets are traversed has an impact on the convergence of the method. We will show the different strategies of the traversal order. A predicate EXHAUSTED controls the number of the solutions that one wants to investigate. In any case, strategy < N,

where N = X ′ .

Once a given strong satisfaction set E S is selected, the remaining set E ′ of those edges that must be satisfied at deeper levels is obtained at step (3.9.2). A check is made at step (3.9.3) whether removing those edges gives a solution L E ′ that could be completed at deeper d + 1. If this is not the case, the next satisfaction set E is checked according to some new strategy ( a strategy counter is increased). Otherwise, a current solution T A is completed with a recursive call to the schedule completion procedure for deeper levels made at the step (3.9.5).

Dependence satisfaction strategies

As mentioned in Subsection 7.6.1, the problem of selecting the dependence edges to be strongly satisfied at a given loop depth d has so far been solved by proposing heuristics [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part II: multidimensional time[END_REF][START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF][START_REF] Bondhugula | Automatic transformations for communication-minimized parallelization and locality optimization in the polyhedral model[END_REF]. Darte [START_REF] Darte | Complexity of multi-dimensional loop alignment[END_REF] has formally shown that this problem is NP-complete. Given a dependence graph G = (V, E) having |E| edges, there are 2 |E| possible subsets of the edges to be strongly satisfied. {If the set E must does not belong to X then break early } One can exhaustively try all the possible subsets -a naive enumeration method. We rather propose to rank those subsets according to some preferable properties. This oracle, named SELECTSATSET, is called from the enumerative search algorithm shown in Figure 7.24.

(3.5) IFE must / ∈ X CONTINUE {Select the family X ′ ⊂ X such that X ′ = {E S ∈ X : E S ⊃ E must } } (3.6) X ′ ← SELECT(X , E must ) (3.7) N ← |X ′ | (3.8) strategy ← 0 (3.9) WHILE (¬ f ound ∧ EXHAUSTED(strategy, N)) (3.9.1) E S ← SELECTSATSET(X ′ , strategy) (3.9.2) E ′ ← E -E S (3.9.3) IF H ∩ L E ′ = / 0 strategy ← strategy + 1 CONTINUE (3.9.5) f ound ← f ound ∨ COMPLETESCHEDULE(L E ′ , E ′ , d + 1) (4
The first and most preferable strategy is the case when E S = / 0. This choice forces the next scheduling level d + 1 to satisfy all the dependence edges E satisfied at previous level d. If this is possible, we will obtain a permutable loop band d . . . d + 1 as in [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF]. If this is not a feasible solution, then at least one edge e ∈ E has to be strongly satisfied. The joint work of Pouchet and Bondhugula [START_REF] Pouchet | Loop transformations: convexity, pruning and optimization[END_REF] shows the simple heuristic that decides to satisfy the maximal number of edges at a given depth if the E S can not be empty. We also believe that following this scheme mostly likely increases the chance of getting the non-empty completion of the schedule at the next level.

Thus, the next preferable strategy is the one that selects {E S ∈ X : max|E S |} As Kennedy has shown in his approach to loop selection problem [START_REF] Allen | Optimizing Compilers for Modern Architectures[END_REF], the greedy heuristic of selecting the solution that covers the maximal number of the edges at once might fail in the case when this maximal subset does not cover those edges that prevent obtaining the solution.

We improve the heuristic of Pouchet and Bondhugula by a shortcutting step since we determine upfront the set of those edges that must be satisfied E must . This set is not determinable in the case of general affine schedules. But in our case we have a restricted search space and we also have an enumeration of the possible solutions built upfront.

The goal is to obtain H

∩ L ′ = / 0
The goal of the procedure MUSTSATISFY(H) called at step (3.3) of the algorithm in Figure 7.24 is to obtain the set E must ⊂ E S of those edges that must be satisfied in order to H ∩ L E-E S = / 0. Obviously if L e i ∩ H = / 0 then the dependence edge e i has to belong to the set E must .

Although the pruning step we have shown can cut down the search space of possible satisfaction sets, it only gives the necessary condition, but not the sufficient condition for H ∩ L E-E S = / 0.

Distribution/fusion scheme

The search method discussed so far does find a complete solution in terms of complete A part of the scheduling matrices. But it does not represent the loop distribution(fusion) schemes since the β part of the scheduling matrices is not constructed. Indeed, the scheme presented implicitly assumes that all the statements are fused (β part of all the schedule matrices being the zero vector).

Forcing all the statements to be fused within one loop might prevent obtaining any legal solution. Indeed, if an original program contains non-perfectly nested loops, some statements do not belong to the same loop -they are distributed. It is not always legal to transform such a program to a completely fused version.

A classical literature on loop fusion and distribution [START_REF] Kennedy | Maximizing loop parallelism and improving data locality via loop fusion and distribution[END_REF][START_REF] Megiddo | Optimal weighted loop fusion for parallel programs[END_REF][START_REF] Darte | On the complexity of loop fusion[END_REF][START_REF] Darte | Loop shifting for loop compaction[END_REF] shows the benefits of both: distribution enables more transformations, including parallelism, while loop fusion is beneficial for optimizing the memory access locality by reducing the distance between producer and consumer pairs. Darte [START_REF] Darte | On the complexity of loop fusion[END_REF] has shown that the problem of determining the optimal loop fusion structure is NP-complete.

The polyhedral model naturally incorporates the loop fusion/distribution decisions, through β part of the scheduling matrix. A systematic way to approach the fusion/distribution problem in the polyhedral model was shown in [START_REF] Bondhugula | A model for fusion and code motion in an automatic parallelizing compiler[END_REF][START_REF] Pouchet | Combined iterative and model-driven optimization in an automatic parallelization framework[END_REF]. While expressing the set of legal fusion/distribution structures in the polyhedral model is not difficult problem [START_REF] Pouchet | Combined iterative and model-driven optimization in an automatic parallelization framework[END_REF], selecting the desired structure out of the combinatorial space is hard problem.

We do not incorporate the search for fusion/distribution structures in our search space. The reason is that we want to decouple this particular search problem from the problem of searching the complete multidimensional schedules.

CONCLUSIONS 123

Our original combinatorial problem is the problem of searching for sets of edges to be satisfied strongly at particular scheduling depth . This problem is modelled as a decision tree that represents the possible subsets of the edges that could be satisfied strongly at the same time.

The combinatorial problem of enumerating the fusion structures is the problem of enumerating the legal partitions of the set of statements S = {S 1 , . . . , S p }. If p is the number of statements, the total number of partitions is the Bell number B p , which gives the worst case complexity of the traversal. We perform the maximal distribution of the statement sets, according to the strongly connected components of the dependence graph G = (V, E). This corresponds to one of the heuristics used by Bondhugula [START_REF] Bondhugula | Effective automatic parallelization and locality optimization using the polyhedral model[END_REF] in his scheduling approach.

The reason for choosing the maximal distribution strategy is the following: one can always maximally distribute the statements into partitions according to the SCCs of a dependence graph. In addition, this choice gives the most freedom in choosing the subsequent multidimensional schedules. This guarantees us that we will get the most expressive space of legal multidimensional schedules within our search space. Once the full space of fully distributed multidimensional schedules is obtained, one might use one of the fusion strategies already developed in [START_REF] Bondhugula | A model for fusion and code motion in an automatic parallelizing compiler[END_REF][START_REF] Pouchet | Combined iterative and model-driven optimization in an automatic parallelization framework[END_REF] to get the desired fusion structure, if it is possible to merge the solution sets into tighter constraints given by fusion.

In order to incorporate the max distribution strategy into the search procedure, we slightly modify the search procedure. At each depth of the search procedure shown in Figure 7.24, after selecting the set of the edges E S to be removed at the next scheduling level, reform the dependence graph G ′ = (V, E -E S ) and compute the set of SCCs for this new graph G ′ ⊂ G. Topologically sort SCCs, and form an ordered sequence (SCC 1 , . . . , SCC m ). For each statement S i ∈ SCC j assign its static scheduling part (β part) to j, where j is the position of the SCCin a topological order. After this step, recursively call the procedure COMPLETESCHEDULE for each SCCindependently.

Conclusions

We have presented the general search strategy used for traversing the legal transformation search space. The current approaches to automatic transformation in the polyhedral model are based either on linear programming (one-shot, best-effort heuristics of Feautrier [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part I: one dimensional time[END_REF] and Bondhugula [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF]) or on the exhaustive search of the space of legal affine transformations (iterative compilation of Pouchet [START_REF] Pouchet | Iterative optimization in the polyhedral model: Part i, one-dimensional time[END_REF]).

We propose to represent the set of legal solutions as a discrete set of selected and legal loop permutations, together with enabling shifting factors that are represented as the convex sets. We use the decision diagram data structure to represent those solution sets. Our search space is a subset of the full space of affine transformations.

Contrary to the iterative approach of Pouchet, which traverses the full space of legal and distinct affine transformations, our search space is restricted, but this enables the control over the complexity of the solution. On one side, the solution set is restricted to the discrete set of legal loop permutations, but on the other side, each legal loop permutation is augmented with a convex set of shifting factors that decides whether this permutation satisfies the dependence strongly or weakly.

The fact that we have an explicit control over the weak/strong satisfaction of the dependence edges, enables us to traverse the space of multidimensional schedules in an exhaustive way (if needed). This is not possible in the current multidimensional affine schedule space traversal approach of Pouchet [START_REF] Pouchet | Iterative optimization in the polyhedral model: part ii, multidimensional time[END_REF]. Pouchet employs a greedy heuristic that forces the unique choice of dependence edge satisfaction at each scheduling level, to keep the search space tractable 9 . Our contribution achieves a degree of flexibility and controllability of the search space that is not possible with the current iterative approach. 9. Though the latest contribution [START_REF] Pouchet | Loop transformations: convexity, pruning and optimization[END_REF] shows (only) a theoretical way to achieve this flexibility Nevertheless, our search space is not as exhaustive as the one of Pouchet. Because we restrict the scheduling coefficients of the dynamic scheduling component ( Section 3.1) to the loop permutation only, we might miss some important transformations such as loop skewing, loop reversal or loop slowing. But this is a necessary price to pay if we cannot accept the running time of the exhaustive iterative search used in the feedback-directed iterative compilation [START_REF] Pouchet | Iterative optimization in the polyhedral model: part ii, multidimensional time[END_REF].

The future work

The search strategy might be improved in several directions. First, the algorithmic improvements are necessary in order to avoid the overhead of recomputing the same solution sets during the construction of multidimensional schedules. This could be achieved by caching the already computed sets to avoid the redundant recomputations.

The representation of the convex part of the solution relies on the operations on the polyhedral sets, that in turn require the Chernikova's algorithm [START_REF] Chernikova | Algorithm for finding a general formula for the non-negative solutions of a system of linear inequalities[END_REF][START_REF] Verge | A note on Chernikova's algorithm[END_REF]. But it could be observed that our convex sets have a special form, that of the two-variable per constraint [START_REF] Cohen | Improved algorithms for linear inequalities with two variables per inequality[END_REF]. Those sets are the subsets of general polyhedra, and they could be operated on with less costly techniques than the general techniques used for polyhedral sets [START_REF] Upadrasta | Potential and challenges of two-variable-per-inequality subpolyhedral compilation[END_REF].

Our search strategy enables the branch-and-bound metaheuristic. But in order to make it effective, we have to provide a systematic way of pruning the search space, cutting the non-profitable branches of the search tree. This could be achieved by having a partially evaluated cost functions that could provide the estimate of the cost, even without having a completed schedule. Our current cost function (Chapter 6) does not have this property.

Chapter 8

Conclusions and perspectives

In this thesis we have presented the theory, design and implementation of the search based program transformation strategy based on top of the three-address code polyhedral compiler. We have shown that in order to obtain the precise execution cost-model function, it is necessary to integrate the polyhedral model framework directly into the three-address code compiler. This is in contrast with the traditional approach for polyhedral compilation based on source-to-source compilers.

Integrating the polyhedral compilation framework directly into the three-address code compiler posed several challenges that were not investigated in the known literature. We have provided some efficient solutions for those problems, concerning the efficient representation of the scalar and memory-based dependences. We have also presented an unifying approach for representing all the transformations, including tiling purely through the scheduling functions.

We have shown a precise execution cost-model function that precisely captures the low-level details, such as SIMD vectorization, of the target architecture. Contrary to simplistic linear cost-model functions, our function is a complex, non-linear characterization. The evaluation of this function requires a search based strategy, contrary to the approaches that use linear programming to directly optimize the simple linear cost-functions. Our cost-model function is evaluated for each point of the constructed search space.

Providing an efficient search strategy requires the construction of the expressive, but size-limited search space of legal program transformations. We took a novel approach: contrary to the current iterative optimizers that explore the huge space of affine schedules, we build a search space of discrete set of transformations. Our search space is a subset of the full affine search space. For an efficient representation of the sets we propose to use the decision diagrams and we provide the basic operators on those sets.

For assessing the feasibility of the approach, we have provided the practical implementation of the techniques in GRAPHITE polyhedral framework that is a part of the production quality compiler GCC. We will summarize the detailed contributions of this thesis:

Three-address code polyhedral compilation framework Traditional polyhedral program transformation frameworks are implemented as the source-to-source program translators. This design is natural and practical for research, but it is somewhat simplistic when it comes to precise modelling of the target machine execution model. We provide the design and implementation of the sophisticated polyhedral transformation framework operating directly on the three-address code in SSA form. We have provided an in-depth investigation of the issues related to the direct polyhedral compilation of the three-address code. To the best of our knowledge, it is one of the first widely available and published polyhedral compiler that operates on the three-address code.

Lazy approach for relaxing scheduling constraints

The problem of effectively handling memory-based data dependences is well-known. This problem is particularly exposed in the three-address code polyhedral compiler, since the low-level code contains many temporary variables for storing intermediate results. Those temporary variables induce many memory-based dependences. The traditional approach to eliminating memory-based dependences is the memory expansion -expanding scalars into arrays, or arrays into higher-dimensional arrays. This approach eliminates all memory-based dependences, but it might have a prohibitive cost. We propose the lazy approach for handling memory expansion. The memory-based data dependences are firstly ignored, and the live range violation analysis is used to compute the minimal set of dependences that have been violated. Only the violated memory-based dependences are removed through expansion. This approach works well with search based strategy, since the cost of the memory expansion might be evaluated for multiple transformation candidates.

Semantical transparency

The traditional source-to-source compilers are sensitive to syntactical details of the input program. The same computational kernel that is written in a slightly different style (by introducing the scalar temporaries for example) might be ignored by an optimizer, simply because it does not conform to the syntactical constraints of the source-to-source compiler. This limitation is one of the most important reasons why the current polyhedral compilers are not able to transparently compile legacy codes or industry-standard benchmarks such as SPECint or SPECfp. Providing the polyhedral framework that is part of the production quality compiler, such as GCC, and operating on the low-level intermediate representation provides the semantical transparency for the end user -an user might still write the programs in the sequential style, while the compiler is responsible for extracting the intrinsic semantics of the written program. This form of the semantical transparency is possible in GRAPHITE due to the fact that the polyhedral framework operates at the stage where the code is transformed into SSA form on which the essential scalar optimizations are performed. A code in such a form captures the essential semantics of the input program and is much less sensitive to the peculiarities of the syntactic form of the program. In addition, the lazy approach to memory-expansion can transparently remove the spurious data dependences that might be inadvertently introduced by the programmer.

Precise performance predicting cost-model

The one-shot, best-effort scheduling heuristics are based on using (integer) linear programming to optimize simple linear objective cost functions. The linear cost-models are based on an abstract computational metrics, and they are not adequate for modelling the complex, low-level aspects of the target machine execution models. One such an aspect is SIMD vectorization, which is highly target specific. We propose a novel approach to this problem by proposing to use machine specific accurate cost-model functions. We demonstrate the feasibility of the approach, by constructing the precise cost-model function for SIMD architectures. Those precise machine-specific cost-models cannot be expressed as linear functions and they could not be handled by linear programming machinery. Instead, they have to be evaluated at each point of the search space. The critical aspect for the effectiveness of our cost function evaluation is the construction of the feasible search space. This relies on our next contribution, which is a construction of effective transformation search space.

Efficient transformation search space

We propose the novel method for constructing the search space of legal polyhedral transformations. Our method is based on providing the discrete sets of legal transformations. Those sets are represented as decision diagrams and they are subset of the full search space of affine transformations. Discrete sets limit the possible affine transformations, but they enable us to control the complexity of the search space. This is crucial contribution in enabling the search-based methodology that is efficient. In the case of multidimensional scheduling, we provide more flexibility: we enable the enumeration of the dependence satisfaction strategies. The current iterative transformation approaches that construct the space of multidimensional schedules take the greedy heuristic for selecting the data dependence satisfaction at different scheduling levels, to keep the search spaces tractable. They could miss some points in the search space that our enumeration strategy enumerates. We show that we can afford this level of flexibility, since we restrict the legal solution sets early in the search process.

Perspectives

We have provided a novel approach for the polyhedral program transformations based on the search strategies that optimize the complex, machine dependent cost-model functions. Still, there are numerous possible improvements and future directions that are worth investigating. We will give brief ideas for the possible improvements and directions for the prospective research topics:

Even more precise cost-model The precise cost model function that we have proposed includes a very rough estimation of the data locality. While it can capture the spatial locality of the given memory access pattern, it cannot precisely model the temporal locality. The model could be extended with cache miss equations [START_REF] Fraguela | Probabilistic miss equations: Evaluating memory hierarchy performance[END_REF] that fit very well with the current model. The architectural interplay between SIMD parallelism and thread-level coarse grained parallelism [START_REF] Dubey | Evaluating performance tradeoffs between finegrained and coarse-grained alternatives[END_REF] could be modelled in the cost function as well.

Parametric cost function

The cost-function that we have presented works on a premise that the loop iteration domains are known at the compile time. This is necessary because the total cost depends on the exact number of iterations of each loop. This also simplifies the comparison of the costs, since the costs are expressed as simple integer numbers. But the polyhedral model allows the parametric loop bounds -bounds that are expressed as symbolic constants that are not known at compile time. We would like to investigate the technique of computing the number of iterations [START_REF] Clauss | Counting solutions to linear and nonlinear constraints through Ehrhart polynomials: applications to analyze and transform scientific programs[END_REF] for parametric iteration domains. This would require representing the costs as symbolic expressions, rather than integer constants. This could however fundamentally change the approach: the best solution will depend on the parameter values, and it could not be obtained automatically. Some assistance from the user, who understands the domain of the problem and can estimate the problem size, would be required.

Machine learning assisted portability of the cost-model

The analytical cost-model function that we have presented relies on the machine specific instruction costs. We have obtained those costs for several architectures, through microbenchmarking and knowledge of the exact instruction latencies. Obtaining the instruction costs for the new architectures is a time-consuming process which hinders the portability of the cost-model based approach. In order to facilitate the porting of the cost-model to the new architectures, it interesting to investigate the possibility of employing machine learning, like in [START_REF] Cavazos | Automatic performance model construction for the fast software exploration of new hardware designs[END_REF], to automatically obtain the instruction costs and other platform-specific coefficients that are used in the cost-function.

Scalability

We propose the search space construction method that can control the explosion of the transformation search space by limiting the number of non-convex solutions. Still, we leave some parametric part (the shifting factor) that is represented as convex polyhedra. This part is currently operated on by the standard polyhedral techniques. But the special form of this parametric part is amenable to representation as the two-variable per constraint system [START_REF] Cohen | Improved algorithms for linear inequalities with two variables per inequality[END_REF], which could be operated on with less costly techniques than the general techniques used for polyhedral sets [START_REF] Upadrasta | Potential and challenges of two-variable-per-inequality subpolyhedral compilation[END_REF]. An interesting research direction is the investigation of the sub-polyhedral domains. Sub-polyhedral domains are less powerful than the traditional polyhedral techniques, but the complexity of operations on those sets is guaranteed to be polynomial or linear.

Extending the scope of the analyzable programs

The traditional polyhedral model representation is restricted to static control programs (SCoP). This is the fundamental restriction of the polyhedral model. Programs with irregular control, non-affine memory accesses and loop bounds are not amenable to analysis within the polyhedral model. The methods [START_REF] Benabderrahmane | The polyhedral model is more widely applicable than you think[END_REF][START_REF] Barthou | Fuzzy array dataflow analysis[END_REF][START_REF] Groesslinger | The Challenges of Non-linear Parameters and Variables in Automatic Loop Parallelisation[END_REF] to overcome those limitations have been proposed, but they are still not widely used, due to the complications in code generation for non-affine programs. An interesting direction is to investigate the static polyhedral techniques, coupled with dynamic techniques [START_REF] Rus | Sensitivity analysis for automatic parallelization on multi-cores[END_REF][START_REF] Mehrara | Parallelizing sequential applications on commodity hardware using a low-cost software transactional memory[END_REF] for program parallelization.

Appendix A

Application of Farkas lemma

P e 1 = {(i, j, i ′ , j ′ )|i ′ = i + 1 ∧ j ′ = j ∧ 1 ≤ i ≤ N -2 ∧ 0 ≤ j ≤ N -1}
The faces of the dependence polyhedron are:

       i -1 ≥ 0 j ≥ 0 N -2 -i ≥ 0 N -1 -j ≥ 0 (a S 2
1,1a S 1 1,1 )i + (a S 2 1,2a S 1 1,2 ) j + ω S 2 1 -ω S 1 1 + a S 2 1,1 -1 ≥ 0 Combining them with Farkas multipliers gives: λ 0 + λ 1 (i -1) + λ 2 ( j) + λ 3 (N -2i) + λ 4 (N -1j) λ 0 ≥ 0, λ 1 ≥ 0, λ 2 ≥ 0, λ 3 ≥ 0, λ 4 ≥ 0 Reorganizing the terms so that the induction variables and parameters could be equated to the legality constraint:

i(λ 1 -λ 3 ) + j(λ 2 -λ 4 ) + N(λ 3 + λ 4 ) + λ 0 -λ 1 -2λ 3 -λ 4

After equating the left hand and right hand sides of the Farkas lemma one obtains:

             a S 2 1,1 -a S 1 1,1 = (λ 1 -λ 3 ) a S 2 1,2 -a S 1 1,2 = (λ 2 -λ 4 ) 0 = (λ 3 + λ 4 ) ω S 2 1 -ω S 1 1 + a S 2 1,1 -1 = λ 0 -λ 1 -2λ 3 -λ 4 λ 0 ≥ 0, λ 1 ≥ 0, λ 2 ≥ 0, λ 3 ≥ 0, λ 4 ≥ 0
One can immediately notice:

λ 3 = 0 ∧ λ 4 = 0          a S 2 1,1 -a S 1 1,1 = (λ 1 ) ≥ 0 a S 2 1,2 -a S 1 1,2 = (λ 2 ) ≥ 0 ω S 2 1 -ω S 1 1 + a S 2 1,1 -1 = λ 0 -λ 1 λ 0 ≥ 0, λ 1 ≥ 0, λ 2 ≥ 0, λ 3 ≥ 0, λ 4 ≥ 0      a S 2
1,1a S

Definition 2 . 1 . 1 (Definition 2 . 1 . 2 (Definition 2 . 1 . 3 (Definition 2 . 1 . 4 (

 211212213214 Affine function). A n-dimensional function f : R m → R n is affine if and only if it can be expressed as follows: f (x) = Ax + b where b = (b 1 , . . . , b n ) T ∈ R n is a column vector and A ∈ R n×m is a matrix with n rows and m columns. Affine relation). A relation F ⊂ R m × R n is affine if and only if: ∀v ∈ R m , ∀w ∈ R n : (v, w) ∈ F ⇐⇒ w = Av + b where b = (b 1 , . . . , b n ) T ∈ R n is a column vector and A ∈ R n×m is a matrix with n rows and m columns. Affine hyperplane). For each vector a ∈ R n and scalar b ∈ R the set X of all vectors x ∈ R n such that: a T x = b defines an affine hyperplane. Affine half-space). For each vector a ∈ R n and scalar b ∈ R the set X of all vectors x ∈ R n satisfying the following affine inequality: a T x ≥ b defines a topologically closed affine halfspace. If we replace ≥ by > we will get a definition of topologically open affine halfspace.

Figure 3 . 2 -

 32 Figure 3.2 -Read-Write classification of the dependences

Figure 3

 3 Figure 3.3 -Matvect -matrix vector multiplication name condition affine constraints execution condition

Figure 3 Figure 3 . 7 -

 337 Figure 3.6 -Original scheduling matrices

Figure 3 Figure 3 .

 33 Figure 3.9 -Legal transformation matrices

2 .Figure 4 . 1 -

 241 Figure 4.1 -A gcc compilation flow

Figure 4 .

 4 2 shows the steps inside GRAPHITE framework: -SCoP outlining. The maximal subgraphs of CFG that have static control property and statically analysable affine data access patterns are extracted from the GIMPLE three-address code. -Mapping three-address code into the polyhedra. For the outlined SCoPs (static control parts) the polyhedral representation is reconstructed and mapped to the GIMPLE intermediate representation.

  (a) iteration count is symbolically determinable (but it might be unknown at compilation time if it depends on the global parameter) (b) loop bounds are constants or affine expressions of the outer loop induction variables and global parameters (c) loop is single exit loop 2. otherwise, all the statements inside a basic block are checked: (a) calls to functions with side effects are not allowed (pure and const function calls are allowed) (b) the only memory references that are allowed are accesses through arrays with affine subscript functions (c) the conditional statements are expressed as affine functions of induction variables and global parameters

Figure 4 . 5 -Figure 4 . 6 -Breaking SCoP 1 .

 45461 Figure 4.3 -CFG graph with region nesting

Figure 4 . 7 -

 47 Figure 4.7 -Matvect kernel -source code and internal representation

Figure 4 . 8 -

 48 Figure 4.8 -Relative nesting and ordering of loops and basic blocks

Figure 4 . 9 -

 49 Figure 4.9 -Polyhedral representation of three-address code

Figure 4 . 10 -

 410 Figure 4.10 -Polyhedral representation the origianl source code

Figure 4 .

 4 Figure 4.11 -Matvect kernel -cross bb dependences, inductions and reductions

Figure 4 .

 4 11: there is cross basic block dependence between a scalar write in B 4 in the statement D.3 = D.2 + pre.3 and its subsequent use in B 6 in the statement b[i_2] = D.4. An explicit single-element array crossu_bb_dep[0] is introduced to represent this dependence.

Figure 4 .

 4 Figure 4.12 -Explicit data dependence graph

Figure 4 . 13 -

 413 Figure 4.13 -Explicit data dependence graph

  5549_6]; D.5554_12 = D.5552_9 / D.5553_11; D.5555_13 = D.5548_5 -D.5554_12; X[i1_37][i2_38] = D.5555_13; D.5556_14 = B[i1_37][i2_38]; D.5557_17 = D.5551_8 * D.5551_8; D.5558_20 = D.5557_17 / D.5553_11; D.5559_21 = D.5556_14 -D.5558_20; B[i1_37][i2_38] = D.5559_21; i2_22 = i2_38 + 1; if (i2_22 != 500) goto <bb 3>; else goto <bb 5>;

Figure 4 .

 4 Figure 4.14 -source code and IR representation for ADI kernel

Figure 4 .

 4 Figure 4.15 -Tiling the iteration space

Figure 5

 5 Figure 5.2 -matrix multiplication -after PRE

  Figure 5.5 -GIMPLE code with CFG

Figure 5 . 7 -

 57 Figure 5.7 -Matmult Data Dependence Graph

Figure 5 .

 5 Figure 5.9 -legal -loop distribution with an expansion

Figure 5 .

 5 Figure 5.12 -A flow of polyhedral compilation based on violation analysis

Figure 6 . 1 -

 61 Figure 6.1 -Main loop kernel in Convolve

  r ( vh = 0 ; vh < K ; vh +=8) { v s [ 0 : 7 ] = { 0 , 0 , . . . , 0 } ; f o r ( j = 0 ; j < K ; j ++) { c = f i l t e r [ i ] [ j ] ; v f i l t e r [ 0 : 7 ] = { c , c , . . . , c } ; v s [ 0 : 7 ] += v f i l t e r [ 0 : 7 ] * image [ v+ i ] [ vh+ j : vh+ j + 7 ] } o u t [ v ] [ vh : vh + 7 ] += v s [ 0 : 7

Figure 6 . 2 -

 62 Figure 6.2 -Convolve Vectorization Examples

Figure 6

 6 Figure 6.3 -compilation flow

Figure 6 . 4 -Figure 6 . 5 -

 6465 Figure 6.4 -Cost model evaluation: comparison of predicted and actual impact of vectorization alternatives on the Cell SPU

Figure 6 . 6 -

 66 Figure 6.6 -Cost model evaluation: comparison of predicted and actual impact for convolve kernel on PPC970 and Cell SPU

Figure 7 . 1 -Figure 7 . 3 -Figure 7 . 4 -Figure 7 . 5 -

 71737475 Figure 7.1 -Uniform dependences and skewed hyperplane

4 .

 4 when comparing the two vectors a and b lexicographically, it is sufficient to find the first position l where a 0 = b 0 , . . . , a l-1 = b l-1 , a l < b l , while the ordering of the remaining positions is irrelevant.

Figure 7 . 6 -

 76 Figure 7.6 -Two statements, uniform dependences

Figure 7 .

 7 Figure 7.8 -n-ary decision tree representing the solution set L e 2

2 Figure 7 . 15 - 2 .Figure 7 . 16 -

 27152716 Figure 7.15 -The final result of PROJECT(L, X 1 )

1 0w w n n n n n n n n n n n n n n n 1 ⊥

 11

Figure 7 . 19 - 1 1 and x e 2 1

 71912 Figure 7.19 -The global solution set L 1 with decision variables x e 1 1 and x e 2 1

A S 1 Figure 7 . 21 -

 1721 Figure 7.20 -Full ranked A matrix

Figure 7 . 22 -

 722 Figure 7.22 -Complete solution obtained

8 .( 1 )( 1 . 1 )( 3 . 1 )( 3 . 2 )( 3 . 3 )

 8111313233 For the implementation convenience, the schedule matrices of those statements whose depth is less than dim max (S ) are padded with zeros Procedure COMPLETESCHEDULE INPUT: L, E, d L -a solution set found so far (for levels 1 ≤ l < d), satisfying the set of edges E E -a subset of edges satisfied at the given or deeper levels d -a scheduling dimension to be completed IF d = dim max (S ) { A complete multidimensional solution is obtained } provide schedule properties (1.2) save the complete solution L ′ (1.3) RETURN True (2) f ound ← False (3) WHILE (¬ f ound) IF empty (T A ← next solution from L) BREAK {compute orthogonal subspace of transformation } H ← ORTHOGONAL(T A ) {Compute the set of edges that must be strongly satisfied in order to give the non-empty solution: H ∩ L ′ = / 0} E must ← MUSTSATISFY(H) { Select the family X of possible edge satisfaction sets E S ∈ X for transformation T A } (3.4) X ← SELECT(L, T A )

)Figure 7 . 24 -

 724 Figure 7.24 -Schedule completion strategy

7. 7 n 111 7. 9 1 . 1 .

 7111911 -ary decision tree representing the solution set L e 1 . . . . . . . . . . . . . . . . . . . . 111 7.8 n-ary decision tree representing the solution set L e 2 . . . . . . . . . . . . . . . . . . . . The result of the operator EMBED(L e1 , {A S 0 1,• }) . . . . . . . . . . . . . . . . . . . . . . 112 7.10 set L ′ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 7.11 set L ′′ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 7.12 INTERSECT(L ′ , L ′′ ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 7.13 set L . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 7.14 Split tree after removing X 1 node . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 7.15 The final result of PROJECT(L, X 1 ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 7.16 Combining the solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 7.17 the solution set L e 1 with decision variable x e 1 . . . . . . . . . . . . . . . . . . . . . . 115 7.18 the solution set L e 2 with decision variable x e 2 . . . . . . . . . . . . . . . . . . . . . . 116 7.19 The global solution set L 1 with decision variables x e 1 1 and x e 2 1 . . . . . . . . . . . . . . . 116 7.20 Full ranked A matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 7.21 Non-full ranked A matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 7.22 Complete solution obtained . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 7.23 The global solution set L 1 with decision variables x e 1 1 and x e 2 1 . . . . . . . . . . . . . . . 118 7.24 Schedule completion strategy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121

Table 3 .

 3 1 -A summary of affine conditions for a dependence edge S 1 → S 2 (S 1 , (i) T ) and (S 2 , (i ′ , j ′ ) T ) respectively.

Table 3 .

 3 2 -A summary of dependences in matvect kernel at S 2 . Figure3.4 shows a full dependence graph for matvect kernel in Figure3.3. The dependences are summarized in Table3.2.

	S1	δ
		δ o
		δ o
	S2	δ
		δ -1
			Figure 3.4 -A RDDG for the matvect kernel
	source sink type
	S 1	S 2	true
	S 1	S 2	output
	S 2	S 2	true
	S 2	S 2	output
	S 2	S 2	antidependence

  do for example.

	f o r ( t t 2 = 0 ; t t 2 <= f l o o r d (N , 6 4 ) ; t t 2 ++) {
	f o r ( t 2 =max ( 6 4 * t t 2 , 0 ) ; t 2 <=min ( 6 4 * t t 2 +63 ,N ) ; t 2 ++) {
	S1 ( i = t 2 ) ;
	}
	}

Table 5 .

 5

	Lazy

1 -The number of the parallelized loops for Polybench

[START_REF]Polybench, the polyhedral benchmark suite[END_REF]

suite

end if end for end for end for 6.3.2 Experimental setup

		N 1	N 2	N 3	N 4	∆ 1	∆ 2	∆ 3	∆ 4
	interp_fp 512	16			1,2	1,0,2
	interp	512	16			1,2	1,0,2
	bkfir	512	32			1,0	1,1
	dct	8,8	8,8	8,8		8,0	0,1	1,8
	convolve	128	128	16	16	128, 1,	128, 1,
						0,	0,	16,	1,
						128	1	0	0
	H264	12,7	12,7			1	1
	dissolve	128	128			1	128
	alvinn	512,32 32,32			1,1	512,512
	MMM	16	16	16		16,0 0,1	1,16
	MMM T	16	16	16		16,0 0,1	1,16

min then cost min ← φ(x); x min ← x

Table 6 .

 6

1 -Benchmarks

  [START_REF] Barthou | Fuzzy array dataflow analysis[END_REF] -the solution set L e 2 with decision variable x e 2

			1
	1 w w p p p p p p p p p p p p p p p	A S 1 1,•	2

  Both edges strongly satisfied at level 1. Enables innermost DOALL loop. Does not require further scheduling dimension.

Table 7 .

 7

7 -Summary of the set of one-dimensional schedules

  1 1,1 ≥ 0 a S2 1,2a S 1 1,2 ≥ 0 ω S2 1 -ω S 1 1 + 2a S 2 1,1a S 1 1,1 -1 ≥ 0 5.9 legal -loop distribution with an expansion . . . . . . . . . . . . . . . . . . . . . . . . . 69 5.10 Read/Write instruction interleaving and variable live ranges . . . . . . . . . . . . . . . . 69 5.11 Speedups on Power platform. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 5.12 A flow of polyhedral compilation based on violation analysis . . . . . . . . . . . . . . . 75 6.1 Main loop kernel in Convolve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 6.2 Convolve Vectorization Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 6.3 compilation flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 6.4 Cost model evaluation: comparison of predicted and actual impact of vectorization alternatives on the Cell SPU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 6.5 Cost model evaluation: comparison of predicted and actual impact of vectorization alternatives on PPC970 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 6.6 Cost model evaluation: comparison of predicted and actual impact for convolve kernel on PPC970 and Cell SPU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 7.1 Uniform dependences and skewed hyperplane . . . . . . . . . . . . . . . . . . . . . . . 104 7.2 Uniform dependences and orthogonal hyperplane . . . . . . . . . . . . . . . . . . . . . 104 7.3 Uniform dependences and outermost parallel hyperplane . . . . . . . . . . . . . . . . . 104 7.4 Simple uniform dependence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 7.5 Two uniform dependences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 7.6 Two statements, uniform dependences . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

INTRODUCTION

A perfect loop nest is a set of nested loops if all statements are contained in an innermost nested loop. Otherwise a loop nest is called imperfect loop nest. In this work we consider a general class of imperfectly nested loops, assuming that perfectly nested loops are contained in this class.

Mainly used for presentation purposes

https:/gcc.gnu.org

so called 'latch'[START_REF] Fsf | GNU Compiler Collection (GCC) Internals Manual[END_REF] basic blocks that jump back to the loop header

If we want to emphasize the different domains of variable definitions, we will use i for the symbol names in IR and the mathematical notation i for the occurrence of the variable in the polyhedral model.

this happens in the PRE (Partial Redundancy Elimination) optimization pass that is scheduled to run before GRAPHITE

GRAPHITE FRAMEWORK based compilers, because of their rich semantical analysis that could be employed.

A parallelized loop is a loop whose iterations could be executed by independent threads, as in a DOALL[START_REF] Allen | Optimizing Compilers for Modern Architectures[END_REF] loops.

For this illustration, we use N=128 and K=16, due to current implementation restrictions.

b is typically not known at compilation time; nevertheless, we are only interested in its alignment modulo the VF, which is generally available.

in our case it is GCC compiler

The low 1.59x measured speedup for alternative

[START_REF] Allen | Automatic translation of fortran programs to vector form[END_REF].2b on the Cell SPU is due to an aliasing bug in GCC that results in bad scheduling. The out-of-order wide-issue (5 slots) PowerPC970 is less sensitive to this, but on the in-order 2-width-issue SPU performance drastically suffers as a result. The cost model obviously cannot (and should not) predict compiler bugs, however it can, as in this case, help reveal them.

The Γ part corresponds to parametric shift, allowing a scheduling function to be dependent on global parameters (usually the problem size).

Table 7.5 -Compressed encoding of the set L e 1 Representing convex rational and integer solution sets is a well studied problem in convex optimization theory [START_REF] Chernikova | Algorithm for finding a general formula for the non-negative solutions of a system of linear inequalities[END_REF][START_REF] Schrijver | Theory of linear and integer programming[END_REF].

Omega [START_REF] Pugh | The omega test: a fast and practical integer programming algorithm for dependence analysis[END_REF] library uses Fourier-Motzkin elimination extended to integer linear programming and integer set operations. Several efficient implementations based on dual representation [START_REF] Schrijver | Theory of linear and integer programming[END_REF] have been proposed. The most popular one was implemented in Polylib, the work of Le Verge and Wilde [START_REF] Wilde | A library for doing polyhedral operations[END_REF]. It is based on Chernikova's algorithm for transforming between constraint representation (the set of affine constraints) and generator representation(set generated by a set of vertices, rays and lines) [START_REF] Chernikova | Algorithm for finding a general formula for the non-negative solutions of a system of linear inequalities[END_REF][START_REF] Verge | A note on Chernikova's algorithm[END_REF].

The most recent developments are PPL [START_REF] Bagnara | The parma polyhedra library: Toward a complete set of numerical abstractions for the analysis and verification of hardware and software systems[END_REF] and ISL [START_REF] Verdoolaege | isl: An integer set library for the polyhedral model[END_REF] libraries for manipulating the convex sets of rational and integer solutions respectively.

Non-convex solution sets

As we have shown in the previous section, we split the solution set into the non-convex and convex part. For representing convex sets, we use techniques that are already well developed.

While the representation of the convex sets has sound mathematical and algorithmic background, the representation of the non-convex sets poses more complexity problems. Each non-convex problem is particular and thus requires dedicated techniques with heuristics well adapted to the problem.

Following is the work that describes our approach to representing the non-convex set of all the legal loop selections for a single dimension. In a Section 7.9 we will extend this to multidimensional solutions.

Representing non-convex solutions

We propose to use decision trees as a representation of non-convex solution sets. We have imposed a constraint on the possible set of coefficient values for the row Θ S k of A part of the scheduling matrix of statement S in Equation 7.9. A row is composed of all zeros except single integer constant 1 at some position s. Instead of representing the full row, we can only represent the position s, 1 ≤ s ≤ dim(S) of the integer constant 1 in a row.

Taking the solution set shown in Table 7.2, we encode the legal solution set as shown in Table 7.5.

The same is done for the solution set L e 2 , as shown in Table 7.6.

After encoding the full matrix row as single integer value, we construct the decision tree representing the solution sets. Even though one might use a simple representation of the sets, such as linked lists or arrays, we propose to use the n-ary decision trees for that purpose. The reason is the fact that we will require the set operations to be performed efficiently. A tree structure enables us to achieve this goal.

The leaf node ⊥ represents a 'no solution' -a given selection of the decision variables gives an illegal schedule. The rest of the leaf nodes represent a legal solutions and they are labelled by the polyhedron

x e

.17 -the solution set L e 1 with decision variable x e 1 1

We have shown in Subsection 7.6.1 the different heuristics considering this problem. Those heuristics try to find the best weak/strong satisfaction choice for each scheduling level at once. We want to take the different approach: we implicitly enumerate all those decisions in the decision tree of legal solutions.

For that purpose we introduce a decision variable x e d for each dependence edge e ∈ E at each scheduling level d. The meaning of the decision variable is the same as used in Feautrier's multidimensional scheduling [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part II: multidimensional time[END_REF]: the decision is made whether a given edge is weakly satisfied(x e d = 0) or strongly satisfied(x e d = 1) Taking those definitions into an account, we build a decision tree for the first scheduling dimension in the same way as we did in 7.8 while adding a decision variable x e d that selects the convex part of the solution to correspond to strict weak/strong satisfaction.

The original decision tree that unconditionally enforces strong satisfaction on the edge e 1 was shown in Figure 7.7. If we now split the dependence satisfaction decision into weak and strong satisfaction (variable x e 1 k ), we get the decision tree as shown in Figure 7.17. The leaf nodes are labelled by convex sets of the legal shifting factors giving the desired edge satisfaction condition (weak/strong). The same splitting is done for edge e 2 , as shown in Figure 7.18.

Combining solutions

In order to obtain the global solution for the first scheduling dimension, one has to intersect the two solution sets for edges e 1 and e 2 . But the decision trees for those two sets have two different decision variables: x e 1 k and x e 1 k respectively. Thus, the result of the intersection operation combines both of them (using embedding operator) into the final solution shown in Figure 7.19.

The global solution set L 1 is the exhaustive enumeration of all the possible one-dimensional schedules within the space defined in Subsection 7.7.2. The decision variables x e 1 1 and x e 2 1 partition the space of convex sets of legal shifting factors giving different edge satisfaction. Some interesting properties of the solutions contained in this set are summarized in Table 7.7.

Orthogonal completion

As shown in Table 7.7, some one-dimensional solution sets do not satisfy all the dependence edges at the first scheduling dimensions. What is more, in order to provide permutable loop band we would like to enforce the satisfaction of all edges at deeper levels, regardless of whether they have been satisfied at earlier levels or not.

A S 1 1,• 1 v v l l l l l l l l l l l l l l l l l l θ S 2 1 = j, θ S 2 2 = i This particular solution enables the outermost loop level parallelism (DOALL loop). This is only one of the possible solutions. The search method presented in the next Section enables the exploration of the full search space of multidimensional schedules.

Towards search methodology

In this section we will propose a search methodology based on decision tree set representation. The goal is to exhaustively explore the space of multidimensional affine transformations within the constraints defined in Subsection 7.7.2.

For that purpose, we will define some more operators on the decision trees and then we will propose the recursive procedure for building the sets of complete multidimensional schedules.

We will also explore the possibility of splitting the dependence graph G = (V, E) into strongly connected components and performing the search independently. This direction is strictly connected with loop fusion/distribution [START_REF] Kennedy | Loop distribution with arbitrary control flow[END_REF][START_REF] Roth | Loop fusion in high performance fortran[END_REF][START_REF] Darte | On the complexity of loop fusion[END_REF][START_REF] Darte | Temporary arrays for distribution of loops with control dependences[END_REF].

Just In Time Compilation

We have provided the polyhedral framework for the direct compilation of the three-address codes. Our technique was investigated in the context of the static compiler, translating the source code to binary machine code.

The fact that the polyhedral technique is implemented directly on the three-address-code representation could be extended and used in the JIT compilers that directly operate on the bytecode. Indeed, Polly project, based on LLVM compiler framework [START_REF] Grosser | Polly -polyhedral optimization in llvm[END_REF] goes in this direction. The crucial problem to be solved though is the scalability of the polyhedral techniques, especially the transformation part.

Going beyond the affine transformations

The affine transformations expressed in the polyhedral model cover a wide range of the classical loop transformations and enable coarse-grained parallelization techniques. Nevertheless, there are programs that could be expressed in the polyhedral model, but for which the coarse-grained parallelism could not be expressed through affine transformations [START_REF] Beletska | Coarse-grained loop parallelization: Iteration space slicing vs affine transformations[END_REF]. There are solutions that go beyond the affine transformations, and they are based on computing the transitive closure [START_REF] Bielecki | Calculating exact transitive closure for a normalized affine integer tuple relation[END_REF][START_REF] Beletska | Computing the transitive closure of a union of affine integer tuple relations[END_REF] of affine relations. Even though the dependence relations are affine, their transitive closure might be non-affine [START_REF] Bielecki | Calculating exact transitive closure for a normalized affine integer tuple relation[END_REF]. Unfortunately, the problem of code generation for non-linear forms is particularly hard and not-so-well studied problem [START_REF] Groesslinger | The Challenges of Non-linear Parameters and Variables in Automatic Loop Parallelisation[END_REF] for which solutions are still awaiting.

List of Figures