Les Fabian

Théo, Chukri, Manu, Florent Deux Nicolas

Julien Langou

Atanas Rountev

Evaluating the complexity of an algorithm is an important step when developing applications, as it impacts both its time and energy performance. Computational complexity, which is the number of arithmetic operations of a program, is easy to characterize for affine programs, that are composed of simple loops with static bounds. Data movement (or, I/O) complexity is more complex to evaluate as it refers, when considering all possible valid schedules, to the minimum required number of I/O between a slow memory with large capacity, and a fast storage location with limited capacity.

This thesis presents IOOpt, a fully automated tool that automatically computes symbolic bounds on the I/O complexity of an affine program. Given a description of an affine program, it automatically computes: 1. a lower bound of the I/O complexity as a symbolic expression of the cache size and program parameters; 2. an upper bound that allows one to assess the tightness of the lower bound; 3. a tiling recommendation (loop permutation and tile sizes) that matches the upper bound. This tool is implemented as an open-source program, and a demonstration is available at https://iocomplexity.corse.inria.fr/. This was made possible by combining advanced mathematical results with powerful program analysis tools such as polyhedral representation. Our data movement model is notably based on the one designed by Hong&Kung (red-blue pebble game), adapted to be able de decompose and recompose complex programs.

We demonstrate the effectiveness of our approach by deriving symbolic bounds for a wide range of affine programs. In many cases, the bounds found by IOOpt can be proven optimal, or close to optimal by a constant factor.

Je tiens en premier lieu à remercier mon directeur de thèse, Fabrice Rastello, qui m'a suivi et soutenu depuis mon stage de L3 en 2015, jusqu'à l'aboutissement de ce doctorat, et m'a donné l'opportunité de faire mon stage de M2 aux États-Unis, stage où ont été posées les bases de ce qui est devenu ma thèse. Au delà du cadre scientifique, c'est également une relation humaine de grande qualité qui s'est tissée au cours de ces années, devant le tableau comme en montagne.

Je remercie les rapporteurs, Laura Grigori et Sebastian Hack, d'avoir pris le temps de se plonger dans mon manuscrit de thèse, et pour leurs retours très positifs. Un grand merci également aux autres membres du jury, en particulier Sven Verdoolaege pour ses retours sur le manuscrit, ainsi que Laure Gonnord et Uday Bondhugula.

Résumé

L'évaluation de la complexité d'un algorithme est une étape importante lors du développement d'une application, de par son impact sur son temps d'exécution et sa consommation énergétique. La complexité arithmétique, qui est le nombre d'opérations effectuées par un programme, est en général facile à caractériser pour un programme affine, composé de boucles simples avec des bornes statiques. La complexité de mouvement de données (ou complexité I/O) est plus complexe à évaluer. En effet, elle se réfère à la quantité minimale de données à transférer entre une mémoire lente de grande capacité, et une mémoire rapide de capacité limitée, et ce en considérant tous les ordonnancements valides des opérations.

Dans cette thèse, nous présentons IOOpt, un outil automatisé capable de calculer des bornes symboliques sur la complexité I/O de programmes affines. Étant donné une description d'un programme affine, il génère : 1. une borne inférieure sur la complexité I/O sous la forme d'une expression symbolique dépendant de la taille de la mémoire rapide et des paramètres du programme ; 2. une borne supérieure qui permet d'évaluer la précision de la borne inférieure 3. une recommandation de tuilage (taille des tuiles et permutation des boucles) qui permet d'atteindre la borne supérieure. Cet outil est disponible en open-source, et peut être testé sur la page https://iocomplexity.corse.inria.fr/. Le développement de cet outil a été possible grâce à la combinaison de résultats mathématiques avancés avec des outils puissants d'analyse de programmes, tels la réprésentation polyédrique. Notre modèle de mouvement de données est notamment basé sur celui proposé par Hong&Kung (« Red-blue pebble game »), adapté pour pouvoir décomposer et recomposer des programmes complexes.

L'applicabilité de cette approche est démontrée en combinant des bornes symboliques sur un large panel de programmes affines. Dans de nombreux cas, les bornes trouvées par IOOpt sont optimales, ou proches de la borne optimale par un facteur constant.

Complexité I/O La complexité I/O est définie sur la base d'un modèle de mémoire à deux niveaux : une mémoire d'accès rapide mais de taille limitée, et une mémoire d'accès plus lent mais de capacité illimitée. Lors de l'exécution d'un programme, les données d'entrée se trouvent initialement dans la mémoire lente. Les données peuvent être déplacées entre les deux niveaux de mémoire, et doivent être dans la mémoire rapide pour être utilisées dans un calcul. Le coût I/O d'un programme, pour un ordre donné des opérations de calcul et de mémoire, est le nombre de transferts depuis la mémoire lente vers la mémoire rapide. La complexité I/O d'un programme est le coût I/O minimum qui peut être obtenu en réordonnant les opérations.

Il est en général impossible de calculer la complexité I/O exacte d'un programme car le nombre d'ordonnancements possible est combinatoire. Ceci motive la recherche de bornes inférieures et supérieures sur cette complexité.

Représentation des programmes Afin de pouvoir représenter les différents ordonnancements possibles, un programme est décrit par ses opérations, et par les dépendances de données entre celles-ci. Cela est naturellement représenté par une structure de graphe orienté, abrégé CDAG en anglais.

Un exemple de tel graphe pour un programme exécutant la multiplication de deux matrices carrées de taille N est présenté sur la figure 2.2b (page 6).

Bornes inférieures Nous utilisons deux approches distinctes pour trouver des bornes inférieures sur la complexité I/O. L'une est basée sur la technique dite du « S-partitionnement », et l'autre sur le front d'onde (wavefront) de graphes. La combinaison de ces deux approches permet de traiter des programmes avec différents motifs de dépendances de données, et de les combiner au sein d'un même programme.

Nous détaillons ici l'idée générale de l'approche par partitionnement, qui est de raisonner sur l'intensité opérationnelle maximale d'une région du CDAG, c'est-à-dire la quantité maximale d'opérations de calcul qui peuvent être effectuées avec un nombre limité et déterminé d'opérations d'I/O. Plus précisément, étant donné un paramètre T (dont la valeur sera fixé ultérieurement), on cherche à borner la taille du plus grand sous-graphe qui puisse être exécuté avec au plus T loads (transferts de la mémoire lente à la mémoire rapide), en supposant que la mémoire rapide contient déjà S données. Par un simple argument de comptage, une telle borne se traduit en une borne inférieure sur le nombre de sous-séquences d'au moins T loads nécessaires pour exécuter le programme complet, ce qui ce traduit enfin en une borne sur le coût I/O du programme.

Pour borner la taille du sous-graphe, il faut remarquer que sa « frontière », composée des prédécesseurs des noeuds du sous-graphe qui n'en font pas partie, est de taille bornée par S + T . Il reste à borner la taille d'un sous-graphe en fonction de la taille de sa frontière, et pour cela on utilise le fait que les graphes de programmes affines sont très réguliers et géométriques. Ainsi, les opérations d'un programme constitué de trois boucles parfaitement imbriquées comme la multiplication de matrices peuvent être naturellement représentées par des points au sein d'un cube. La dernière pièce du raisonnement est fournie par une famille d'inégalités mathématiques (l'inégalité de Loomis-Whitney et sa généralisation, l'inégalité de Brascamp-Lieb), qui permet de borner la taille d'un ensemble en fonction de la taille de ses projections en dimension inférieure (voir figure 2.4, page 8).

Pour la multiplication de matrices, la borne obtenue par cette méthode est 2N 3 √ S

. Cette borne est optimale : il est possible de trouver un ordonnancement des opérations qui atteint cette borne (en ignorant les termes négligeables).

Bornes supérieures La question qui se pose, une fois une borne inférieure établie, est de savoir si celle-ci est atteignable. Il s'agit pour cela d'exhiber un ordonnancement des opérations du programme le plus efficace possible en terme d'I/O. Autrement dit, une borne supérieure sur la complexité I/O d'un programme est simplement donnée par un ordonnancement valide. Ce qui rend ce problème difficile est qu'on cherche a obtenir des bornes paramétriques, valables pour toutes les valeurs des paramètres (dimensions des tableaux, capacité de la mémoire rapide), et qu'il est nécessaire de trouver un « bon » ordonnancement des opérations pour que la borne soit la plus proche possible de l'optimum.

La principale opération de réordonnancement qui permet d'optimiser les transferts de données est le tuilage. Il s'agit de partitionner le domaine d'itération d'un nid de boucles en sous-domaines appelés tuiles, qui sont ensuite exécutées les unes après les autres. Des exemples de codes tuilés pour la multiplication de matrices se trouvent page 10, figure 2.5.

Notre approche pour trouver des bornes inférieures peut se résumer comme suit. En premier lieu, le programme est analysé pour isoler les parties où la transformation de tuilage peut être appliquée. Une première étape sélectionne ensuite un ensemble réduit de permutations des boucles de tuilage. Ceci procure pour chaque permutation un programme tuilé avec des tailles de tuiles paramétriques. Ensuite, un algorithme calcule une expression symbolique du coût I/O pour chaque permutation, ainsi que des contraintes sur la taille des tuiles. Enfin, ces expressions peuvent être utilisées pour optimiser numériquement la taille des tuiles pour des valeurs données des paramètres, ou pour en déduire une borne symbolique qui ne dépend plus des tailles de tuile mais uniquement de S et des paramètres du programme, et peut donc être directement comparée aux bornes inférieures calculées précédemment.

Sur la multiplication de matrices, notre algorithme génère un tuilage de taille T i × T j × 1, pour un coût I/O de N 3 /T i +N 3 /T j +N 2 . La minimisation symbolique de cette expression sous la condition que toutes les données d'une tuile tiennent dans la mémoire rapide de taille S donnent une borne supérieure sur la complexité I/O de 2N 3 √ S+1-1

+ N 3 . Cela correspond bien à la borne inférieure, à des termes asymptotiquement négligeables près.

Chapter 1

Introduction

Several factors can impact the execution time of an application on a computer: the most obvious one is the number of arithmetic operations, but the cost of moving data between CPU registers and memory can also be critical. For several decades now, CPUs have seen their computational power (FLOPS) improve at a higher rate than memory latencies and bandwidth. This means that, in order to use a processor to its full capacity, the minimization of data transfers between registers, CPU caches and RAM is getting increasingly important.

While latency costs can usually be hidden by hardware mechanisms like prefetching, the mismatch between computational power and memory bandwidth poses a fundamental limit. The balance between computation and data transfers can be expressed through a program's operational intensity, which is the ratio between the number of arithmetic operations and the volume of data transfers. This quantity can be compared with the machine balance of a CPU, which is the ratio between peak arithmetic performance and memory bandwidth. Only if a given program has an operational intensity that exceeds the machine balance can it exploit the full potential of the CPU; otherwise its execution speed will be limited by data transfers.

This motivates the need to study the data movement complexity of programs, also called I/O complexity. Similarly to computational complexity, which characterizes the number of operations made by a program, I/O complexity quantifies the volume of data movement operations to and from memory that are needed to perform a computation, for any valid schedule of operations. This quantity can be very challenging to compute precisely, and has been the subject of a long line of work, starting with the seminal paper of Hong & Kung [START_REF] Hong | I/O complexity: The Red-Blue Pebble Game[END_REF] that first presented a method for deriving lower bounds on the I/O of programs. Concurrently, advanced methods in program analysis such as polyhedral compilation have been developed, enabling the analysis and representation of data dependencies of complex programs (precisely affine programs, that are composed of a combination of loop nests accessing multidimensional arrays, where loop bounds and array accesses are affine function of loop counters and program parameters). Elango et al. [START_REF] Elango | On Characterizing the Data Access Complexity of Programs[END_REF] first demonstrated how these two lines of work could be combined to automate the data movement analysis of affine programs, including stencils and linear algebra kernels.

In this thesis, we present a fully automated method for computing both lower and upper bounds on the I/O complexity of affine programs. It has been implemented as two opensource tools called IOLB and IOUB.

IOLB (for I/O Lower Bounds) takes as input an affine program and outputs a parametric lower bound with scaling constant on its I/O. This means that the bound is a symbolic expres-sion, depending on program parameters (matrix dimensions for example) and on the size of the fast memory S, and that it does not rely on asymptotic notation. For instance, for the basic multiplication of square matrices of dimensions N ×N , the bound output by IOLB is 2N 3 √ S -3S. As such, this bound can be evaluated for any value of N and S, and compared for example to a cache miss count for a specific execution. Our algorithm for lower bounds is able to decompose a program into sub-regions and to combine different methods on each sub-program to derive a bound as tight as possible. IOLB was extensively evaluated on 30 algorithms specified in PolyBench [START_REF] Pouchet | PolyBench/C 4.2[END_REF], with several first-time I/O lower bounds demonstrations on these programs.

Similarly, IOUB (for I/O Upper Bounds) automatically computes a symbolic expression of an upper bound on the data movement cost of a program, from an abstract representation of its iteration space and array accesses. To do so, it generates a set of candidate tiling transformations, using some pruning techniques to avoid combinatorial explosion, and associates to each of them a symbolic expression of the I/O cost (which is an upper bound in the general case but is actually exact in many cases). It then compares these expressions using numerical evaluation to single out the best ones. For concrete values of fast memory size and program parameters, the tool can also suggest a loop permutation and numerical tile sizes that minimize I/O. Finally, it is able to use computer algebra to determine tile sizes as a function of the fast memory size S, and to get an upper bound that is directly comparable to the lower bound expression from IOLB. For instance, on square matrix-matrix multiplication, the upper bound from IOUB is 2N 3

√ S+1-1
+ 1, which matches the lower bound up to lower order terms, and effectively proves that the I/O complexity for this problem is 2 N 3 √ S

(1 + o( 1)). IOUB was evaluated on several tensor contraction and convolution programs, giving optimal bounds in many cases.

The rest of this manuscript is organized as follows. Chapter 2 gives a high-level view of the concepts and methods used to derive lower and upper bounds on I/O complexity, using examples to provide intuition to the reader. Chapter 3 delves into the details of automated I/O lower bound derivation, starting with the formal framework needed to reason about them, then presenting the two main methods (K-partitioning and wavefront), as well as how to combine them. We provide detailed algorithms for the automation of I/O lower bound analysis, and an evaluation of our implementation on a large selection of affine programs. Chapter 4 similarly presents the I/O upper bound analysis method, first describing the model and theoretical framework, then explaining how to build an automated program analysis, and showing its application on tensor contraction and convolution kernels. Finally, conclusions and suggestions for future work are provided in Chapter 5.

Chapter 2

Overview and Background

In this chapter, we give a high-level view of the concept of I/O complexity, and illustrate the techniques for deriving lower and upper bounds on an example. The goal is to provide intuitions on the representations and reasoning methods that will be more thoroughly detailed in Chapters 3 and 4.

I/O Complexity

Our framework is a simple memory model with two levels: a small and fast memory (of fixed size S) and a slow but infinite memory. Initially, all program data resides in the slow memory. Data can be transferred back and forth between the two memory levels, and a piece of data has to be in the fast memory in order to be used in a computation. The I/O cost of a program, for a given schedule of arithmetic and I/O operations, is the number of data transfers from slow to fast memory (or loads) 1 . Now, for a given program, we are interested in the lowest I/O cost that can be achieved by reordering operations: this quantity is called the I/O complexity of a program.

Computing the I/O complexity exactly is usually intractable, as one would have to simulate every possible valid schedule of operations. Hence the need to find lower and upper bounds on this complexity. To get an precise estimation of the I/O complexity of a program, 1. We choose in this work to count only loads (transfers from slow to fast memory) and not stores (transfers from fast to slow memory). This simplifies the reasoning while still providing meaningful values.

CPU

Fast memory capacity S Slow memory unbounded capacity the lower and upper bounds should be as close as possible. In some cases (such as matrix multiplication that will be used in this section), it is possible to find matching lower and upper bounds (up to lower order additive terms), but even for simple programs this requires complex reasoning. It is important to note that the bounds that we derive are valid for a particular implementation of an algorithm, and not for any program computing a given problem.

I/O cost Figure 2.1 -Memory model f o r ( i = 0 ; i < N ; i ++) f o r ( j = 0 ; j < N ; j ++) f o r ( k = 0 ; k < N ; k ++) S : C [ i ] [ j ] += A [ i ] [ k ] * B [ k ] [ j ] ; (

Program representation

Since we need to represent whether a given schedule of operations is valid or not, a program is represented by its operations, and data dependencies between them. This is naturally abstracted as a directed graph, called a CDAG (computational directed acyclic graph, see Def. 1). For programs that are composed of nested loops, vertices of this graph can be represented as points in a multi-dimensional space.

Example: Matrix Multiplication Let us take the example of matrix-matrix multiplication. The corresponding C code and the graph representation are shown on Fig. 2.2. On Fig. 2.2b, vertices in green, red and light gray respectively represent input data from arrays A, B and C. Vertices in dark gray are instances of statement S (i.e., each vertex corresponds to a multiplication and an addition). Edges represent data dependencies: all direct predecessors of a vertex must be in fast memory to perform the computation.

The arithmetic complexity of this program (counting the multiplication and accumulation as a single operation) is C = N 3 .

Considering the original schedule given by Fig. 2.2a, and an optimal memory management (evicting the piece of data that will be needed latest from the small memory when it is full), and supposing the fast memory size S to be much smaller than N 2 , the I/O cost for this schedule is on the order of N 3 -intuitively, the condition S << N 2 means that all data (matrices A, B and C) does not fit in the fast memory. Indeed, to compute a single row of C, the totality of matrix B needs to be loaded in memory, so at most S -1 words of matrix B can stay in fast memory between two consecutive iterations of the outer loop on i, and each iteration demands at least N 2 -S + 1 load operations. The outer loop iterates N times, so the total number of operations when S << N 2 is indeed of the order of N 3 . The question is to know whether this can be improved, and up to which point.

Lower bound derivation

We employ two distinct approaches to find lower bounds on data movement: one based on the S-Partitioning approach [START_REF] Hong | I/O complexity: The Red-Blue Pebble Game[END_REF], and one based on graph wavefronts [START_REF] Elango | On characterizing the data movement complexity of computational DAGs for parallel execution[END_REF]. Combining these two approaches is essential for handling of a large class of programs, as they are complementary and work on different data dependence patterns. A lower bound for a program exhibiting a combination of both kinds of patterns can combine results from both approaches. In this introductory section, we only present a high-level overview of the S-Partitioning approach, and apply it to matrix multiplication, to familiarize the reader with the reasoning and terminology used.

The general idea of the S-partitioning approach for proving lower bounds is to reason on the maximal operational intensity of a sub-CDAG, i.e. the maximum number of computations that can be done with a fixed number of memory operations (loads). More precisely, for a given parameter T (whose value will be chosen later), we are looking for an upper bound on the size of the largest subgraph that can be executed doing only T loads, and supposing the fast memory is full, i.e. S pieces of data are already loaded. Supposing we know such a bound U (T ), this implies that any valid schedule must contain at least ⌊|V |/U (T )⌋ sub-sequences with T loads, where V is the set of computational vertices in the CDAG. This leads to a lower bound on the number of loads of

Q low = T • |V | U (T ) . ( 2

.1)

Let us now explain how the upper bound U (T ) can be computed. The idea is to use information on the "boundary" of a set of vertices P that can be executed with T load instructions, to derive a bound on the cardinality of P . Intuitively, the "boundary" corresponds to values that need to be in fast memory to execute vertices in P , so its size is bounded. When the dependence pattern in the CDAG is regular, geometric inequalities that relate the cardinality of a set of points in a multi-dimensional space to cardinalities of lower-dimensional projections of those points can be used to derive a bound on |P |.

More formally, the In-set In (P ) of P is the set of all predecessors of the vertices in P that do not belong to P (see Fig. 2.3a). In (P ) represents values that are dependencies of operations in P , but are not computed within P . Since all values in In (P ) must be in fast memory in order to execute the operations in P , they must either already be in fast memory at the beginning of the execution, or be explicitly loaded. At most S values from In (P ) can be initially present in the fast memory, and by definition T values can be loaded to execute P . Thus the size of In (P ) must be less than (S + T ).

In many programs, CDAGs have a natural embedding in a multi-dimensional space, where each dimension corresponds to a loop counter, and dependences between statements are regular and orthogonal. When this is the case, the size of the In-set In (P ) can be bounded by the size of its projections on hyperplanes. Fig. 2.3b illustrates this relation in a 2-dimensional setting. This relation is one of the key insights to the automation of lower bound computation.

Let us go back to matrix multiplication. In this example, vertices corresponding to statement S are naturally represented as points in a three-dimensional lattice, each dimension P In(P) ). With that representation, we observe that the size of the In-set of a vertex set of this particular graph must be greater than or equal to the cardinality of the orthogonal projections of P onto the three elementary planes:

|In (P )| ≥ proj x (P ) + proj y (P ) + proj z (P ) . (2.2)
Conveniently, a mathematical result known as the Loomis-Whitney inequality [START_REF] Loomis | An inequality related to the isoperimetric inequality[END_REF] relates the size of a vertex set in a three-dimensional space to the sizes of its three 2D projections (see Fig. 2.4):

|P | ≤ proj x (P ) • proj y (P ) • proj z (P ) 3 2 . (2.3)
This result can be generalized to arbitrary dimensions and any set of (not necessarily orthogonal) projections, and is called the Brascamp-Lieb inequality.

Optimizing the conjunction of inequalities (2.2) and (2.3) leads to the following bound on |P |: (mathematical details can be found in Sec. 3.3.1)

U (T ) = (S + T ) 3 3 2 . (2.4)
Finally, setting T = 2S and plugging the expression into equation (2.1) gives the optimal bound

IO mm ≥ Q low = 2S • N 3 S 3 2 ≈ 2N 3 √ S . ( 2.5)
A key contribution of this work is the automation of this process.

Upper bound derivation

Once we found a lower bound, the question is whether it can be attained. An upper bound on I/O complexity is simply the cost of a valid program schedule.

Thus, for a program with fixed parameter values, it is possible to generate a full execution trace, and count the number of data transfers, using an optimal cache policy (always evicting the piece of data that will be used the latest). However, this is costly and only valid for fixed parameter values, while we would like to derive parametric expressions. Moreover, the goal is to find an upper bound as low as possible, so we need to use a "good" schedule.

The main rescheduling transformation to optimize data transfers is tiling. It consists in partitioning the complete iteration domain of a loop nest into smaller parts, usually (hyper-)rectangles, called tiles, and executing each tile atomically. Fig. 2.5 shows examples of tiled codes for matrix multiplication.

Our approach to finding upper bounds on I/O complexity can be summarized as follows. Beforehand, the program is analyzed to find the part on which the tiling transformation can be applied. A first step selects a subset of tiling loop permutations, using analytical reasoning to reduce the search space. For each permutation, this provides a parametrically tiled program. Then, an algorithm based on polyhedral calculus computes a symbolic expression of the I/O cost for each tiling permutation, as well as constraints on tile sizes. Finally, these expressions can be used to provide actual tile sizes for numerical parameter values, or to provide a symbolic bound depending only on S but not on tile sizes, that can be directly compared with the lower bound.

f o r ( i 1 = 0 ; i 1 < N ; i 1 += T i ) f o r ( j 1 = 0 ; j 1 < N ; j 1 += T j ) f o r ( k1 = 0 ; k1 < N ; k += Tk ) f o r ( i = i 1 ; i < i 1 + T i ; i ++) f o r ( j = 0 ; j < j 1 + T j ; j ++) f o r ( k = 0 ; k < k1 + Tk ; k ++) C [ i ] [ j ] += A [ i ] [ k ] * B [ k ] [ j ] ;
(a)

f o r ( j 1 = 0 ; j 1 < N ; j 1 += T j ) f o r ( i 1 = 0 ; i 1 < N ; i 1 += T i ) f o r ( k1 = 0 ; k1 < N ; k += Tk ) f o r ( k = 0 ; k < k1 + T i ; k ++) f o r ( i = i 1 ; i < i 1 + T j ; i ++) f o r ( j = 0 ; j < j 1 + Tk ; j ++) C [ i ] [ j ] += A [ i ] [ k ] * B [ k ] [ j ] ; (b)
Figure 2.5 -Two valid tiled schedules for Matrix-matrix multiplication.

f o r ( i 1 = 0 ; i 1 < N ; i 1 += T i ) f o r ( j 1 = 0 ; j 1 < N ; j 1 += T j ) f o r ( k1 = 0 ; k1 < N ; k += Tk ) f o r ( i = i 1 ; i < i 1 + T i ; i ++) f o r ( j = 0 ; j < j 1 + T j ; j ++) f o r ( k = 0 ; k < k1 + Tk ; k ++) C [ i ] [ j ] += A [ i ] [ k ] * B [ k ] [ j ] ; A B C N 2 N 2 N 2 Ti • N N 2 Ti • N Ti • N Tj • N Ti • Tj Ti • Tk Tj • Tk Ti • Tj Tk Tk • Tj Tj Tk Tk 1 1 1 1 Figure 2.

-Array footprints

Let us illustrate it on the matrix multiplication example. A tiled schedule for matrix multiplication consists of six loops, and the permutation of the three outer loops and the three inner loops is free. Figure 2.5 shows two valid tiled schedules (supposing for clarity that Ti, Tj and Tk are divisors of N). We use a heuristic to select a subset of relevant permutations, which is detailed in Sec. 4.3. For now, let us focus on code (a) in Fig. 2.5.

Figure 2.6 shows the symbolic expression of the memory footprint of each array for each subloopnest. For each array, the last level for which the footprint is supposed to fit in the fast memory is highlighted. This means that, for some fixed values of i1 and j1, at the beginning of the third loop iterating of k1, all Ti•Tj values C[i][j] for i = i1...i1 + Ti, j = j1...j1 + Tj are loaded and kept in fast memory for the whole iteration. The same goes at the beginning of the fourth loop (iterating on i), for values

A[i][k] and B[k][j], i = i1...i1 + Ti, j = j1...j1 + Tj, k = k1...k1 + Tk.
This implies that all these pieces of data fit in the fast memory, i.e.

Ti • Tk + Tj • Tk + Ti • Tj ≤ S (2.6)
The outer loops iterating on i1, j1 and k1 iterate respectively N/Ti, N/Tj and N/Tk times. In

f o r ( i 1 = 0 ; i 1 < N ; i 1 += T ) f o r ( j 1 = 0 ; j 1 < N ; j 1 += T ) f o r ( k = 0 ; k < N ; k ++) f o r ( i = i 1 ; i < i 1 + T ; i ++) f o r ( j = 0 ; j < j 1 + T ; j ++) C [ i ] [ j ] += A [ i ] [ k ] * B [ k ] [ j ] ;
IO = N 3 Ti + N 3 Tj + N 2 (2.7)
Finally, to be able to compare this expression with the lower bound, Ti and Tj need to be expressed as functions of S and N. As the goal is to minimize I/O, this amounts to finding values of Ti, Tj and Tk that minimize expression (2.7) under condition (2.6). This can be solved manually or using computer algebra, and the optimal value is Ti = Tj = √ S + 1-1, and Tk = 1, leading to

IO = 2N 3 √ S + 1 -1 + N 2 . ( 2.8)
This expression matches the lower bound up to lower order terms. Since the optimal value of T k is 1, one loop of the tiled code can be removed, and the final code attaining this bound is shown on Fig. 2.7.

Our tool, IOOpt, is also able to find a numerical solution for fixed values of N and S. For instance, for N = 1500 and S = 1024, IO is minimized for Ti = Tj = 31, and gives a total number of loads IO = 65516129.

On an actual machine, several levels of memory are present, so multiple levels of tiling are needed to produce efficient code. IOOpt is able to model such multi-level memory hierarchies and optimize the I/O for all levels at once. The result can then be used in conjunction with additional optimization techniques (vectorization, efficient register use...) to produce highly efficient code for modern CPUs.

Affine programs and automation

This chapter gave a high-level view of the models and methods to produce I/O bounds, but did not detail so far the automation process, which is a key part of this thesis.

To do so, we have to choose a class of programs to work with: the class of affine programs stands in a nice spot combining expressiveness and applicability of our techniques, and comes with a convenient analysis tool ecosystem. Affine programs cover a wide set of key algorithms, as exemplified with the 30 algorithms in PolyBench/C [START_REF] Pouchet | PolyBench/C 4.2[END_REF] that span popular dense linear algebra, stencils/convolutions, and dynamic programming techniques.

An affine program is composed of (potentially nested) for loops and statements, operating on variables and multi-dimensional arrays. Loop bounds are affine expressions of surrounding loop indices and program parameters (a program parameter is a symbolic constant during the compilation, like array dimensions). Similarly, access functions made by statements to arrays are affine expressions of surrounding loop indices and program parameters. As such, the control flow of such programs is statically analyzable.

This class corresponds to the class of programs that can be analyzed by polyhedral compilation techniques [START_REF] Feautrier | Polyhedron model[END_REF]. Powerful analysis tools such as ISL [START_REF] Verdoolaege | ISL: An Integer Set Library for the Polyhedral Model[END_REF] and barvinok [START_REF] Verdoolaege | Counting Integer Points in Parametric Polytopes Using Barvinok's Rational Functions[END_REF] can represent these programs abstractly as points in a multidimensional space, and do some operations on these points such as analyzing data flow dependences, counting the number of points satisfying some property as a symbolic expression... This is particularly interesting for our I/O analysis, as it allows us for instance to uncover regular data dependencies that correspond to projections in our lower bound derivation.

For example, in the matrix multiplication code of figure 2.2a, all iterations of statement S are represented as follows in ISL notation:

[N ] → {S[i, j, k] : 0 ≤ i < N ∧ 0 ≤ j < N ∧ 0 ≤ k < N } (2.9)
Here the matrix dimension N is a program parameter (a symbolic constant), and every instance of statement S is represented as an integer point in a 3-dimensional space.

In the same way, data dependencies can be represented parametrically. The following ISL formula represents all dependencies from data in array A:

[N ] → {A[i, k] → S[i, j, k] : 0 ≤ i < N ∧ 0 ≤ j < N ∧ 0 ≤ k < N } (2.10)
The corresponding vertices and edges are highlighted on the CDAG from Fig.

Lower Bounds

In this chapter, we present our results and tools for the derivation of I/O lower bounds. Formal definitions and theoretical background needed for the main analysis are introduced in Section 3.1. Then, Section 3.2 provides insights on how complex programs can be decomposed to derive tighter bounds. Two proof techniques, namely the K-partition and the wavefront based proofs are respectively described in Section 3.3 and Section 3.4. An overview of the complete framework is provided in Sec. 3.5. We demonstrate the power of our approach by running it on a full benchmark suite of affine programs: Sec. 3.6 reports the data movement complexities for the 30 algorithms benchmarked in PolyBench, as well as tensor contraction and convolution programs. Finally, related work is discussed in Sec. 3.7.

Foundations

Before delving into the main matter, we present some background and discuss prior results needed for the developments in this chapter.

CDAG

The formalism and methodology we use is strongly inspired by the foundational work of Hong & Kung [START_REF] Hong | I/O complexity: The Red-Blue Pebble Game[END_REF]. It was introduced in the context of data movement lower bounds and models the execution of an algorithm on a processor with a two-level memory hierarchy.

In this formalism, an algorithm is abstracted by a graph -called a CDAG -, where vertices model execution instances of arithmetic operations and edges model data dependencies among the operations. The data movement (or I/O) complexity of a CDAG is formalized via the red-white pebble game (a variation of Hong & Kung's red-blue pebble game). In this game, a vertex of a CDAG can hold red and white pebbles. Red pebbles represent values in the fast memory (typically a cache or scratchpad), and their total number is limited. White pebbles represent computed values, that can be loaded into the fast memory. A value can be computed only when all its operands reside in the fast memory: a red pebble can be placed on a vertex in the CDAG if all its predecessors hold a red pebble, a white pebble is placed alongside the red one. Values that have been computed can be loaded in and discarded from the fast memory at any time: a red pebble can be placed or removed from a vertex holding a white pebble. However a value can only be computed once: once a vertex holds a white Parameters: N, M;

Input: A[N], C[M]; Output: A[N]; f o r (t=0;t<M;t++) f o r (i=0; i<N; i++) A[i] = A[i] * C[t]; (a)
Parameters: N, M; pebble, it cannot be removed. The I/O cost of an execution of the game is the number of loads into the fast memory: the number of times a red pebble is placed alongside a white one. Contrary to Hong & Kung's original model, this formalism does not allow recomputation of the value at a vertex. This follows many previous efforts [4,3,[START_REF] Bilardi | A characterization of temporal locality and its portability across memory hierarchies[END_REF][START_REF] Bilardi | A Lower Bound Technique for Communication on BSP with Application to the FFT[END_REF][START_REF] Demmel | Communicationoptimal Parallel and Sequential QR and LU Factorizations[END_REF][START_REF] Elango | On characterizing the data movement complexity of computational DAGs for parallel execution[END_REF][START_REF] Elango | On Characterizing the Data Access Complexity of Programs[END_REF][START_REF] Irony | Communication Lower Bounds for Distributed-Memory Matrix Multiplication[END_REF][START_REF] Savage | Extending the Hong-Kung Model to Memory Hierarchies[END_REF]. This assumption is necessary to be able to derive lower bounds for complex CDAGs by decomposing them into subregions. Another slight difference with prior work is that it only models loads and not stores -this means the generated bounds are clearly also valid lower bounds for a model that counts both loads and stores. Since the number of loads dominates stores for most computations, the tightness of the lower bounds is not significantly affected.

Input: A[N], C[M]; Output: S M-1 [N]; f o r (0 ≤ t < M and 0 ≤ i < N ) i f (t==0): S 0,i = A[i] * C[0]; e l s e : S t,i = S t-1,i * C[t]; (b) t i (c)
Let us define formally the notion of CDAG, and our pebbling game.

Definition 1 (Computational Directed Acyclic Graph (CDAG)). A Computational Directed Acyclic Graph (CDAG) is a tuple G = (V , E, I) of finite sets such that (V , E) is a directed acyclic graph, I ⊆ V is called the input set and every v ∈ I has no incoming edges.

Definition 2 (Red-White Pebble Game). Given a CDAG G = (V , E, I), we define a complete Sred-white pebble game (S-RW game for short) as follows: In the initial state, there is a white pebble on every input vertex v ∈ I, S red pebbles and an unlimited number of white pebbles. Starting from this state, a complete game is a sequence of steps using the following rules, resulting in a final state with white pebbles on every vertex.

(R1) A red pebble may be placed on any vertex that has a white pebble.

(R2) If a vertex v does not have a white pebble and all its immediate predecessors have red pebbles on them, a red pebble may be placed on v. A white pebble is placed alongside the red pebble. (R3) A red pebble may be removed from any vertex.

The cost of a S-RW game is the number of applications of rule (R1), corresponding to the number of transfers from slow to fast memory. A CDAG (see Fig. 3.1c) represents a single dynamic execution of a program, and can be very large. To be able to analyze programs of realistic size with reasonable resources, we use a compressed representation called Data-flow graph (DFG). Another advantage of such a representation is that it is parametric, i.e. a single DFG can represent CDAGs of different sizes, depending on program parameters. A DFG represents an affine computation, which is the class of programs that can be handled by the polyhedral model [START_REF] Feautrier | Polyhedron model[END_REF]. We use terminology and syntax inspired on the ISL library [START_REF] Verdoolaege | ISL: An Integer Set Library for the Polyhedral Model[END_REF], and illustrate them with the example of Fig. 3.1. Formal definitions can be found in the tutorial [START_REF] Verdoolaege | Presburger Formulas and Polyhedral Compilation[END_REF].

D A = [N ] → {A[i] : 0 ≤ i < N } D C = [N ] → {C[t] : 0 ≤ t < M} D S = [M, N ] → {S[t, i] : 0 ≤ t < M ∧ 0 ≤ i < N } |D S | = MN (c) Node domains R e 1 = [N ] → {A[i] → S[0, i] : 1 ≤ i < N } R e 2 = [M, N ] → {C[t] → S[t, i] : 0 ≤ t < M ∧ 0 ≤ i < N } R e 3 = [M, N ] → {S[t, i] → S[t + 1, i] : 0 ≤ t < M -1 ∧ 0 ≤ i < N } (d) Edge relations
Vertex domains As one can see on Fig. 3.1c, to each loop is associated a "geometric" space dimension (t and i here) so that each vertex of the CDAG lives in a multidimensional iteration space, its domain, that can be algebraically represented as a union of parametric Z-polyhedra bounded by affine inequalities. Definition 4. A domain is a Z-polyhedron, i.e. a set of points in the lattice Z d for some integer d, bounded by affine inequalities. In ISL notation, it is called a set and is written:

S = [p 1 , . . . , p m ] → {X[i 1 , . .

. , i d ] : P }

where:

p 1 , . . . , p m are (integer) parameter variables, -X is an identifier for the space to which points in S belong (two points with the same coordinates but different identifiers are treated as distinct points), i 1 , . . . , i d are coordinate variables, -P is a Presburger formula on parameter and coordinate variables (in most cases a conjunction of affine inequalities I 1 ∧ . . . ∧ I t ).

It is also often useful to manipulate unions of such domains, and ISL provides appropriate objects for this. Standard operations (union, intersection, difference,. . . ) are available, as well as a cardinality operation (denoted |D|), which computes the cardinality of a given domain as a (quasi-)polynomial function of the parameter variables (this is provided by the barvinok library [START_REF] Barvinok | A Polynomial Time Algorithm for Counting Integral Points in Polyhedra when the Dimension is Fixed[END_REF][START_REF] Verdoolaege | Counting Integer Points in Parametric Polytopes Using Barvinok's Rational Functions[END_REF]). As an example (see Fig. 3.2c), the domain D S of statement S is a Z-polyhedron with parameters M and N made up of all integer points (t, i) such that 0 ≤ t < M and 0 ≤ i < N . The number of points in this set (cardinality) is |D S | = MN . Note that the space within which all the points of a statement (S here) live is identified with the name of the statement, using the notation S[t, i].

Edge relations A set of edges of the CDAG is represented using a relation (ISL map), which is a set of pairs between two spaces, from the domain space to the image space.

Definition 5. A relation is a binary relation between points in two domains, subject to affine constraints. It is called a map in ISL, and written as follows:

R = [p 1 , . . . , p m ] → {X[i 1 , . . . , i d ] → Y [j 1 , . . . , j ′ d ] : I 1 ∧ • • • ∧ I t } The domain Dom (R) of R is the set of points X[i 1 , . . . , i d ] having an image through R. The image Im (R) of R is the set of points Y [j 1 , . . . , j ′ d ] that have a predecessor through R. As an example (see Fig. 3.2d), the data flow from statement S[t, i] (definition of A[i] in S) to statement S[t + 1, i] (use of A[i] in S) is represented using the relation R e 3 .
In addition to standard set operations, ISL can compute (an over-approximation of) the transitive closure of a relation, denoted R * . The following operations are also supported: image of a domain D through a relation R (denoted R(D)), and composition of two relations R 1 and R 2 , denoted R 1 •R 2 (this is left composition, going the opposite way from usual functional notation). Composition restricts the image space of the resulting relation to points where the composition relation makes sense:

Dom (R 1 • R 2 ) = R -1 1 (Im (R 1 ) ∩ Dom (R 2 )), Im (R 1 • R 2 ) = R 2 (Im (R 1 ) ∩ Dom (R 2 )
). As with domains, we will sometimes manipulate unions of such relations.

A Data-flow graph (DFG) A DFG is a graph G = (S, D). Each vertex S ∈ S of the graph represents a (static) statement or an input array of the program. Each vertex S is associated with a parametric iteration domain D S and a list of enclosing loops (empty for input arrays). Each edge d = (S a , S b ) ∈ D represents a flow dependency between statements or input arrays. Each edge is associated with an affine relation R d between the coordinates of the source and sink vertices. The DFG is a compact (exact) representation of the dynamic CDAG where a single vertex/edge of the DFG represents several vertices/edges of the dynamic CDAG. While all the reasoning and proofs can be done by visualizing a CDAG, the actual heuristic described in this chapter manipulates its compact representation, allowing us to translate graph methods [START_REF] Elango | On Characterizing the Data Access Complexity of Programs[END_REF] into geometric reasoning. Fig. 3.2b shows the DFG for our simple stencil code.

Partitioning

One key idea from Hong & Kung was the design of a mapping between any valid sequence of moves in the red-blue pebble game and a convex partitioning of the vertices of a CDAG and thereby the assertion of an I/O lower bound for any valid schedule in terms of the minimum possible count of the disjoint vertex-sets in any valid 2S-partition (see below) of the CDAG.

The argument is the following: any execution can be decomposed into consecutive segments doing exactly (but for the last one) S loads. There are at most S vertices in fast memory before the start of each segment. Considering the set of computed vertices in one of these segments, we can bound the size of its "frontier" by 2S: there can be at most S vertices in fast memory before the execution of the segment, and by construction there are exactly S loads.

Smith et al. [START_REF] Michael | A Tight I/O Lower Bound for Matrix Multiplication[END_REF] introduced a generalization of this argument, leading to tighter bounds in many cases. The idea is to decompose the execution into segments with T loads. This leads to a (S + T )-partitioning lemma instead of the original 2S. We provide formal definitions below.

Definition 6 (In-set). Let G = (V , E) be a DAG, P ⊆ V be a vertex set in G. The In-set of P is the set of vertices outside P with a successor inside P . Formally,

In (P ) = {v ∈ V ∖ P , ∃(v, w) ∈ E ∧ w ∈ P } Definition 7 (K-bounded set). Let G = (V , E) be a DAG. A vertex set P ⊆ V is called K-bounded if In (P ) ≤ K. Definition 8 (K-partition of a CDAG). Let G = (V , E, I) be a CDAG. A K-partition of G is a collection of subsets V 1 , V 2 , . . . , V m of V \ I such that: (1) {V 1 , V 2 , . . . , V m } is a partition of V \ I, i.e. ∀i j V i ∩ V j = ∅ and m i=1 V i = V \ I. (2) There is no cyclic dependence between V i 's. (3) Every V i is K-bounded, i.e. ∀i In (V i ) ≤ K.
Lemma 1 ((S + T )-Partitioning [START_REF] Hong | I/O complexity: The Red-Blue Pebble Game[END_REF]). Let T > 0. Any complete calculation R of the red-white pebble game on a CDAG G using at most S red pebbles is associated with a (S + T )-partition of the CDAG such that

Q R ≥ T • (h -1),
where Q R is the number of applications of rule (R1) in the game and h is the number of subsets in the partition.

In particular, an upper bound U on the size of a (S + T )-bounded set directly translates into a data movement lower bound:

Q(G) ≥ T • |V \I| U -1 .
This bound can actually be improved to:

Q(G) ≥ T • |V \ I| U . (3.1)
Indeed, the proof of Lemma 1 establishes a correspondence between the number of sets in a (S + T )-partition of G and a partition of an execution of the game in segments containing exactly S loads, except maybe the last one. The subtraction by one is due to this last segment.

When |V \ I| /U is an integer, the segment contains exactly S loads and thus the subtraction is not necessary. We actually want to be able to compute lower bounds for CDAGs in which no vertices are tagged as input (this is particularly useful when doing decomposition, see Sec. 3.2) The following lemma (Lemma 2) establishes such a result. The main idea is as follows: we tag some vertices as input, getting a new CDAG on which Lemma 1 applies and gives some bound. The additional I/O cost is at most the number of input vertices that were added, so we get a lower bound for the original CDAG by subtracting this number from the bound. See [START_REF] Elango | On characterizing the data movement complexity of computational DAGs for parallel execution[END_REF] for a complete proof.

Definition 9 (Sources). Let G = (V , E) be a DAG, P ⊆ V be a vertex set in G. The sources of P are the vertices of P with no predecessors in P . Formally, Sources (P ) = {v ∈ P , ∄u ∈ P , (u, v) ∈ E} Lemma 2 ((S + T )-Partitioning I/O lower bound, no input case [START_REF] Elango | On characterizing the data movement complexity of computational DAGs for parallel execution[END_REF]). Let S be the capacity of the fast memory, let G = (V , E, ∅) be a CDAG, and let h be the minimum number of subsets in a (S + T )-partition of G I = (V , E, I = Sources (V )) for some T > 0. Then, the minimum I/O for G satisfies:

Q(G) ≥ T • (h -1) -|Sources (V )| .

Using projection to bound the cardinality of K-bounded sets

The key idea behind the automation of data movement lower bound computation is the use of geometric inequalities through an appropriate program representation. Vertices of a CDAG are mapped to points in a multidimensional geometric space E ≃ Z d through some mapping ρ (where dimensions are typically loop indices), and regular data dependencies in the CDAG are represented as projections on a lower-dimensional space.

The condition "set of vertices P ⊂ V is K-bounded" in the CDAG corresponds to a condition of the form "the size of the projections of ρ(P ) in E is bounded by K". Finding a bound on the size of a K-bounded set in a CDAG can thus be reduced to: finding a bound on the size of a set E in a geometric space, given cardinality bounds on some of its projections. This correspondence is developed in Sec. 3.3. In this section, we take it for granted and only introduce the mathematical notations and results.

There exist inequalities for doing exactly what we need, namely the discrete Brascamp-Lieb inequality, introduced by Christ et al. [START_REF] Christ | Communication Lower Bounds and Optimal Algorithms for Programs that Reference Arrays -Part 1[END_REF] as a discrete analogue to the one established by Brascamp and Lieb for metric spaces [START_REF] Herm | Best constants in Young's inequality, its converse, and its generalization to more than three functions[END_REF].

Theorem 1 (Brascamp-Lieb inequality, discrete case [START_REF] Christ | Communication Lower Bounds and Optimal Algorithms for Programs that Reference Arrays -Part 1[END_REF]). Let d and d j be nonnegative integers and φ j :

Z d → Z d j be group homomorphisms for 1 ≤ j ≤ m. Let 0 ≤ s 1 , s 2 , . . . , s m ≤ 1. Suppose that: rank (H) ≤ m j=1 s j • rank φ j (H) for all subgroups H of Z d (3.2)
Then:

|E| ≤ m j=1 φ j (E) s j for all nonempty finite sets E ⊆ Z d . (3.3)
A special case is when the φ j are the canonical projections on Z d-1 along basis vectors,

and 1) . Fig. 2.4 illustrates this special case in three dimensions.

s 1 = • • • = s m = 1 d-1 , giving a bound of the form |E| ≤ d i=1 φ j (E) 1/(d-
The issue, when trying to apply Theorem 1, is that (3.2) has to be true for all subgroups, which can obviously be quite difficult to establish. However, as all the coefficients are integers and bounded by d, the number of distinct inequalities in (3.2) is bounded. The set of admissible s j is thus a (convex) polyhedron and it has been shown [START_REF] Christ | On Holder-Brascamp-Lieb inequalities for torsion-free discrete Abelian groups[END_REF] that the polyhedron defined by these inequalities is computable. The algorithm is actually combinatorial and quite complex, so we do not use it in our present work. Instead, we use the following result, which restricts the set of subgroup for which (3.2) has to be verified.

Definition 10 (Lattice of subgroups). The lattice of subgroups generated by subgroups H

1 , H 2 , . . . , H m of a group G is the closure of {H 1 , H 2 , . .

. , H m } under group sum and intersection.

Lemma 3 (Lattice of subgroups in Brascamp-Lieb [START_REF] Valdimarsson | The Brascamp-Lieb polyhedron[END_REF]). Theorem 1 holds with the weaker condition:

rank (H) ≤ m j=1 s j • rank φ j (H) for all subgroups H ∈ L φ 1 ,φ 2 ,...,φ m (3.2b)
where L φ 1 ,φ 2 ,...,φ m is the lattice of subgroups generated by Ker(φ 1 ), Ker(φ 2 ), . . . , Ker(φ m ).

The subgroup lattice is not necessarily finite, so this does not give a tractable algorithm, but this will be sufficient in most cases, as loop nests are usually of quite small dimensions, and data dependencies are not too complex. In our practical implementation, we use a timeout. We add projections: the more projections the tighter the bound. Each time we add a projection we update the lattice of subgroups and check for Conditions 3.2b to be satisfied. The process (adding projections) stops if we reach the time-out. This does not mean that the algorithm fails, rather that the bound will potentially be less tight (cf. Sec. 3.5).

A situation that often occurs is when the subgroups Ker(φ j ) are linearly independent: in this case there is no need to compute the generated lattice of subgroups, simply testing on each kernel individually for Conditions 3.2 is sufficient (see proof in [START_REF] Christ | Communication Lower Bounds and Optimal Algorithms for Programs that Reference Arrays -Part 1[END_REF], Sec. 6.3).

Choosing s j 's The goal is to solve the following problem: "Given a set of projections (group

homomorphisms in Z d ) φ 1 , . .

. , φ m and a constant K, find an upper bound (as tight as possible) on the cardinality of a set

E ⊂ Z d satisfying φ j (E) ≤ K."
For any coefficients s 1 , . . . , s j satisfying (3.2), Theorem 1 gives the following bound on |E|:

|E| ≤ m j=1 φ j (E) s j ≤ m j=1 K s j = K j s j .
To get a bound as tight as possible on |E|, we want to minimize the right-hand side of this inequality. This amounts to minimizing j s j while satisfying the constraints in (3.2). Since these constraints are linear inequalities, the optimal choice for s j 's can be obtained by a linear solver.

In the special case where the φ j 's are orthogonal projections along basis vectors (and m = d), kernels are linearly independent and the linear program is:

Minimize j s j s.t. ∀1 ≤ i ≤ d, 1 ≤ j i s j
and its solution is, as expected,

s 1 = • • • = s d = 1 d-1 .

DFG Paths

A fundamental object in our lower bound analysis is a DFG-path, which is simply a directed path in a DFG. The relation R p of a DFG-path p = (e 1 , . . . , e k ) is the composition of the relations of its edges:

R p = R e 1 • • • • • R e k . We
S b [⃗ x] → S a [A • ⃗ x + ⃗ b],
where A is not full-rank. Intuitively, a chain circuit corresponds in the CDAG to "iterative" dependencies, for instance every statement S i,j in a 2-dimensional loop depending on the result of statement S i-1,j . Broadcast paths correspond to some data being reused multiple times, for instance a variable x being used by every statement S i in a one-dimensional loop. The dimension of the kernel of A in the definition above corresponds to the dimension of the set of statements that use a single piece of data: it is of dimension d if it is used in every iteration of a d-dimensional loop. In both cases, these are regular data reuse patterns that can be exploited by our geometric approach.

In Fig. 3.2, path p = (e 3 ) is a chain circuit, going from S to itself with translation vector

⃗ b = (1, 0). Path p ′ = (e 2 ) is a broadcast path, with relation R p ′ = R e 2 = {C[t] → S[t, i] : 0 ≤ t < M ∧ 0 ≤ i < N }. The inverse relation is the linear function ⃗ I → A• ⃗ I + ⃗ b, with A = (1 0), ⃗ b = (0). The kernel of A is {(0, i), i ∈ R}.

CDAG decomposition

To derive data movement lower bounds for a complex program, it is essential to be able to decompose it into subregions for which we can compute lower bounds, and then sum the complexity for each subregion. The no recomputation condition is necessary for such a decomposition. Under this hypothesis, it is quite straightforward to see that a decomposition into disjoint subregions is sufficient. In this section, we provide a more general decomposition lemma, using the fact that vertices of a subregion that will not be counted as loads can also be part of another subregion. We then explain how it is applied on the DFG representation, distinguishing two cases: combining a fixed number of program regions (see example in Fig. 3.4); and summing over all iterations of a loop (see example in Fig. 3.3), which amounts to combining an unbounded (parametric) number of program regions. We stress that the CDAG partitioning method (Sections 3.1.3 and 3.1.4) and the CDAG decomposition method (this section) are two distinct things, used at different stages in the global algorithm.

Non-disjoint Decomposition Lemma

Definition 11 (sub-CDAG, no-spill set). Let G = (V , E, I) be a CDAG, and

V i ⊂ V . The sub- CDAG G |V i of G is the CDAG with vertices V i , edges E i = E ∩ (V i × V i ) and input vertices I i = I ∩ V i .
The no-spill set of G |V i is the subset of vertices of V i \ I i with either:

1. no outgoing edges in E i , or

no incoming edges in E i and at most one outgoing edge in E i

The may-spill set of

G |V i is the complement of its no-spill set in V i . Lemma 4 (CDAG decomposition). Let G = (V , E, I) be a CDAG. Let V 1 , V 2 , . . . ,
Q(G) ≥ k i=1 Q(G |V i ).
Proof. Let R be an optimal S-RW-game on G, with cost Q = Q(G). For all i, we denote Q i the cost of R restricted to V i , that is the number of applications of rule (R1) on vertices in the may-spill set of V i . Since the may-spill sets are pairwise disjoint, clearly

Q ≥ k i=1 Q i . For all i, we will build from R a valid game R i for G |V i with cost Q i . This will show that Q(G |V i ) ≤ Q i , from which follows k i=1 Q(G |V i ) ≤ k i=1 Q i ≤ Q = Q(G)
, establishing the result. To build the game R i from R, we proceed as follows:

1. Remove every move involving vertices outside V i .

2. For no-spill vertices in V i without successors, remove every application of rule (R1) and (R3).

3. For no-spill vertices v in V i with no predecessors and one successor w, move the single application of rule (R2) on v just before the application of (R2) on w, add a (R3) move for v just after this point and remove all subsequent (R1) and (R3) moves on v.

It is clear that after step 1, we have a valid game for G |V i . Indeed conditions for (R1) and (R3) are trivially preserved, and since we kept all (R2) moves on vertices in V i , there will always be the necessary red pebbles to apply (R2).

Step 2 also gives a valid game, since no (R2) moves can depend on such a vertex having a red pebble.

Step 3 is also a valid transformation because since v does not have any predecessor in V i (and is not an input vertex), it can be activated via (R2) and any given point. Therefore activating it just when it is needed, and applying (R3) just after is valid. This preserves all (R1) moves on the may-spill set of G |V i and removes all (R1) moves on its no-spill set, thus the cost of

R i is indeed Q i .
IOLB implements two different mechanisms that make use of the non-disjoint decomposition lemma. The basic one (bounded combination -Sec. 3.2.2) simply decomposes the CDAG into a bounded number of sub-CDAGs (e.g.,corresponding to different sub-regions of the code), computes the corresponding I/O complexities, and combines them. The more complex one (loop parametrization -Sec. 3.2.3), decomposes the CDAG into an unbounded number of sub-CDAGs by "slicing" the iteration space of a loop nest. IOLB combines the two mechanisms. The following example illustrates the decomposition lemma for loop parametrization.

Illustrating example Consider Example 2 on Fig. 3.3. The CDAG can be decomposed into M -1 identical subgraphs, as shown on Fig. 3.3c (each subgraph G |V t , t = 1, . . . , M -1 corresponds to iteration t of the loop enclosing S 1 , and iterations t -1 and t of the loop enclosing S 2 ). On each of these subgraphs, the may-spill set contains the two "bottom" rows (because vertices in the "top" row have no successor in the sub-CDAG). Thus the may-spill sets of these subgraphs are pairwise disjoint and the I/O for the whole CDAG is greater than the sum of the individual I/O for each subgraph by Lemma 4.

On each subgraph G |V t , the wavefront method (Sec. 3.4) can be applied, giving a lower bound on intersect, the individual complexities can be summed over t = 1, . . . , M -1, providing a lower bound for the whole CDAG:

I/O of Q(G |V t ) ≥ N -S.
Q(G) ≥ (M -1)(N -S).

Bounded combination

The main procedure of IOLB (Sec. 3.5.3) selects a bounded set of (possibly overlapping) sub-CDAGs and computes their individual complexities. The objective of Alg. 1 is to combine (sum) as many non-interfering (disjoint may-spill sets) complexities as possible. It does so using a greedy approach: Assume there are two sub-CDAGs both with a "high" complexity but with non-disjoint may-spill sets. Alg. 1 will select the one with the highest complexity, recompute the complexity of the second after removing the intersecting part, and then sum them up. The overall set of sub-CDAGs is iteratively processed this way (and the complexities summed-up) until empty or negligible complexities remain. The comparison (what is "higher") is done using instances of parameter values, simply evaluating the corresponding symbolic expressions. It should be emphasized that the final bound is a valid lower bound for any parameter values, the instances of parameter values are only used for heuristics.

Let us have a look at the example on Fig. 3.4. In the original code (3.4a), notice that k is the outer loop index, meaning that A[k] will have been modified either in the current loop iteration or the previous one depending on the order between i and k (Floyd-Warshall exhibits the same pattern, with three loops instead of two). This is made clear in the single-assignment form (3.4b), and can be visualized in the CDAG representation (3.4c). The dependences on input values are grayed out in (3.4b) and omitted in (3.4c), and we will ignore them in the discussion to keep the explanations simple.

1 function combine_subQ input : A DFG G, an instance I, a set of complexities Q output: A combined complexity Q I 2 Q I = 0; 3 Let G ′ be a copy of G; 4 while Q ∅ do 5 let Q such that Q(I) = max Q Q(I); 6 Q = Q -{Q}; 7 if Q(I) = 0 then return Q I ; 8 if G ′ ∩ Q.may-spill ∅ then 9 Recompute Q assuming CDAG G ′ ; 10 Q = Q ∪ {Q}; 11 else 12 Q I := Q I + Q; 13 G ′ = G ′ -Q.may-spill 14 return Q I Algorithm 1:
Summing lower bound expressions by removing interferences Considering only the statement vertex S in the DFG, the dependency analysis gives the following relations:

R 1 = {S[k -1, i] → S[k, i] : 1 ≤ k < N ∧ 0 ≤ i < N } R 2 = {S[k -1, k] → S[k, i] : 1 ≤ k < N ∧ 0 ≤ i < k} R 3 = {S[k, k] → S[k, i] : 0 ≤ k < N ∧ k < i < N }
The image spaces of R 2 and R 3 provide a natural decomposition of the CDAG into two non-interfering sub-CDAGs, as shown in (3.4d). On each part, the pattern is similar to that of Example 1 on page 14, and the geometric approach gives a lower bound (omitting lower order terms)

Q(G i ) ≥ N 2
2S . Since they do not interfere, Alg. 1 will return their sum Q(G) ≥ N 2 S , independently of the parameter instance.

Loop parametrization

As done on the example above, IOLB can compute the I/O complexity of some inner loop nests of a bigger enclosing loop nest and sum them. To this end, our scheme performs what we call loop parameterization. Loop parameterization considers each individual sub-CDAG where the outermost indices are fixed (our algebraic formulation allows us to consider such indices as parameters without the need to explicitly enumerate them) enriched by their input vertices. Taking the notations

V i = V i ∪ In (V i )
parametrizing the outer "t" loop with t = Ω (with Ω a parameter) allows us to compute Q Ω , a (parametric) lower bound for each individual value of Ω = 1, . . . , M -1 (The sub-CDAG for

Parameters: N; Input: A[N]; Output: A[N]; f o r (k=0;k<N;k++) f o r (i=0; i<N; i++) A[i] = f(A[i],A[k]);
(a) C-like code

Parameters: N; Input: A[N]; Output: S N -1 [N]; f o r (0 ≤ k < N and 0 ≤ i < N ) i f (k==i==0): S 0,i = f(A[0],A[0]); e l s e i f (k==0): S 0,i = f(A[i],S 0,0 ); e l s e i f (i<=k): S k,i = f(S k-1,i , S k-1,k ); e l s e i f (i>k): S k,i = f(S k-1,i , S k,k ); (b) Corresponding single assignment form k i R 1 R 2 R 3 (
Q = 1≤Ω<M Q Ω = 1≤Ω<M Q |{v∈V , t=Ω} = 1≤Ω<M N -S = (M -1) • (N -S).
In more complex cases, the parametric bound can depend on the outer loop parameter Ω, and we use formulas for sum of polynomials.

In ISL terms, this is done by making the outer loop index a parameter. Here the original domain

D S 1 = [M, N ] → {S 1 [t, i] : 0 < t < M ∧ 0 < i < N } becomes D Ω S 1 = [M, N , Ω] → {S 1 [t, i] : t = Ω ∧ 0 < t < M ∧ 0 < i < N }.
The corresponding parts of the algorithms are highlighted in Algorithm 6 on page 38.

K-partition bound derivation

In this section, we explain how to apply the geometrical reasoning of Sec. 3.1.4 on a CDAG G = (V , E), using its compact representation as a DFG. We also present, in 3.3.1, a generaliza-tion of one of the techniques introduced in [START_REF] Dongarra | Revisiting matrix product on master-worker platforms[END_REF][START_REF] Lowery | Improving the Communication Lower Bounds for Matrix-Matrix Multiplication[END_REF][START_REF] Michael | A Tight I/O Lower Bound for Matrix Multiplication[END_REF]] that these authors used to derive a tighter lower bound for matrix multiplication.

To apply Lemma 2 on G, we need to find a lower bound on the minimum number of subsets h in any K-partition of G. The general reasoning is as follows:

1. Embed V in a geometric space through a map ρ : P ⊆ V → E ⊆ Z d , such that two disjoint subsets of V are mapped to disjoint subsets of Z d . We have

ρ(P ) ≤ |P | .
2. Use the DFG representation to find a subset V ′ ⊆ V and a set of projections (group homomorphisms) φ 1 , . . . , φ m with the property that:

Any K-bounded set P ⊆ V ′ \ Sources (V ′ ) satisfies φ j (ρ(P )) ≤ K. (3.4)
3. Using Theorem 1, derive an upper bound U on ρ(P ) for any K-bounded P . This provides a lower bound

|V ′ \Sources(V )| U on the number h of disjoint K-bounded sets in V ′ \ Sources (V ′ ).

Geometric embedding, DFG-paths and projections

Let S k be some fixed DFG-vertex (corresponding to one program statement). Let Q 1 , . . . , Q m be DFG-paths all ending in S k , with a common image space

D k = {S k [i 1 , . . . , i d ] : . . . }.
The embedding ρ is defined as:

ρ(P ) = {(i 1 , . . . , i d ) | S k [i 1 , . . . , i d ] ∈ P }.
Vertices corresponding to statement S k are mapped to their corresponding d-dimensional point, and other vertices are ignored.

Definition 12 (embedded projections). For a given path Q with relation R Q , the geometric projection φ Q is defined as follows:

-If the path is a broadcast path with

R Q = {S j [j 1 , . . . , j d ′ ] → S k [i 1 , . . . , i d ] :
. . . } for some statement S j (not necessarily S k ), then the projection is directly given by the path relation

φ Q (i 1 , . . . , i d ) = (j 1 , . . . , j d ′ ). -If the path is a chain circuit with R Q = {S k [i 1 , . . . , i d ] → S k [i 1 + δ 1 , . . . , i d + δ d ] : . . . }, then the
projection is the orthogonal projection on the hyperplane in Z d defined by orthogonal vector δ = (δ 1 , . . . , δ d ). Its explicit formulation can be computed but is not needed here.

In the case of a broadcast, it is straightforward that φ Q satisfies (3.4), because R -1 Q (P ) is included in In (P ) for any P ⊆ V , so φ Q (ρ(P )) ≤ |In (P )| ≤ K for any K-bounded P .

In the case of a chain circuit, let us call I Q (P ) = R -1 Q (P ) \ P . This is basically the In-set of P restricted to edges corresponding to DFG-path Q, so I Q (P ) ⊂ In (P ). The projection φ Q associates one point to each straight line directed by δ. Since P ⊆ V \ Sources (V ), there is at least one point in φ Q (ρ(P )) for every nonempty chain in P , and φ Q (ρ(P )) ≤ |In (P )| ≤ K.

Example Consider paths p 1 = (e 2 ) and p 2 = (e 3 ) in Fig. 3

.2. p 1 is a broadcast path with relation {C[t] → S[t, i]}, so the corresponding projection is φ 1 (t, i) = (t). p 2 is a chain path with relation {S[t, i] → S[t + 1, i]}, so the corresponding projection is φ 2 (t, i) = proj (1,0) (t, i) = (0, i) (see Fig. 2.3b).

Summing projections

In some cases, the parts of the In-set of a vertex set associated with two given path relations are actually disjoint. Let Q 1 and Q 2 be two such paths, such that R -1 Q 1 (P ) ∩ R -1 Q 2 (P ) = ∅ for any P ⊆ V \ Sources (V ). If these are two broadcast paths, then since R -1 Q i (P ) ⊂ In (P ), any K-bounded set P satisfies the stronger inequality:

φ Q 1 (ρ(P )) + φ Q 2 (ρ(P )) ≤ K The same holds if Q 1 is a chain circuit and R -1 Q 1 (P ) ∩ R -1 Q 2 (P ) = ∅
, by a similar argument. We say two paths Q 1 and Q 2 are independent for domain

D S if R -1 Q 1 (D S ) ∩ R -1 Q 2 (D S ) = ∅.
We can build the DFG-path interference graph: vertices are paths Q 1 , . . . , Q m and there is an edge between any independent pair of paths. In this graph, if vertices

Q i 1 , . . . , Q i t form a clique, then Any K-bounded set P ⊆ D k satisfies s φ i s (P ) ≤ K.
Example Looking again at Example 1, it is straightforward to check that paths p 1 and p 2 are independent, so a K-bounded set P actually satisfies φ 1 (P ) + φ 2 (P ) ≤ K.

Combining several such inequalities, such that every projection occurs at least once, leads to a general constraint of the form:

m j=1 β j φ j (E) ≤ K,
for some positive coefficients β j . This is achieved by finding a set of maximal cliques covering all vertices in the interference graph, and summing the corresponding inequalities. Computing these β j 's is the role of function coeffInterf in Algorithm 4.

In this case, a tighter bound can be derived:

|E| ≤ m j=1 φ j (E) s j ≤ K j s j j s j m j=1 s j β j s j .
The following lemma establishes this result. Proof. We use Lagrange multipliers to find the constrained maximum of the function

ψ : x ∈ R m → m j=1 x s j j . L(x, λ) = m j=1 x s j j -λ         m j=1 β j x j -C        
Partial derivatives are:

∂L ∂x j = s j x s j -1 j k j x s k k -λβ j , 1 ≤ j ≤ m ∂L ∂λ = C - m j=1 β j x j
Setting them to be 0, we get:

λβ j x j = s j m k=1 x s k k , 1 ≤ j ≤ m
Summing for 1 ≤ j ≤ m gives:

λ m j=1 β j x j = λC =        m k=1 x s k k        • m j=1 s j
From which we derive:

x j = Cs j β j j s j , 1 ≤ j ≤ m
And finally:

ψ(x) ≤ m j=1 x s j j = m j=1 Cs j β j i s i s j = C j s j j s j m j=1 s j β j s j .
With this more general formulation, the choice of the s j coefficients is more involved than the linear optimization problem of Sec.3.1.4, and is developed in Sec. 3.3.3.

Kernel subgroup lattice

As already mentioned, the subgroup lattice (Def. 10 used in Lemma 3) is not necessarily finite, so it is better to build it step-by-step, updating it each time we add a new path. We set a time limit for the computation to converge, and do not add the path if this limit is reached. Function subspace_closure in Algorithm 2 tentatively updates the current subgroup lattice with a new one, returning the original lattice in case of a timeout.

Finding paths

The function that generates the set of paths P = {Q 1 , . . . , Q m } for a DFG-vertex S is named genpaths (Alg. 3). Starting from S, it uses a simple backward traversal (backward DFS) that favors walking through predecessors with largest domain. As the number of paths can be combinatorial, IOLB sets a timeout to avoid a computational blow-up. For a path P = (S 1 , . . . , S t = S), we store sub-path relations R S i →S for every intermediate statement S i . This is necessary to "remember" exactly which CDAG vertices are included in the computation.

if ∃H ∈ L ′ , H ∩ K L ′ then L ′ = L ′ ∪ {H ∩ K}; 5 else if ∃H ∈ L ′ , H + K L ′ then L ′ = L ′ ∪

Computing the lower bound

Once we have found a path combination, it is quite straightforward to apply the theoretical results introduced above. This is detailed in function sub_paramQ_bypartition in Alg. 4. Here, the role of the function call to coeffInterf is to compute the coefficients β j . It does so by finding a clique cover of the DFG-path independence graph and summing the constraints formed by each clique as explained in Sec. 3.3.1. The values for coefficients s j that satisfy inequalities 3.2b are then determined using convex optimization (IPOPT [START_REF] Wächter | On the Implementation of an Interior-Point Filter Line-Search Algorithm for large-scale nonlinear programming[END_REF] in our case) so as to minimize as much as possible the quantity

U = (S + T ) j s j j s j m j=1 s j β j s j . ( 3.6)
Indeed this expression being an upper bound on ρ(P ) (see Lemma 5), minimizing it has the effect of tightening 1 the computed I/O complexity. The constraints in (3.2b) describe a convex polyhedron, but the objective function (3.6) is not convex. In the basic case when all projections are simply bounded by (S + T ), the objective is U := (S + T ) j s j , so a natural objective is to minimize j s j in this generalized case. It can be easily checked that (3.6) is indeed equal to this when β j = 1 m for all j and s 1 = • • • = s j . Notice that the first factor in the expression of U depends only on the value of the sum.

So once s j is fixed, it is natural to minimize the second factor m j=1 s j β j s j , which is convex as objective: minimize j s j and then j s j β j s j constraints: ∀H ∈ L, j s j rank φ j (H) ≥ rank (H) where φ j = proj ⊥ Ker(P j ) } ;

T := 1 j s j -1 S ; foreach broadcast path P i = (S 0 , S 1 , . . . , S t = S) ∈ P do

U := m j=1 (S+T )s j β j i s i s j ; Q := max |D| U × T -
D ms := D ms ∪ t j=0 R -1 S j →S (D) ;
foreach chain circuit P i = (S 0 , S 1 , . . . , S t = S) ∈ P do

D ms := D ms ∪ t j=1 R -1 S j →S (D) ∪ R -1 S 0 →S (D) ∩ k j R -1 P k (D) ;
function Ker input : path P j output: linear space K such that the orthogonal projection φ j = proj ⊥ Ker(P )

if isBroadcast(P j ) then return Ker(j 1 . . . j d ′ ) where

R P j = {T [j 1 , . . . , j d ′ ] → [i 1 , . . . , i d ] : . . . } ;
if isChain(P j ) then return (δ 1 , . . . , δ d ) where a function of (s 1 , . . . , s m ). This convex optimization problem can be solved with an appropriate tool, such as IPOPT [START_REF] Wächter | On the Implementation of an Interior-Point Filter Line-Search Algorithm for large-scale nonlinear programming[END_REF]. The last step amounts to setting an appropriate value for T that provides a lower bound of |D S | U ×T -|I| (see Lemma 2) as big/tight as possible. Here, D S corresponds to V \Sources (V ) in the CDAG view, and I is the frontier of the domain, corresponding to Sources (V ). T is chosen as 1 j s j -1 S, because it maximizes the first term asymptotically. Finally, we store the may-spill set corresponding to the CDAG for which this lower bound is valid.

R P j = {S[i 1 , . . . , i d ] → S[i 1 + δ 1 , . . . , i d + δ d ] : . . . } ;
Taking the example given in Figure 3.1, it is sufficient to check condition (3.2b) on

H 1 = {(0, i)}, H 2 = {(t, 0)},

Extensions: Small Dimensions and Reduction Detection

We introduced two extensions to the method described above, that allow IOLB to derive much tighter bounds for some programs. For example, on a convolution (Fig. 3.5), the lower bound derived by the initial algorithm is several asymptotic orders below the upper bound. This lower bound can be improved by taking into account (i) the reductions across multiple dimensions, which refines the dependence analysis performed in the first steps of the algorithm, and (ii) small dimensions. Small dimensions are dimensions with a much smaller number of iterations than the size of the fast memory S. This property can be used to improve the power factor of the asymptotic bound.

Using Reductions

A reduction is the successive application of an associative and commutative binary operator over a set of values. For example, in the convolution kernel described in Figure 3.5, there is a reduction over 3 dimensions (c, h and w) with the addition operation. We only consider reductions that sum all the values along some dimensions, called reduced dimensions. This property can be used to sum the values of a set in any order, instead of having a fixed sequence of summation.

Because the program is expressed as a loop nest, its reduction is sequential. Thus, when the DFG is extracted from a code, there is a chain of dependencies linking nodes along reduced dimensions. When there are more than 2 reduced dimensions, this negatively impacts our study of DFG-paths, thus the quality of the extracted projections φ j and the upper bound found on the size of the K-bounded set E.

We detect a simple class of reductions by performing a pattern-matching on affine dependencies in the DFG. In this class, the reduced dimensions are defined over a rectangular domain and that the summation is performed in a lexicographic order over a permutation of these reduced dimensions. Once a reduction is detected, we replace the sequential chain of dependencies (between the nodes of the reduction) by two sets of broadcast dependencies: (i) starting from the computations using the final result of the reduction to all the nodes of the reduction, and (ii) starting from the initialization of the reduction to the nodes of the reduction.

On the convolution example of Figure 3.5, if we do not detect the reduction, the dependence on Out is along dimension s, and the corresponding projection found by IOLB is φ 1 (b, c, f , x, y, h, w) = (b, c, f , x, y, h, 0). But if we detect it, the projection of the new path becomes φ 1 (b, c, f , x, y, h, w) = (b, 0, f , x, y, 0, 0), which leads to a bigger kernel and better constraints on s j .

Adapting the Method to Small Dimensions

The geometrical argument of the K-partition method consists, at a high level, in bounding the size of the largest sub-CDAG P with a given operational intensity, by bounding its extent along each dimension. The bound obtained in the end is some power of the cache size S, which makes sense when loop limits along all dimensions are large (i.e., of the same magnitude and with their product >> S). However, in some programs, a subset of dimensions has small sizes in most real life implementations: in convolutions, H and W are almost always less than 10. When this is the case, the extent of P along small dimensions is more constrained by the total extent of the CDAG along this dimensions, than by the geometrically derived bound depending on S.

Formally, we assume that some of the dimensions of the iteration domain are small, i.e., their size is orders of magnitude smaller than S. For notation simplicity, let us assume that they correspond to the first q dimensions in statement domain D k = {S k [i 1 , . . . , i d ] : . . . }. We adapt the partitioning method by considering an additional group homomorphism φ sd to the application of Theorem 1. This group homomorphism is a projection of the space on the small dimensions:

φ sd = (i 1 , . . . , i q ).
Let N sd = φ sd (ρ(D k )) be the size of the image of D k through φ sd (in most cases, this will be equal to the product of small dimension sizes, e.g., N sd = H • W for convolution). Introducing a coefficient s sd associated with φ sd , we can rewrite Theorem 1 as follows.

Assuming: Note that the constraints on the s j are relaxed due to the contribution of s sd , giving us an opportunity to find smaller values of s j , compared to the usual method. Moreover, in order to get a bound as tight as possible, we should minimize s j first, s sd next, and finally

s j β j s j .
Let us illustrate this method on the convolution example described in Figure 3.5. Let us start with the "standard" derivation, i.e. without taking small dimensions into account. By examining the affine dependencies of the program, we find three projections φ 1 , φ 2 and φ 3 , described in Figure 3.7a.

As stated by Theorem 1, we consider the subgroups H generated from the combination of kernel vectors of the φ j . This yields constraints shown in Figure 3.7b, ignoring the rightmost column. Their corresponding subgroups are shown in the leftmost column. The three projections are independent so β 1 = β 2 = β 3 . The optimization problem can be written as:

φ 1 (b, c, f , x, y, h, w) = (b, 0, f , x, y, 0, 0) (from reduction) φ 2 (b, c, f , x, y, h, w) = (x + h, y + w, c, b) (from Image) φ 3 (b, c, f , x, y, h, w) = (f , h, w, c) (from Filter) φ sd (b, c, f , x, y, h, w) = (h, w) (small dimensions) (a) Projections (group homomorphisms) used in Brascamp-Lieb inequality. ( ⃗ f ) ≤ s 1 + s 3 ( ⃗ b) ≤ s 1 + s 2 (⃗ c) ≤ s 2 + s 3 (⃗ x) ≤ s 1 + s 2 ( ⃗ w) ≤ s 2 + s 3 +s sd ( ⃗ w -⃗ x) ≤ s 1 + s 3 +s sd ( ⃗ w, ⃗ x) 2 ≤ s 1 + s 2 + s 3 +s sd (⃗ y) ≤ s 1 + s 2 ( ⃗ h) ≤ s 2 + s 3 +s sd ( ⃗ h -⃗ y) ≤ s 1 + s 3 +s sd ( ⃗ h, ⃗ y) 2 ≤ s 1 + s 2 +
Minimize s 1 + s 2 + s 3 under constraints:

s 1 + s 2 ≥ 1 s 1 + s 3 ≥ 1 s 2 + s 3 ≥ 1 s 1 + s 2 + s 3 ≥ 2
The solution is s 1 = s 2 = s 3 = 2 3 and gives the following upper bound on ρ(P ) :

U = S + T 3 2 .
Finally we set T = 1 s 1 +s 2 +s 3 -1 S = S, and get the final bound (ignoring lower-order terms corresponding to the input):

Q = BCXY W H U • T = 9 4
BCXY W H S .

Now, we will see that taking small dimensions into account leads to a tighter bound. Here small dimensions are h and w, and N sd = W • H. Compared to the first derivation, we add a projection φ sd , shown in Figure 3.7a. The presence of this new projection adds a new coefficient s sd to the constraints on s j s, as shown in Figure 3.7b.

The optimization problem becomes:

Minimize s 1 + s 2 + s 3 , then s sd under constraints:

s 1 + s 2 ≥ 1 s 1 + s 3 ≥ 1 s 2 + s 3 ≥ 1 s 1 + s 2 + s 3 + s sd ≥ 2
The solution is now s 1 = s 2 = s 3 = s sd = 1 2 , and we get

U = S + T 3 3 2 • (W H) 3 2 .
Finally we set T = 1 s 1 +s 2 +s 3 -1 S = 2S, and get the final lower bound on I/O (ignoring again lower-order terms):

Q = BCXY W H U • T = 2BCXY √ W H √ S .
This is a tighter lower bound for most real-life configurations.

Wavefront bound derivation 3.4.1 Theoretical results

An alternative way to derive data movement lower bounds in the no-recomputation model is the wavefront abstraction. At any point in an execution of a RW-game, the wavefront is the set of vertices that have been computed but whose result is still needed by some successor (sometimes called the set of live vertices). If the size of the wavefront at some point in the execution is greater than the size of the fast memory, then necessarily some vertices have to be spilled to the slow memory and thus loaded using rule (R1). This is formalized in the definition and lemma below: Definition 13 (Wavefront). Let R be an execution of the RW-game on a CDAG G, and v a vertex of G. Consider the time t in the execution just before v has been computed (i.e., just after a white pebble has been placed on v using rule (R2)). The wavefront W R (v) in execution R is the set of vertices that, a time t, have been computed (i.e., have a white pebble) but have some successor that does not. Lemma 6 (Min-wavefront [START_REF] Elango | On characterizing the data movement complexity of computational DAGs for parallel execution[END_REF]). Let S be the capacity of the fast memory, and G = (V , E) be a CDAG. Let w min G = min R (max v∈V |W R (v)|), so that any valid RW-game on G has a wavefront of size at least w min G . Then,

Q ≥ w min G -S.
This lemma is quite general and cannot be applied directly, as w min G is not usually computable. We thus provide the following result, which uses a condition on the structure of the CDAG to get a lower bound on the size of a wavefront in any execution. By Lemma 6, this implies:

Corollary 1. Let G = (V , E) be a CDAG, and V 1 , V 2 be disjoint subsets of V such that every vertex in V 2 is reachable from every vertex in V 1 through some path in G. Let L 1 , . . . , L m be disjoint paths in G, starting in V 1 and ending in V 2 . More formally:        ∀L j = (v j 1 , . . . , v j l j ), v j 1 ∈ V 1 and v j l j ∈ V 2 V 1 × V 2 ⊂ E * Then, w min G ≥ m. V 1 V 2
Q ≥ m -S.
Proof. Let v be the first vertex to be computed among vertices of V 2 in some fixed RW-game on G. By the second condition, every vertex in V 1 must have been computed since v depends on all of them. Just before a red pebble is placed on v, no vertex in V 2 has a red pebble. Therefore there is a live vertex (that has a red pebble and a successor without one) in every path L j . Thus w min G ≥ m.

The common case to use this technique to get a strong data movement lower bound is to combine it with the parametric CDAG decomposition (Sec. 3.2.3).

Example

In the example of Fig. 3.3, Corollary 1 can be applied on each subgraph, with V 1 and V 2 chosen as shown on Fig. 3.8. It can be easily checked that every vertex in V 1 can reach every vertex in V 2 , and there are N disjoint paths L j , so this gives a lower bound

Q(G |V t ) ≥ N -S.

Implementation

To apply this technique, our algorithm tries to uncover a set of disjoint paths satisfying the hypotheses of Corollary 1. To reduce the search space, Algorithm 5 looks for a much more constrained pattern, in which all disjoint paths L j begin and end in different instances of some statement S, with an increment in the innermost parametrized loop index (see Sec. 3.2.3). This amounts to finding an injective circuit in the DFG with a relation of the form

{S[I 1 . . . I d , i d+1 . . . i D ] → S[I 1 . . . I d + 1, i d+1 . . . i D ]}.
Intuitively, we look at two "slices" of the CDAG, each of them of dimension D -d, representing two successive iterations of the body of the loop iterating over dimension d (with index I d and I d + 1, respectively). We try to find subsets of these two slices (corresponding to V 1 and V 2 ) such that every vertex in the first one can reach every vertex in the second one.

In Algorithm 5, the first loop computes the following relations: -R S→S is the union of all path relation of elementary circuits from S to itself, -RL S→S is the union of all affine path relations (where every subpath relation is affine) of elementary circuits S → S, -R S→ * is the union of all path relations of elementary paths from S to any other DFGvertex.

This is then used to compute the relation R I d . Here I d is the index of the innermost parametrized loop, and R I d is the restriction of RL S→S to paths that do exactly one step along dimension d, not changing any other index. This represents a set of disjoint paths going from instances of statement S in "slice" I d to instances of S in "slice" I d + 1.

To apply Corollary 1, we need to restrict R I d to a domain W where every starting vertex reaches every ending vertex. To do so, we first compute X = R complete -(R S→S ) * Dom R I d , that is, all vertices in "slice" I d + 1 that are not reachable from "slice" I d . We then take

W = Dom R I d -R -1 I d (X)
, that is, we only keep vertices in "slice" I d that can reach every vertex in slice Our front end (PET [START_REF] Verdoolaege | Polyhedral Extraction Tool[END_REF]) takes as input a program in C where the to-be analyzed regions (SCoPs -Static Control Parts) are delimited by #pragma scop and #pragma endscop annotations. For PET, all array accesses are supposed not to alias with one another. Any scalar data is assumed to be atomic and all of the same size: our CDAG is not weighted (which is a limitation of our implementation and not a conceptual limitation of the approach). As illustrated by the example of Fig. 3.1 and 3.2 (multidimensional-)array accesses are affine expressions of static parameters and loop indices. A static parameter can be the result of any complex calculation but has to be a fixed value for the entire execution of the region. Loop bounds and more generally control tests follow the same rules (affine expressions). As a consequence, the iteration space is a union of (parametric) polyhedra, and memory accesses (read and writes) are piecewise affine functions. This representation of the region execution that fits into the polyhedral framework [START_REF] Feautrier | Polyhedron model[END_REF] allows IOLB to compute data dependencies using standard data-flow analyses.

I d + 1. Application of Corollary 1 with V 1 = W , V 2 = R I d (W ) and {L j } j = R I d gives Q ≥ |W | -S.
A := {(S i , S j ) ∈ E, R S i →S j is affine and R -1 S i →S j injective}; 6 RL S→S j := (S i ,S j )∈A RL S→S i • R S i →S j ; 7 R S→S j := (S i ,S j )∈E R S→S i • R S i →S j ; 8 R S→ * := R S→ * ∪ R S→S j ;
PET outputs a polyhedral representation of the input C program, from which we extract a Data-flow graph (DFG) G = (S, D) (see Sec. 3.1.2).

Instances of parameter values

As briefly explained in Sec. 3.2, to generate bounds that are as tight as possible, our heuristic needs to take decisions. Such decisions are based on our ability to compare the size of two different domains sizes or even the complexity of two different sub-CDAGs. The overall framework being parametric (it provides complexities that are functions of parameter values and cache size), a total order is obtained by considering a specific instance of parameter values, taken as an additional input alongside the C program. One needs to outline that a specific instance of parameter values is not considered by the algorithm as a precondition: For a given instance, the computed lower bound expression is universal, i.e., is correct for any parameter values. For completeness, several instances are considered, and to each instance I is associated a complexity Q I . As we have Q ≥ Q I for any instance, denoting I the set of all considered instances, they are simply combined as:

Q = max I∈I Q I .

Main algorithm

Alg. 6 contains the skeleton of the main part of IOLB, with links to corresponding subsections.

To make it concrete, we show a step-by-step execution of the algorithm on the cholesky kernel. The pseudo-code and associated DFG for cholesky are reported in Fig. 3.9.

In this example, the K-partition method is the method that yields the strongest bound. To keep things tractable, we will detail only the parts of the algorithm that contribute to this bound: the iteration of the outer loops (lines 4 and 5) for which d = 0 and S = S 3 , and only the K-partition part (lines 8 to 18, corresponding to Sec. 3.3).

The DFG contains three statement vertices {S 1 , S 2 , S 3 } (the vertex corresponding to input array A and the corresponding dependences are omitted as they do not play a role in the lower bound derivation). The main loop of Alg. 6 iterates on those statements and computes some lower bound complexities for each of them.

Procedure genpaths (called at line 11 in Alg. 6) traverses the DFG, searching for chain or broadcast paths ending in S (cf. Sec. 3.1.2). Here, it finds three "interesting paths" for statement S 3 . These are the three paths pointing to S 3 , namely: P 1 = (e 1 ) (chain path), P 2 = (e 2 ) and P 3 = (e 3 ) (broadcast paths). Their relations are: R P 1 = R e 1 , R P 2 = R e 2 , R P 3 = R e 3 (cf. Fig.

3.9c).

The corresponding projections and kernels are: 18

φ 1 (k, i, j) = proj (1,0,0) (k, i, j) = (0, i, j) K 1 = Ker(φ 1 ) = ⟨(1, 0, 0)⟩ φ 2 (k, i, j) = (k, j) K 2 = Ker(φ 2 ) = ⟨(0, 1, 0)⟩ φ 3 (k, i, j) = (k, i) K 3 = Ker(φ 3 ) = ⟨(0,
(Q, G ′ ) = combine_paramQ(Q, G ′ , sub_paramQ_bypartition(P , D S , L, Ω d )); 19 (Q, G ′ ) = combine_paramQ(Q, G ′ , sub_paramQ_bywavefront(S, Ω d )); Q low = input_size(G) + max(0, combine_subQ(Q)); function combine_paramQ input : set of global bounds Q, DFG G ′ , parametrized bound Q(Ω) output: updated Q, G ′ if [Ω Ω ′ ⇒ Q.interf(Ω) ∩ Q.interf(Ω ′ ) = ∅] then Q := Ω Q(Ω); Q.may-spill := Ω Q.may-spill(Ω); Q = Q ∪ {Q}; G ′ := G ′ \ Q.may
f o r ( k = 0 ; k < n ; k ++) A [ k ] [ k ] = s q r t ( A [ k ] [ k ] ) ; // S1 f o r ( i = k +1; i < n ; i ++) A [ i ] [ k ] /= A [ k ] [ k ] ; // S2 f o r ( i = k +1; i < n ; i ++) f o r ( j = k +1; j <= i ; j ++) A [ i ] [ j ] -= A [ i ] [ k ] * A [ j ] [ k ] ; // S3 (
R e 1 = {S 3 [k -1, i, j] → S 3 [k, i, j] : 1 ≤ k < N ∧ k + 1 ≤ i < N ∧ k + 1 ≤ j ≤ i} R e 2 = {S 2 [k, j] → S 3 [k, i, j] : 0 ≤ k < N ∧ k + 1 ≤ i < N ∧ k + 1 ≤ j ≤ i} R e 3 = {S 2 [k, i] → S 3 [k, i, j] : 0 ≤ k < N ∧ k + 1 ≤ i < N ∧ k + 1 ≤ j ≤ i} R e 4 = {S 3 [k -1, i, k] → S 2 [k, i] : 1 ≤ k < N ∧ k + 1 ≤ i < N } R e 5 = {S 1 [k] → S 2 [k, i] : 0 ≤ k < N ∧ k + 1 ≤ i < N } R e 6 = {S 1 [k -1, k, k] → S 1 [k] : 1 ≤ k < N ∧ k + 1 ≤ i < N } (c) Edge relations
D S := D S 3 = {S 3 [k, i, j] : 0 ≤ k < N ∧ k + 1 ≤ i < N ∧ k + 1 ≤ j ≤ i} .
The foreach loop then iterates over {P 1 , P 2 , P 3 }. Here the three paths have domains that almost span the entire domain of S 3 , so the condition at line 12 is always true (γ is a constant between 0 and 1). They also have pairwise orthogonal kernels, so the condition at line 14 is also true at each iteration, and at the end of the foreach loop:

D S = ((D S 3 ∩ Im (P 1 )) ∩ Im (P 2 )) ∩ Im (P 3 ) = {S 3 [k, i, j] : 1 ≤ k < N ∧ k + 1 ≤ i < N ∧ k + 1 ≤ j ≤ i} P = {P 1 , P 2 , P 3 }
At line 18, function sub_paramQ_bypartition derives a lower bound from the set of paths P , which is then added to the set of lower bounds Q by function combine_paramQ.

We can check that P 1 is independent from P 2 and P 3 , but P 2 and P 3 are not (that is R -1

P 1 (D) ∩ R -1 P 2 (D) = ∅, R -1 P 1 (D) ∩ R -1 P 3 (D) = ∅ and R -1 P 2 (D) ∩ R -1 P 3 (D) ∅).
Thus the following inequality holds for any K-bounded set E in the CDAG:

φ 1 (E) + 1 2 φ 2 (E) + 1 2 φ 3 (E) ≤ K. (3.7)
Let s 1 , s 2 , s 3 be the solutions to the following optimization problem:

Minimize σ := j s j s.t. 0 ≤ s 1 , s 2 , s 3 ≤ 1 1 ≤ s 2 + s 3 1 ≤ s 1 + s 3 1 ≤ s 1 + s 2 (3.8)
The discrete Brascamp-Lieb theorem [START_REF] Christ | Communication Lower Bounds and Optimal Algorithms for Programs that Reference Arrays -Part 1[END_REF], combined with Lemma 5 applied on projections φ i , guarantee that, for any K-bounded set E:

|E| ≤ K σ 2s 1 σ s 1 2s 2 σ s 2 s 3 σ s 3 . (3.9)
The solution to (3.8) is

s 1 = s 2 = s 3 = 1 2 , so |E| ≤ 2 • (K/3) 3/2 .
Lemma 2 tells us that, if U is an upper bound on the size of a (S + T )-bounded-set in G, then:

Q(G) ≥ T • |V \ Sources (V )| U -|Sources (V )| .
Here 2 Taking for U our upper bound on |E| provides the following inequality for which the objective is to set a value for T that maximizes its right hand side:

V = D ∪ R -1 P 1 (D) ∪ R -1 P 2 (D) ∪ R -1 P 3 (D), giving: V \ Sources (V ) = {S 3 [k, i, j] : 1 ≤ k < N ∧ k + 1 ≤ i < N ∧ k + 1 ≤ j ≤ i} Sources (V ) = {S 3 [0, i, j] : 1 ≤ i < N ∧ 1 ≤ j ≤ i; S 2 [k, i] : 1 ≤ k < N ∧ k + 1 ≤ i < N } So |V \ Sources (V )| = N 3 6 and |Sources (V )| = N 2 .
Q ≥ T × N 3 /6 2 • (K/3) 3/2 -N 2 ≈ T (S + T ) 3/2 × N 3 /6 2 • (1/3) 3/2 .
Setting T = 2S (so K = S + T = 3S) leads to the following lower bound on Q:

Q ∞ low = (2S) × N 3 /6 2S 3/2 = N 3 6 √ S .
Concerning the parts of the algorithm that were not detailed here: the outermost loop (Line 4) corresponds to the loop parametrization detailed in Sec. 3.2.3: for each loop depth d, outermost indices are fixed (as parameter Ω d -Line 6), and parametrically computed lower bounds are summed (when not interfering -Line 22) over all iterations (Line 23 in com-bine_paramQ). The loop on statements S (Line 5) decomposes the full CDAG into as many "S-centric" sub-CDAGs. The so-obtained bounded set of lower bounds Q are combined using procedure combine_subQ (Line 20) as described in Sec. 3.2.2. To take compulsory misses into account, the size of the input data of the program is added to the expression.

For each statement S, both techniques (K-partition and wavefront, resp. Line 18 and Line 19) generate lower bounds. As opposed to the implicitly considered "S-centric" sub-CDAGs for the wavefront reasoning, an "S-centric" sub-CDAG for the K-partition reasoning (which is built by finding a set P of DFG-paths that terminate at S -Lines 10-17 through function genpaths) does not necessarily span all the S-vertices (D S ) of the CDAG. So several (non-intersecting) sub-CDAGs can be built until no more interesting lower bound can be derived (Line 17).

Experimental evaluation

Implementation

IOLB was implemented in C, using ISL-0.13 [START_REF] Verdoolaege | ISL: An Integer Set Library for the Polyhedral Model[END_REF], barvinok-0.37 [START_REF] Verdoolaege | Counting Integer Points in Parametric Polytopes Using Barvinok's Rational Functions[END_REF] and PET-0.05 [START_REF] Verdoolaege | Polyhedral Extraction Tool[END_REF]. We also used GiNaC-1.7.4 [START_REF] Bauer | Introduction to the GiNaC Framework for Symbolic Computation within the C++ Programming Language[END_REF] for the manipulation of symbolic expressions, and PIP-1.4.0 [START_REF] Feautrier | Parametric Integer Programming[END_REF] for linear programs. IOLB takes as input an affine C program and outputs a symbolic expression for a lower bound on I/O complexity as a function of the problem size parameters of the program and capacity of fast memory. Figure 3.10 shows an example of annotated C code (jacobi-1d from PolyBench), with the output from IOLB.

Evaluation on PolyBench

IOLB was applied to all programs in the PolyBench/C-4.2.1 benchmark suite [START_REF] Pouchet | PolyBench/C 4.2[END_REF]. For each kernel, our tool outputs an I/O lower bound expression Q low , from which we derive an upper bound on operational intensity OI up by forming the ratio of the number of operations and Q low . To evaluate the quality of the results produced by IOLB, we manually generate tiled versions of each kernel, then manually compute parametric data-movement costs as a function of tile sizes and cache size, then manually find the optimal tile sizes and thereby, finally, derive a manually optimized data-movement cost for this kernel. By forming the ratio of the total number of operations and the data-movement cost, we then generate OI manual . In this derivation, we assume that we have explicit control of the cache. Then OI manual is compared with an operational intensity upper-bound obtained by forming the ratio of the number of operations and the data movement lower bound generated by IOLB: OI up .

IOLB runs in less than one second on each of these kernels on a standard computer. Let us use jacobi-1d as an example to illustrate this. IOLB computes a lower bound expression Q low on the number of loads needed for any schedule of the jacobi-1d kernel:

Q low = 2 + N + max 0, T N 4S -N -T - 1 4 N S - 3 4
T S -S + 5 .
void k e r n e l _ j a c o b i _ 1 d ( i n t t s t e p s , i n t n , double

A [ n ] , double B [ n ] ) { i n t t , i ; #pragma scop f o r ( t = 0 ; t < t s t e p s ; t ++) { f o r ( i = 1 ; i < n -1 ; i ++) B [ i ] = 0.33333 * ( A [ i -1] + A [ i ] + A [ i + 1 ] ) ; f o r ( i = 1 ; i < n -1 ; i ++) A [ i ] = 0.33333 * ( B [ i -1] + B [ i ] + B [ i + 1 ] ) ; } #pragma endscop }    * (B[i-1] + );
small LB tool -small dimensions The first term is the input data size, and the second term is obtained by the partitioning technique. Since the expression of Q low can be quite large, we automatically simplify to Q ∞ low by only retaining the asymptotically dominant terms, assuming all parameters N , M . . . and cache size S tend to infinity, and S = o(N , M, . . . ). Finally, from Q ∞ low and the fact that the jacobi-1d kernel performs 6T N operations, we compute an upper bound for the OI of any schedule of the jacobi-1d kernel:

Inputs: n

Q ∞ low = T N 4S OI up = 6T N Q ∞ low = 24S
3.6.3 Parametric bounds for OI # input data is high, and we manually find that high OI is achievable through tiling. IOLB gives a non-trivial OI upper bound that is within a constant of the manually obtained OI lower bound OI manual . The bound is asymptotically tight for 9 of them, and within a factor of 2 for an additional 5. Except for matrix multiplication (gemm), where it matches the best published bound, these are all improvements over previously published results.

(7 kernels) The ratio # ops

# input data is a constant: clearly, these cases do not provide enough operations to enable data reuse. The reported lower bound by IOLB is # input data, which is asymptotically tight for 5 of them, and within a factor of 2 for 1 more.

(2 kernels)

The ratio # ops # input data is high which could indicate potential for tiling and high OI. However our best manual schedule leads to a constant OI which is arbitrarily far from this optimistic ratio. IOLB proves that the code is not tileable, so the best achievable OI is a constant. IOLB finds this upper bound on OI thanks to the wavefront technique. This is better by at least a factor of √ S than any bound that could be obtained by geometric reasoning.

(2 kernels)

There is an arbitrarily large discrepancy between OI up and OI manual . Visual examination shows that, for these cases, IOLB is too optimistic. These codes are actually not tileable in all dimensions, and we believe that it is possible, using more advanced techniques that are currently out of the scope of IOLB, to prove a matching OI upper bound.

Complete lower bound formulae

Table 3.2 shows the complete formulae produced by IOLB. Next to the complete formulae are asymptotic formulae as presented in Table 3

.1. kernel # input data # ops Q ∞ low OI manual ≤ OI ≤OI up ratio Published OI up 2mm N i N k + N k N j 2(N iN j N k 2(N i N j N k √ S ≤ OI ≤ √ S - +N j N l + N i N l +N i N j N l ) +N i N j N l )/ √ S 3mm N i N k + N k N j 2(N iN j N k + N j N l N m 2(N i N j N k + N i N j N l √ S ≤ OI ≤ √ S - +N j N m + N m N l +N i N j N l ) +N j N l N m )/ √ S cholesky 1 2 N 2 1 3 N 3 1 6 N 3 / √ S √ S ≤ OI ≤2 √ S 8 √ S [2] correlation MN M 2 N 1 2 M 2 N / √ S √ S ≤ OI ≤2 √ S - covariance MN M 2 N 1 2 M 2 N / √ S √ S ≤ OI ≤2 √ S - doitgen N p N q N r + N 2 p 2N 2 p N q N r 2N 2 p N q N r / √ S √ S ≤ OI ≤ √ S - fdtd-2d 3N x N y + T 11N x N y T 2 3 √ 3 N x N y T / √ S 11 24 √ 3 √ S ≤ OI ≤ 33 2 √ 3 √ S 36 - floyd-warshall N 2 2N 3 2N 3 / √ S √ S ≤ OI ≤ √ S 8 √ S [2] gemm N i N j + N j N k + N i N k 2N i N j N k 2N i N j N k / √ S √ S ≤ OI ≤ √ S √ S [48] heat-3d N 3 30N 3 T 3 8 3 √ 2N 3 T / 3 √ S 5 2 3 √ S ≤ OI ≤40 • 2 2/3 3 √ S 16 • 2 2/3 - jacobi-1d N 6N T 1 4 N T /S 3 2 S ≤ OI ≤24S 16 48S [22] jacobi-2d N 2 10N 2 T 2 3 √ 3 N 2 T / √ S 5 4 √ S ≤ OI ≤15 √ 3 √ S 12 √ 3 40 √ 2 √ S [22] lu N 2 2 3 N 3 2 3 N 3 / √ S √ S ≤ OI ≤ √ S 8 √ S [2] ludcmp N 2 2 3 N 3 2 3 N 3 / √ S √ S ≤ OI ≤ √ S 8 √ S [2] seidel-2d N 2 9N 2 T 2 3 √ 3 N 2 T / √ S 9 4 √ S ≤ OI ≤ 27 √ 3 2 √ S 6 √ - symm 1 2 M 2 + 2MN 2M 2 N 2M 2 N / √ S √ S ≤ OI ≤ √ S 8 √ S [2] syr2k 1 2 N 2 + 2MN 2MN 2 MN 2 / √ S √ S ≤ OI ≤2 √ S 8 √ S [2] syrk 1 2 N 2 + MN MN 2 1 2 MN 2 / √ S √ S ≤ OI ≤2 √ S 8 √ S [2] trmm 1 2 M 2 + MN M 2 N M 2 N / √ S √ S ≤ OI ≤ √ S 8 √ S [2] atax MN 4MN MN 4≤ OI ≤4 - bicg MN 4MN MN 4≤ OI ≤4 - deriche HW 32HW HW 16 3 ≤ OI ≤32 - gemver N 2 10N 2 N 2 5≤ OI ≤10 - gesummv 2N 2 4N 2 2N 2 2≤ OI ≤2 - mvt N 2 4N 2 N 2 4≤ OI ≤4 - trisolv 1 2 N 2 N 2 1 2 N 2 2≤ OI ≤2 - adi N 2 30N 2 T N 2 T 5≤ OI ≤30 - durbin N 2N 2 1 2 N 2 2 3 ≤ OI ≤4 - gramschmidt MN 2MN 2 MN 2 / √ S 1≤ OI ≤2 √ S 2 √ S - nussinov 1 2 N 2 1 3 N 3 1 6 N 3 / √ S 1≤ OI ≤2 √ S 2 √ S - Table 3.1 -Results on PolyBench benchmarks
We present the complete formulae produced by IOLB for a few reasons. While the lower bounds on I/O obtained by IOLB are lower bounds for any values of the parameters (M, N , S, etc), the asymptotic formulae have to be used with care. Indeed, the asymptotic reasoning, while providing simpler and easier to understand formulae, unfortunately removes the lower bound property. (Negligible negative terms are removed during the asymptotic reasoning.) Furthermore, if the asymptotic assumptions are violated, then the asymptotic formula becomes really off. For example if S is not negligible with respect to M and N , or if one dimension in gemm is small. Also, the asymptotic reasoning entails some assumptions. We chose to assume all parameters N , M . . . and cache size S tend to infinity, and S = o(N , M, . . . )). We can imagine other reasonable asymptotic reasoning. With the complete formula, it is possible to derive them at will. Finally, it can be instructive to understand the form of the complete formulae returned by IOLB, and showing the complete next to the asymptotic expansion explains our asymptotic assumption.

Tensor contractions and convolutions This table also shows parametric lower bounds for tensor contraction and 2D convolution kernels, as an application of the extensions presented in 3.3.4.

Because the lower bound derived is still sound, even if the small dimension hypothesis is not satisfied, we have to combine various lower bounds for different small dimension scenarios together.

For convolution kernels, we considered 5 small dimension scenarios, based on array sizes: (i) no small parameters, (ii) H and W small parameters, (iii) H, W and B small parameters, (iv) H, W , X, Y and B small parameters, and (v) C, H, W and B small parameters. Among those bounds, only the first three scenarios lead to interesting bounds. They correspond to the last three rows in the expression. The previous version of the algorithm (without reduction management or small dimensions) fails to find an interesting bound, and returns the sum of array sizes (first row in the expression).

Related work

The seminal work of Hong & Kung [START_REF] Hong | I/O complexity: The Red-Blue Pebble Game[END_REF] was the first to present an approach to derive lower bounds on data movement for any valid execution schedule of operations in a computational DAG. Their work modeled data movement in a two-level memory hierarchy and presented manually derived decomposability factors (asymptotic order complexity, without scaling constants) for a few algorithms like matrix multiplication and FFT. Several efforts have sought to build on the fundamental lower bounding approach devised by Hong & Kung, usually targeting one of two objectives: i) generalizing the cost model to more realistic architecture hierarchies [START_REF] Savage | Extending the Hong-Kung Model to Memory Hierarchies[END_REF][START_REF] Bilardi | A characterization of temporal locality and its portability across memory hierarchies[END_REF][START_REF] Bilardi | A Lower Bound Technique for Communication on BSP with Application to the FFT[END_REF], or ii) providing an I/O complexity with (tight) constant for some specific class of algorithms (sorting/FFT [START_REF] Aggarwal | The Input/Output Complexity of Sorting and Related Problems[END_REF][START_REF] Ranjan | Strong I/O Lower Bounds for Binomial and FFT Computation Graphs[END_REF], relaxation [START_REF] Ranjan | Upper and Lower I/O Bounds for Pebbling r-Pyramids[END_REF], or linear algebra [START_REF] Irony | Communication Lower Bounds for Distributed-Memory Matrix Multiplication[END_REF]4,3,[START_REF] Demmel | Communicationoptimal Parallel and Sequential QR and LU Factorizations[END_REF]).

In the context of linear algebra, Irony et al. [START_REF] Irony | Communication Lower Bounds for Distributed-Memory Matrix Multiplication[END_REF] were the first to use the Loomis-Whitney inequality [START_REF] Loomis | An inequality related to the isoperimetric inequality[END_REF] to find a lower bound on data movement. This was in the context of gemm (one of the kernels of PolyBench). The asymptotic upper bound on OI from this paper is 4

√ 2 √ kernel Complete lower bound Asymptotic simplified formula 2mm max N i N j + N j N k + N i N k + N j N l + 2, 2 √ S N i N j (N k -1) + 2N i + 2N j + N k -4 √ 2S 2 √ S N i N j N k + 2 √ S N i N l (N j -1) + 2N i + 2N l + N j -4 √ 2S + 2 √ S N i N l N j -2N i N j -2 3mm max N i N k + N j N k + N j N m + N l N m , 2 √ S N i N j (N k -1) + 2N i + 2N j + N k -4 √ 2S 2 √ S N i N j N k + 2 √ S N i N l (N j -1) + 2N i + 2N l + N j -4 √ 2S + 2 √ S N i N l N j + 2 √ S N j N l (N m -1) + 2N j + 2N l + N m -4 √ 2S + 2 √ S N j N l N m -2N j N i -2N j N l -N i N l -6 adi 4N 2 + max 0, (N 2 -4N -S + 5)(T -2) N 2 T atax MN + N + max 0, 1 8 1 S ((2M -1 -8S)(2N -1 -8S) -1) -10S + 2 MN bicg MN + M + N + max 0, 1 8 1 S ((2M -1 -8S)(2N -1 -8S) -1) -10S + 2 MN cholesky max 1 2 N (N + 1), 1 6 1 √ S (N -1)(N -2)(N -3) + 1 2 √ 2 1 S (N -1)(N -2) 1 6 1 √ S N 3 -(N -2)(N -7) -4 √ 2S correlation max MN + 2, 1 2 1 √ S M(M -1)(N -1 + √ 2 2 1 √ S ) -1 2 (M -3)(M + 2N -2) + 2 -4S √ 2 1 2 1 √ S M 2 N covariance max MN + 2, 1 2 1 √ S M(M -1)(N -1 + √ 2 2 1 √ S ) -1 2 (M -3)(M + 2N -2) + 1 -4S √ 2 1 √ S M 2 N deriche HW + 1 HW doitgen max N 2 p + N p N q N r , 2 √ S N q N r N p (N p -1 + 1 √ 2 1 √ S ) -N q N r (N p -1) + 2N p -8 √ 2S -1 2 1 √ S N q N r N 2 p durbin 2N + max 0, 1 2 (N -3)(N -2 -2S) N 2 fdtd-2d max 3N x N y -N y + T -1, 1 2 √ 2 1 √ S (N x -2)(N y -2)(T -1) + 2(N x + 2)(N y + 2) 1 √ 2 1 √ S N x N y T -T (N x + N y -6) -N y -S -23 floyd-warshall max N 2 , 1 √ S (N -1) 3 -(6N -19)(N -2) -8 √ 2S 1 √ S N 3 gemm max N i N j + N j N k + N i N k + 2, 2 √ S N i N j (N k -1) + 2N i + 2N j + N k -4 √ 2S 2 1 √ S N i N j N k gemver N 2 + 8N + 2 + max 0, 1 4 1 S (3N -2)(N -8S) -3S + 1 N 2 gesummv 2N 2 + N + 2 + max 0, 1 2 1 S (N -1)(N -8S) -2S 2N 2 gramschmidt max MN , 1 √ S MN (N -3) -M(N -5 -2 √ S ) -1 2 (N -1)(N -6) -4 √ 2S -3 1 √ S MN 2 heat-3d max (N -10)(N + 2) 2 , 9 3 √ 3 16 1 3 √ S (T -1)(N -3) 3 -3(T -7)(N -3)(N -4) 3 √ 3 16 1 3 √ S N 3 T +42N -T -9 3 √ 3 4 3 √ 4 S -111 jacobi-1d max 2 + N , 1 4
1

S (T -1)(N -3) -T -S + 7 1 S N T jacobi-2d max (N -2)(N + 6), 2 3
√ 3 1 √ S (N -3) 2 (T -1) -4 √ 2 3 √ 3 S -(T -7)(2N -7) + 14 2 √ 3 1 √ S N 2 T lu max N 2 , 2 3 1 √ S (N -2)(N 2 -4N + 6) -2(N 2 -10N + 18) -8 √ 2S 1 √ S N 3 ludcmp max N 2 + N , 1 3 1 √ S (2N -3)(N -1)(N -2) √ 2 1 S (N -1)(N -2) -(2N 2 -15N + 19) -16 √ 2S mvt N 2 + 4N + max 0, 1 6 1 S N (N -1) -2S -4N + 4 N 2 nussinov 1 2 N 2 + 5 2 N -1 + max 0, 1 6 1 √ S (N -3)(N -4)(N -5) + 1 4 1 S √ 2(3N 2 -19N + 6) 1 6 1 √ S N 3 -(N 2 -13N + 22) -8 √ 2S seidel-2d max N 2 , 2 3
√ 3 1 √ S (N -3) 2 (T -1) -(2N -7)(T -5) -4 √ 2 3 √ 3 S + 12 2 3 √ 3 1 √ S N 2 T symm max 1 2 M(M + 1) + 2MN + 2, 2 1 √ S (M -1)(M -2)N -1 2 ((4N + M)(M -5)) 2 1 √ S M 2 N +5(M -2) -8 √ 2S syr2k max 2 + 2MN + 1 2 N (N + 1), 1 √ S (M -1)(N + 1)N + M + 4N -4 √ 2S 1 √ S MN 2 syrk max MN + 1 2 (N + 1)N + 2, 1 2 1 √ S (M -1)(N + 1)N -(M -4)(N -1) -2 √ 2S + 4 1 2 1 √ S MN 2 trisolv 1 2 N (N + 1) + N + max 0, 1 8 1 S (N -1)(N -2) -2N -S + 5 1 2 N 2 trmm max 1 2 M(M -1) + MN + 1, 1 √ S (M -2 + √ 2 √ S )(M -1)N -(M -4)(N -2) -8 √ 2S + 5 1 √ S M 2 N TC abcde-efbad-cf max ABCDE + EFBAD + CF, -2 + F -S + 2C + 2ABDE + 2ABCDE(F-1) √ S 2 √ S ABCDEF TC abcd-dbea-ec max ABCD + DBEA + EC, -2 + E -S + 2C + 2ABD + 2ABCD(E-1) √ S 2 √ S ABCDE TC abc-bda-dc max ABC + BDA + DC, -2 + D -S + 2C + 2AB + 2ABC(D-1) √ S 2 √ S ABCD TC abcdef-dega-gfbc max ABCDEF + DEGA + GFBC, -2 + G -S + 2ADE + 2BCF + 2ABCDEF(G-1) √ S 2 √ S ABCDEFG TC abc-adec-ebd max ABC + ADEC + EBD, -3 + DE -S + 3AC + 3B -ABC + TC abcd-aebf-fdec max ABCD + AEBF + FDEC, -3 + EF -S + 3AB + 3CD -ABCD + 2ABCD(EF-1) √ S 2 √ S ABCDEF conv-2d-s1 max BC(Y + H -1)(X + W -1) + BFXY + FCHW , -2 -S + C + 4F + BY + BX + 2BXY -2BXY F + BFXY (W HC-1) S , -2 -S + C + 4F + BY + BX + 2BXY -2BXFY + 2BXY CF √ HW √ S -2BXY F √ HW S , -2 -S + C + 4F + BY + BX + 2BXY -2BXFY + 2XY CF √ BHW √ S -2XY F √ B √ HW S 2 √ S
BXY CF √ HW Table 3.2 -Complete Lower Bound Formulae for PolyBench, TCCG and convolutions obtained with IOLB method presented in [START_REF] Ballard | Communication lower bounds and optimal algorithms for numerical linear algebra[END_REF] is limited to a few algorithms. Kwasniewski et al. [START_REF] Kwasniewski | Red-blue Pebbling Revisited: Near Optimal Parallel Matrix-Matrix Multiplication[END_REF] implemented an algorithm for parallel matrix-matrix multiplication that matches the communication lower bound for any combinations of matrix dimensions, processors counts and memory sizes. See discussion on [START_REF] Christ | Communication Lower Bounds and Optimal Algorithms for Programs that Reference Arrays -Part 1[END_REF] for more details on these limitations. The studies that are the most related to this chapter are those from Christ et al. [START_REF] Christ | Communication Lower Bounds and Optimal Algorithms for Programs that Reference Arrays -Part 1[END_REF], and Elango et al. [START_REF] Elango | On characterizing the data movement complexity of computational DAGs for parallel execution[END_REF][START_REF] Elango | On Characterizing the Data Access Complexity of Programs[END_REF].

The idea of using a variant of the Brascamp-Lieb inequality to derive bounds for arbitrary affine programs comes from Christ et al. [START_REF] Christ | Communication Lower Bounds and Optimal Algorithms for Programs that Reference Arrays -Part 1[END_REF]. However, the approach they propose suffers from several limitations: 1. The model is based on association of operations with data elements and does not capture data dependencies in a computational DAG. Consequently, it can lead to very weak lower bounds on data movement for computations such as Jacobi stencils. 2. There is no way to (de-)compose the CDAG, and they view all the statements of the loop body (that has to be perfectly nested) as an atomic statement. As a consequence, it is incorrect to use this approach for loop computations where loop fission is possible. 3. The lower bounds modeling is restricted to S-partitioning, leading to very weak lower bounds for algorithms such as adi or durbin. 4. Obtaining scaling constants, in particular with nonorthogonal reuse directions, is difficult, and only asymptotic order complexity bounds can be derived. 5. No automation of the lower bounding process was proposed, but manually worked out examples of asymptotic complexity as a function of fast memory capacity (without scaling constants) were presented.

Elango et al. [START_REF] Elango | On characterizing the data movement complexity of computational DAGs for parallel execution[END_REF] used a variant of the red-blue pebble game without recomputation, enabling the composition of several sub-CDAGs, and the use of a lower-bounding approach based on wavefronts in the DAG. Manual application of the approach for parallel execution was done on specific examples, but no approach to automation was proposed.

The later work of Elango et al. [START_REF] Elango | On Characterizing the Data Access Complexity of Programs[END_REF] was the first to make the connection between paths in the data-flow graph and regular data reuse patterns and to propose an automated compiler algorithm for affine programs. However, their proposed approach suffers from several limitations: 1. Only asymptotic Ω(. . .) data movement bounds were obtainable, without any scaling constants. In contrast, IOLB generates meaningful non-asymptotic parametric I/O lower bound expressions. From these expressions, we can derive asymptotic lower bounds with scaling constants, critical for use in deducing upper limits on OI for a roofline model. 2. Since they were only trying to provide asymptotic bounds without constants, they did not address (de-)composition (asymptotic bounds can be safely summed up even if they interfere). Also, they only considered enumerative decomposition, and not dimension decomposition through loop parameterization that is necessary to obtain a tight bound for their Matmult-Seidel illustrative example. They also only considered the simple non-overlapping notion of interference, and did not allow decomposition of the same statement, required in order to obtain a tight bound for computations like floyd-warshall. 3. Finally, their approach only used the S-partitioning paradigm for lower bounds but not the wavefront-based paradigm, thus leading to very weak bounds for benchmarks such as adi or durbin.

Chapter 4 Upper Bounds

We have presented a method for automatically computing data movement lower bounds for affine programs. This can be very useful to assess the optimization potential of a program, but a first step is to evaluate its actual data movement cost. This raises two questions: 1. how can we estimate the I/O cost of a program under a given schedule, as a parametric expression of parameters and cache size? 2. how do we find an I/O-efficient schedule for a program? Answering these questions can help us assess the tightness of a data movement lower bound, with the best scenario being matching lower and upper bounds: when this is the case, we have found an I/O-optimal schedule. This chapter tries to tackle this challenge for a large class of programs, and presents a method for automating the process, that was implemented in a tool called IOUB.

It includes the first algorithm for computing a symbolic over-approximation of the data movement for a parametric (multi-dimensional) tiled version of an affine code, as well as the first fully automated scheme for expressing as an operations research problem the minimization of this data movement expression.

The chapter is organized as follows. Sec. 4.1 provides some useful formalism along with the required background. Sec. 4.2 defines the tiling transformation, and the associated domains and quantities that are used in the analysis. Sec. 4.3 describes the I/O cost model, and the algorithm that we designed to compute a symbolic expression of the I/O cost of a tiled program. Sec. 4.4 explains how we reduce the search space for tiling loop permutations. Sec. 4.5 connects the previous sections into a complete framework. Sec. 4.6 reports how it was implemented into a tool called IOUB, and reports the results of our experiments on tensor contractions and convolutions. Finally, Sec. 4.7 discusses related work.

f o r ( i = 0 ; i < Ni ; i ++) f o r ( j = 0 ; j < Nj ; j ++) f o r ( k = 0 ; k < Nk ; k ++) C [ i ] [ j ] += A [ i ] [ k ] * B [ k ] [ j ] ;
(a)

f o r ( i 1 = 0 ; i 1 < Ni ; i 1 += T i ) f o r ( j 1 = 0 ; j 1 < Nj ; j 1 += T j ) f o r ( k = 0 ; k < Nk ; k ++) f o r ( i = i 1 ; i < i 1 + T i ; i ++) f o r ( j = j 1 ; j < j 1 + T j ; j ++) C [ i ] [ j ] += A [ i ] [ k ] * B [ k ] [ j ] ; (b)

Background

In this section, we present the class of programs, its representations, and the memory model we consider in our analysis.

Class of programs

Imperfectly nested affine loop programs We consider imperfectly nested loop programs, which can be recursively defined as a sequence of statements and for loops, each being itself composed of imperfectly nested loops. A program parameter is a symbolic constant during the compilation, whose value will be known during the execution of the program (for example, array sizes). The access functions made by statements to arrays are affine expressions of surrounding loop indices and program parameters. We also assume that conditions on loop indices are affine expressions of surrounding loop indices and program parameters [START_REF] Bastoul | Putting Polyhedral Loop Transformations to Work[END_REF].

For example, the 1D-convolution described in Figure 4.2 matches these criteria. This program is even perfectly nested: all loops are nested with a single statement inside. Constants Nc, Nf, Nx and Nw are the program parameters. All loop indices are bounded by affine constraints (e.g. 0 ≤ c < Nc), and all array subscripts are affine expressions (e.g. (x + w, c) for array Image).

Fully-tilable program

The tiling transformation [START_REF] Xue | Loop Tiling for Parallelism[END_REF] is a key loop transformation to improve the data locality of a program. Given a list of consecutive dimensions (called loop band), this transformation groups their iterations into tiles, which are executed atomically. Figure 4.1 shows an example of tiling transformation for the matrix multiplication example: the tiled dimensions are i and j and the iterations are grouped into rectangles of size Ti × Tj (tile shape). Note that when the tiling shape is rectangular, each tiled dimension (e.g. i) is stripmined, giving rise to the tile dimension that will iterate over the tiles (e.g. i1), and the local dimension that iterates inside a tile (e.g. i in the transformed program).

A tiling is legal when there is no cycle of dependencies between the computation of different tiles, that is, when the atomicity condition between tiles can be respected by the schedule of the program. The algorithms presented in this chapter are described on the assumption that the input program is fully tilable using rectangular tiles, which means that all its dimensions can be legally tiled by using rectangular tiles.

In theory, any affine program could be pre-processed using a polyhedral compiler to provide a fully permutable (that is, tilable using rectangular tiles) loop band, but the reality is more complex, as analyzing a non-regular iteration domain involves being able to cope with the simplification and solving of a complex system of symbolic expressions. Nevertheless, several important classes of computation fit our simplified hypothesis, such as a convolution, all kinds of tensor contractions and many linear algebraic kernels, including matrix mul-tiplication. Figure 4.2 shows another example of such computation, which is a simplified convolution. It will be used as a running example for this chapter.

Program representation

For theoretical reasoning, we still use the CDAG abstraction for programs. However representing a full graph is not a practical way of handling programs, so our algorithms manipulate a more convenient intermediate representation. The formalism we use in this chapter is slightly different from the DFG we use for lower bounds, as the focus is more on arrays than on statements, but the underlying framework is the same. Let us formally define this representation.

The iteration domain I S of a statement S is the set of integral values that the surrounding loop indices take during execution. Each point of this space is associated with an execution instance of the statement. Due to the hypotheses made on the loop bounds, the iteration domain of a statement is a Z-polyhedron, that is, a set of integral points whose coordinates satisfy a set of affine constraints. The computation of the cardinality of a set E (denoted |E|symbolic expression as a function of the program parameters), can be efficiently done using the Barvinok algorithm1 [START_REF] Barvinok | A Polynomial Time Algorithm for Counting Integral Points in Polyhedra when the Dimension is Fixed[END_REF].

Each array A is associated with a memory domain M A , which is the multidimensional set of valid array indices for the array A. Given a statement S containing a read or a write to an array A, this occurrence is associated with a memory access function f A : D S → M A , which is an affine function. For any given sub-domain D S ⊂ I S defined as a Z-polyhedron, the associated data footprint f A (D S ) is also a Z-polyhedron that can be easily computed by polyhedral compilation tools.

For the matrix-matrix multiplication example (Figure 4.1, untiled version), the dimensions are D = {i, j, k} and the iteration domain is Tile limit contains five loops on dimensions i, j, k, i, j, from outer to inner. The three outermost loops span the entire iteration domain of size Ni × Nj × Nk, while the two innermost ones span a polyhedron of size Ti × Tj. We refer to the innermost (resp. outermost) part as the intra-tile (resp. inter-tile) dimensions. As we will see later, unlike the permutation of the inter-tile dimensions, that of the intra-tile dimensions is not relevant in our model. In other words, while a different schedule with loops ordered as (from outer to inner) (i1,k,j1,i,j) could lead to a different I/O estimation than for the schedule with loops ordered as (i1,j1,k,i,j) (as in Figure 4.1), our cost model will not be affected if we permute i and j as in loop order (i1,j1,k,j,i). This motivates the following notation to represent the tiling schedule of

I = {(i, j, k) | 0 ≤ i < N i ∧ 0 ≤ j < N j ∧ 0 ≤ k < N k }.

Sub-domains and reuse

To model the I/O cost of a given tiling schedule, we first need some formal definitions.

For an inter-tile loop dimension d j , we define its sub-domain as the set of points in the iteration domain such that the indices of the enclosing loops d k (|D| ≥ k ≥ j) have fixed values. Definition 15. Let d j be some inter-tile loop dimension. We denote by ik the fixed index value at level k. The sub-domain SD d j at level j is defined as:

SD d j (iD, . . . , ij) = {i |D| , . . . , i 1 ∈ I | ∀k ≥ j, ik ≤ i k < ik + T d k }
used in the previous sub-domain and the corresponding I/O can thus be saved. This means that the optimal inverse operational intensity for any sub-domain in the front (resp. in the back) is no more than ID front A,l (resp. ID back A,l ).

Multiple arrays

Now when there are multiple arrays A 1 , . . . , A s , we "cut" the cache into array-specific regions of sizes S 1 , . . . , S s , such that S 1 + • • • + S s = S. The idea is then simply to do the same as before, with each region holding only data from its array. This is slightly restrictive, as one could imagine the optimal I/O schedule using some memory space alternatively for different arrays. But it is necessary to make the problem tractable, and does not prevent IOUB from giving optimal or close to optimal bounds in many cases.

The outermost reuse dimension for array A i is defined is the same way as in the one array case: Definition 19. The outermost reuse dimension l i for array A i is the leftmost index in P such that SDF A,l i ≤ S i .

The result is straightforwardly derived from the one array case, and is summed up in the following lemma:

SDF A i ,1 ≤ S i . (4.1)
Then the following expression is an upper bound on the cost of a red-white pebble game for this program:

IO = IO A 1 + • • • + IO A s
where

IO A i = ID front A i ,l i × I front + ID back A i ,l i × |I back | (4.2)
Conditions (4.1) can be replaced by their sum, avoiding the need to explicitly define S i s:

i SDF A,l i ≤ S.
However this notation makes the assumptions more explicit, and makes it possible the cache partition when solving the equations numerically.

Example On our running example, for array Image, supposing the outermost reuse dimension is x:

IO Image = ID front Image,1 × |I front | + ID back
SDF Image,1 (c1,f1,x) = Nw × Tc ≤ S Image SDF Out,1 (c1,f1,x) = Tf ≤ S Out SDF Filter,2 (c1,f1) = Tf × Nw × Tc ≤ S Filter
This leads to the following bound:

IO = Nc.Nf.(Nx + Nw -1) Tf + Nc.Nf.Nx Tc + Nc.Nf.Nw Assuming Nw × Tc + Tf + Tf × Nw × Tc ≤ S.

Extension to multiple memory levels

The red-white pebble game and the corresponding lower bound reasoning frameworks can be extended to multi-level memory hierarchies [START_REF] Savage | Extending the Hong-Kung Model to Memory Hierarchies[END_REF]. The above computation of I/O can also be extended by simply considering one tiling band per cache level and independently applying the previous reasoning to each level. One only needs to add constraints on tile sizes so that tiles get bigger at each level.

The question that arises is: what is the I/O cost of a multi-level tiling? Our model can provide a cost for each memory level, but how to aggregate these costs is not straightforward. In a real-life CPU with multiple levels of cache, it is likely that the I/O cost is dominated by the level where the cost is the closest to the bandwidth. In this work, we chose to use a weighted sum, where the cost at each level is multiplied by the inverse bandwidth.

Optimization problem

Having a symbolic expression as a function of the tile sizes makes it possible to express the problem of choosing tile sizes as a non-linear optimization problem (NLP). Thus, for fixed values of S and program parameters, an NLP solver such as IPOPT [START_REF] Wächter | On the Implementation of an Interior-Point Filter Line-Search Algorithm for large-scale nonlinear programming[END_REF] can be used to find optimizing tile sizes for this particular permutation. One still needs to be cautious, as the tile sizes found by such tools are not guaranteed to be integers, and rounding can have an impact on the I/O cost. This is only valid for a fixed permutation of tiling loops, as well as fixed outermost reuse dimensions for each array. Indeed, since tile sizes are variables, conditions (4.1) come as constraints that are verified afterwards. The next section explains how to choose these parameters.

Loop permutation selection

We have so far only focused on computing the required I/O cost for a given permutation of the inter-tile loops. Unfortunately, the search space of all possible permutations grows exponentially with the number of dimensions and with the memory depth. However, it is straightforward to see that many permutations are equivalent in terms of I/O cost: for example, switching the two outer loops c1 and f1 in the conv-1D example in Figure 4.3 would not change the I/O costs. Moreover, some permutations are better than others: for example, it is not worth exploring the permutation whose innermost dimensions do not allow any data reuse between two successive tiles. In this section, we provide some insights into how to select a subset of relevant permutations for a given code.

Reuse for an array along a dimension

First, we define a notion of reuse: for a given dimension d and an array A, there is reuse for array A on dimension d if, when putting dimension d innermost (that is, setting d 1 = d), the sub-domain data footprint at level 2 for A does not increase much compared to level 1:

SDF A,2 -SDF A,1 ≪ SDF A,1
This allows us to decide which dimensions are worth putting at the innermost levels.

In simple cases, this notion is not ambiguous: for array Out in the conv-1D example, c is a reuse dimension as setting it innermost would lead to SDF Out 2 = SDF Out 1 = Tf × Tx, while f is not a reuse dimension as setting it innermost would lead to SDF Out 2 -SDF Out 1 = (Nf -Tf) × Tx (unless Nf = Tf but in that case it means that f iterates only once). For more complex subscript expressions, such as for array Image, the criterion is less clear: setting w as innermost and assuming it is not part of the tile would lead to comparing Tw -1 with Tx (as SDF Image 2 = (Tx + Tw -1) × Tc, SDF Image 1 = Tx × Tc). In other words, a reuse criterion requires an "oracle" (such as the user) to provide information about small and large dimensions.

This also provides information on how to choose outermost reuse dimensions, defined in the previous section. For an array A, this dimension will be the highest one such that all lower dimensions are reuse dimensions for array A.

Algorithm for permutation selection

Assuming such an oracle and a reuse criterion, Algorithm 7 selects a set of representative permutations. This algorithm builds a permutation from the innermost to the outermost dimensions while keeping track of the set of remaining dimensions (variable D ′ ), and of the set S of arrays having a potential reuse along each dimension d (variable R). When a dimension is chosen, the set of potential reuse is updated. When there is no more potential reuse, an arbitrary permutation of the remaining dimensions is chosen.

Figure 4.4 illustrates the steps of the algorithm on the 1D convolution (Figure 4.2), and how it narrows the loop permutation space to three permutations (one of them is the one used for the running example in Figure 4.3). Note that some permutations are pruned during the selection process (red cross in Figure 4.4), because they have strictly less reuse potential than others. For example, selecting dimension c as the innermost dimension allows us to have reuse on array Out, but dimension w allows us to have reuse on one more array (Out and Image).

Putting it all together

The IOUB tool combines the permutation selection and I/O cost computation as follows: it first selects a list of permutations and generates the corresponding set of tiled versions (with

P := ∅; forall (d, S) ∈ R such that ∄(d ′ , S ′ ) ∈ R, S ⊊ S ′ do D ′ d := D ′ \ {d}; R d := {(d ′ , S ′ ∩ S), (d ′ , S ′ ) ∈ R \ {(d,

Parameter values Operations research

Computer algebra

f o r ( i = 0 ; i < Ni ; i ++) f o r ( j = 0 ; j < Nj ; j ++) f o r ( k = 0 ; k < Nk ; k ++) C [ i ] [ j ] += A [ i ] [ k ] * B [ k ] [ j ] ; {(i, j, k), (i, k, j), (k, j, i)} f o r ( i 1 = 0 ; i 1 < Ni ; i 1 += T i ) f o r ( j 1 = 0 ; j 1 < Nj ; j 1 += T j ) f o r ( k = 0 ; k < Nk ; k ++) f o r ( i = i 1 ; i < i 1 + T i ; i ++) f o r ( j = j 1 ; j < j 1 + T j ; j ++) C [ i ] [ j ] += A [ i ] [ k ] * B [ k ] [ j ] ; IO = N i N j N k 1 T i + 1 T j + 1 N k T i T j + T i + T j ≤ S U B = N i N j 2N k √ S+1-1 + 1
Figure 4.5 -Summary of the full framework of IOUB, illustrated on matrix-matrix multiplication tile sizes as parameters) using Algorithm 7. For each permutation, it computes a symbolic expression (as a function of program parameters, tile sizes, and cache size) that represents the required I/O, using the method described in Sec. 4.2 and 4.3, as well as constraints on tile sizes. Then, this can be fed to an NLP solver to find the best tile sizes for each permutation. The final I/O cost is the minimum over all versions, and IOUB also provides a basic tiled code that implements the corresponding tiling scheme. It can also be useful to have a symbolic bound that does not depend on tile sizes, but only on program parameters and cache size. A discussion on how to derive such an expression for the examples we used in experiments can be found in Sec. 4.6.2. This is summarized in Fig. 4.5.

Experiments

Using our new I/O cost computation method, and our improvements to the lower bound algorithm for computing I/O complexity, IOUB is able to provide tight bounds for both tensor contractions (even with small dimensions) and convolutions. The general code for convolution is shown on Fig.

Benchmarks

We illustrate this by running it on representative benchmarks: For convolutions, we considered 11 different layers of Yolo9000 [START_REF] Redmon | YOLO9000: Better, Faster, Stronger[END_REF]. The parameter values for each layer are shown in Fig. 4.8. For tensor contractions, we considered the ones from TCCG [START_REF] Springer | Design of a High-Performance GEMM-like Tensor-Tensor Multiplication[END_REF].

The benchmark Python script from TCCG source code gathers 73 tensor contraction kernels (originating from various sources), that can be reduced to 49 different kernels, once duplicates are identified. For each kernel, TCCG selects the sizes of every dimension so that they are multiple of 8, roughly equal, and so that the product of all sizes is around 200 × 2 20 .

We can further reduce the number of relevant kernels that need to be considered in our analysis, by grouping them according to the number of dimensions of each array, and the number of dimensions shared between them. Indeed, the array layouts (dimension order) do not impact our analysis. This yields eight classes of tensor contraction kernels, described in Figure 4.9.

Symbolic upper bound expressions

The method presented in this chapter provides a symbolic expression for the I/O of a program, as a function of program parameters and tile sizes. To compare them with the parametric lower bound expressions, we would like to remove tile sizes from the expression, and express them as functions of the cache size instead. In general, this is too complicated to solve. However, for the benchmarks we consider, this is possible by using only a few assumptions.

Let us start with the matrix multiplication example from

+ CF TC abcd-dbea-ec UB = 2ABCDE √ S+1-1 + CE TC abc-bda-dc UB = 2ABCD √ S+1-1 + CD TC abcdef-dega-gfbc UB = 2ABCDEFG √ S+1-1 + BCFG TC abc-adec-ebd UB = 2ABCDE √ S+1-1 + BDE TC ab-cad-dcb UB = 2ABCD √ S+1-1 + AB TC ab-ac-cb UB = 2ABC √ S+1-1 + BC TC abcd-aebf-fdec UB = 2ABCDEF √ S+1-1 + CDEF 2D Convolution UB = CFHW XY 1 XY + 1 H∆W + (H+∆-1)(W +X-1) H∆ 2 W X where ∆ = -HW +W + √ H 2 W 2 +4HSW -2HW 2 +4SW +4S+W 2 2(HW +W +1)

64

The IO cost expression and the constraints on tile sizes are:

IO = IO A + IO B + IO C = Ni • Nj • Nk • 1 Ti + 1 Tj + 1 Nk (4.3) SDF A,1 + SDF B,1 + SDF C,1 = Ti + Tj + Ti • Tj ≤ S (4.4)
Then, we consider square tiles, by assuming that Ti and Tj are equal to the same value T . We also assume that the tile completely fills the cache, which means that we consider that inequality (4.4) is actually an equality. While this may not be the best solution, these hypotheses will still provide a valid bound. We obtain the following equality:

2T + T 2 = S.
It has a unique positive solution:

T = √ S + 1 -1.
We can then plug this value back into expression (4.3) to get our symbolic bound:

IO = Ni • Nj • 2Nk √ S + 1 -1 + 1
Note that the dominant term matches the dominant term of the I/O lower bound for matrix multiplication.

Here, we chose the tiling scheme manually, but this can be automated using the NLP solver with some relevant numerical values of parameters. It will return a permutation and numerical tile sizes, and we can use the information on which dimensions have T d = 1 or T d = N d to guide our symbolic solving. If the problem sizes are significantly bigger than √ S, the best tiling found by the solver will precisely be this one.

To generalize this reasoning to tensor contraction kernels, we only need to add some extra information on the expected order of magnitude of tile sizes.

Dimensions in a tensor contraction computation can be divided into three groups such that after "merging" the dimensions in each group, the computation is equivalent to a matrix multiplication. This corresponds to the "shared dimensions" in Fig. 4.9. The condition we impose is that the products of tile sizes inside each group are equal. For example abc-adecebd, the three groups are {a, c}, {b} and {d, e}, and the condition is

T a T c = T b = T d T e .
The expressions are shown in Figure 4.10. They correspond to the "general case" where the input size is not the bottleneck, i.e., when parameters are sufficiently large compared to S.

The expression for convolution is quite complex, but an asymptotic analysis shows that the highest order term is 2CFXY

√ HW √ S
when C, F, X and Y are sufficiently large, which matches the third term in the lower bound (the B factor does not appear since it is equal to 1 in all our benchmarks). As a sanity check, we confirm that the upper bound found is always above the lower bound, which was not trivial for the convolution bounds, and that both curves are decreasing or are constant along S. Both bounds are close to each other, ranging from at most a factor of 3 between them for Yolo9000-2 on the smallest value of S, to almost the same values for the ab-ac-cb TC kernel, which is also known as matrix multiplication. Moreover, the upper bound and the lower bound come closer for the largest values of S, which shows a match between the asymptotic dominant term of both bounds. When the cache size becomes very large, the dominant cost becomes the initial loading of the input data, hence the matching bounds and the horizontal line for some benchmarks.

Comparison of upper and lower bounds for different cache sizes

Implementation

IOUB is implemented as a command-line tool, written mostly in Python and using the ISL [START_REF] Verdoolaege | ISL: An Integer Set Library for the Polyhedral Model[END_REF] and Sympy [START_REF] Meurer | SymPy: Symbolic Computing in Python[END_REF] libraries, as well as the IPOPT [START_REF] Wächter | On the Implementation of an Interior-Point Filter Line-Search Algorithm for large-scale nonlinear programming[END_REF] NLP solver. Figure 4.12 show an example of input file and output from IOUB for a convolution kernel and a 3-level memory hierarchy. In its current version, the input format is a yaml file describing loops and arrays manually, but IOUB could be combined with polyhedral tools like PET [START_REF] Verdoolaege | Polyhedral Extraction Tool[END_REF] and PluTo [START_REF] Bondhugula | A Practical Automatic Polyhedral Program Optimization System[END_REF] to generate this intermediate representation directly from C code.

Related work

Computing an I/O complexity upper bound for an algorithm is the most reasonable way to assess the tightness of a lower bound. While this computation is usually done by hand using ad hoc techniques specific to each studied algorithm [START_REF] Olivry | Automated Derivation of Parametric Data Movement Lower Bounds for Affine Programs[END_REF][START_REF] Demmel | Communication-Optimal Convolutional Neural Nets[END_REF][START_REF] Aggarwal | The Input/Output Complexity of Sorting and Related Problems[END_REF][START_REF] Michael | A Tight I/O Lower Bound for Matrix Multiplication[END_REF][START_REF] Kwasniewski | Red-blue Pebbling Revisited: Near Optimal Parallel Matrix-Matrix Multiplication[END_REF][START_REF] Ranjan | Upper and Lower I/O Bounds for Pebbling r-Pyramids[END_REF], Fauzia et al. [START_REF] Fauzia | Beyond Reuse Distance Analysis: Dynamic Analysis for Characterization of Data Locality Potential[END_REF] proposed a heuristic that directly reasons on the CDAG, which unfortunately does not scale to real programs. Finding an upper bound for a fixed architecture can also be viewed as finding an optimized program transformation that minimizes data movement costs, which also implies being able to evaluate this cost. Thus, restricting the analysis to affine programs and using the polyhedral framework appears to be appropriate for this problem. However, while the seminal scheduling algorithm from PluTo [START_REF] Bondhugula | A Practical Automatic Polyhedral Program Optimization System[END_REF] is able to expose tilable loops and generate tiled code, it can only handle fixed tile sizes. This is because parametric tiling is not an affine transformation: a tiled affine code with parametric tile sizes is no longer affine. A consequence is that existing polyhedral tools cannot be used to evaluate the I/O cost of such a code. These tools include PolyFeat [START_REF] Bao | Static and Dynamic Frequency Scaling on Multicore CPUs[END_REF], which computes an approximation of the number of capacity misses for arbitrary affine programs, as well as other algorithms [START_REF] Bao | Analytical modeling of cache behavior for affine programs[END_REF][START_REF] Chatterjee | Exact Analysis of the Cache Behavior of Nested Loops[END_REF][START_REF] Ghosh | Cache Miss Equations: a Compiler Framework for Analyzing and Tuning Memory Behavior[END_REF] which focus on precisely modeling conflict misses. Some of these algorithms are restricted to a small class of programs, and most of them can only model a one-level cache. They all need to consider fixed parameters and fixed tile sizes. Other works [START_REF] Yuki | Automatic Creation of Tile Size Selection Models[END_REF][START_REF] Shirako | Analytical Bounds for Optimal Tile Size Selection[END_REF][START_REF] Li | Analytical Cache Modeling and Tilesize Optimization for Tensor Contractions[END_REF] implemented ad hoc computation of this cost function for very restricted sub-classes of programs, but our algorithm is the first to be able to automatically generate a symbolic expression of it for arbitrary parametric tiled affine programs, in a multi-level cache setting.

We showed how we use this cost function to find the best loop permutations and tile sizes for a given architecture, using operations research. Once again, this is out of the scope of polyhedral analysis, as cost functions are not affine. So the usual optimization strategies used by polyhedral compilers [START_REF] Grosser | Polly -Performing Polyhedral Optimizations on a Low-Level Intermediate Representation[END_REF][START_REF] Bondhugula | A Practical Automatic Polyhedral Program Optimization System[END_REF][START_REF] Verdoolaege | Polyhedral Parallel Code Generation for CUDA[END_REF], mostly based on parametric integer programming [START_REF] Feautrier | Parametric Integer Programming[END_REF], cannot be used. Some literature exists on tile size selection, but it mostly deals with performance model design and is either handcrafted for specific kernels and access function patterns [START_REF] Liu | Tile Size Selection for Optimized Memory Reuse in High-level Synthesis[END_REF][START_REF] Prajapati | Simple, Accurate, Analytical Time Modeling and Optimal Tile Size Selection for GPGPU Stencils[END_REF] or uses machine learning [START_REF] Yuki | Automatic Creation of Tile Size Selection Models[END_REF] and even cache simulators [START_REF] Mehta | Tile Size Selection Revisited[END_REF]. Polyhedral analysis is only used for the tiling transformation. Renganarayana and Rajopadhye [START_REF] Renganarayanan | Positivity, Posynomials and Tile Size Selection[END_REF] showed that most performance models for tile size selection used in the literature are polynomials that are operations research-friendly. Our automatically generated cost function, which matches the distinct-access model promoted by Ferrante et al. [START_REF] Ferrante | On Estimating and Enhancing Cache Effectiveness[END_REF], fits into this category.

Chapter 5

Conclusion

Summary of Results

This thesis presents a set of methods and a theoretical framework to automatically analyze the data movement complexity of affine programs. More precisely, for a two-level memory hierarchy with a fast memory of limited size S, our algorithms are able to find:

-a symbolic, non-asymptotic lower bound on the number of loads (transfers from slow to fast memory) that are necessary to execute this program under any valid schedule, as a function of S and program parameters; -a symbolic upper bound on the number of loads, obtained by comparing several possible tiling transformations; -a concrete tiled code with numerical tile sizes for fixed values of program parameters and fast memory size. These two algorithms have been implemented as open-source tools, and evaluated on multiple examples. On tensor contractions and convolution kernels, the expressions output by both tools match up to lower-order terms, effectively closing the theoretical problem of their I/O complexity. The lower bound algorithm was also applied on all programs in PolyBench, and is able to match or even improve over the state-of-the-art in many cases.

Limitations

Despite these promising results, IOLB and IOUB still suffer from limitations, due to the inherent complexity of these problems.

Lower bounds One category of programs for which the bounds we get are not tight are stencils (e.g. jacobi and seidel in PolyBench). This type of computation exhibits complex data dependencies (array subscripts like A[i+1][j-1]) that are not found in linear algebra kernels for instance. While some careful manual analysis can be used to slightly improve the constant factor of the expression compared to the automatic analysis, bridging the gap between lower and upper bounds is still an open problem. Moreover, the techniques presented here (projection and wavefront-based) do not cover all possible data dependency patterns, even inside the class of affine programs, and some other methods would have to be developed in order to prove tight lower bounds on gramschmidt for instance.

Upper bounds Here the limitations come more from the rapidly increasing complexity of expressions when trying to tackle more involved programs. When iteration domains get complex, for instance when loop counters depend on one another, or loop bounds are more involved functions of parameters, cardinalities computed by Barvinok's algorithm also get more complex, with several possible expressions depending on relations between parameters. This can greatly increase the size of the optimization problem, making it currently intractable. Another limit is the method to generate symbolic upper bounds: currently we rely on solving a polynomial equation, which is not generally doable when the degree exceeds four.

Finally, such automatic analyses can only be applied on a specific class of programs. It seems realistic to extend this work to get tight bounds on a large part of practical polyhedral programs, and even on some non-polyhedral ones, but many programs will certainly stay out of the scope of this type of automated analysis anyway.

Future work

The work presented here could be continued in many directions. Some of the limitations we pointed out could be overcome: for instance the affine requirement could be relaxed by using under-or over-approximations, and more complex wavefront patterns could be included in IOLB. The technique for generating upper bounds that do not depend on tile sizes in IOUB could also be generalized to cases when the polynomial equation is not solvable exactly, by relaxing the problem of finding the precise expression of a tile size that completely fills the cache to only finding a size that does not exceed the cache capacity.

Applications As we explained earlier IOUB is able to provide numerical tile sizes suggestions that minimize I/O in our model. The next step is to integrate these suggestions into a code generation framework, and get running programs that effectively use a CPU as close to its peak performance as possible. As of today, producing high-performance code is a difficult problem, for which a lot of manual tuning and machine-specific expertise is often needed. More and more solutions for automatic code generation are developed (TVM [START_REF] Chen | TVM: An Automated End-to-End Optimizing Compiler for Deep Learning[END_REF], Halide [START_REF] Ragan-Kelley | Halide: a language and compiler for optimizing parallelism, locality, and recomputation in image processing pipelines[END_REF]), with various compromises between expressiveness and code performance. Multilevel tiling is a key transformation to achieve this goal, but it must be combined with finetuned basic blocks that are able to exploit modern CPU mechanisms such as vectorization, memory prefetching, instruction-level parallelism, and adapt to machine-specific details. Choosing tiling loop sizes and permutations is a crucial step, and an analytical model such as the one we developed can be very useful to guide the search in a highly combinatorial space. It should be coupled with some autotuning, in order to take into account various aspects that we do not model, such as data layout, cache line sizes, partial tiles handling, control-flow overhead, etc.

One notable project is BLIS [START_REF] Field | BLIS: A Framework for Rapidly Instantiating BLAS Functionality[END_REF], which aims at proving a portable high-performance implementation of BLAS-like linear algebra. They use a static tiling scheme, manually derived from an analytical model, for all configurations. This scheme is the one automatically derived for matrix-matrix multiplication by IOUB in a 4-level memory setting with ordinary matrix dimensions. In the context of tensor computations and convolutions, great efforts have been done by Nicolas Tollenaere et al. [START_REF] Tollenaere | Efficient convolution optimisation by composing micro-kernels[END_REF] They achieve performance very close to native libraries like oneDNN, with higher flexibility and automation.

More general extensions Let us suggest two interesting directions. The first one concerns parallelism: the work presented here mostly deals with sequential programs, however most applications in today's computing exhibit some level of parallelism. While some data-movement lower bounds techniques can be transposed in a parallel setting to prove bounds on the volume of communication, this is not always generalizable, and the challenge of generating and analyzing highly communication-efficient parallel programs remains vastly unsolved. A first major step would be to extend the analytical model from IOUB to distributed memory communication volume, and be able to provide data distribution and communication scheme suggestions. Manual exploration shows that it should be possible to derive communicationefficient algorithms like 2.5D matrix multiplication [START_REF] Solomonik | Communication-Optimal Parallel 2.5D Matrix Multiplication and LU Factorization Algorithms[END_REF] with a similar approach.

Another direction of research would be to extend our analytical I/O cost model to include more complex aspects, such as data layout: on an actual machine data is grouped in cache lines, and conflicts can appear, so evaluating the costs and benefits of reorganizing the way data is stored (packing) could be valuable.

  Graph representation (CDAG) for N = 2.

Figure 2 . 2 -

 22 Figure 2.2 -Matrix-matrix multiplication.

Figure 2 . 3 - 2 Figure 2 . 4 -

 23224 Figure 2.3 -In-set and projections. When data dependencies are regular and orthogonal, In (P ) is bounded by the size of projections of P .

Figure 2 . 7 -

 27 Figure 2.7 -Best tiled schedule for matrix multiplication.

Figure 2 . 8 -

 28 Figure 2.8 -Ilustration of equations (2.9) and (2.10) on the CDAG for N = 2.

Figure 3 . 1 -

 31 Figure 3.1 -Example 1. (a) C-like code. (b) Corresponding single assignment form. (c) Corresponding CDAG. Input nodes A[N] (resp. C[N]) are in grey (resp. white), compute nodes are in black.

Definition 3 (

 3 I/O complexity). The I/O (or data movement) complexity of a CDAG G for a fast memory capacity S, denoted Q(G), is the minimum cost of a complete S-RW game on G.

3. 1 . 2 A

 12 compact representation of the CDAG: the data-flow graph f o r (0 ≤ t < M and 0 ≤ i < N ) i f (t==0): S[0,i]=A[i] * C[0]; e l s e : S[t,i]=S[t-1,i] * C[t];

Figure 3 . 2 -

 32 Figure 3.2 -DFG for Example 1

  are only interested in two specific types of DFG-paths, depending on their relation: -chain circuits, which are cycles from one DFG-vertex S to itself, such that the path relation R p is a translation S[⃗ x] → S[⃗ x + ⃗ b]. -broadcast S a , S b -paths, which are elementary paths (from a S a to S b -S b possibly equal to S a ) in which all DFG-edges but the first one are injective edges, such that the inverse of the corresponding relation R p is an affine function

  Figure 3.3 -Example 2

  c) Corresponding CDAG for N=5. Input nodes A[N] are omited. (d) Decomposition into two non-interfering sub-CDAGs. Sources are in gray. May-spill sets are encircled.

Figure 3 . 4 -Example 3 Ω

 343 Figure 3.4 -Example 3

Lemma 5 .

 5 Let 0 ≤ s 1 , s 2 , . . . , s m ≤ 1 and C > 0. Let x j be nonnegative integers and β j > 0 for 1 ≤ j ≤ m such that m j=1 β j x j ≤ C. Then

2 L ′ = L; 3 while not timeout do 4

 234 Lattice of subgroups L, subgroup to add K output: updated set of subspaces L ′

1 .

 1 Observe that any values of s j s lead to a correct bound function sub_paramQ_bypartition input : paths P = {P 1 , . . . , P m } with domain D and lattice L output: complexity Q I := P i ∈P R -1 P i (D); d := dim (D); (β 1 , . . . , β m ) := coeffInterf(P , D); (s 1 , . . . , s m ) := convex-opt { variables: {s 1 , . . . , s m ∈ Q + }

  |I|, 0 ; Q.may-spill := may-spill(P , D); function coeffInterf input : paths P = {P 1 , . . . , P m }, domain D output: coefficients (β 1 , . . . , β m ) such that β j φ j (E) ≤ K for any K-bounded set E G := graph with V := P 1 , . . . , P m and E := (P i , P j ), R -1 P i (D) ∩ R -1 P i (D) ∅; I := set of maximal independent sets of G such that every node belongs to at least one set (greedy construction); β j := #{I ∈ I , P j ∈ I}/ |I |; function may-spill input : paths P = {P 1 , . . . , P m }, domain D output: may-spill set D ms D ms := ∅;

Algorithm 4 :Figure 3 . 5 -

 435 Figure 3.5 -Convolution kernel.

  and the optimization problem is: Minimize s 1 + s 2 and then s . s 1 ≥ 1, s 2 ≥ 1 The solution is s 1 = s 2 = 1, and U = ((S + T ) /2) 2 = S 2 (because T = 1 s 1 +s 2 -1 S = 1 • S). |D S | = MN and |I| = N + M, so Q ≥ MN S 2 × S -N -M.

Figure 3 .Figure 3 . 6 -

 336 Figure 3.6 -Lower bounds for convolution, without and with IOLB extensions (only dominant terms are shown).

  rank φ j (H) + s sd • rank (φ sd (H)) for all subgroups H ∈ L φ 1 ,φ 2 ,...,φ m (3for all nonempty finite sets E ⊆ Z d .(3.3s)The upper bound on ρ(P ) becomes:

  s 3 +s sd (b) Constraints on s j from Brascamp-Lieb inequality.

Figure 3 . 7 -

 37 Figure 3.7 -Partitioning method -Brascamp-Lieb inequality applied to a convolution, without and with small dimensions (H and W).

Figure 3 . 8 -

 38 Figure 3.8 -Application of Cor. 1 on a sub-CDAG. Paths L j (single edge) are shown in bold.

2 D 3 R

 23 1 function sub_paramQ_bywavefront input : DFG G = (S, D), vertex S ∈ S, parametrized dimensions Ω d output: lower bound Q := dim (S); S→S := ∅; RL S→S := ∅; R S→ * := ∅; 4 foreach S j ∈ S in topological order, from S (excluded) to S (included) do

5

 5

9 R 10 R 11 W 12 Q 13 QAlgorithm 5 :

 9101112135 I d := RL S→S ∩ {S[I 1 . . . I d , i d+1 . . . i D ] → S[I 1 . . . I d + 1, i d+1 . . . i D ]}; complete := {S[I 1 . . . I d , i d+1 . . . i D ] → S[I 1 . . . I d + 1, i ′ d+1 . . . i ′ D ]}; := Dom R I d -R -1 I d R complete -(R S→S ) * (Dom R I d ) ; := max (|W | -S, 0); .may-spill := W ∪ R S→ * (W ) ∩ R -1 S→ * R I d (W ) Derivation of a lower bound with the wavefront method 3.5 Complete framework 3.5.1 DFG construction

6 Ω 8 while elapsedTime < timeout do 9 Let 10 P 16 P

 6891016 0, 1)⟩ function program_Q input : Data-flow graph G = (S, D), an instance I output: lower bound Q low Q = ∅; Let D be the max loop depth; foreach loop level 0 ≤ d < D do foreach S ∈ S surrounded by at least d + 1 loops do d := [I 1 , . . . , I d ] → {S[i 1 , . . . , i D ] : i 1 = I 1 ∧ . . . i d = I d }; 7 Let G ′ be a copy of G; D S be the parametrized domain of S in G ′ ; := ∅, L := ∅; 11 foreach P i ∈ genpaths(G ′ , S, Ω d ) do 12 if |D S ∩ Im (P i )| ≥ γ • |D S | then 13 K i := Ker(P i ); 14 if L := subspace_closure(B, K i ) changed then 15 D S := D S ∩ Im (P i ); := P ∪ P i ; 17 if P = ∅ then exit while loop;

6 (

 6 b) DFG (input nodes are omitted)

Figure 3 . 9 -

 39 Figure 3.9 -Cholesky decomposition

Figure 3 . 10 -

 310 Figure 3.10 -Input file and output from IOLB

Figure 4 . 1 -

 41 Figure 4.1 -Matrix-matrix multiplication: (a) untiled version (b) tiled version.

  f o r ( c = 0 ; c < Nc ; c ++) f o r ( f = 0 ; f < Nf ; f ++) f o r ( x = 0 ; x < Nx ; x ++) f o r (w = 0 ; w < Nw ; w++) Out [ f ] [ x ] += Image [ x+w ] [ c ] * F i l t e r [ f ] [ w ] [ c ] ;

Figure 4 . 2 -

 42 Figure 4.2 -Running example -1D Convolution

A 3 -

 3 dimensional rectangular subset of this domain will be defined as[l i , u i ] × [l j , u j ] × [l k , u k ],where, for each dimension d ∈ D, 0 ≤ l d ≤ u d < N d . The example code also contains three array accesses:-for array A, with domainM A = {x, y | 0 ≤ x < N i ∧ 0 ≤ y < N k } and access function f A (i, j, k) = (i, k); -for array B, with domain M B = {x, y | 0 ≤ x < N k ∧ 0 ≤ y < N j }and access function f B (i, j, k) = (k, j); -and for array C, with domain M C = {x, y | 0 ≤ x < N i ∧ 0 ≤ y < N j } and access function f C (i, j, k) = (i, j).

4. 2

 2 Loop permutation and tiling 4.2.1 Tiling transformation Let us consider the example of Figure 4.1 that reports a tiled and non-tiled code for matrixmatrix multiplication. The original code contains three dimensions i, j, and k. The tiled code f o r ( c1 = 0 ; c1 < Nc ; c+=Tc ) f o r ( f1 = 0 ; f1 < Nf ; f +=Tf ) f o r ( x = 0 ; x < Nx ; x ++) f o r ( c = c1 ; c < c1+Tc ; c ++) f o r (w = 0 ; w < Nw ; w++) f o r ( f = f1 ; f < f1+Tf ; f ++) Out [ f ] [ x ] += Image [ x+w ] [ c ] * F i l t e r [ f ] [ w ] [ c ] ;

Figure 4 . 3 -

 43 Figure 4.3 -Tiled code for 1D-convolution for tiling schedule (w, c, f , x) , T c = Tc, T f = Tf, T x = 1, T w = Nw

Figure 4 . 1 :

 41 ((i, j, k), T i = Ti, T j = Tj, T k = 1 where the ordered list P = (i, j, k) describes the permutation of the inter-tile loop dimensions, and T = {T i = Ti, T j = Tj, T k = 1} describes the tile sizes. Similarly, for the example in Figure4.2, the tiling shown in Figure4.3 is represented as ((w, c, f , x), T c = Tc, T f = Tf , T x = 1, T w = Nw}). Note that the loop on w in the inter-tile loops (outermost loops), and the loop on x in the intra-tile loops are both omitted in the code as they would have only one iteration. The permutation of intra-tile loops in the code is arbitrary.Let us provide a more formal definition. Definition 14. A tiling schedule is defined as a tuple (P , T ) where P = d j |D|≥j≥1 is a permutation of dimensions D representing the inter-tile loop order, and T = {T d } d∈D represent the tile dimensions. In P , d |D| represents the outermost loop dimension, while d 1 represents the innermost loop dimension that encloses the tile: the permutation order is outer (leftmost) to inner (rightmost). The tile size for dimension d is T d .

Lemma 7 .

 7 Let us consider a program with s arrays A 1 , . . . , A s , and a tiling (P = (d |D| , . . . , d 1 ), T = {T d } d∈D ) for this program. Let S 1 , . . . , S s be integers such that S 1 + • • • + S s = S, and for all i = 1, . . . , s:

Image, 1 × |I back | = 1

 11 Tf × Nw.Nc.Nf + 1 Nw.Tf × Nw.Nc.Nf.(Nx -1) = Nc.Nf.(Nx + Nw -1) TfThe same computation for arrays Out and Filter yields (assuming the outermost reuse dimensions are respectively x and f ):IO Out = Nc.Nf.Nx/Tc IO Filter = Nc.Nf.NwConditions on outermost reuse dimensions are:

  input : dimensions D, arrays A, reuse oracle output: list of permutations P function genPerm(D ′ , R) if D ′ = ∅ then return {()}; if S = ∅ for all (d, S) ∈ R then return {P } where P is an arbitrary permutation of D ′ ;

7 :Figure 4 . 4 -

 744 Figure 4.4 -Application of Algorithm 7 to the 1D convolution kernel

4 . 7 .

 47 An example of tensor contraction is given on Fig. 4.6.

  f o r ( a = 0 ; a < A ; a++) f o r ( b = 0 ; b < B ; b++) f o r ( c = 0 ; c < C ; c ++) f o r ( d = 0 ; d < D ; d++) f o r ( e = 0 ; e < E ; e++) Out [ a , b , c ] += I n 1 [ a , d , e , c ] * I n 2 [ e , b , d ] ;

Figure 4 . 6 -

 46 Figure 4.6 -Tensor contraction kernel (abc-adec-ebd)

Figure 4 . 7 -

 47 Figure 4.7 -Convolution kernel

Figure 4 .

 4 1, with tiling (i, j, k), T i = Ti, T j = Tj, T k = 1 . Kernel I/O upper bound TC abcde-efbad-cf UB = 2ABCDEF √ S+1-1

Figure 4 . 10 -

 410 Figure 4.10 -Combined parametric I/O bounds of tensor contraction (TC) and 2D convolution kernels. S is the small memory size and other uppercase letters are problem sizes (for TC kernels A. . . F are tensor dimensions, for convolution see Fig. 4.7).

Figure 4 .

 4 Figure 4.11 shows a comparison between lower and upper bounds for the considered benchmarks. For several cache sizes (from 16 kB to 4MB), the lower bound is computed by plugging the actual values in the symbolic expressions shown above, and the upper bound is obtained by running the NLP solver on the optimization problems generated by IOUB.

Figure 4 .

 4 Figure 4.11 -I/O bounds (lower is blue and upper is orange) of tensor contractions (TCCG) and convolutions (Yolo) for different cache sizes

  {H + K};

	6	else return L ′ ;
	7	return L if timeout
		Algorithm 2: Update the subgroup lattice with a new kernel

1 function genpaths input : a Data-flow graph G = (S, D), a statement S ∈ D output: set of paths P 2 start from S and backward traverse to build any possible path that reaches S; 3 drop paths for which R S ′ →S (D S ′ ) has lower dimensionality than D S ; 4 drop paths if ¬(isBroadcast(P ) ∨ isChain(P )); Algorithm 3: Generate paths

Table 3 .

 3 1 reports, for each kernel in PolyBench: -the size of the input data and the number of operations; -the simplified I/O lower bound Q ∞ low from IOLB; -the parametric lower and upper bounds on operational intensity OI manual and OI up =

	# ops Q ∞ low -the best known published OI up , when it exists. , and their ratio;
	The 30 reported benchmarks can be divided into four categories, corresponding to table
	divisions:
	1. (19 kernels) The ratio	# ops

Programmes affines et automatisation L'automatisation des techniques présentées ci-dessus est une contribution clé de cette thèse. Pour que cela soit possible, il a fallu choisir une classe de programme sur laquelle travailler, à savoir la classe des programmes affines. Cette classe constitue un bon compromis entre expressivité et applicabilité des techniques de calcul d'I/O, et dispose d'un riche écosystème d'outils d'analyse. Les programmes affines comprennent un large éventail d'algorithmes, notamment l'algèbre linéaire dense, les stencils, les convolutions, et de nombreux programmes utilisant la programmation dynamique. Un programme affine se compose de boucles for (possiblement imbriquées) et d'expressions, qui opèrent sur des variables et des tableaux multidimensionnels. Les bornes des boucles sont des expressions affines des paramètres du programme et des indices de boucles qui les entourent. De même les accès aux tableaux sont des fonctions affines des indices de boucles et des paramètres. Ainsi, le flot de contrôle d'un programme affine est analysable statiquement.Les programmes affines peuvent être analysés par les outils de compilation polyédrique, comme ISL et barvinok. Ces outils représentent les programmes par des points dans un espace multi-dimensionel, et peuvent effectuer des opérations sur des ensembles de points, analyser les dépendances de données entre instructions, compter le nombre de points dans un ensemble défini paramétriquement, entre autres. Cela est particulièrement utile pour notre calcul de bornes inférieures, pour mettre en évidence des dépendances régulières qui correspondent à des projections.Expériences et résultats Ces algorithmes ont été implémentés dans des outils disponibles librement, et évalués sur un certain nombre de programmes, comprenant les 30 programmes de PolyBench, une suite d'évaluation de référence de programmes affines. Notamment, sur les contractions de tenseur et les convolutions, ainsi que sur plusieurs algorithmes d'algèbre linéaire, les bornes inférieures et supérieures sont égales asymptotiquement, ce qui permet de clore la question de l'évaluation de leur complexité I/O.

From here on, we omit lower-order additive terms to keep things readable. The full formula output by IOLB can be found in Sec

3.6.4.

S for all of them. IOLB returns √ S for 4 of these kernels, and 2 √ S for the other 3. The

√ SN3

2ABC(DE-1)√ S 2 √ S

2AB(CD-1)√ S 2 √ S

2AB(C-1) √ S 2 √ S ABC

In theory, Barvinok's algorithm has an exponential complexity in the number of dimensions of the set. In practice, its execution time is usually less than a few seconds, when applied to the polyhedral sets commonly encountered during program analysis.

Remerciements

 Definition 16. The sub-domain data footprint for an array A at level j is defined as: SDF A,j (iD, . . . , ij) = f A SD d j (iD, . . . , ij)

Taking the example of the 1D-convolution of Figure 4.3, the sub-domain SD d 2 (c1,f1) for loop dimension

The sub-domain data footprint at level 2 for array Image for this program is:

Whenever there is a non-empty overlap between the data used by two consecutive subdomains, estimating it allows us to refine our cost model. This leads us to define the intersub-domain reuse as follows:

Definition 17. The inter-sub-domain reuse for an array A at level j is:

A sharp eye would have observed that the notion of reuse from a "previous" sub-domain is meaningful for all sub-domains but the first one in the loop. To handle this subtlety, for a given loop dimension d j (assuming its loop index starts at 0), we split the iteration space into two sub-domains: the front domain

and the back domain

This allows us to define, for a given sub-domain and an array A, its inverse density as the ratio:

ID A,j (iD,...,0) = SDF A,j (iD, . . . , ij)

As we will see later, the inverse density is an over-approximation of the best attainable inverse operational intensity (data movement per computation unit). As an example, for 1Dconvolution:

To avoid complicated expressions for non-rectangular domains, we consider the maximum value for the inverse density that we specialize for the front and the back:

ID A,j (iD, . . . , ij)

In many cases, the sub-domain sizes and intersections do not depend on loop indices, and there is no need to take the maximum.

Cost model

Now, we want to evaluate the I/O cost of a given tiling schedule. A tiling schedule defines the order of "computation operations" (rule (R2) in red-white pebble game terminology), but the order of memory operations still needs to be defined.

Single array

Let us start with the case when there is a single array A. We consider a tiling schedule

The basic idea for scheduling I/O operations is simple: for each tile, first bring all input data in fast memory, then execute all computations inside the tile.

For this to be possible, the data footprint of a single tile has to be smaller than the cache capacity S:

For simplicity, we restrict our I/O cost model to situations where this condition holds. It would be possible to compute it when is it not the case, but the I/O cost would be high. Since our goal is to derive an upper bound as close as possible to the optimal, this is not necessary.

The basic idea described above can be improved: sometimes, consecutive tiles actually use the same data. To be more accurate, we want to find the largest loop nest such that the input data for all tiles in the loop nest still fits in fast memory. The data then only needs to be brought in fast memory once for each such "macro-tile". Finally, if there is some overlap in the input data for two consecutive tiles, only the additional data needs to be fetched.

Formally, we define the outermost reuse dimension for array A, as the "first" (smallest possible l) dimension d l in P such that the sub-domain data footprint of A is less than the cache capacity S. Definition 18. The outermost reuse dimension l for array A is the leftmost index in P such that SDF A,l ≤ S.

The total I/O cost for array A derived from the inverse density as follows constitutes an upper bound on the cost of an optimal red-white pebble game: