Chong Li

Zhen Zhang

Sheng Yang

Yu Xia

Jie Zhao

Keywords: Parallel programming, polyhedral compilation, automatic parallelization

Pour être honnête, j'ai toujours douté des performances impressionnantes et de la magie surprenante du modèle polyédrique avant de rejoindre le groupe PARKAS. Les obscures abstractions mathématiques du modèle m'ont toujours empêché d'atteindre le fond du monde polyédrique. J'ai donc pris la décision d'entrer dans ce monde en poursuivant un doctorat dans le domaine de la recherche. Aujourd'hui, J'arrive finalement au bout de ce chemin, avec les conseils, les soutiens et les assistances d'un grand nombre de personnes.

Je tiens à exprimer ma plus profonde gratitude à mon superviseur, Albert Cohen, pour m'avoir conduit au monde de la compilation polyédrique et pour sa patience, sa volonté et sa passion pour m'aider à approfondir des idées de recherche pratiques et intéressantes. Albert m'a toujours fourni des retours très utiles via des discussions en personne et des échanges de vues en ligne, me laissant ainsi beaucoup de liberté et d'espace pour explorer les problèmes et rechercher des solutions de manière autonome. Il m'a également donné beaucoup de conseils et une foule de suggestions grâce à sa vision prospective, m'inspirant de suivre les orientations de recherche prometteuses dans le monde académique et les questions d'actualité. Ce fut une excellente expérience conseillée par Albert et je suis impatient de retravailler avec lui.

Un grand merci à Michael Kruse, à Oleksandr Zinenko, à Chandan Reddy et à Riyadh Baghdadi pour les conversations intéressantes et utiles que j'ai eues au cours de mon doctorat. Parler avec eux m'aide à mieux comprendre le contexte mathématique de la compilation polyédrique et leurs suggestions constructives sur mon travail m'aident toujours à progresser. Merci à Marc Pouzet, Francesco Zappa Nardelli et Timothy Bourke d'avoir dirigé l'équipe de PARKAS et créé un environnement harmonieux pour le bon déroulement de ma thèse. Je tiens également à remercier Ulysse Beaugnon, Lélio Brun et Guillaume Iooss, qui m'ont beaucoup aidé lors de mon inscription à l'école doctorale et à l'université ainsi que de nombreux autres travaux quotidiens

Many thanks to Michael Kruse, Oleksandr Zinenko, Chandan Reddy and Riyadh Baghdadi for the interesting and helpful conversations during my doctoral study. Talking with them helps me better understand the mathematical background of polyhedral compilation and their constructive suggestions on my work always help me make progress. Thanks to Marc Pouzet, Francesco Zappa Nardelli and Timothy Bourke for leading the PARKAS team and establishing a harmonious environment for the smooth progress of my thesis. I also acknowledge Ulysse Beaugnon, Lélio Brun, and Guillaume Iooss, who helped me a lot on my registration in the doctoral school and the university as well as many other daily activities.

I would also like to acknowledge Bill McColl and his team, including Chong Li, Zhen Zhang, Sheng Yang, and Yu Xia, for hosting me during my visit at Huawei Research at Paris. We have evolved good friendships and they will last for my whole life. The collaboration work at Huawei for building a connected, intelligent world is an unforgettable experience. Thanks to your support and grant for my visit to the Huawei Research Center at Beijing, China. I am expecting to cooperate with you in the future. I am very grateful to the valuable time and efforts from the reviewers of my PhD dissertation, iii I also feel very thankful to the contributors of the PPCG code generator and the isl integer solving library, Sven Verdoolaege and Tobias Grosser. This work would not be possible without the contributions of their work. I would also like to acknowledge the contributors and developers of a large number of underlying components of polyhedral compilers and optimizers, including pet, Pluto, CLooG, Omega, iscc, etc., without which the understanding of the polyhedral model would be an unimaginable task.

I have to express my deep gratefulness to some important persons who walked with me the path of undergraduate and graduate study, including Baoliang Li, Jinchen Xu, Zhifeng Chen, Guanghui Liang and Xu Zhao, without some of whom my programming knowledge and skills would be impossible.

Last but not least, I want to acknowledge my family for their endless support to my education from the childhood. Everyone needs a house to live in, but a supportive family is what builds a home.

Paris, 10 December 2018

Jie Zhao iv

Introduction

In the early stage of high-performance computing, increasing clock frequencies was the main source of performance gain. Since the breakdown of Dennard scaling a dozen years ago, most CPU manufacturers have been focusing on multicore processors as an alternative of raising clock frequencies from one generation to next. Figure 1.1 shows the semi-centennial trend of clock frequency and number of cores per chip.

Multicore processors are nowadays ubiquitous on almost all platforms, ranging from supercomputers ranked on the TOP500 list to personal laptops and mobile devices. In addition, their pervasiveness in the embedding computing domain of multimedia and image processing due to the recent process of neural networks also validates the dominance of multicore processors in all realms of computing. Apart from the resulting improvements in performance due to the introduction of multicore processors, one of the most challenging issues is the difficulty to effectively exploit parallelism on such devices. On the one hand, a multicore processor is allowed to implement multiprocessing freely by coupling the cores on the device either tightly or loosely, leading to a variety of memory hierarchies and resulting in the diversity in architectures. On the other hand, the evolution of parallel architectures also calls for the innovation of programming languages amenable to the memory hierarchy, giving rise to the design of both general-purpose and domain-specific parallel programming languages and complicating the programmability issue further.

Even though end users may be equipped with the knowledge of the high-level programming language of a platform after a long-term study or training, it is still a complex and error-prone task to deploy the code written by end users on the target architecture. An optimizing compiler is not only responsible for translating the code implemented by a high-level general-purpose/domain-specific language into a low-level executable program, but is also expected to automatically apply both high-level and low-level transformations, especially those performance-critical loop transformations, for exploiting parallelism and improving locality, thus releasing the burden of end users from taking the hardware information into consideration at the beginning of programming.

In the domain of scientific and engineering applications, a large number of computationally intensive applications spend most of the execution time on nested loops, making the polyhedral model [START_REF] Feautrier | Polyhedron Model[END_REF] a very competitive and promising approach to solving the above problems. The polyhedral model is a powerful mathematical abstraction of loop nests, providing a way to reason about loop transformations by abstracting each iteration of a statement as an integer point in a "polyhedron" and mapping a multi-dimensional logical execution date [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem. part ii. multidimensional time[END_REF] for defining its lexicographic execution order. As a role of bridging the gap between high-level programming interfaces and underlying hardware, the polyhedral model has made a great deal of progress in the past few decades, but it is now facing new emergent challenges brought by both modern architectures and programming languages.

From General-purpose Languages to Domain-specific Languages

Thanks to the significant advances in dependence analysis [START_REF] Feautrier | Dataflow analysis of array and scalar references[END_REF][START_REF] Pugh | The omega test: A fast and practical integer programming algorithm for dependence analysis[END_REF][START_REF] Vasilache | Violated dependence analysis[END_REF][START_REF] Baghdadi | Improved loop tiling based on the removal of spurious false dependences[END_REF], schedule transformation [Fea92a, Fea92b, LL97, BHRS08, BAC16, UC13, ABC18] and code generation [AI91, Che, QRW00, Bas04, VBC06, GVC15], the polyhedral model has been brought to the front scene in automatic parallelization and locality optimization. There exist a large number of mature polyhedral compilation frameworks and loop optimizers, including both research projects [BHRS08, CCH08, VCJC + 13] and commercial productions [TCE + 10, GGL12, CSG + 05, BGDR10, LLS06]. Such compilers usually take a general-purpose (intermediate) language as input and generate optimized high-level/low-level code amenable to the target architecture as demand. Despite that, the optimality of the code generated by such polyhedral 1.2. Architecture Diversity compilers still remains elusive, falling behind the performance of heavily hand-tuned codes written by an expert.

Part of the reason of performance gap between the generated codes of optimizing compilers and hand-written programs is due to the conservativeness that a compiler has to possess in nature as system software, reducing the optimization space of automatic transformations. Worse yet, the absence of the ability to reason about the domain knowledge about the implementation strategies from a piece of code also constraints such compilers, missing aggressive and/or global optimizations that can be performed by hand. Domain-specific languages (DSLs) are proposed to obtain high performance and now very prevalent in many application domains. The polyhedral model was successfully integrated with DSLs, such as those for optimizing DSLs for graphical dataflow language [START_REF] Somashekaracharya | Polyglot: A polyhedral loop transformation framework for a graphical dataflow language[END_REF][START_REF] Sbîrlea | Polyhedral optimizations for a data-flow graph language[END_REF], stencil computations [HVF + 13], etc. Recently, due to the revolution in machine learning caused by the success of deep learning, the polyhedral community is also expected to resolve the problem of bridging neural network applications and high-performance hardware accelerators. A DSL may be a standalone language or more often embedded in a general-purpose language, like Halide [RKBA + 13] in C++, TensorComprehensions [VZT + 18] and TVM [CMJ + 18] in Python1 , DeepDSL [START_REF] Zhao | Deepdsl: A compilation-based domainspecific language for deep learning[END_REF] in Scala, etc. A domain-specific compiler leverages specialized internal representations for expressing domain-specific knowledge, extending its optimization space by enabling such domain-specific high-level transformations. Representative DSL compilers for such applications include the TensorComprehension framework for automating the deployment of neural network applications on multicore platforms and the PolyMage compiler [MVB15, MAS + 16] for Halide [RKBA + 13], a DSL for writing high-performance image processing code.

While the polyhedral model eases the translation of both general-purpose languages and DSLs on modern architectures, it often suffers from scalability challenges to various input languages. Even though some internal representations like Hailde IR and PENCIL [BBC + 15] were proposed as the solution to this problem, the polyhedral model still faces many painful problems due to its incompetence for dynamic control and non-affine applications.

Architecture Diversity

Generally speaking, a multicore system is supposed as homogeneous if the system includes only identical cores, or heterogeneous otherwise. The Pluto optimizer [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF] provides a systematic, end-to-end way for automatic parallelization and locality optimization on homogeneous multicore systems, taking into consideration the memory hierarchy problem by automating simply/complex tile shapes [START_REF] Bondhugula | Diamond tiling: Tiling techniques to maximize parallelism for stencil computations[END_REF]. The emergence of Graphics Processing Units (GPUs) brought new challenges not found in homogeneous systems to the polyhedral model, calling for source-to-source polyhedral compilers capable of generating correct codes for both host processors and device accelerators, further complicating the code generation issue.

Unlike CPUs that can run efficiently when data is resident in Caches, GPUs have a variety of different kinds of processing units, leading to a more complicated memory hierarchy. For instance, Table 1.1 lists the memory hierarchy of fastest NVIDIA Tesla GPUs in the past five years. PPCG [VCJC + 13] is considered as one of the most successful polyhedral compilers for heterogeneous systems, exploiting parallelism and locality optimization as in traditional homogeneous systems but also automating the management of memory system on devices and communications between host and device. Like the diamond tiling technique in Pluto, a hybrid/hexagonal tiling approach [GCH + 14] is also implemented in PPCG for further improving the performance of generated code. Some follow-up PPCG-based researches focus on parametric tiling [START_REF] Juega | Adaptive mapping and parameter selection scheme to improve automatic code generation for gpus[END_REF] and the mapping and separation of multi-level parallelisms in the accelerators of a heterogeneous system [START_REF] Shirako | Optimized two-level parallelization for gpu accelerators using the polyhedral model[END_REF].

Besides shared memory strategy, message passing is also used as the inter-core communication method in distributed systems and heterogeneous systems, requiring the code generator of an optimizing compiler to express the explicit communication with libraries like Message Passing Interface (MPI). Polyhedral compilation frameworks targeting on minimizing communication volume [START_REF] Bondhugula | Compiling affine loop nests for distributed-memory parallel architectures[END_REF][START_REF] Reddy | Effective automatic computation placement and data allocation for parallelization of regular programs[END_REF] or handling the mixture of regular/irregular loop nests [RDE + 15] were proposed for such multicore systems.

Similarly, accelerators in heterogeneous systems are not restricted to GPUs. For example, configurable devices like Field-Programmable Gate Array (FPGA) [START_REF] Uday Bondhugula | Automatic mapping of nested loops to fpgas[END_REF][START_REF] Pouchet | Polyhedralbased data reuse optimization for configurable computing[END_REF] can also be the target of an optimizing compiler, followed by some researches in high-level synthesis area [ZLC + 13, WLC14]. These together with the above mentioned architectures are calling for a strict portability of the polyhedral model to multiple platforms. Recent work integrating multicore parallelism and Single Instruction Multiple Data (SIMD) vectorization [TNC + 09, KVS + 13] not only addressed the code generation issue but also implemented a different

Beyond Parallelization and Locality Optimization

scheduling strategy in a tile.

Even though the advances made by the polyhedral community on so many accelerators, there still exits a long way to achieving architectural portability. The emergence of artificial intelligence applications has brought new challenges for this issue. For instance, Table 1.2 summarizes the latest xPUs used for modern artificial intelligence accelerators.

Beyond Parallelization and Locality Optimization

On modern multicore processors, parallelism due to the increased core numbers on a single chip and locality caused by memory hierarchy are the two main objectives considered by compiler designers. As a result, optimizing compilers like a polyhedral optimizer are usually expected to be capable of automatic parallelization and locality optimization. Unfortunately, parallelization and locality optimization are sometimes contradictory with each other by putting conflict constraints on the objective function of scheduling algorithms, forcing them to make a tradeoff between parallelism and locality for achieving optimality of performance.

As a loop transformation aiming at improving locality while preserving the parallelism that has been exploited by a scheduling algorithm, loop tiling [START_REF] Irigoin | Supernode partitioning[END_REF] has been long considered as foreign to optimizing compilers; even for the polyhedral model, it could not be easily expressed using an affine function a decade ago. Thanks to its recent advances, the polyhedral model has been proved to be promising in automating loop tiling. A cost-model-based scheduling algorithm like the Pluto scheduler [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF] or its variants put into practice the automation of simple tile shapes in the polyhedral model. A follow-up trend on the automatic tiling technique focuses on more complex tile shapes like diamond [START_REF] Bondhugula | Diamond tiling: Tiling techniques to maximize parallelism for stencil computations[END_REF] and hexagonal [GCH + 14] working with arbitrary affine dependences, and overlapped and split shapes [KBB + 07] restricted to constant dependence vectors.

Loop fusion [START_REF] Kennedy | Optimizing Compilers for Modern Architectures: A Dependence-based Approach[END_REF] is another loop transformation to enhance locality and reduce synchronizations across multiple loop nests. There have also been successful advances on loop fusion in the polyhedral model [BGDR10, MLY14, JB18], providing a variety of fusion heuristics to modern optimizers.

Storage optimization is also a research direction of polyhedral compilation. Array contraction, for example, is a long considered automatic memory footprint optimization in the polyhedral world. The applicable domain is still constrained to special cases like stencil computations although researchers made a lot of efforts in this direction, including the universal-occupancyvector-based [START_REF] Strout | Scheduleindependent storage mapping for loops[END_REF], lattice-based [START_REF] Darte | Lattice-based memory allocation[END_REF] and storage-hyperplane-based [START_REF] Somashekaracharya | Smo: An integrated approach to intra-array and inter-array storage optimization[END_REF] techniques, etc.

In spite of the exciting progresses made by the polyhedral community on automatic parallelization and locality optimization, there still exist a large number of opening issues awaiting supports and efforts. The latest research trend also tried to integrate the polyhedral model with dynamic/runtime techniques [KPP + 15, SRC15, BKP + 16, SPR17] for extending the scope of the tool, leaving much room for the extensions in this field.

Combining Languages and the Polyhedral Model

The polyhedral model so far is successful in so-called "static control parts" (SCoPs) where loop nests satisfy certain statically predictable restrictions. There is an increasingly emergent demand on its applicability to non-affine domains to cope with the complexity of modern multicore architectures. A notable direction among the open challenges is the incompetence of the polyhedral model to handle non-affine applications and transformations. Such non-affine applications usually involve dynamic data-dependent control flow and/or non-affine expressions that go beyond the scope of the polyhedral model, while non-affine transformations2 are usually not expressible using existing techniques.

A representative polyhedral-based approach on non-affine applications are the work of handling while loops [START_REF] Benabderrahmane | The polyhedral model is more widely applicable than you think[END_REF], along with a great deal of work with special focus on sparse matrix computations [SGO13, VSHS14, VHS15, SLC + 16, VMP + 16]. The former misses more aggressive optimizations when handling less expressive dynamic conditions than a general while loop, while an inspector/executor scheme is usually constrained to a subset of sparse matrix computations.

Combining Languages and the Polyhedral Model

With regard to non-affine transformations, overlapped tiling [KBB + 07] is a representative technique gaining much attentions recently due to its compatibility with other optimizations like fusion, scratchpad memory allocation, etc. when optimizing image processing pipelines. Unlike a standard tile shape exploited by current polyhedral compilers, additional overlapped regions are introduced for exploiting inter-tile parallelism in such transformation by jointing consecutive tiles. Unfortunately, no implementations of overlapped tiling in a general-purpose polyhedral framework have been reported except PolyMage [START_REF] Teja Mullapudi | Polymage: Automatic optimization for image processing pipelines[END_REF] , a DSL compiler for image processing pipelines.

This dissertation describes a combined language and polyhedral approach to extend the application domain of the polyhedral model in non-affine applications and express non-affine transformations in the model. On the one hand, we study the parallelizing compilation and loop nest optimization of an important class of programs where counted loops have a dynamic data-dependent upper bound. Such loops are amenable to a wider set of transformations than general while loops with inductively defined termination conditions: for example, the substitution of closed forms for induction variables remains applicable, removing the loopcarried data dependences induced by termination conditions; such loops can also be viewed a generalization of sparse matrix computations using compressed data layout stores nonzero elements only as the latter can be easily generalized by subtracting the lower bound from the upper bound.

On the other hand, we revisit overlapped tiling in polyhedral compilation and demonstrate how to derive tighter tile shapes with less redundant computations, by enabling overlapped tiles based on a well-defined general-purpose intermediate representation. It releases the overlapped tiling in polyhedral model from being restricted to a domain-specific language while not introducing sophisticated rescheduling and custom code generation in a polyhedral framework.

Given the diversity of multicore architectures and the difficulty of programming on these platforms, a polyhedral compilation approach has become a compelling alternative for writing parallel code on these targets. Our approach is driven by combining an intermediate language and the polyhedral model, not only removing the conservativeness caused by using a general-purpose language hindered by the difficulty of static analysis but also avoiding the implementation of a DSL compiler for the portability to different architectures. By coupling with such an intermediate language, one may define coding rules predominantly related to restricting the non-statically predictable manners, allowing for better optimizations when translating such programs into the code on target machines using a polyhedral framework. More importantly, leveraging such an intermediate language also eases the code generation for different architectures, making the portability issue a straightforward task.

Our method on counted loops with a dynamic data-dependent upper bound goes beyond the state of the art in fully automating the process, specializing the code generation algorithm to the case of dynamic counted loops and avoiding the introduction of spurious loop-carried dependences. The experimental results on representative irregular computations ranging from dynamic programming, computer vision and finite element methods to sparse matrix linear algebra validate that the method is applicable to general affine transformations for locality optimization, vectorization and parallelization.

Our algorithm to generalize overlapped tiling allows for tighter overlapped tile shapes than the state of the art, further improving the performance of image pipelines on both general-purpose multicores and heterogeneous accelerators by integrating with transformations including alignment and scaling of stages in the pipeline, loop fusion, scratchpad allocation, hybrid tiling, etc. The experimental evaluation on the PolyMage benchmarks and representative iterated stencils validates the effectiveness and general applicability of our technique on both general-purpose multicores and accelerators.

The organization of the dissertation is as follows. In the first part, we described in this chapter an introduction to the problem we aim at in this dissertation, followed by Chapter 2 providing the technical background on the polyhedral model and the intermediate language used in our approach. The second part presents our method to handle non-affine application, i.e., parallelization and optimization of counted loops with a dynamic data-dependent upper bound, including the motivation and related work in Chapter 3, scheduling algorithm in Chapter 4, code generation method in Chapter 5 and experimental results in Chapter 6. Generalizing overlapped tiling in the polyhedral model regarding the non-affine transformations is introduced in the third part, comprising Chapter 7 describing the motivation and related work, Chapter 8 explaining the polyhedral implementation of the method and Chapter 9 evaluating the proposed technique on both homogeneous and heterogeneous architectures. We finally conclude the dissertation in the last part, Chapter 10, by summarizing the topics studied in the dissertation and discussing directions for further research.

Background

After a thirty-year evolution, the polyhedral model has become a powerful optimizer in the domain of automatic parallelization and optimization. There have been a great number of open-source and/or commercial implementations of polyhedral compilation in both research and industry worlds. Comparing with unimodular matrices [START_REF] Banerjee | Unimodular Matrices[END_REF][START_REF] Wolf | A loop transformation theory and an algorithm to maximize parallelism[END_REF] used in parallelizing compilers, the polyhedral model is equipped with (1) wider range of applications due to the capability to transform imperfect loop nests, (2) more powerful expressiveness by modeling almost each kind of loop transformations and (3) greater optimization space by compositing more transformations at one time.

As a consequence, polyhedral compilation nowadays is gradually becoming the state of the art of almost each domain of parallelizing compilers. In this chapter, we would first introduce the background of polyhedral compilation for a better understanding of the underlying principle of the polyhedral model. To cope with polyhedral compilation for non-affine applications and transformations, we would next present the intermediate language used in the dissertation.

Polyhedral Compilation

As we introduced in the previous section, the polyhedral community so far has made a great deal of progress in all realms of computing. Nonetheless, polyhedral compilation is long considered as too abstract for those people outside the polyhedral world. Part of the reason is due to the painfully theoretical descriptions in existing polyhedral publications; more importantly, the underlying principle of the polyhedral model involves a variety of concepts from linear algebra, static analysis, etc., making the use of the tools elusive for end users.

To make it easier to understand the polyhedral model, we introduce the background of polyhedral compilation in this section. One may refer to [START_REF] Feautrier | Polyhedron Model[END_REF] for a much detailed description on fundamental concepts and definitions. In general, we would first give an overview of modern polyhedral compilation and then explain how programs are represented in the model. Next, the transformations that can be modeled in the polyhedral model are presented by comparing

Polyhedral Compilation

Dependence analysis. First, polyhedral extraction is serving as the frontend, parsing the code fragment of input languages and extracting polyhedral representations with regard to statement instances and array elements for the program. The duty of a parser is checking whether the input is a "SCoP", meaning the code fragment is statically predictable and polyhedral compilation would represent the input program with finite internal representations if true. pet [START_REF] Verdoolaege | Polyhedral extraction tool[END_REF] and clan [cla] are two representative, practical parsers for polyhedral compilation, generally used in a variety of mature polyhedral compilers. Such representations of statement instances and their relationships with the array elements they access are used to compute dependence relations by solving an integer linear programming (ILP) problem.

Polyhedral compilation differentiates dependence analysis from conventional methods by refining the analysis from statement-wise to instance-wise [START_REF] Feautrier | Dataflow analysis of array and scalar references[END_REF]. Dependence relations can be further separated into value-based dependences, the result of data flow analysis [START_REF] Dror | Array-data flow analysis and its use in array privatization[END_REF][START_REF] Maslov | Lazy array data-flow dependence analysis[END_REF], used for preserving the semantic of programs, and memory-based dependences, studied for the purpose of improving data reuse [PW92, VBCG06, BCVT13].

Schedule transformation.

Secondly, schedule transformation is the core of polyhedral compilation, producing a new schedule by taking into consideration target architectures. In other words, schedule transformation is the process of mapping a new logical execution date for each integer point in a polyhedron, accomplished by invoking the underlying ILP solver. The process of schedule transformation could also be considered as a composition of different loop transformations with the purpose to fully exploit parallelism and data locality.

We take the 1D iterated stencil shown in Figure 2.2(a) as an illustrative example. Iterated stencils are a class of computations updating an array element using its neighbors, commonly found in computational fluid dynamics, image processing, partial differential equations, etc. The original iteration space of the 1D stencil code is shown in Figure 2.2(b), indicating the computation proceed first along t axis and then i axis. Instead, the transformed iteration space after schedule transformation in Figure 2.2(c) implies the computation should first follow t axis and then t + i axis. The instances of the statement are represented by integer points in iteration space, coordinated with each other by a blue arrow denoting a dependence relation.

Tiling along t axis and i axis in Figure 2.2(b) is illegal since such tiling may produce dependence cycles between tiles along i axis, prohibiting the data locality along t axis of the original iteration space. On the contrary, one may benefit from the data locality along both axes on the transformed iteration space as tiling along the axes may not result in dependence conflicts.

In fact, schedule transformation could be understood as the reconstruction of the basis of iteration space, attained by a scheduling algorithm like [Fea92a, Fea92b, LL97, BHRS08, BAC16, UC13, ABC18] and their variants in libraries. While the scheduling algorithm proposed by Feautrier [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem. i. onedimensional time[END_REF][START_REF] Feautrier | Some efficient solutions to the affine scheduling problem. part ii. multidimensional time[END_REF] Figure 2.2 -1D iterated stencil and its iteration space before and after scheduling overhead, the hyperplane1 partitioning technique [LL97] also failed to minimize the order of synchronization even though it takes into account communications. The cost-model-based scheduler [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF] developed in the Pluto compiler was designed to overcome such flaws and has been demonstrated as effective in practice by a variety of implementations.

Schedule transformation is considered as the most difficult component of polyhedral compilation. A scheduling algorithm is tightly coupled with dependence relations produced by the frontend: while it should preserve the semantic of the program by being constrained to dependence relations, it is also expected to minimize dependence distances for reuse purposes, thereby improving data locality. A scheduling algorithm is not only responsible for exploiting different compositions of loop transformations, e.g., the composition of loop tiling and skewing could be triggered in the example of Figure 2.2, but also obligated to exploit both fine-grained and coarse-grained parallelisms. Moreover, a scheduling algorithm should also take into consideration the requirements from follow-up code generation: for example, schedule transformation should allow the insert of thread-/warp2 -level synchronizations when generating CUDA code on GPUs.

Code generation. Finally, code generation constitutes the backend of polyhedral compilation, made up by two phases with one building abstract syntax trees (ASTs) from the results of polyhedral extraction and schedule transformation while the other generating expected high-performance code executable on target platforms. As code generation would not change the execution order of programs, it is allowed to not take into consideration dependence relations, scanning the iteration space and generating code according to execution dates defined by the transformed schedule. As a consequence, code generation is also referred to as polyhedral scanning.

When building ASTs from the results of previous steps, a code generator manages to determine loop bounds and conditionals of control flow by seeking solutions for an optimization problem subject to integers. The generated ASTs could then be passed to emit instructions amenable to different programming models on target machines, facilitating the portability to different architectures.

One representative implementation of code generation in polyhedral compilation is the convex-based algorithm [START_REF] Ancourt | Scanning polyhedra with do loops[END_REF], generating code by first constructing a convex for all polyhedra in the iteration space. Figure 2.3(a) shows an illustrative iteration space composed of two polyhedra with one comprising all red square points and the other made of blue circle points, followed by a diagram of convex-based algorithm in Figure 2.3(b). One of the flaws of this algorithm is the generated code may include multiple nested if conditionals governing the correct execution of each polyhedron in the iteration space, promoting some code generators like Codegen+ working on hoisting if conditionals [Che]. The other code generation technique was proposed by Quilleré et al. [START_REF] Quilleré | Generation of efficient nested loops from polyhedra[END_REF] and implemented in the CLooG generator and its variants [START_REF] Bastoul | Code generation in the polyhedral model is easier than you think[END_REF][START_REF] Vasilache | Polyhedral code generation in the real world[END_REF][START_REF] Grosser | Polyhedral ast generation is more than scanning polyhedra[END_REF]. Unlike the convex-based algorithm, the method used in CLooG may first split the polyhedra into distinct regions as shown in Figure 2.3(c), producing code by scanning each of such regions individually.

As one may also find in Figure 2.1, an ILP solver is at the core of polyhedral compilation, providing each step with minimal flexible integer solutions, thereby achieving the manipulation of polyhedral transformations. There exist a variety of libraries for solving ILP, including isl [START_REF] Verdoolaege | Isl: An integer set library for the polyhedral model[END_REF], Omega [KMP + 96], PIP [START_REF] Feautrier | Parametric integer programming[END_REF], PolyLib [START_REF] Loechner | Polylib: A library for manipulating parameterized polyhedra[END_REF] and PPL [START_REF] Bagnara | The parma polyhedra library: Toward a complete set of numerical abstractions for the analysis and verification of hardware and software systems[END_REF], etc., differing each other by using different algorithms and data structures. some widely used optimizers like LLVM/Polly, PPCG, etc. One is without doubt free to choose any other combinations of representations for their implementation, as such representations could always be mutually transformed.

Before the introduction to polyhedral representations, we prefer to first present some mathematical concepts. This is because such mathematical concepts are the underlying expressions in integer manipulation libraries. More importantly, it would be a painful task to explain polyhedral representations if bypassing such mathematical concepts. If the readers are interested in a more detailed description or some more concentrated examples on such mathematical descriptions, we would suggest to refer to the work of Grosser [START_REF] Grosser | A decoupled approach to high-level loop optimization: tile shapes, polyhedral building blocks and low-level compilers[END_REF].

Integer sets. An integer set is a set of n-tuple integers subject to a group of affine constraints relating such n-tuple integers with m-tuple constant parameters, with n representing the dimensionality of the set, m the dimensionality of the parameters in constraints. Mathematically, an integer set can be written as

S = {(i 1 , i 2 , ..., i n ): f ((i 1 , i 2 , ..., i n ), (p 1 , p 2 , ..., p m ))} (2.1)
An integer set is called named integer set when assigning a name to the integer tuple. In practice, one may think an n-tuple integer set as the collection of loop iterators, m-tuple constant parameters the parameters of programs. A constraint function usually comes as the conjunction of multiple inequalities.

Integer maps. An integer map is a binary relation mapping an n 1 -dimensional integer set, i.e., the domain of the map, to another n 2 -dimensional integer set, i.e., the range of the map, subject to a set of affine constraints on the integer sets and constant parameters. An integer map can be generalized as

M = {(i 1 , i 2 , ..., i n 1 ) ! (i 1 , i 2 , ..., i n 2 ): f ((i 1 , i 2 , ..., i n 1 ), (i 1 , i 2 , ..., i n 2 ), (p 1 , p 2 , ..., p m ))} (2.2)
An integer map can be interpreted as the relation between statement instances with their accessed data locations or dependence relations. In the same way, a named integer map represents an integer map between two named integer tuples.

Named union sets and Named union maps. We use named union sets to refer to the union of different named integer sets, and named union maps for the union of different named integer maps. Named union sets can be used to express all statement instances of a SCoP, while named union maps may either define execution order on statement instances or relate statements and the data they access.

Given the above mathematical concepts, we can explain polyhedral representations in a much easier way. Following the previous subsection, we still use an illustrative example for explanation. One may obtain the iteration space of the loop nest listed in Figure 2.4(a) as shown in Figure 2.4(b), with one polyhedron for statement S1, composed of blue points, and the other for statement S2, depicted with the collection of all red points. f and g are affine functions of their indices. In our work, we would use the following representations.

for (i=1; i<=4; i++) { S1:

a[i]=f(i); for (j=1; j<=i; j++) S2: b[i][j]=g(a[i]); } (a) An illustrative loop nest 1 2 3 4 1 2 3 4 j i (b)
The iteration space of the loop nest shown in Figure 2.4(a)

Figure 2.4 -An illustrative loop nest and its iteration space in polyhedral compilation Iteration Domain. Iteration domain is the collection of all statement instances, represented using a named union set with each named component covering all instances of one statement followed by a set of inequalities for bounds. The iteration domain of the code above can be expressed with (2.3).

Domain = {S 1 (i ):1∑ i ∑ 4; S 2 (i , j ):1∑ i ∑ 4 ^1 ∑ j ∑ i } (2.3)
Access Relations. Access relations are a set of relations coordinating statement instances with the data locations they access, modeled by a set of named union maps together with some inequalities for bounds. (2.4) describes the access relations of the example, consisting of a Write relation and a Read relation. By refining access relations with read and write relations, polyhedral compilation is free to compute dependence relations easily. Furthermore, a write relation can also be split into may-write and must-write relations for the purpose of aggressive optimizations.

Polyhedral Compilation

Write = {S 1 (i ) ! a(i ):1∑ i ∑ 4; S 2 (i , j ) ! b(i , j ):1∑ i ∑ 4 ^1 ∑ j ∑ i } Read = {S 2 (i , j ) ! a(i ):1∑ i ∑ 4 ^1 ∑ j ∑ i } (2.4)
Schedule. As we already mentioned in previous context, polyhedral compilation would map a multi-dimensional logical execution date to each statement instance. Schedule is such a multi-dimensional execution date assigned to a statement instance, expressed using a binary relation between two different multi-dimensional integer tuples. A lexicographically smaller schedule implies an earlier execution of the statement instance. An original schedule is the execution date assigned to a statement instance before schedule transformation. For example, the original schedule of the code in Figure 2.4(a) can be written as (2.5). A new schedule would be computed after schedule transformation if the scheduling algorithm may find a better execution date with respect to parallelism and data locality.

Schedule = {S 1 (i ) ! (i , 0); S 2 (i , j ) ! (i ,1, j )} (2.5)
Dependences. Dependences represent the access conflicts between statement instances, written as named union maps and used for guaranteeing the execution date of a producer be lexicographically smaller than that of the consumer, therefore enforcing the correctness of any transformations enabled by polyhedral compilation. A refinement from traditional statement-level dependences to instance-level dependences in the polyhedral model makes the expression of dependences more complicated, like the constraints after the named maps shown in (2.6) indicating the dependences are described with regard to statement instances.

Dependence = {S 1 (i ) ! S 2 (i , j ):1∑ i ∑ 4 ^1 ∑ j ∑ i } (2.6)

Loop Transformations

Given the polyhedral representations, one may apply any loop transformations and/or their compositions. Loop transformations could be attained by schedule transformations, i.e., reordering the statement instances. Considering the example shown in Figure 2.2, schedule transformation triggers loop skewing by specifying a new execution date to the statement instances of 1D iterated stencil.

If we use the polyhedral representations introduced in the last subsection to express the transformation, we may obtain the original schedule as {S 1 (t , i ) ! (t , i )} and {S 1 (t , i ) ! (t , t + i )} for the new schedule after applying loop skewing. As a result, a transformation could be expressed as (t , i ) ! (t , t + i ) for representing loop skewing. Suppose t and i as variables, the underlying principle of schedule transformation could be interpreted as seeking a coefficient matrix and a constant vector such that

" c 11 c 12 c 21 c 22 #" t i # + " c 10 c 20 # = " t t + i # (2.7)
It is straightforward to solve the above system of linear equations by hand as the problem is heavily simplified for the sake of illustration. However, polyhedral compilation may have to solve it automatically by resorting to an ILP solver; worse yet, the practical problems faced by the polyhedral model would be much more complicated. We would not go further into the underlying structure of schedule transformation but invent the readers to refer to the work of Bondhugula [START_REF] Kumar | Effective Automatic Parallelization and Locality Optimization Using the Polyhedral Model[END_REF] for a detailed mathematical explanation.

To generalize, schedule transformation can be understood as solving a coefficient matrix and a constant vector such that a transformation between two integer tuples can be accomplished. Each row of the coefficient matrix can be interpreted as a hyperplane. Note that the two integer tuples could differ with regard to dimensionality: for example, a scalar dimension could be introduced to achieve loop fusion. Besides, a more complex example of loop transformation would be loop tiling, increasing the dimension of the input tuple by doubling those components requiring tiling.

Schedule transformation by manipulating integer sets broadens the optimization space of polyhedral compilation and simplifies the composition of different loop transformations compared with traditional compilation models like unimodular matrices [Ban93, WL91], while the latter applies loop transformations by means of elementary matrix operations.

Besides, the loop transformations covered by unimodular matrices are also very restricted, including loop interchange, skewing and reversal; the polyhedral model is rather capable for automating a wider set of loop transformations, widening the optimization space by enabling loop fission, fusion, index set splitting [START_REF] Griebl | Index set splitting[END_REF], peeling, strip-mining [START_REF] Kelly | A unifying framework for iteration reordering transformations[END_REF], tiling, unroll and jam [START_REF] Kumar | Effective Automatic Parallelization and Locality Optimization Using the Polyhedral Model[END_REF][START_REF] Bastoul | Improving data locality by chunking[END_REF], unrolling, unswitching, etc.3 A recent work [YGK + 13] also makes it possible to model algorithmic changes which could not be achieved by other techniques, further enriching the transformations of polyhedral compilation.

Schedule Trees

Apart from index set splitting, all the loop transformations modeled by polyhedral compilation could be facilitated by operating on the schedule representation, i.e., named union maps.

This, however, does not mean the operations on named union maps would be ease of use.

A typical drawback of such method can be found when comparing the lexicographic order of two integer sets, since such operations can only be applied on integer sets with the same dimension.

More importantly, the above mentioned schedule representation could not be easily extended to handle non-affine transformations as named union maps hold the "injectivity" and "singlevaluedness" properties. "injectivity" indicates a schedule representation allows different statement instances to be assigned the same logical execution date for expressing inner parallelism; "single-valuedness" refers to a schedule representation would only assign a single execution date to a statement instance, preventing the statement instance from being executed more than once. Clearly, the latter would not allow the implementation of non-affine transformations like overlapped tiling requiring multiple executions of a statement instance.

In this dissertation, we rely on a well-defined schedule representation that would make the expression of non-affine transformations possible in polyhedral compilation. As the schedule construction may decompose a dependence graph recursively and compute a partial schedule for each component independently, a schedule representation would naturally have the form of a tree [START_REF] Grosser | Polyhedral ast generation is more than scanning polyhedra[END_REF]. The schedule representation is thus called "schedule tree".

There have been some schedule representations proposed in the past, including the Kelly's abstraction [START_REF] Kelly | Code generation for multiple mappings[END_REF], "2d+1"-schedules [GVB + 06], etc., that can be viewed as an encoding of schedule trees. Like named union maps, such encoding methods are usually restricted due to missing the ability to facilitate non-affine applications and transformations.

To give an intuitive impression on schedule trees, we depict the schedule tree representation of the code shown in Figure 2.4(a) in Figure 2.5. A schedule tree is constructed by recursively building partial schedule trees which in turn constructed by schedule nodes. For example, the schedule tree in the figure is composed of two sub-trees rooted at one filter node representing statement S 1 (i ) and the other for statement S 2 (i , j ). A partial schedule tree comprises one or more schedule nodes for expressing different semantics. We would next introduce the basic node types in schedule trees. For a complete description of schedule trees and nodes, please refer to [START_REF] Grosser | Polyhedral ast generation is more than scanning polyhedra[END_REF].

Domain.

A domain node in schedule trees is a named union set, appearing as the root of a schedule tree and covering the collection of all statement instances that should be scheduled by the schedule tree. For the sake of simplicity, we would sometimes represent a domain node as "domain" like what we have done in Figure 2.5 rather than writing in the form like (2.3).

Context.

A context node is used to introduce constraints on symbolic constants of the schedule tree. Symbolic constants introduced by a context node could serve as parameters of programs, usually omitted when they are only referenced by the domain node. In practice, parameters passed to compilers like tiling sizes are usually introduced in a context node for determining bounds of new schedules. domain

[S 0 (i ) ! (i ); S 1 (i , j ) ! (i )] sequence {S 0 (i )} {S 1 (i , j )} [S 1 (i , j ) ! ( j )]
Figure 2.5 -The schedule tree of the code in Figure 2.4(a)

Filter. A filter node can be the root of a partial schedule tree, expressed with a named union set and representing the statement instances to be scheduled by its descendants. The partial schedule tree rooted at a filter node has to be retained to another rooted at a domain node, constructing the final schedule tree recursively. A filter node is guarded by a pair of braces.

Band. A band node is used to express the partial schedule on its parent node, written as a named union map. A band node in schedule trees is also integrated with more operations governing code generation, allowing for flexibility for more complicated cases. Named union maps in a band node could be piecewise quasi-affine for expressing schedules like tiling. We use square brackets to denote band nodes in schedule trees.

Sequence/Set. A sequence/set node always appears as the parent node of a group of filter nodes, forcing the children to be executed in a given/arbitrary order. An explicit support for a sequence/set representation makes it possible to break up the instances of a statement into separated parts.

Mark. A mark node can introduce any kinds of information to schedule trees. The use of mark nodes provides a great compatibility to schedule trees with other intermediate representations. For example, one may use a mark node to attain the information about representation mismatching, informing the follow-up code generator to handle this mark node with a custom implementation.

We introduce the basic node types for schedule trees here because we believe they are sufficient to understand the principle of schedule tree representation. Some other node types would be introduced in the following context, together with their uses in our work. The readers may also find more examples of schedule trees throughout the thesis.

A Platform-Neutral Compute Intermediate Language

In this subsection, we would introduce a platform-neutral compute intermediate language, PENCIL [BBC + 15], that we rely on to facilitate non-affine extensions. Combining languages with polyhedral approaches have been proved effective in many applications as we introduced in Section 1. We choose PENCIL as an intermediate representation in addition to polyhedral representations by taking into consideration the following properties.

PENCIL provides a sequential semantic in accordance with the philosophy of widely used programming languages by hiding target-specific hardware information, allowing programmers to follow existing way of programming with such languages by only adding lightweight annotations. This allows for the compatibility of integrating with both general-purpose and domain-specific languages.

PENCIL simplifies the analysis of non-affine expressions and eases the implementation of such extensions. With PENCIL, arrays must be declared through the C99 variable-length array syntax [START_REF] Iso | The ansi c standard (c99)[END_REF]. The C99 type qualifiers/keywords static const restrict or a macro pencil_attributes expending to these type qualifiers/keywords must be used to declare the array function arguments. This allows the polyhedral model know about the length of arrays, and that arrays do not overlap during optimizations.

Pointer declarations and definitions are allowed in PENCIL, but pointer manipulation (including arithmetic) is forbidden except that C99 array references are allowed as arguments in function calls. Pointer dereferencing is neither allowed except for accessing C99 arrays. The restricted use of pointers can essentially eliminate aliasing problems for moving data between different address spaces of hardware accelerators.

A PENCIL for loop must have a single iterator, an invariant start value, an invariant stop value and a constant increment (step). Invariant here requires the value must not change in the loop body. To some extent, such structured for loops may simplify the polyhedral transformations. Considering the fact that recursive calls are not supported by accelerator programming languages like CUDA and OpenCL, recursive calls are excluded from PENCIL. However, we are allowed to extend the semantic of PENCIL for such extensions as long as they are needed.

As shown in Figure 2.6 is the high level overview of PENCIL compilation flow. PENCIL can be the target of a domain-specific compiler, followed by a polyhedral framework, therefore delivering information between a domain-specific language and polyhedral compilation. A typical representative application in DSLs of PENCIL is its use in the early stage of TensorComprehensions [VZT + 18]4 . One is also allowed to write a general-purpose language with PENCIL specifications to model non-affine extensions, extending the polyhedral approaches to handle more complex cases. For example, a combined polyhedral technique with PENCIL was used to handle user-define reductions [START_REF] Reddy | Reduction drawing: Language constructs and polyhedral compilation for reductions on gpu[END_REF] which would not be possible without PENCIL.

With regard to the code generation of the PENCIL compilation flow, there have already been

PART II HANDLING NON-AFFINE APPLICATIONS

Dynamic Counted Loops

While a large number of computationally intensive applications spend most of their time in static control loop nests-with affine conditional expressions and array subscripts, several important algorithms do not meet such statically predictable requirements, going beyond the scope of polyhedral compilation. A well-known non-affine extension to remove the limitation of the polyhedral model in this direction is the work of Benabderrahmane et al. [START_REF] Benabderrahmane | The polyhedral model is more widely applicable than you think[END_REF] for handling general while loops. We are interested in the class of computational kernels involving dynamic counted loops. These are regular counted loops with numerical constant strides, iterating until a dynamically computed, data-dependent upper bound. Such bounds are loop invariants, but often recomputed in the immediate vicinity of the loop they control; for example, their definition may take place in the immediately enclosing loop.

Dynamic counted loops play an important role in numerical solvers, media processing applications, and data analytics, as we will see in the experimental evaluation. They can be seen as a special case of while loop that does not involve an arbitrary, inductively defined termination condition. The ability to substitute their counter with a closed form-an affine induction variable-makes them amenable to a wider set of transformations than while loops. Dynamic counted loops are commonly found in sparse matrix computations, but not restricted to this class of algorithms. They are also found together with statically unpredictable, non-affine array subscripts.

The purpose of this part is to further extend the ability of polyhedral compilation for handling such non-affine applications by enabling a wider set of loop transformations. We will first present the background along this research direction and then introduce our solution in the next two chapters, followed by some experimental results and discussions.

Background and Motivation

The polyhedral framework of compilation unifies a wide variety of loop and array transformations using affine (linear) transformations. The availability of a general-purpose method to generate imperative code after the application of such affine transformations [QRW00, Bas04, GVC15] brought polyhedral compilers to the front scene, in the well-behaved case of static control loops.

Limitations of Previous Work

While significant amount of work targeted the affine transformation and parallelization of while loops [GL94, Col94, GC95, Col95, GGL98, GGL99, BPCB10, JCD + 14], these techniques face a painful problem: the lack of a robust method to generate imperative code from the polyhedral representation. One representative approach to model while loops in a polyhedral framework, and in the code generator in particular, is the work of Benabderrahmane et al. [START_REF] Benabderrahmane | The polyhedral model is more widely applicable than you think[END_REF].

This work uses over-approximations to translate a while loop into a static control loop iterating from 0 to infinity that can be represented and optimized in the polyhedral model. It introduces exit predicates and the associated data dependences to preserve the computation of the original termination condition, and to enforce the proper termination of the generated loops the first time this condition holds. These data dependences severely restrict the application of loop transformations involving a while loop, since reordering of the iterations of the latter is not permitted, and loop interchange is also restricted.

The framework was also not fully automated at the time of its publication, leaving much room for the interpretation of its applicable cases and the space of legal transformations it effectively models. Speculative approaches like the work of Jimborean et al. also addressed the issue [JCD + 14], but a general "while loop polyhedral framework" compatible with arbitrary affine transformations has yet to emerge. In this dissertation, we make a more pragmatic, short term step: we focus on the special case of dynamic counted loops where the most difficult of these problems do not occur.

There has also been a significant body of research specializing on high-performance implementations of sparse matrix computations. Manually-tuned libraries [BAA + 14, BG09, BG11, LBG + 12, MCG04, VDY05] are a commonly used approach, but it is tedious to implement and tune for each representation and target architecture. A polyhedral framework that can handle non-affine subscripts has a greater potential to achieve transformations and optimizations on sparse matrix computations, as illustrated by Venkat et al. [START_REF] Venkat | Loop and data transformations for sparse matrix code[END_REF].

As a result, we would like to propose an automatic polyhedral compilation approach to parallelize and optimize dynamic counted loops that can express arbitrary affine transformations and achieve performance portability. We are allowed to make full use of systems of affine inequalities as implemented in state-of-the-art polyhedral libraries [START_REF] Verdoolaege | Isl: An integer set library for the polyhedral model[END_REF] for our purpose. Moreover, following what has been implemented in the work [SCF03, SLC + 16], we expect not to resort to more expressive first-order logic with non-interpreted functions/predicates such as the advanced analyses and code generation techniques of Wonnacott et al. [START_REF] Pugh | Nonlinear array dependence analysis[END_REF], while

Background and Motivation

avoiding the complexity and overhead of speculative execution.

Static Control Parts

The polyhedral compilation framework was traditionally limited to static control loop nests. It represents a program and its semantics using iteration domains, access relations, dependences and schedules. The statement instances are included in iteration domains. Access relations map statement instances to the array elements they access. Dependences capture the partial order on statement instances accessing the same array element (one of which being a write).

The schedule implements a (partial or total) execution order on statement instances that is compatible with dependences.

Consider the running example in Figure 3.1. The upper bounds, m and n, of the j-loop and k-loop are computed in their common enclosing loop and updated dynamically as the iloop iterates. As a result, it is not possible to classify the whole loop nest as a SCoP, and traditional polyhedral techniques do not directly apply. Tools aiming at a greater coverage of benchmarks-such as PPCG or LLVM/Polly-will abstract the offending inner loops into a black box, greatly limiting the potential for locality-enhancing and parallelizing optimizations.

#pragma scop //begin of our scop for (i=0; i<100; i++) { S0: m = f(i); S1: n = g(i); //begin of the scop of traditional techniques for (j=0; j<m; j++) for (k=0; k<n; k++) S2: S(i, j, k); //end of the scop of traditional techniques } #pragma endscop //end of our scop As an alternative, one may narrow the SCoP by only considering the j-/k-loop nest and treating the dynamic upper bounds as symbolic parameters, enabling polyhedral transformations without problems. This, however, either introduces more frequent synchronizations by exploiting fine-grained parallelism when targeting on CPU targets, or misses the data locality along the outermost loop dimension and the opportunity to exploit full-dimensional parallelism on GPU platforms.

Our Solution

To extend the polyhedral framework to dynamic computed loops, we may need to address the following problems.

Modeling Control Dependences. Undeniably, the polyhedral model in its current form cannot handle dynamic counted loops. We would first derive a static upper bound for such dynamic conditions to make dynamic counted loops amenable to polyhedral compilation. To solve this problem, we may rely on the computation of an affine upper bound for all dynamic trip counts that a given loop may reach, using a combination of additional static analysis and dynamic inspection. Revisiting the polyhedral compilation framework [START_REF] Benabderrahmane | The polyhedral model is more widely applicable than you think[END_REF] of arbitrary while loops, we introduce exit predicates for dynamic counted loops, modeling the control dependence of the original loop through additional data dependences from the definition of these exit predicates to every statement in the loop body.

Achieving Exact Dependence Analysis. Dynamic counted loops are commonly found in sparse matrix computations involving indirect array subscripts, preventing polyhedral model from achieving exact dependence analysis. We leverage the PENCIL language for eliminating alias suspicion and ambiguous analysis in the polyhedral model, allowing for the exact dependence analysis even in the presence of indirect array subscripts.

Eliminating the Effect of Over-approximations. Due to the over-approximation caused by deriving a static upper bound, we need to eliminate the introduced empty iterations for performance improvement. We extend the schedule-tree-based algorithm [START_REF] Grosser | Polyhedral ast generation is more than scanning polyhedra[END_REF] to enable the full automation of imperative code generation after the application of affine transformations, targeting both CPU and GPU architectures.

Our method goes beyond the state of the art [BPCB10, JCD + 14, VHS15] in fully automating the process, specializing the code generation algorithm to the case of dynamic counted loops, and avoiding the introduction of spurious loop-carried dependences or resorting to speculative execution. We conduct experiments on representative irregular computations, including dynamic programming, computer vision, finite element methods, and sparse matrix linear algebra. We validate that the method is applicable to general affine transformations for locality optimization, vectorization and parallelization.

Extension Nodes in Schedule Trees

Our work follows the idea behind the work of Benabderrahmane et al. [START_REF] Benabderrahmane | The polyhedral model is more widely applicable than you think[END_REF] by using overapproximations and modeling control dependences, the latter, however, misses a systematic code generation algorithm. The difficulty to generate early exits for over-approximations is 3.3. An Overview of Our Approach because such statements are not included in the iteration domain modeled by polyhedral compilation. It is difficult to model such statements at the time of the publication of the approach on general while loops. Fortunately, we may now leverage the schedule tree representation for such purpose.

The schedule tree representation has the same expressiveness with traditional polyhedral representations but it allows for the modeling of non-affine extensions in polyhedral compilation. In the case of dynamic counted loops, we may rely on the extension node of schedule trees to introduce additional domain elements to be scheduled.

Recall the mathematical concepts and polyhedral representations we introduced in Subsection 2.1.2, an extension node can be expressed using a named union map relating the outer schedule dimensions with the set of array elements accessed by the statement. In our case, we may abstract an early exit statement as a virtual statement accessing a scalar data named "exit". As a result, we may express a general extension node for such early exit statements as the following

{(d 1 , d 2 , ..., d n ) ! exit()} (3.1)
with (d 1 , d 2 , ..., d n ) representing the outer schedule dimensions and exit() for the data accessed by the statements.

A similar use of extension nodes in PPCG [VCJC + 13] is the creation of data copying statements for locality optimization and the introduction of thread-level synchronization instructions. A statement introduced by an extension node may be scheduled even it is originally excluded by the iteration domain of schedule trees.

An Overview of Our Approach

We may explain our approach by starting with dependence analysis. As shown in Figure 3 We base our formalism and experiments on the schedule tree representation [START_REF] Grosser | Polyhedral ast generation is more than scanning polyhedra[END_REF]. Schedule trees can be flattened into a union of relations form, with each relation mapping the iteration domain of individual statements to a unified logical execution time space.

Since dynamic counted loops cannot be appropriately represented in the iteration domain, a state-of-the-art polyhedral compiler like PPCG may only model the outer loop, abstracting away the j-loop and k-loop, as the schedule tree of Figure 3.3. Following Benabderrahmane's work [START_REF] Benabderrahmane | The polyhedral model is more widely applicable than you think[END_REF], we can derive two static upper bounds, u 1 and u 2 , that are greater than or equal to m and n. The domain and access relations of statement S 2 can be over-approximated accordingly, and represented parametrically in u 1 and u 2 . This representation can be used to compute a conservative approximation of the dependence relation for the whole schedule tree.

Based on this dependence information, one may derive a correct schedule using the Pluto algorithm or one of its variants [BHRS08, VCJC + 13], to optimize locality and extract parallelism. The resulting schedule tree may indeed be seen as a one-dimensional external domain and schedule enclosing a two-dimensional inner domain and schedule controlled by two additional parameters, u 1 and u 2 , as will be seen in Figure 4.4.

The final step is to generate code from the schedule tree to a high level program. The generation of the abstract syntax tree (AST) follows the approach implemented in isl [Ver10], traversing the schedule tree and specializing the code generation algorithm to integrate target-specific constraints, e.g., nested data parallelism and constant bounds. Before encountering a filter node associated with a dynamic counted loop, the exit predicate and its controlled loop body is seen as a single black-box statement by the AST generation algorithm. When passing the filter node constraining the dynamic upper bound, it is necessary to complement the standard code generation procedure with dedicated "dynamic counted loop control flow". This involves either (on GPU targets) the reconstruction of the exit predicate and the introduction of an early exit (goto) instruction guarded by the predicate or (on CPU targets) the replacing the over-approximated static upper bound with the dynamic condition and the removing of the introduced control flow. Our algorithm generates code in one single traversal of the schedule tree1 . domain

[S 0 (i ) ! (i ); S 1 (i ) ! (i ); S 2 (i ) ! (i )] sequence {S 0 (i )} {S 1 (i )} {S 2 (i )}
Figure 3.3 -Original schedule tree of the example

Related Work

The polyhedral framework is a powerful compilation technique to parallelize and optimize loops. It has become one of the main approaches for the construction of modern parallelizing compilers. Its application domain used to be constrained to static control, regular loop nests. But the extension of the polyhedral framework to handle irregular applications is increasingly important given the growing adoption of the technique. The polyhedral community invested significant efforts to make progress in this direction.

A representative application of irregular polyhedral techniques is the parallelization of while loops. The polyhedral model is expected to handle loop structures with arbitrary bounds that are typically regarded as while loops. Collard [START_REF] Collard | Space-time transformation of while-loops using speculative execution[END_REF][START_REF] Collard | Automatic parallelization of while-loops using speculative execution[END_REF] proposed a speculative approach based on the polyhedral model that extends the iteration domain of the original program and performs speculative execution on the new iteration domain. Parallelism is exposed at the expense of an invalid space-time mapping that needs to be corrected at run time.

Beyond polyhedral techniques, Rauchwerger [START_REF] Rauchwerger | Parallelizing while loops for multiprocessor systems[END_REF] proposed a speculative code transformation and hybrid static-dynamic parallelization method for while loops. An alternative, conservative technique, consists in enumerating a super-set of the target execution space [GL94, GC95, GGL98, GGL99], and then eliminating invalid iterations by determining termination detection on the fly. The authors present solutions for both distributed and shared memory architectures.

Benabderrahmane et al. [START_REF] Benabderrahmane | The polyhedral model is more widely applicable than you think[END_REF] introduce a general framework to parallelize and optimize arbitrary while loops by modeling control-flow predicates. They transform a while loop as a for loop iterating from 0 to +1. Compared to these approaches to parallelizing while loops in the polyhedral model, our technique relies on systems of affine inequalities only, as implemented in state-of-the-art polyhedral libraries. It does not need to resort to the first-order logic such as non-interpreted functions/predicates, it does not involve speculative execution features, and it makes dynamic counted loops amenable to a wider set of transformations than general while loops.

A significant body of work addressed the transformation and optimization of sparse matrix computations. The implementation of manually tuned libraries [BAA + 14, BG09, BG11, LBG + 12, MCG04, VDY05] is the common approach to achieve high-performance, but it is difficult to port to each new representation and to different architectures.

Sparse matrix compilers based on polyhedral techniques have been proposed [START_REF] Venkat | Loop and data transformations for sparse matrix code[END_REF], abstracting the indirect array subscripts and complex loop-bounds in a domain-specific fashion, and leveraging conventional Pluto-based optimizers on an abstracted form of the sparse matrix computation kernel. We ought to extend the applicability of polyhedral techniques one step further, considering general PENCIL code as input, and leveraging the semantical annotations expressible in PENCIL to improve the generated code efficiency and to abstract non-affine expressions.

Scheduling Dynamic Counted Loops in the Polyhedral Model

Preparation

A dynamic counted loop with a dynamic counted upper bound and a static lower bound is referred to as the normalized format of dynamic counted loops. For example, the code shown in Figure 3.1 is such a normalized format. However, these loops can be easily normalized by subtracting the lower bound from the upper bound, as shown in Figure 4.1(c). This transformation may introduce non-affine array subscripts since the lower bound may not be affine; we assume the dependence analysis will conservatively handle such subscripts, leveraging PENCIL annotations to refine its precision [CBF95, BBC + 15]; we may also symbolically eliminate identical non-affine expressions on the left and right-hand side.

Deriving a Static Upper Bound

Some forms of while loops may also be modeled, as long as an affine induction variable can be identified and assuming the variant part of the exit condition reduces to this induction variable. For example, the while loop shown in Figure 4.2(a) is extracted from the equake program of SPEC2000 benchmarks. It could be normalized to the format shown in Figure 4.2(b) without changing the semantic of the program.

Deriving a Static Upper Bound

To make a dynamic counted loop amenable to a polyhedral representation, our approach assumes that a static control upper bound u on the dynamic number of iterations is available.

The general idea is that a dynamic counted loop can always be converted into a static for loop enclosing an if statement whose condition checks the dynamic bound. 1 One may determine the u parameter statically or dynamically.

Static Approaches

The u parameter can be approximated statically, as the dynamic upper bounds are functions of outer enclosing loop variables: a typical solution relies on Fourier-Motzkin elimination, projecting out enclosing dimensions and eliminating non-affine constraints.

For instance, the following set of dynamic conditions is extracted from the HOG benchmark of the PENCIL benchmark suite, which we may use Fourier-Motzkin elimination for eliminating the max/min operations and finally deriving a static upper bound.

8 > > > > < > > > > : lb x = max( f x (x), 1) lb y = max( f y (y), 1) ub x = min(g x (x), 1) ub y = min(g y (y), 1) (4.1)
The u parameter can also be determined in other ways, from array size declarations or additional user-defined predicates in PENCIL [BBC + 15]. We use the C99 type qualifiers/keywords static const restrict when declaring an array argument of a PENCIL function, guaranteeing the array argument do not alias and thereby allowing for the static derivation of the u parameter. When such static methods fail, MAXINT or any type-dependent bound remains a valid approximation, but a tighter bound is preferable to avoid lifting induction variables to a wider integral type.

Dynamic Approaches

Besides static analysis, dynamic inspection prior ahead of the loop nest of interest may be practical in some cases. For example, in sparse matrix computations, u may be computed by inspecting the maximum number of non-zero entries in a CSR format. We may infer that the static upper bound of the sparse matrix shown in All in all, affine bounds on the u parameter can generally be derived automatically, at compilation or run time, and the tightness of the approximation does not have an immediate impact on performance.

Modeling Control Dependences

To model control dependences on dynamic conditions, we introduce additional data dependences associated with exit predicates and their definition statements.

An exit predicate definition and check is inserted at the beginning of each iteration of a dynamic counted loop. At code generation time, all statements in the body of the counted loop will have to be dominated by an early exit instruction conditioned by its predicate. This follows Benabderrahmane's method for while loops [START_REF] Benabderrahmane | The polyhedral model is more widely applicable than you think[END_REF], but without the inductive computation and loop-carried dependence on the exit predicate. Of course, we delay the introduction of goto instructions/changing back to the dynamic conditions until code generation, to keep the control flow in a statically manageable form for a polyhedral compiler. For example, the code in The control dependences are therefore converted into data dependences between definition statements and the body of dynamic counted loops. Each statement in a dynamic counted loop is associated with a list of exit predicates. These predicates should be attached to the band node dominating the dynamic counted loop, and will be used to guard or terminate the execution within the over-approximation iteration domain bounded by the u parameters.

Scheduling

The u parameter and conversion of control dependences make it possible to approximate dynamic counted loops in the polyhedral model, at the expense of traversing a larger iteration space. We may thus apply any affine scheduling on this "approximated static control program", to safely compute a correct schedule tree preserving all dependences. Based on the result of scheduling, we may leverage the mark node and extension node for accomplishing transformations for dynamic counted loops while preserving the correctness of the program.

Schedule Construction

The original domain node, as shown in Figure 3.3 can be expressed as

Domain = {S 0 (i ); S 1 (i ); S 2 (i ):0∑ i < 100} (4.2)
Note that statement S 2 is modeled as an one-dimensional statement without our technique.

With the introduction of the u parameter and conversion of control dependences, we may obtain a new domain node as

Domain = {S 0 (i ); S 1 (i ):0∑ i < 100; S 2 (i , j , k):0∑ i < 100 ^0 ∑ j ∑ u 1 ^0 ∑ k ∑ u 2 } (4.3)
and a correct schedule tree representation shown as below.

domain

[S 0 (i ) ! (i ); S 1 (i ) ! (i ); S 2 (i , j , k) ! (i )] sequence {S 0 (i )} {S 1 (i )} {S 2 (i , j , k)} [S 2 (i , j , k) ! ( j ); S 2 (i , j , k) ! (k)]
Figure 4.4 -A correct schedule tree of the example Note that this schedule tree has not been applied any transformations and can be viewed as the original schedule tree representation. The dynamic parameters are assigned at their definition statements, and then virtually read by statement S 2 implicitly guarded by the negation of the exit predicates. This can be modeled as read and write (affine) access relations:

Read = {S 2 (i , j , k) ! m[] | 0 ∑ i < 100 ^0 ∑ j < u 1 ^0 ∑ k < u 2 ; S 2 (i , j , k) ! n[] | 0 ∑ i < 100 ^0 ∑ j < u 1 ^0 ∑ k < u 2 } (4.4) Write = {S 0 (i ) ! n[]; S 1 (i ) ! m[] | 0 ∑ i < 100} (4.5)
According to the variant of the Pluto algorithm implemented in isl [START_REF] Verdoolaege | Isl: An integer set library for the polyhedral model[END_REF], one may set the validity dependences, associated with semantics preservation, to

Validity = (Read °1 ± Write 0 + Write °1 ± Read 0 + Write °1 ± Write 0 ) \(Schedule ¡ Schedule 0 ) (4.6)
and the proximity dependences, associated with locality enhancement, to

Proximity = (Read °1 ± Write 0 + Write °1 ± Read 0 + Write °1 ± Write 0 +Read °1 ± Read 0 ) \ (Schedule ¡ Schedule 0 ) (4.7)
where Schedule represents the original schedule constructed from the original code according to the procedure above, and the 0 (primed) maps distinguish iterations in dependence.

We can then compute a new schedule that applies the variant of the Pluto algorithm using In other words, the scheduling algorithm may safely compute a new schedule, starting from the original one shown in Figure 4.4, preserving all dependences and attempting to minimize the reuse distance.

Schedule Transformation

Once a correct schedule tree representation can be obtained, we are allowed to leverage any types of schedule nodes to apply schedule transformations. Indeed, it is possible to apply transformations and generate code without any special handlings on the current schedule tree. However, the generated code would waste a great number of iterations due to the overapproximations caused by the u parameters and conversion of control dependences, thereby inhibiting performance improvements.

Marking Dynamic Counted Loops. As we introduced before, a mark node can be used for retaining any pieces of information to schedule trees. We therefore are allowed to insert mark nodes above the nodes representing the dynamic counted loops, implying the child node would be considered as dynamic counted loops. As a consequence, one may obtain a schedule tree with a mark node shown in Figure 4.5. The only information we need to retain to a mark node is a string, i.e., "dynamic_counted_loops 00 . The reason behind our intention to leverage a mark node comes from the strategy used in PPCG [VCJC + 13] for generating thread-level synchronization instructions.

One may now use the schedule tree with mark nodes shown in Figure 4.5 for schedule transformation. As there are two separated, tilable band nodes in the schedule tree, a scheduling algorithm would identify the outer band node, corresponding the outermost i -loop in Figure 3.1, leading to a coarser parallelism. Without such abstraction, it would be impossible to model the whole program as a SCoP, as the polyhedral model could only obtain a schedule domain

[S 0 (i ) ! (i ); S 1 (i ) ! (i ); S 2 (i , j , k) ! (i )] sequence {S 0 (i )} {S 1 (i )} {S 2 (i , j , k)} mark: "dynamic_counted_loops 00 [S 2 (i , j , k) ! ( j ); S 2 (i , j , k) ! (k)]
Figure 4.5 -Inserting mark nodes in the schedule tree tree like Figure 3.3 shows, missing dependence information about S 2 and therefore failing to identify the outer tilable band.

Recall that we have mentioned in the previous section, an alternative way is to narrow the SCoP by only considering the j-/k-loop nest and treating the dynamic upper bounds as symbolic parameters, but such modeling strategy may introduce more frequent synchronizations by exploiting fine-grained parallelism or misses the data locality along the outermost loop dimension.

After applying schedule transformation, we need to handle the introduced mark node. One may notice that there exist two dynamic counted loops in the example, but we only introduce one mark node since the j -loop and k-loop are combined into one band node. The Pluto algorithm or its variants would always try to combine loops into one tilable band because such combination would exploit nested parallelism and hereby creating opportunities for more transformations. As transformations have been applied (for example, the scheduling algorithm may strip-mine the outer band or tile the inner band), we are free to split such combined band node. As a result, the schedule tree would be transformed into the following format (we do not represent the transformations that may have been applied in the figure).

In other words, a mark node would be broadcast to each dimension of a combined band node after splitting no matter whether such dimension is a dynamic counted loop. In case where a normal loop is also marked with such mark node, one can determine by checking whether the loop iterator appears in the predicates introduced for conversion of dynamic control. We will explain this issue further in the next section.

Non-affine Extensions. Extension nodes can now be used to replace each occurrence of the mark node. A mark node can only be used to attach additional information but not for custom implementations. As can be found from Figure 4.6, the loop dimension information is not present in a mark node as such information cannot be determined when introducing mark domain

[S 0 (i ) ! (i ); S 1 (i ) ! (i ); S 2 (i , j , k) ! (i )] sequence {S 0 (i )} {S 1 (i )} {S 2 (i , j , k)} mark: "dynamic_counted_loops 00 [S 2 (i , j , k) ! ( j )] mark: "dynamic_counted_loops 00 [S 2 (i , j , k) ! (k)]
Figure 4.6 -Mark nodes with split band nodes nodes. Besides, a mark node may also be broadcast after schedule transformation, making the retaining of further information like loop dimensions impractical at that stage. domain

[S 0 (i ) ! (i ); S 1 (i ) ! (i ); S 2 (i , j , k) ! (i )] sequence {S 0 (i )} {S 1 (i )} {S 2 (i , j , k)} extension: [i , j ] ! exit() [S 2 (i , j , k) ! ( j )] extension: [i , j , k] ! exit() [S 2 (i , j , k) ! (k)]
Figure 4.7 -Replace each mark node with an extension node Figure 4.7 shows the schedule tree with each mark node being substituted by an extension node. Unlike mark nodes, an extension node can be used for recording loop nest information, calling for a custom implementation during code generation. We use an explicit expression in each extension node for the illustrative purpose, followed by the context used in the implementation of schedule trees. An extension node can be represented with a map, relating the loop nest information to an early exit statement which implies an early exit instruction should be emitted.

Emitting early exit statements could be implemented differently depending on the target architecture. Generally, we force the code generator to change the introduced u parameters back to the original dynamic conditions when targeting CPUs, while a goto statement would be introduced for GPUs. The reason behind the different implementations is due to the different programming models used on different targets, as we will explain in detail next section.

Generation of Imperative Code

Once a new schedule is produced, additional transformations can be applied on band nodes, to implement loop tiling or additional permutations, strip-mining for vectorization, etc. Eventually, one needs to return to imperative code through a so-called code or AST generation algorithm. AST generation is a performance-critical step in any polyhedral framework. We extend the code generation scheme of Grosser et al. [START_REF] Grosser | Polyhedral ast generation is more than scanning polyhedra[END_REF], itself derived from the algorithm by Quilleré et al. [START_REF] Quilleré | Generation of efficient nested loops from polyhedra[END_REF] and its CLooG enhancements and implementation [START_REF] Bastoul | Code generation in the polyhedral model is easier than you think[END_REF].

When the Grosser et al. algorithm traverses the band nodes in a schedule, it projects out the local schedule constraints from the domain node. As the dynamic upper bounds are not modeled in the iteration domain (the domain node in the schedule tree and subsequent filter nodes), the generated loops will iterate from 0 to u. It is thus necessary to emit an early exit statement (for GPU architectures) or change the over-approximated static upper bound back to the original dynamic condition (for CPU architectures). Besides, the introduced control flow can also be removed when generating code for CPU targets, reducing the control overhead.

Extending the Schedule Tree

Let us first recall the extensions we made to the schedule tree in the last section. The Grosser et al. algorithm is not able in its original form to generate semantically correct code for our extended schedule tree. However, it can be easily modified to handle the special case of exit predicates that are homogeneous over all statements in a sequence or set node of the schedule tree (e.g., all statements in a band of permutable loops). This is facilitated through the syntactic annotation of dynamic counted loops using so-called mark nodes in the schedule tree. A mark node may attach any kind of information to a subtree; we used it here to specify which band nodes and which dimensions in those bands involve dynamic counted loops. To account for affine transformations combining static and dynamic counted loops (e.g., loop skewing), mark nodes are inserted at every dimension.

One may insert an extension node in a schedule tree to extend its iteration domain, e.g., to

insert a new statement with a specific iteration domain. In our case, we replace each mark node with an extension node, inserting a guard statement with the proper exit predicate. In a first pass, all exit predicates are attached to the band node; a follow-up traversal through the predicate list lets the AST generator detect whether a dimension of the band node is a dynamic counted loop, and position early exits at the right level.

In fact, the code generation scheme is designed to complement the standard code generation procedure with a dedicated "dynamic counted loop template". This template serves as a post-processing step after code generation, involving the reconstruction of the exit predicate.

To this stage, we may have to consider the features of target architectures and branch the code generation template.

We generate CUDA code when targeting on GPU architectures. However, it does not mean our technique is only restricted to such programming APIs. When launching a CUDA kernel, a fix-length loop bound should be specified for execution configuration. We may therefore emit an exit statement every time when an extension node is encountered to force the overapproximated loops terminate correctly. When we implement a custom code generation scheme for CPU architectures, we may generate loops amenable to OpenMP directives, implying an exit statement is not allowed to jump out of the parallel region but a dynamic loop bound is permitted. As a result, we choose to change the introduced u parameters back the original dynamic computed bounds and remove the associated predicates for eliminating control overheads.

Generating Early Exits

When scanning the schedule tree to generate early exits for GPU targets, the AST generator creates a goto AST node for each of the above-mentioned extension nodes. A goto statement can be generated from the AST node using the following steps.

Positioning a goto Statement. As shown in (3.1) and Figure 4.7, an extension node is ex- pressed with a named union map, relating the outer schedule dimensions with the introduced exit statement. Such outer schedule dimensions could be used for positioning the goto of each exit statement.

Generating the Guard. During modeling control dependences, a predicate was introduced at the beginning of each iteration of a dynamic counted loop. Such predicates would be recorded and passed over to the AST generator. A goto statement should also be guarded by a predicate by negating one of the conditions of the introduced predicate dominating the body of the corresponding dynamic counted loop, as the conditions may be a conjunction of multiple dynamic counted loops, e.g., the code shown in Figure 4.3(b). This can be attained by checking the appearance of the loop iterator of current dimension in the conditions.

Changing Back to Dynamic Conditions

Determining Whether the Loop is a Dynamic Counted Loop. As we explained in subsection 4.4.2, the band node with a mark node specifying dynamic counted loop would be split after schedule transformation, with each dimension marked with the same mark node by broadcasting the latter. This broadcast may also mark a static loop as dynamic. Of course, all of the mark nodes would be replaced with extension nodes. When generating a goto statement by checking each occurrence of such extension nodes, the AST generator should determine whether the current dimension is a dynamic counted loop. Following the idea of generating the guard, this can also be achieved by inspecting the appearance of the loop iterator of current dimension in the conditions of the predicate. One may determine the current dimension is a static, normal loop provided the loop iterator is not present in the conditions.

Counting Labels. A label destination is required when using a goto statement. As a result, the AST generator should maintain a global label counter, enforcing the exit statements jump to the corresponding destinations. The label counter is incremented each time a dynamic counted loop is encountered, enforcing uniqueness.

Changing Back to Dynamic Conditions

When targeting on CPU architectures, it may not be allowed to jump in or out of the parallel region using an early exit statement like goto, but one may change the over-approximated static upper bound u back to the original dynamic condition. The information to facilitate such replacement can be attached to an AST annotation node and be the same with those of the goto AST node in GPU case except the label counter. The code generation is similar to the case of GPU case, being accomplished by the following steps.

Positioning Dynamic Counted Loops. Like what we have explained in GPU case, the AST generator has to pick out dynamic counted loops from the band node due to the broadcast of mark nodes. The method has been introduced above, and the AST generator can follow the same scheme used in GPU case for recognizing such dynamic counted loops.

Substituting the u Parameters. The AST code generator may look up the predicate list and extract the condition corresponding to the current dimension. The right-hand side of the condition can be taken out to substitute the introduced u parameter of a dynamic counted loop.

Removing Control Overheads. Similarly, each occurrence of the earlier introduced dynamic conditions at the beginning of each iteration can now be degenerated for eliminating such control overheads.

Changing back to the dynamic condition for a dynamic counted loop is straightforward, but special cares have to be taken to handle cases with multiple associated predicates. One may construct a max operation comprising all the associated predicates as the upper bound of a dynamic counted loop, without removing these introduced control flow since they have to be there to preserve the semantic of the code.

This schedule-tree-based code generation algorithm enables all kinds of loop transformations, with the most challenging one being loop fusion. When fusing two dynamic counted loops, the two sets of predicates are considered, and the early exit statements/max-operation-based dynamic upper bounds are guarded by/composed of their statementwise conjunction/them. As shown in Figure 5.1 is the original and fusion result of two dynamic counted loops. One may conclude from the figure for CPU architectures easily. A normal loop can be treated as a specific case of dynamic counted loop by reasoning on its static upper bound as a predicate.

Unfortunately this scheme efficiently supports a single dynamic counted loop only, and does not deal with the expression of parallelism in these loops.

Flat and Nested Parallelisms

As shown in Figure 4.4, the canonically constructed schedule tree isolates two nested band nodes to represent different levels of the loop nest. This works fine when the target architecture is a shared memory multiprocessor. As an illustrative example, Figure 5.2 is the generated code for a shared memory multiprocessor after the application of loop tiling on the code in Figure 3.1 with the outermost i-loop being parallelized. We also depict the corresponding schedule tree representation in Figure 5.3 for reference.

#pragma omp parallel for for (i=0; i<100; i++) { m = f(i); n = g(i); for (jj=0; jj<m/BB+1; jj++) for (kk=0; kk<n/CC+1; kk++) for (j=0; j<min(m, jj*BB+BB); j++) for (k=0; k<min(n, kk*CC+CC); k++) S(i, jj, kk, j, k); } Figure 5.2 -Code generation with loop tiling for CPU However, when targeting GPU accelerators or producing fix-length vector code, we usually expect to combine nested bands to express parallelism at multiple levels, and a constant iteration count may also be required for data-parallel dimensions. We therefore consider two cases depending on the need to extract parallelism across more than one band. domain

[S 0 (i ) ! (i ); S 1 (i ) ! (i ); S 2 (i , j , k) ! (i )] sequence {S 0 (i )} {S 1 (i )} {S 2 (i , j , k)} extension: [i , j /BB] ! exit() [S 2 (i , j , k) ! ( j /BB)] extension: [i , j /BB, k/KK] ! exit() [S 2 (i , j , k) ! (k/CC)] extension: [i , j /BB, k/CC, j ] ! exit() [S 2 (i , j , k) ! ( j )] extension: [i , j /BB, k/CC, j , k] ! exit() [S 2 (i , j , k) ! (k)]
Figure 5.3 -The schedule tree representation of code shown in Figure 5.2

Flat parallelism within a band

Let us first discuss the case of regenerating imperative code for one or more nested dynamic counted loops within a single band. As a first step, one may systematically generate conditional statements on exit predicates at the innermost level. Figure 4.3(b) shows an example illustrating this approach. The predicates of both loops are included in a single conditional, and generated under the inner loop. Notice that this approach is compatible with affine loop transformations such as loop interchange, not expressible in [START_REF] Benabderrahmane | The polyhedral model is more widely applicable than you think[END_REF] due to the presence of spurious loop-carried dependences.

Our approach is generally applicable in the context of loop interchange except when one attempts to permute a dynamic counted loop with its enclosing affine loop governing the dynamic condition. As we introduce redundant, empty iterations to dynamic counted loops, we inject early exits for eliminating the effect of such over-approximation. In this special interchange case, the introduction of early exits may not guarantee the semantic of original programs. We therefore leave out the introduction of such early exits in this special interchange case. As an illustrative example, Figure 5.4(b) shows the result of permuting the dynamic counted loop with its governing affine loop in the example of Figure 5.4(a). Such treatment may not remove the empty iterations introduced by the over-approximation but the semantic of the program could be preserved. In summary, our approach always tries to remove or minimize the introduced empty iterations but puts the correctness of the program in the first place, implying that we may have to sacrifice the performance by skipping the injection of early exits in some rare cases.

for (i=0; i<N; i++) { m = f(i); for (j=0; j<m; j++) S(i,j); } Yet one still needs to generate early exits in order to avoid traversing a potentially large number of empty iterations. We may extract the iterators one by one from the predicate list and generate the corresponding exit statements from the innermost outwards. The exit predicates are generated in the form of multiple conditionals rather than else branches, as shown in Figure 5.4 and 5.6. Unlike Jimborean et al. [JCD + 14], we do not need speculation on the number of iterations, since we do not deal with general while loops; our technique always executes the same number of iterations as the original programs.

Loop tiling is a special case that should be taken into account. Loop tiling involves the insertion of one or more additional schedule dimensions through strip-mining. When strip-mining a dynamic counted loop, there should be an exit statement at both levels. For the point loop-iterating within a tile-the common case above applies. For the tile loop-iterating among tiles-we align its bounds and strides to follow the structure of the inner loop, so that its counter can also be compared systematically with the same bound.

Nested parallelism across bands

Targeting GPU accelerators or producing fix-length vector code motivates the exploitation of data parallelism within dynamic counted loops, in combination with other nested loops. Since dynamic counted loops result in nested bands in the schedule tree, the combined exploitation of multiple levels of parallelism including one or more dynamic counted loops requires special treatment that is not directly modeled by affine sets and relations. The constraints on the grid of multi-level data parallelism require the collection of bound information across nested bands: when launching a kernel, the parameters of the grid must be known and may not evolve during the whole run of the kernel.

Unfortunately, the statements between nested bands that occur in dynamic counted loops are used to initialize dynamic upper bounds. Statements in the body of these dynamic counted loops depend on those definition statements, through the added dependences modeling the original dependence of the dynamic loop. Still, one can sink these definition statements inside, within the dynamic counted loops, as a preprocessing step. Figure 5.5(a) shows the code after sinking the definition statements of the example in Figure 3.1, followed by a depiction on its schedule tree in Figure 5.5(b).

Note that both dynamic definition statements change into 3-dimensional statements due to the inward movement. In the schedule representation, we use a "band" to represent the band node after such operation, which can be expressed as a piecewise schedule, [{S 0 (i , j , k)

! (i ); S 1 (i , j , k) ! (i ); S 2 (i , j , k) ! (i )}, {S 0 (i , j , k) ! ( j ); S 1 (i , j , k) ! ( j ); S 2 (i , j , k) ! ( j )}, {S 0 (i , j , k) ! (k); S 1 (i , j , k) ! (k); S 2 (i , j , k) ! (k)}].
As a result, the nested bands can be combined again, with no intervening computation or control flow.

The inward movement of these definition statements is safe with the introduction of the upper bound u-parameter. Yet as a side-effect of this movement, each definition will be redundantly evaluated as many times as the number of iterations of the dynamic counted loop itself. This is the price to pay for a fixed upper bound on the iterations.

Once again, this overhead may be mitigated with additional strip-mining of the outer loops, to better control the value of u, effectively partitioning the loop nest into coarse-grain sub computations amenable to execution on a heterogeneous target. Figure 5.6 shows an example after the application of loop tiling on the code in Figure 3.1, and one may also refer to Figure 5.7 for the schedule representation.

As the nested bands are combined into a single one, a polyhedral framework would identify it with multiple dimensions of parallelism, partitioning the loop nest into coarse sub-problems that can be solved independently on heterogeneous platforms. for (i=0; i<100; i++) for (j=0; j<u1; j++) for (k=0; k<u2; k++){ m=f(i); n=g(i); if (j<m && k<n) S(i, j, k); } (a) Sinking the dynamic definitions of the example in Figure 3.

1 domain band sequence {S 0 (i , j , k)} {S 1 (i , j , k)} {S 2 (i , j , k)} (b)
The schedule tree of the code shown in Figure 5.5(a)

Figure 5.5 -Sinking the dynamic definition and its schedule tree representation

General Applicability to Loop Transformations

One of the benefits of our approach with respect to Benabderrahmane et al.'s work [START_REF] Benabderrahmane | The polyhedral model is more widely applicable than you think[END_REF] is its compatibility to various loop transformations. We analyze in this subsection the general applicability of our approach for each loop transformation presented in Subsection 2.1.3 and their combinations.

Loop Transformations of Unimodular Matrices

We first analyze the case of loop transformations covered by unimodular matrices [START_REF] Banerjee | Unimodular Matrices[END_REF][START_REF] Wolf | A loop transformation theory and an algorithm to maximize parallelism[END_REF], i.e., loop interchange, skewing and reversal. One may view unimodular matrices as loop transformations for a single statement because the loop body may always be abstracted as a black box and the structure stays unchanged under such transformations.

Our approach is compatible to loop interchange, as we explained in Subsection 5.5.1 with a special treatment designed for permuting dynamic counted loops with the enclosing affine loop. Similarly, such strategy can also be applicable to loop reversal since the iterations of a dynamic counted loop would be traversed in a reversed order, guaranteeing the correctness of our technique by introducing over-approximations. Fortunately, such cases rarely happen in practice as the control dependences caused by the dynamic conditions prevent such permutation and they are not seen in our experiments.

[S 0 (i , j , k) ! (i /AA); S 1 (i , j , k) ! (i /AA); S 2 (i , j , k) ! (i /AA)]
extension: [i /AA, j /BB] ! exit()

[S 0 (i , j , k) ! ( j /BB); S 1 (i , j , k) ! ( j /BB); S 2 (i , j , k) ! ( j /BB)] extension: [i /AA, j /BB, k/CC] ! exit() [S 0 (i , j , k) ! (k/CC); S 1 (i , j , k) ! (k/CC); S 2 (i , j , k) ! (k/CC)] extension: [i /AA, j /BB, k/CC, i ] ! exit() [S 0 (i , j , k) ! (i ); S 1 (i , j , k) ! (i ); S 2 (i , j , k) ! (i )] extension: [i /AA, j /BB, k/CC, i , j ] ! exit() [S 0 (i , j , k) ! ( j ); S 1 (i , j , k) ! ( j ); S 2 (i , j , k) ! ( j )] extension: [ii/AA, j /BB, k/CC, i , j , k] ! exit() [S 0 (i , j , k) ! (k); S 1 (i , j , k) ! (k); S 2 (i , j , k) ! (k)] sequence {S 0 (i , j , k)} {S 1 (i , j , k)} {S 2 (i , j , k)}

Loop Transformations in Code Generation

As we explained in Subsection 2.1.3, some loop transformations including loop peeling, unrolling, unswitching, are achieved by code generation. Our approach is valid for such cases since the code generator may only change the loop structure instead of reordering statement instances. For example, we are allowed to apply our method on each version after loop peeling and/or unswitching, and the predicate may be introduced before each instance after unrolling.

Similar to the case of unimodular matrices, loop transformations achieved by code generation can also be viewed as transformations applied on loop nests with a single statement, since the loop body of each version after such transformations stays unchanged.

Other Loop Transformations

The analysis on index set splitting [START_REF] Griebl | Index set splitting[END_REF] can follow the case of loop transformations achieved by code generation, as we always split the index set by introducing an affine parameter, implying the loop structure after splitting may stay unchanged.

Strip-mining [START_REF] Kelly | A unifying framework for iteration reordering transformations[END_REF], unrolling-and-jam [START_REF] Kumar | Effective Automatic Parallelization and Locality Optimization Using the Polyhedral Model[END_REF][START_REF] Bastoul | Improving data locality by chunking[END_REF] and loop tiling can be put together as the first two transformations can be viewed as special case of loop tiling. The interchange involved in tiling will not change the order of a dynamic counted loop and its governing loop, neither the iterations of the dynamic counted loop, meaning the introduction of early exits should always be correct.

The solution to loop fusion are discussed in Subsection 5.4 and the validation of the correctness is therefore straightforward. One may see loop fission as a reverse transformation of fusion, and the general applicability of our method for fission is also validated. Loop fusion and fission are transformations that apply on multiple statements since they change the body of the loop nest.

As our method is correct on each loop transformation, it should also be correct on all combinations of these transformations.

Experimental Results

Our framework takes a C program as input, and resorts to PENCIL [BBC + 15] extensions only when dealing with indirect accesses (subscripts of subscripts), implying that all arrays are declared through the C99 variable-length array syntax with the static const restrict qualifiers, allowing PPCG to derive the size of the arrays offloaded on the accelerator despite the presence of indirect accesses, and telling that these arrays do not alias.

We use PPCG [VCJC + 13] to generate target codes, a polyhedral compiler that performs loop nest transformations, parallelization, data locality optimization, and generates OpenCL or CUDA code. The version ppcg-0.05-197-ge774645-pencilcc is used in our work. In a follow-up auto-tuning step, we look for optimal parameter values for tile sizes, block sizes, grid sizes, etc. for a given application and target architecture.

The experiments are conducted on a 12-core, two-socket workstation with an NVIDIA Quadro K4000 GPU. Each CPU is a 6-core Intel Xeon E5-2630 (Ivy Bridge). Sequential and OpenMP code are compiled with the icc compiler from Intel Parallel Studio XE 2017, with the flags -Ofast -fstrict-aliasing (-qopenmp). CUDA code is compiled with the NVIDIA CUDA 7.5 toolkit with the -O3 optimization flag. We run each benchmark 9 times and retain the median value. Median rather than the mean, for more stability. Long discussion there, this is not idea either in general, but more suitable here. Note that the median pushes for an odd number of runs.

Dynamic Programming

Dynamic programming is an alternative method of greedy algorithms to guarantee an optimal solution. In computer science, dynamic programming implies the optimal solution of the given optimization problem can be obtained by the combination of optimal solutions of its subproblems, by solving the same sub-problems recursively rather generating new ones. Dynamic counted loops are usually involved in these problems. We investigate two representative dynamic programming problems-change-making and bucket sort.

Typically, the change-making problem is used to find the minimum number of coins that can add up to a certain amount W and to count how often a certain denomination is used, but it has a much wider application than just currency. The algorithm is also used to count how often a certain denomination is used.

Suppose N denominations are provided, each of which is d i (0 ∑ i < N ). As long as the given amount W > d i , the frequency of the i -th denomination will be incremented by 1. As a result, d i appears as a bound of the inner dynamic counted loop, enclosed by an outer loop iterating over the total number of denominations. Our technique successfully parallelizes the inner dynamic counted loop and generates the CUDA code in conjunction with a loop interchange optimization. We show the performance with different number of denominations N under different amount constraints W in Figure 6.1. It can be concluded from the figure that the performance improvement grows with the rise of the the number of denominations. Bucket sort is a generalization of counting sort, sorting by first scattering the N elements of a given array into a set of M buckets, sorting each bucket individually, and finally gathering the sorted elements in each bucket in order. Due to the comparison operations, a sorting algorithm is inherently not the candidate for parallelization. However, it is possible to parallelize and optimize the gathering step of bucket sort.

We consider a uniform random distribution of elements of the input array. The algorithm has to gather size[i ] elements in the i -th bucket, whose static upper bound can be set as N . The dynamic counted loop controlled by the bucket size is captured by our method and parallelized in the form of CUDA code on GPUs. The performance with different array sizes N and different bucket numbers M is shown in Figure 6.2, indicating the speedup rises along with the increase of the number of buckets involved. 4 8 16 32 64 128 256 512 1024

1 3 5 7 N Speedup M = 16K M = 32K M = 64K M = 128K
Figure 6.2 -Performance of the bucket sort on GPU

HOG Benchmark

The HOG benchmark is extracted from the PENCIL benchmark suite1 , a collection of applications and kernels for evaluating PENCIL compilers. The distribution of intensity gradients or edge directions describe the local object appearance and shape within an image. When processing an image, the HOG descriptor divides it into small connected regions called cells. A histogram of gradient directions is then compiled for the pixels within each cell. The descriptor finally concatenates these histograms together. The descriptor also contrast-normalize local histograms by calculating an intensity measure across a block, a larger region of the image, and then using this value to normalize all cells within the block to improve accuracy, resulting in better invariance to changes in illumination and shadowing.

The kernel of the HOG descriptor contains two nested, dynamic counted loops. The upper bounds of these inner loops are defined and vary as the outermost loop iterates. The dynamic parameter is an expression of max and min functions of the outer loop iterator and an array of constants. We derive the static upper bound parameter u from the BLOCK_SIZE constant, a global parameter of the program to declare the size of an image block.

Since we target a GPU architecture, we ought to extract large degrees of parallelism from multiple nested loops. As explained in the previous section, we sink the definition statements of dynamic parameters within inner dynamic counted loops and apply our AST generation scheme for a combined band for GPU architecture. We may then generate the CUDA code with parameter values for tile sizes, block sizes, grid sizes, etc. We show performance results with and without host-device data transfer time, in Figure 6.3, considering multiple block sizes.

The detection accuracy improves with the increase of the block size. Our algorithm achieves a promising performance improvement for each block size, and our technique can obtain a speedup ranging from 4.4£ to 23.3£ while the PENCIL code suffers from a degradation by we apply a non-affine shift by the dynamic lower bound as discussed earlier. The maximum number of non-zero entries in a row is the static upper bound and may be set as the u parameter. It can be derived through an inspection. As a result, the references of indirect array subscripts can be sunk under the inner dynamic counted loop, exposing a combined band in the schedule tree.

Venkat et al. [START_REF] Venkat | Nonaffine extensions to polyhedral code generation[END_REF] generate two optimized versions for the CSR format on GPU target. The first version parallelizes the outer loop by strip-mining the latter into two dimensions with the outer dimension mapped to thread blocks and the inner dimension to threads, mapping each row of the sparse matrix to a thread; the inner loop is executed sequentially. This version is referred to as "scalar". The second version, i.e., "vector", also parallelizes the outer loop but assigns multiple threads to each row of the sparse matrix; in addition, the inner loop is further strip-mined to allow for intra-wrap reduction. In both cases, only the outer loop is identified as a permutable band while the "vector" version exploits intra-wrap reduction by assigning multiple thread to a row, parallelizing the inner loop in without combining it with its outer loop and missing loop transformations across loop nest.

Our technique can identify the inner dynamic counted loop and parallelize both loops, exposing a higher degree of parallelism and allowing for loop tiling rather than (nested) strip-mining when comparing with Venkat et al.'s scalar/vector version. Tiling on CSR format may reduce the impact of thread imbalance issue when the input sparse matrix is unstructured, i.e., the number of non-zero entries in each varies greatly. We introduce an atomic operation in the generated code for preserving the correctness of reduction of sparse matrix computations.

We show the performance in Figure 6.5, using the matrices obtained from the University of Florida sparse matrix collection [START_REF] Davis | The university of florida sparse matrix collection[END_REF] as input, and the properties of the input matrices are listed in Table 6.1. We also show the performance of a manually-tuned library-CUSP [START_REF] Bell | Implementing sparse matrix-vector multiplication on throughput-oriented processors[END_REF] in the figure for a comparison with hand-written implementations.

Overlapped Tiling

A large number of efforts have been taken to improve data locality and parallelism for computationally intensive applications especially for iterated stencils, resulting in various loop tiling transformations including simple shapes like rectangular, parallelogram [BHRS08, VCJC + 13] and complex shapes like overlapped, split [KBB + 07], diamond [START_REF] Bondhugula | Diamond tiling: Tiling techniques to maximize parallelism for stencil computations[END_REF], hexagonal [GCH + 14], etc.

The polyhedral framework has been brought to the front scene by its ability to analyze and optimize general-purpose loop nests. Its main scheduling and code generation algorithm remain limited to classical tile shapes however, leading to inefficient wavefront inter-tile parallelism with pipelined startup, while custom solutions with more complex tile shapes exist to exploit inter-tile parallelism along the axes of the iteration domain and data while improving data locality.

There have been several loop tiling techniques specializing on time-iterated stencils, either constrained to constant dependence vectors [KBB + 07] or being difficult to extend to other areas like image processing pipelines [BBP17, GCH + 14]. Image processing pipelines are a class of computations arising in computer vision and computational photography, consisting of a directed acyclic graph of filtering stages and edges representing dependences between stages. The stages in a pipeline usually exhibit abundant data parallelism but require locality optimization to achieve high performance, making manually exploiting a difficult task.

PolyMage [START_REF] Teja Mullapudi | Polymage: Automatic optimization for image processing pipelines[END_REF] is the state-of-the-art polyhedral compilation framework automatically generating high-performance schedules for such image processing pipelines, benefiting from the full inter-tile parallelism enabled by overlapped tiling [KBB + 07]. It takes a DSL inspired by Halide [RKBA + 13] as input, computes manually-written-competitive schedules for image processing pipelines, and generates high-performance imperative code. The PolyMage framework implements overlapped tiling by finding bounding hyperplanes of a tile, implying these bounding hyperplanes have to be as inclined as possible to preserve the dependence vectors at all levels of the pipeline, leading to a looser tile shape than expected.

In this part, we intend to implement a schedule-tree-based overlapped tiling technique. Compared with the "DSL+compiler" approach of the PolyhMage framework, we still leverage the PENCIL language and the well-defined polyhedral representation for achieving such nonaffine transformation. Due to the introduction of overlapped regions, a statement instance may be assigned multiple execution dates, violating the "single-valuedness" property hold by many previous polyhedral representations. We would show how such violation can be preserved and implemented in the schedule tree representation.

Background and Motivation

Loop Tiling

Loop tiling has been integrated into polyhedral compilation frameworks, being implemented as a post-scheduling transformation for exploiting data locality and parallelism. Typically, a polyhedral compiler, e.g., PPCG, abstracts the input program as sets and maps defined by systems of affine inequalities, before constructing an affine schedule respecting all dependences carried by statement instances to better exploit data locality and expose parallelism. A follow-up tiling transformation is performed automatically, embedding the computation into a higher dimensional space of tile and point dimensions, but currently limited to classical, rectangular or parallelogram tile shapes.

Decoupling loop tiling from the scheduler may prohibit tile-level concurrent start. An alternative way allowing full inter-tile parallelism involves a tighter coupling of loop tiling and affine scheduling, like the Pluto compiler does by introducing diamond tiling [START_REF] Bondhugula | Diamond tiling: Tiling techniques to maximize parallelism for stencil computations[END_REF]. Tile-level concurrent start along a face of the iteration space is possible if there are no inter-tile dependences parallel to this face, forcing the face to be evicted or linearly independent from the candidate bounding hyperplanes of the scheduler1 .

The schedule found by such evictions could not be better than those found by the standard scheduling algorithm with respect to the dependence distance minimizing cost function, implying a different schedule has to be used for intra-tile parallelism and/or vectorization, complicating the scheduling process and follow-up code generation.

An alternative way to eliminate pipelined startup and drain is to modify the tile shape obtained by existing polyhedral compilation frameworks [KBB + 07]. Overlapped tiling is constructed by adding an additional (shaded) region to the left of the tile obtained by existing frameworks, jointing consecutive tiles for exploiting inter-tile parallelism. The shaded regions between consecutive tiles have to be recomputed. Split tiling is obtained by splitting the overlapped tile shape into two sub-tiles, with one being the shaded region and the other consisting of all the remaining points executed in order. Each sub-tile can be executed concurrently with those of neighboring tiles. As an illustrative example, Figure 7.1 shows the comparison between

Expansion Nodes in Schedule Trees

tile. In the example shown in Figure 7.3, the (dashed) slopes caused by dependences between stages A and B are extended, constructing a looser shaded region between tiles, leading to more redundant computations than one expects.

Our solution

We propose an overlapped tiling technique for eliminating the redundant (circled) points in shaded regions. To construct a tighter overlapped tile, we first fuse the stages in this example and let a general polyhedral framework perform a rectangular tiling on the iteration space regardless of the correctness. We then expand the bounding faces of a tile by taking into consideration the inter-tile dependences, without necessitating further extending the slopes between stages A and B . This can be implemented by modifying the schedule tree intermediate representation [START_REF] Grosser | Polyhedral ast generation is more than scanning polyhedra[END_REF] in general polyhedral frameworks. As a result, we can construct a bounding face like the dashed slopes shown in Figure 7.3, eliminating the redundant points from shaded regions. Such overlapped tile is referred to as acute trapezoid tile.

One may also construct an overlapped tile by first resorting to the scheduler of a polyhedral compiler, transforming the iteration space into the form shown in Figure 7.4. In this case, a PolyMage-like technique may construct an overlapped tile shape like the solid slopes show, but we can still achieve a tighter shape by (1) resorting to a scheduler to shift the iteration space, (2) performing a rectangular tiling and (3) extending the left bounding face of a tile by operating on the schedule tree representation, still eliminating the redundant points from shaded regions. This shape is referred to as right trapezoid tile.

We design and implement our approach in a source-to-source polyhedral compiler targeting on image processing pipelines written in a general-purpose language. As a result, our method would be neither restricted to a domain-specific language nor does it introduce sophisticated rescheduling and custom code generation in a polyhedral framework. We leverage the schedule tree representation instead, and this allows us to construct tighter tile shapes than the state of the art, minimizing the shaded region of redundant computations across overlapped tiles.

Our technique also goes beyond the state of the art by generating code for both generalpurpose multicores and heterogeneous accelerators. We validate its general applicability by conducting additional experiments on iterated stencils, providing a comparison between overlapped tiling and other state-of-the-art techniques.

Expansion Nodes in Schedule Trees

As we introduced before, the "single-valuedness" property is usually hold by many existing polyhedral representations, preventing such encoding methods from expressing non-affine transformations like overlapped tiling. The reason is the overlapped regions introduced by overlapped tiling require the polyhedral representation may map more than one execution by the Pluto scheduling algorithm or its variants is missing concurrent start, limiting the performance of the generated code especially those stencil-based applications.

Overlapped [START_REF] Pananilath | An optimizing code generator for a class of lattice-boltzmann computations[END_REF]. Unlike overlapped tiling and split tiling, diamond tiling may work with arbitrary affine dependences. The introduction of scheduling to find tiling dimensions does not only complicate the scheduling but also increase the code generation time in practice. We show in our evaluation section that our technique could achieve competitive performance on iterated stencils with diamond tiling by carefully selecting tile sizes. Our technique is also applicable to image processing pipelines.

Hybrid hexagonal/classical tiling was proposed by Grosser et al. [GCH + 14] to exploit full inter-tile parallelism of iterated stencils on GPU architectures. It can be seen as a generalization of diamond tiling, allowing for partial concurrent start by constructing a hexagonal tile shape along the time and the first space dimension and classical tiling along the other space dimensions. Grosser et al. [START_REF] Grosser | The relation between diamond tiling and hexagonal tiling[END_REF] also show a comparison between diamond tiling and hexagonal tiling. We compare our technique with hexagonal tiling in our experiments on GPU architectures.

By revisiting overlapped tiling in polyhedral compilation frameworks especially for image processing pipelines, we construct a much tighter overlapped tile shape for improving the performance of such applications. Halide [RKAP + 12, RKBA + 13] is a domain-specific language for image processing pipelines, decoupling algorithms from schedules for easy optimizations for such benchmarks, allowing users to experiment with schedules without touching the algorithms. Manually or autotuning approaches [AKV + 14] to finding schedules usually takes a long time to facilitate fascinating performance. The polyhedral model is a promising solution to automatically search schedules for image processing pipelines by integrating with transformations like overlapped tiling, fusion, scratchpad allocation, etc.

Acute Trapezoid Tiling and Right Trapezoid Tiling

Overlapped tiling is an efficient tiling technique allowing for tile-level concurrent start. As we described above, one may choose to obtain an acute or a right trapezoid tile by fusing or shifting first, with both achievable by operating on the schedule tree representation.

Representations in Schedule Trees

As we would leverage schedule trees to express both acute trapezoid tiles and right trapezoid tiles, let us recall the schedule tree representation and explain how we can make use of it. In polyhedral compilation, schedules in the polyhedral model are used to define the execution order of programs, including both the original and those generated by scheduling algorithms like the Pluto scheduler or its variants. A schedule has in nature the form of a tree, making an explicit tree representation have the same expressiveness with previous encoding methods but simplify the implementation of non-polyhedral operations.

A statement instance is expressed by a named multi-dimensional vector with the name identifying the statement and the coordinates corresponding to iteration variables of the enclosing loops. The collection of all statement instances, i.e., the iteration domain, is expressed using Presburger formulas [START_REF] Pugh | Static analysis of upper and lower bounds on dependences and parallelism[END_REF], retained in a domain node. For example, the iteration domain of the code shown in Figure 7.2 can be expressed as

{S A (i ):1∑ i < N ; S B (i ):2∑ i < N °1; S C (i ):4∑ i < N °3}
with loop boundaries included.

A statement instance is also mapped to a multi-dimensional logical execution date [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem. part ii. multidimensional time[END_REF] for defining its lexicographic execution order. Such mapping is referred to as a schedule, expressed by a piecewise multi-dimensional quasi-affine function over the iteration domain and contained in a band node. A band node is derived from tilable band in the Pluto framework [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF], defining permutability and/or parallelism properties on a group of statements as well. A rectangular tiling regardless of the correctness of the example in Figure 7.2 is written as

[{S A (i ) ! (i /T ); S B (i ) ! (i /T ); S C (i ) ! (i /T )}, {S A (i ) ! (i ); S B (i ) ! (i ); S C (i ) ! (i )}]
with the former piece representing tile loops (iterating among tiles) and the latter representing point loops (iterating within tiles).

A filter node selects a subset of statement instances introduced by an outer domain/filter node. Filter nodes usually appear as children of a sequence/set node expressing a given/arbitrary order on its children. As an illustrative example, Figure 8.1 is the original schedule tree of the example shown in Figure 7.2, indicating each filter node selects a subset (a stage) of the domain node. The sequence node defines the three stages should be executed in order, followed by a loop iterating over the iteration space of each stage.

domain sequence {S A (i )} [S A (i ) ! (i )] {S B (i )} [S B (i ) ! (i )] {S C (i )} [S C (i ) ! (i )]
Figure 8.1 -The original schedule tree of the code in Figure 7.2

There have been many other node types existing in schedule trees [START_REF] Grosser | Polyhedral ast generation is more than scanning polyhedra[END_REF], among which we focus on the expansion node as introduced in subsection 7.2. An expansion node can expand a statement instance to one or more instances, constituting a new set of statement instances to be scheduled by the schedule tree. A standard loop tiling would partition the iteration space into smaller blocks, each of which is disjoint with each other, making it difficult to construct an overlapped tile without an expansion node. With an expansion node, we are free to choose one statement instance of a stage in a tile and expand it to as many instances as we expect. These expanded statement instances would joint with those of neighboring tiles, resulting in overlapped tiles over the whole iteration space.

Acute Trapezoid Tiling

Let us first consider implementing overlapped tiling by modifying both sides of a tile. Figure 8.2 shows the result of rectangular tiling on the iteration space of the example listed in Figure 7.2 regardless of the correctness. Such rectangular tile can be obtained by first strip-mining

Schedule Generation

We can now generate schedules for both acute and right trapezoid tiling. Considering acute trapezoid tiling, strip-mining is performed on each stage, followed by fusing the tile loops iterating among tiles. Correctness is enforced by introducing expansion nodes below the point loops of stage A and B iterating within a tile, with boundaries updated as explained in the last subsection. We finally get a schedule tree for acute trapezoid tiling as shown in Figure 8.6. domain

[S A (i ) ! (i /T ); S B (i ) ! (i /T ); S C (i ) ! (i /T )] sequence {S A (i )} expansion: [S A (i ) ! S A ( j )] [S A (i ) ! (i )] {S B (i )} expansion: [S B (i ) ! S B ( j )] [S B (i ) ! (i )] {S C (i )} [S C (i ) ! (i )]
Figure 8.6 -Schedule tree of acute trapezoid tiling If we let a scheduling algorithm transform the iteration space of the pipeline, we may get a shifted schedule before tiling, leading to a two-dimensional tilable band. However, we need to separate the point loops from the band, making it possible for introducing expansion nodes to stages A and B , as the generated schedule tree shown in Figure 8.7. The code generator of a polyhedral compilation framework can take these schedule trees for code generation, with overlapped tiling enabled in the generated code.

domain [S A (i ) ! (i /T ); S B (i ) ! ((i + 1)/T ); S C (i ) ! ((i + 3)/T )] sequence {S A (i )} expansion: [S A (i ) ! S A ( j )] [S A (i ) ! (i )] {S B (i )} expansion: [S B (i ) ! S B ( j )] [S B (i ) ! (i + 1)] {S C (i )} [S C (i ) ! (i + 3)]

Removal of Control Overheads

Let us still consider the stage B of the code in Figure 7.2 and suppose the tile size be T. An expansion node in schedule trees is used to map each domain element of a domain/filter node to one or more images, forming a wider set of elements scheduled by the domain/filter node. Expanding the point loop may change a rectangular/parallelogram tile obtained from existing polyhedral compilation frameworks into an overlapped tile, avoiding necessitating a scheduling algorithm like the PolyMage framework does. We may therefore insert an expansion node below a filter node representing the point loop as shown in Figure 8.6 and 8.7. An expansion node here is represented as a map expanding the original domain S B (i ) to its image set S B ( j ).

Introduced for iterating the images of the expansion node, j is an unbounded parameter in the schedule tree, producing an unbounded domain for the subtree of stage B . One may take into consideration the (upper and lower) bounds of j by adding the original bounds on all elements of stage B , guaranteeing the images of the expansion node would not exceed the original bounds. We can write it formally as

lb ∑ j ∑ ub (8.1)
where lb and ub represent the lower and upper bounds.

This condition alone cannot produce an overlapped tile, as the boundaries of a tile are still remain unchanged. To expand the boundaries of a rectangular/parallelogram tile as explained in previous subsections, we let b(i /T )c denote the greatest integer less than or equal to its parameter, representing the tile number in our case. The original points executed by the (b(i /T )c + 1)-th tile (b(i /T )c is 0 for the first tile) can be expressed by

T £b(i /T )c∑ j ∑ T £ (b(i /T )c + 1) °1 (8.2)
We can change it into

T £b(i /T )c°l ∑ j ∑ T £ (b(i /T )c + 1) °1 + r (8.3)
where l represents the number of expanded points introduced by expansion nodes to the left boundary, and r to the right boundary. In practice, l or r may be 0, e.g., right trapezoid tiling.

Expanding the first point of a stage in a tile is straightforward, simplifying schedule transformations when changing a rectangular/parallelogram tile into an overlapped one. We find, however, that the selection of the point from which the images of an expansion node should be generated may have heavy impact on the control overheads in generated code.

One may obtain a bounded map representing an expansion node by conjuncting constraints 8.6. Comparing the Two Trapezoid Tile Shapes (8.1) and (8.3), implying the point at T + 2forL-tile and the point at T for R-tile are selected for expansion, isolating a partial tile (L-tile) from full tiles (R-tile and all the remaining if they exist)1 and introducing more control overheads in generated code. This isolation may lead to a performance degradation when there are much more stages in a group.

A better solution to remove control overheads in generated code is integrating partial tiles with full tiles. Each starting point of stage B executed by a full (other than the first) tile can be expressed as T £b(i /T )c, while T £b(i /T )c + 2(T £b(i /T )c is equal to 0) being used for the partial (the first) tile. One is free to choose any point of a stage covered by a tile for generating images of expansion nodes without changing the meaning of the generated code, implying one may choose any point other than T £b(i /T )c in a full tile. We thus choose the point at T £b(i /T )c + 2 as starting point for full tiles, enforcing the uniqueness with the partial tile and removing control overheads.

Finally, one may write another condition constraining the selection of starting point of full tiles as

i = T £b(i /T )c + s (8.4)
where s represents the shifting of the starting point in each tile due to the bounds on the whole domain. We may finally obtain an expansion node by constraining the map of an expansion node with a set of conditions consisting of (8.1), (8.3) and (8.4) as follows.

{S B (i ) ! S B ( j ):(8.1) ^(8.3) ^(8.4)} (8.5)

Similarly, we may also obtain the map of expansion nodes like (8.5) and the schedule tree as shown in Figure 8.7, with r in (8.3) being equal to 0, removing control overheads for right trapezoid tiling.

Comparing the Two Trapezoid Tile Shapes

Given a schedule tree with expansion nodes written as (8.5), one may obtain an acute trapezoid tile by first fusing the stages of an image processing pipeline or a right trapezoid tile by first shifting the iteration space. One may distinguish the difference between the two trapezoid tile shapes by comparing Figure 8.3 and Figure 8.5.

Apart from the shape, the two tile shapes also differ from each other with respect to data locality. For the sake of simplicity, we first show the generated code with different tile shapes of the example in Figure 7.2 in Figure 8.8 and Figure 8.9.

Complementary Transformations

One of the main purposes of our technique is to enhance the performance of image processing pipelines by optimizing overlapped tiling in polyhedral compilation frameworks. Following the domain-specific code generator PolyMage, we may also need some other transformations for further improving the performance of such benchmarks.

On the other hand, to cover a wider coverage of benchmarks, we also target on iterated stencil code. We design our technique to handle all cases of iterated stencil code, including multidimensional cases and multiple statements as explained in the last subsection. Exploiting full-dimensional parallelism in practice may not be necessary, especially when targeting multi-dimensional cases on general-purpose multicores for single-thread performance. We therefore also need to take it into consideration.

Alignment and Scaling

Constructing overlapped tiling is allowed only in the case of constant dependence vectors, making it not straightforward to exploit inter-tile parallelism in practice when heterogeneous stages exhibiting different dimensions and/or complex access patterns are grouped together.

Alignment and scaling of stages can be introduced, following the PolyMage framework, to achieve near-neighbor dependences, changing the dependence vectors into constant.

Alignment can be achieved by introducing a scalar dimension in dependence vectors, followed by shifting among dimensions for eliminating non-constant dependences. Up-sampling and down-sampling are introduced for scaling schedules of stages appropriately, obtained by multiplying a scaling factor for each stage in a group.

Fusion

Loop fusion is an important transformation implemented by polyhedral compilers for exploiting data locality. Benefiting from alignment and scaling, some stages of the pipeline may be fused together, creating opportunities for exploiting overlapped tiling across more stages.

One can make full use of the fusion heuristic adopted by a polyhedral compiler by setting compilation options, but it may not be good enough for image processing pipelines even with the aggressive fusion heuristic. The criteria the PolyMage code generator provides is fusing a successor stage with its only child when it has only one child by viewing the pipeline as a directed acyclic graph consisting of nodes representing stages and edges denoting dependences between stages, followed by iteratively attempts for fusing opportunities until no fusion can be found. We reimplement this heuristic in our technique.

Reducing Memory Footprint

Loop fusion transforms the pipeline into several groups, each of which consists of a set of intermediate stages and an output stage, requiring storage allocation optimization for improving performance.

Those values produced by intermediate stages are only used within a tile, implying they can be discarded when they are not live after the computation of the tile. These intermediate values can therefore be allocated in small scratchpad memory rather than full buffers, leading to better locality and improving the performance when integrating with overlapped tiling and the transformations mentioned above. Indexing expressions generated for such scratchpads can be determined according to the conditions defined in expansion nodes.

Hybrid Tiling

When targeting on multi-dimensional iterated stencils, we are allowed to restrict overlapped tiling only to the time dimension and a subset of space dimensions, leveraging existing polyhedral frameworks to perform rectangular/parallelogram tiling on the remaining space dimensions, just like diamond tiling [START_REF] Bondhugula | Diamond tiling: Tiling techniques to maximize parallelism for stencil computations[END_REF] and hybrid hexagonal/classical tiling [GCH + 14]. This lower dimensional overlapped tiling may change the tile shapes for a multi-dimensional case, as we show in Figure 8.11 comparing the difference between a full and partial dimensional overlapped tile shapes for 2d stencil code. A full dimensional overlapped tiling would form a base of a pyramid, extending along both dimensions of the space, while a partial one only overlaps along one dimension. We compare the performance of our generated code with both the naïve and optimized version generated by PolyMage. Table 9.1 also shows the speedup of our code over the optimized version of PolyMage by only differentiating the overlapped tiling shapes, validating the effectiveness of our tiling techniques. We also show the performance of the Halide1 manual schedule written by experts and automatic scheduling algorithm [MAS + 16] inspired from PolyMage. The OpenCV version we use is 2.4.9.1. We set the overlapped tile sizes the same as the optimized version of the PolyMage framework. Acute trapezoid tiling is used in our experiment, as the data locality exploited by right trapezoid tiling does not make sense in this case, as we have already explained.

Bilateral Grid

Bilateral Grid [CPD07, PKT + 09] is a data structure enabling fast edge-aware image processing, localizing operations involved including bilateral filtering, edge-aware painting, local histogram equalization, etc. This benchmark smoothes images while preserving their edges by first constructing a bilateral grid and then sampling the grid along each dimension, producing

Iterated Stencils

To validate the general applicability of our technique, we also conduct experiments on three representative iterated stencils. The detailed information about the examples and tiling sizes we use in this subsection is list in Table 9.2, with the execution times of each sequential code shown in the last column. One can obtain the execution time of each CPU/GPU case by combining with Figure 9.8, 9.9, 9.10 and 9.12. We use right trapezoid tiling for experimenting because it has a better locality for iterated stencils than acute trapezoid tiling. We first run the sequential code of the stencils and record the execution time as a baseline reference. We compare the performance with state-of-the-art diamond tiling enabled by the Pluto compiler and parallelogram tiling enabled default by PPCG, showing the speedups of different techniques.

The 1/2/3d-heat benchmarks we use in this subsection are iterated stencils solving the heat equations, iteratively updating data element using three-point and five-point stencils respectively. When selecting tile sizes, we follow the sizes chosen by diamond tiling [START_REF] Bondhugula | Diamond tiling: Tiling techniques to maximize parallelism for stencil computations[END_REF]. However, one may have to choose the sizes of parallelogram tiling carefully. The diamond tiling paper selects 1024£1024 for parallelogram tiling when comparing the performance, preventing the tiles from executing in wavefront parallelization because the tile size along time dimension is greater than the iteration time. We therefore choose to sacrifice locality to enable wavefront parallelization and ensure enough tiles on the wavefront.

One flaw of overlapped tiling is the redundant computation caused by shaded regions between neighboring tiles. One may thus have to select tile sizes for constructing a sharp overlapped tile, minimizing redundant computations as much as possible. We find 32£8192, 4 £ 64 £ 64 and 4 £ 16 2 £ 256 in practice for these stencils as the best tile sizes. The selection of diamond tile sizes follows its publication. We show the performance comparisons of these stencils in Figure Figure 9.8, 9.9 and 9.10, demonstrating overlapped tiling may achieve similar performance to diamond tiling by enabling inter-tile parallelization but without introducing complex rescheduling step and compilation overhead penalty.

Note that the slight performance gap between overlapped tiling and diamond tiling is due to the recomputation nature of overlapped tiling although an implementation of the latter can save compilation time. The purpose of our experiment by comparing with diamond tiling PART IV

CONCLUSIONS 10 Conclusions and Perspectives

In this chapter, we conclude the thesis with an overview of our contributions and a prospect of future work.

Conclusions

The polyhedral model has become a powerful candidate for compilation tools achieving automatic parallelization and locality optimization. Its success is due to its impressive effectiveness on a variety of computational problems, and to its encouraging integration efforts into both research and industry compilers. Its applicability has also witnessed several extensions, resolving various challenges in domain specific areas and at different scales. However, polyhedral compilation was long criticized for its restrictions to affine applications and transformations.

In this dissertation, we overcome some static control limitations of the polyhedral model from the following angles.

Handling Non-affine Applications. We first studied the parallelizing compilation and optimization of an important class of loop nests where counted loops have a dynamically computed, data-dependent upper bound. Such loops are amenable to a wider set of transformations than general while loops. To achieve this, we introduce a static upper bound and model control dependences on data-dependent predicates by revisiting a state-of-the-art framework to parallelize arbitrary while loops. We specialize this framework to facilitate its integration in schedule-tree-based affine scheduling and code generation algorithms, covering all scenarios from a single dynamic counted loop to nested parallelism across bands mapped to GPUs with fixed-size data-parallel grids.

Our method relies on systems of affine inequalities, as implemented in state-of-the-art polyhedral libraries. It takes a C program with PENCIL functions as input, covering a wide range of non-static control application encompassing the well studied class of sparse matrix computations. The experimental evaluation using the PPCG source-to-source compiler on representa-tive irregular computations, from dynamic programming, computer vision and finite element methods to sparse matrix linear algebra, validated the general applicability of the method and its benefits over black-box approximations of the control flow.

Automating Non-affine Transformations. As the second contribution of the dissertation, we revisited overlapped tiling in a polyhedral compilation framework for optimizing image processing pipelines and iterated stencils. These classes of computations exhibit abundant data parallelism but require locality optimizations for improving performance. Our technique allows for tighter overlapped tile shapes than the state of the art, further improving the performance of such pipelines on both general-purpose multicores and heterogeneous accelerators by integrating with transformations including alignment and scaling of stages in the pipeline, loop fusion, scratchpad allocation, hybrid tiling, etc. Our technique can generate both acute and right trapezoid tile shapes, and has been implemented in the PPCG source-to-source compiler running on a general-purpose C input. We validated the general applicability of the approach and its benefits over a state-of-the-art framework.

Future Work

The thesis extends the application domain of polyhedral compilation to non-affine cases by combining the polyhedral model with a well-defined intermediate language, allowing for aggressive program transformations and automatic code generation strategies without resorting to speculative optimizations nor custom time-consuming rescheduling algorithms.

There still exist a large number of remaining challenges along the non-affine extensions.

Dynamic Control Extensions to Deep Learning. In this dissertation, we studied the dynamic control issue in polyhedral compilation. The polyhedral model nowadays has been integrated with deep learning by automatically converting a high-level description of convolutional/recurrent networks into high-performance implementations amenable to the target processors and accelerators [VZT + 18, BRR + 18]. Such polyhedral-model-based software stacks are widely welcomed as a promising solution to optimizing custom operators that do not fit existing library calls. What is missing is the support of dynamic control extensions in such compiler stacks, restricting the applicable cases to dense tensor operations only.

The taco library [KKC + 17] models sparse tensor algebra and allows for the automatic code generation of such operations. However, the library still remains being a prototype implementation for sparse tensors, leaving a variety of performance-crucial optimizations including tiling, autotuning, etc. that have already been implemented in or integrated with the polyhedral-model-based tools outside their frameworks. More importantly, the missing of code generation for heterogeneous platforms is still an opening issue. The parallelization and optimization of dynamic counted loops may inspire existing polyhedral tensor compiler stacks for supporting dynamic control cases, achieving an end-to-end compilation flow of 10.2. Future Work sparse tensor operations for both homogeneous and heterogeneous architectures.

Trade-off between Memory Locality and Redundant Computation. Tiling and fusion are two effective transformations for exploiting locality optimization. While a much larger space of different compositions of loop transformations could be exploited by the polyhedral model, the difficulty to reason about trade-offs between different criteria is also heavily exacerbated.

In spite of the redundant computations, overlapped tiling is chosen to benefit from inter-tile parallelism while preserving locality in image processing pipelines due to its ability to allow for aggressive storage optimization. However, such redundant computations would be heavily exacerbated when fusing excessive stages in an image pipeline.

By constructing tighter trapezoid shapes, we successfully minimized the introduced recomputations in overlapped tiling. A fusion heuristic based on dynamic programming [START_REF] Jangda | An effective fusion and tile size model for optimizing image processing pipelines[END_REF] was also proposed for optimizing the fusion strategy. Overlapped tiling would not be effective without storage reduction, while the latter in turn could be exploited by fusing as many stages as possible but result in more recomputations between tiles. However, we may still expect for a well-defined cost model to maximize locality and parallelism while minimizing redundant computations by integrating overlapped tiling and fusion heuristic.

Support to Dynamic High Performance Languages. In a similar way, the evolution of programming languages in the field of scientific computing also has to make a compromise with regard to ease of optimization, generality of coverage, difficulty in development, support of dynamism, etc. The scope of polyhedral compilation has been extensively widened by increasing the kinds of supported languages, from general-purpose languages to domainspecific languages with regard to the generality of coverage or from high-level descriptions to intermediate representations in respect of abstraction level of programming philosophy.

An interesting possible direction to further broaden the scope of polyhedral compilation is to model dynamic features of languages like C++ or the recently released Julia [START_REF] Bezanson | Julia: A fresh approach to numerical computing[END_REF].

Particularly, the invention Julia integrates the merits of numerous programming languages, raising a variety of new challenges to optimizing system software. The dynamic control extensions made in this dissertation open the door to addressing new research problems in such dynamic languages by combining an intermediate language with the polyhedral model, providing a flexible compilation flow for lowering high-level, dynamic features into static controls amenable to polyhedral compilation and expanding the support to runtime dynamic checks.

Scalability of Polyhedral Compilation. The polyhedral compilation was long considered outside the domain of real-world applications due to its missing to handle dynamic control and non-affine programs. We provided a systematic way for handling dynamic, non-affine applications and generalizing non-affine transformations, evaluating the scalability of the idea behind our techniques to a great number of real-world applications. Still, the performance of our generated code falls behind manually written libraries in some cases, leaving further room to optimize existing techniques. For example, we may further extend the framework to leverage the dynamic parallelism APIs provided by CUDA programming; we may also consider the compositions of overlapped tiling with other transformations like fusion and unrolling as mentioned before.

Finally, compilation time complexity also limits the scalability due to the underlying principles of polyhedral methods. Recently, some efforts attempted to relax the complexity of the underlying integer linear programming problem [START_REF] Acharya | Polyhedral autotransformation with no integer linear programming[END_REF] or leveraging statement clustering methods [START_REF] Mehta | Improving compiler scalability: Optimizing large programs at small price[END_REF]. Nonetheless, such approaches rely heavily on exact dependence analyses which is still a very time-consuming step. A well-defined frontend capable of relaxing the time complexity of both dependence analysis and schedule transformation might be a promising solution, possibly in conjunction with a helper intermediate language for optimization purposes.

  Figure 1.1 -The trend of clock frequency and number of cores per chip in the past 50 years

Figure 3

 3 Figure 3.1 -Example with dynamic counted loops

  Figure 4.1 -A sparse matrix computation and its normalized format

  Figure 4.1(a) shows an example of such computations using Compressed Sparse Row (CSR) format. One may represent the sparse matrix as in Figure 4.1(b), with the additional arrays for storing the information about dynamic conditions.

Figure 4

 4 Figure 4.2 -A while loop and its normalized format

  Figure 4.1(b) is 2 with an inspection. Alternatively, one may think about the transformation changes the sparse matrix into the form of

  Figure 4.3(a) is preprocessed as the version in Figure 4.3(b) before constructing the affine representation.

Figure

  Figure 4.3 -Conditional abstraction

New_Schedule = schedule

 schedule Domain under Schedule respecting Validity and minimizing Proximity (4.8)

  Figure 5.4 -An interchange example

  domain extension: [i /AA] ! exit()

Figure 5 . 7 -

 57 Figure 5.7 -The schedule tree representation of the code shown in Figure 5.6

W

  Figure 6.1 -Performance of change-making on GPU

Figure 8 . 7 -

 87 Figure 8.7 -Schedule tree of right trapezoid tiling

Figure 8

 8 Figure 8.11 -Full and partial dimensional overlapped tiling

Table 1 .

 1 1 -The memory hierarchy of NVIDIA Tesla GPUs in the past five years

	Tesla Product	Tesla K40	Tesla M40	Tesla P100	Tesla V100
	GPU	GK180 (Kepler) GM200 (Maxwell) GP100 (Pascal) GV100 (Volta)
	Memory Size	Up to 12 GB	Up to 24 GB	16 GB	16 GB
	L2 Cache Size	1536 KB	3072 KB	4096 KB	6144 KB
	Shared Memory Size	16/32/48 KB	96 KB	64 KB	96 KB *
	Register File Size/SM	256 KB	256 KB	256 KB	256 KB
	Register File Size/GPU	3840 KB	6144 KB	14336 KB	20480 Kb

* The shared memory size of GV100 (Volta) is confinable up to 96 KB.

Table 1

 1

		.2 -List of xPUs used for AI Accelerators
	Abbreviation	Full name	Manufacturer Released year
	APU	Accelerated Processing Unit	AMD	2011
	BPU	Brain Processing Unit	Horizon Robotics	2017
	DPU	Deep Learning Processing Unit Dataflow Processing Unit	Deephi Tech Wave Computing	2016 2017
	EPU	Emotion Processing Unit	Emoshape	2017
	HPU	Holographic Processing Unit	Microsoft	2017
		Intelligence Processing Unit	GraphCore	2017
	IPU	Intelligence Processing Unit	Mythic	2018
		Image Processing Unit	Google	2017
	NPU	Neural Network Processing Unit	Vimicro	2016
	SPU	Stream Processing Unit	AMD	2006
	TPU	Tensor Processing Unit	Google	2016
	VPU	Vision Processing Unit	Intel	2016
	ZPU	Zylin CPU	Zylin AS	2015

  was complained due to the missing of considering communication

	Chapter 2. Background
	for (t=0; t<T; t++)
	for (i=1; i<N-1; i++)
	A[t+1][i]=0.25*(A[t][i+1]+2.0*A[t][i]+A[t][i-1]);
	(a) 1d iterated stencil
	t
	i
	(b) Original iteration space of the code in Figure 2.2(a)
	t
	t + i
	(c) Transformed iteration space of the code in Figure 2.2(a)

  .1, statement S 2 does not have data dependences on other statements. However, there are output dependences among definition statements of dynamic parameters m and n. To faithfully capture the scheduling constraints, one should also model the control dependences of S 2 over both headers of the enclosing dynamic counted loops. Such control dependences can be represented as data dependences between the definition statements of dynamic upper bounds and S 2 .To establish such a dependence relation, an exit predicate may be introduced before each statement of the loop body, like in the framework of Benabderrahmane et al.[START_REF] Benabderrahmane | The polyhedral model is more widely applicable than you think[END_REF]. The resulting dependence graph is shown in Figure3.2. The solid arrows represent the original (output) dependences between definition statements of dynamic parameters, and the dashed arrows represent the data dependences converted from the exit conditions' control dependences.

	e s 0 !s 0		e s 1 !s 1
	S 0		S 1
	e s 0 !s 2	S 2	e s 1 !s 2
	Figure 3.2 -Dependence graph of the example
	By capturing control dependences as affine relations from the definition of exit predicates to
	dominated statements in loop bodies, one may build a sound abstraction of the scheduling

constraints for the loop nest. This technique is applicable to arbitrary while loops, in conjunction with a suitable code generation strategy to recover the exact control flow protected by the exit predicate, and by over-approximating the loop upper bound as +1. This is the approach explored by Benabderrahmane et al., but the resulting polyhedral representation is plagued by additional spurious loop-carried dependences to update the exit predicate, removing many useful loop nest transformations from the affine scheduling space. In the more restricted context of dynamic counted loops, it is possible to eliminate those loop-carried dependences as the exit predicate only depends on loop-invariant data.

Table 6

 6

	.1 -Summary of the input sparse matrices
	Matrices	Symmetric # of Nonzero Entries Rows£Columns
	cant	yes	4007383	62451£62451
	consph	yes	6010480	83334£83334
	cop20k_A	yes	2624331	121192£121192
	mac_econ_fwd500	no	1273389	206500£206500
	mc2depi	no	2100225	525825£525825
	pdb1HYS	yes	4344765	36417£36417
	Pres_Poisson	yes	715804	14822£14822
	pwtk	yes	11634424	217918£217918
	rma10	no	2374001	46835£46835
	tomographic1	no	647495	73159£59498

Our method beats the scalar version of Venkat et al. in most cases due to the higher degree of parallelism. Benefiting from the intra-wrap parallel reduction, the vector version of Venkat PART III AUTOMATING NON-AFFINE TRANSFORMATIONS

  tiling and split tiling [KBB + 07] were proposed targeting on concurrent start by modifying the tile shape obtained by a Pluto-like scheduler. The latter should be implemented by either splitting the shaded region of overlapped tiles or introducing more difficult scheduling process for finding different bounding faces of a tile. No implementations of overlapped tiling on general-purpose platforms have been reported though it was implemented in domain-specific compilers for image processing pipelines [MVB15, RKBA + 13] or stencil code generator[START_REF] Holewinski | High-performance code generation for stencil computations on gpu architectures[END_REF] with OpenCL [ZGG + 12]. There was also an implementation of split tiling for iterated stencils on GPU architectures [GCK + 13]. Comparing with these approaches, our technique covers a wider application domain, improving performance over the state of the art by constructing a much tighter overlapped tile.Davis et al.[START_REF] Davis | Transforming loop chains via macro dataflow graphs[END_REF] proposed to construct the fusion-based and shifting-based overlapped tiles for improving performance of applications involving stencil computations. The image processing pipeline used in their work only involves stencil operations but not sampling nor histogram operations. Our work covers all the basic operations involved in image processing pipelines. The work of Davis et al. was neither trying to minimize the redundant computation of overlapped tiling. Bondhugula et al.[START_REF] Bandishti | Tiling stencil computations to maximize parallelism[END_REF][START_REF] Bondhugula | Diamond tiling: Tiling techniques to maximize parallelism for stencil computations[END_REF] proposed a generalizing formalism for diamond tiling in the polyhedral model by introducing a rescheduling step in the Pluto compiler. There has been a great amount of work [EM90, GV15, MHLK17, OG09, SGM + 15, SSPS11] reported on the evaluation of diamond tiling. It was also generalized to handle iterated stencils defined over periodic data domain with index set splitting technique [BBC + 14] and Lattice-Boltzmann method

Table 9

 9

		.1 -Summary of the PolyMage benchmark
	Benchmark	Stages	CPU execution time (ms)	Speedup over PolyMage
			naïve (1 core)	PolyMage (32 cores)	Our work (32 cores)	(32 cores)
	Bilateral Grid	7		66.01		5.57	5.41	1.03
	Camera Pipeline	32	116.32		5.95	5.87	1.01
	Harris Corner Detection	11	246.88		5.10	5.10	1.00
	Local Laplacian Filter	99	480.48		35.35	27.08	1.31
	Multiscale Interpolation	49	209.10		20.07	16.44	1.22
	Pyramid Blending	44	350.49		17.18	15.41	1.11
	Unsharp Mask	4	142.16		5.01	3.68	1.36
	Benchmark	Stages GPU execution time (ms) Speedup over PPCG
				PPCG	Our work
	Bilateral Grid		7	4.67	4.25	1.10
	Camera Pipeline		32	4.29	2.18	1.97
	Harris Corner Detection	11	1.89	1.07	1.77
	Local Laplacian Filter	99	16.73	11.12	1.50
	Multiscale Interpolation	49	12.86	7.76	1.66
	Pyramid Blending		44	8.33	5.78	1.44
	Unsharp Mask		4	2.17	1.29	1.68

Table 9

 9 .2 -Problem sizes and tile sizes of the iterated stencils

		Problem sizes		Tile sizes		Baselinetime (s)
			standard diamond overlapped
	heat-1d 1000£160000 128 £ 1024 2048 2	32£8192	6.99
	heat-2d 1000£4000 2	16£64 2	64 3	4 £ 64 2	7.17
	heat-3d	100£150 3	16 3 £ 256 16 3 £ 256 4 £ 16 2 £ 256	1.21

TensorComprehensions and TVM here are used to refer to the DSLs rather than the compiler stacks.

An affine transformation should be "signle-valued", i.e., an one-to-one mapping function of the integer points on iteration space.

A hyperplane is the projection of an n dimensional space on its n °1 dimensional sub-space.

A warp is a set of threads arranged lengthwise on a loom and crossed by the woof.

Some of these loop transformations, i.e., loop peeling, unrolling, unswitching, are achieved by code generation rather than schedule transformation, since these transformations change the loop structure rather than reorder statement instances.

PENCIL was introduced for bridging the Halide IR and polyhedral representations in the prototype implementation but was later removed from the framework due to simplification considerations. However, it helps the developers construct the early prototype implementations and lays a solid foundation for follow-up development of the framework.

Another difference with[START_REF] Benabderrahmane | The polyhedral model is more widely applicable than you think[END_REF] where multiple traversals were needed.

This is easier than a general while loop, since the dynamic bound check remains continuously false after its first falsification.

https://github.com/pencil-language/pencil-benchmark

A hyperplane can also be understood as the higher dimensional analog of a face in 3D space.

A full tile is completely contained in the iteration space while a partial one is not but has a non-empty intersection with the iteration space [KRR + 07].

Version: commit: 8c23a1970faba9b06bf7145d2653618fb978479e.

Acknowledgements

To be honest, I had always been doubting about the impressive performance and the surprising magic of the polyhedral model before I joined the PARKAS group. The obscure mathematical abstractions of the model always prevented me from getting to the bottom of the polyhedral world. I therefore came to a decision to go into this world by pursuing a PhD degree majoring this research direction. Today, I eventually reach the end of this path, with the guidance, advices, supports and assistances from a large number of persons to name.

Acknowledgements

Mary W. Hall and Cédric Bastol. Thanks to Mary W. Hall for the help to a better comprehension about sparse computations, with an in-depth look into your approach and the NVIDIA implementations of the applications. Thanks to the practical comments and criticisms from Cédric Bastol for helping improve the quality of this work.

Remerciements Acknowledgements

Je suis très reconnaissant au temps précieux et aux efforts des rapporteurs de ma thèse, Mary W. Hall et Cédric Bastol. Merci à Mary W. Hall pour l'aide à une meilleure compréhension des calculs clairsemés, avec un examen approfondi de votre approche et des implémentations NVIDIA des applications. Merci aux commentaires pratiques et aux critiques de Cédric Bastol pour avoir aidé à améliorer la qualité de ce travail.

Code Generation for a Single Loop

The final step is converting the AST to a high level program. When a goto AST node of a dynamic counted loop is captured, a goto statement conditioned by its predicates is enforced after the loop body, as well as a label destination after the loop itself. The associated predicates are gathered in a conjunction and wrapped as one conditional, with loop iterators instantiated according to the loop level. A label is inserted after each dynamic loop as a target for a goto statement.

for (i =0; i<N; i++) { for (j =0; j<u1; j ++) { m=f(i); if(j<m) S1(i,j); } } for (i =0; i<N; i++) { for (j =0; j<u2; j ++) { n=g(i); if(j<n) S2(i,j); } } (a) Before fusion for (i =0; i<N; i++) { for (j =0; j< max (u1 ,u2); j ++) { m=f(i); n=g(i); if(j<m) S1(i,j); if(j<n); S2(i,j); if(j >=m && j >=n) goto label0 ; } label0 :; } In fact, the injection of early exits may only be affected by the iteration reordering of a dynamic counted loop, i.e., loop reversal, and/or the change of the dynamic condition, i.e., the interchange with the governing affine loop. As a result, the validation of the correctness of our method on loop skewing is straightforward, as no such transformations happen in skewing and the introduced predicates before each iteration of the dynamic counted loop would also be updated with respect to the result of skewing.

To conclude, our method on dynamic counted loops are always correct for loop transformations covered by unimodular matrices, and any combinations of these transformations.

Chapter 8. Acute Trapezoid Tiling and Right Trapezoid Tiling

In Figure 8.8, the point loops of individual stages are distributed. In other words, the tiles of individual stages are executed one after another in an acute trapezoid tile shape. On the contrary, the point loops of these stages are fused in Figure 8.9, minimizing the intra-tile producer-consumer relation distances.

for (c0=0; c0<N/T+1; c0++){ for (c1=max(T*c0-3,1);c1<min(T*(c0+1)+3,N); c1++)

-Code generated by right trapezoid tiling

The right trapezoid tile should be preferred as it holds a better intra-tile data locality than the acute one. However, such data locality has a very little impact on the performance of the generated code when targeting on image processing pipelines since the tile height is usually not large enough for benefiting from the former. Iterated stencils may exploit such data locality when the tile size along the time dimension is large enough, as we will explain in the experiment.

Handling Multi-statement/-dimensional Cases

Multiple Statements

The ability to handle multiple statements is important since such cases may happen in practice, as we will show in the evaluation. Suppose we have a similar schedule tree as shown in Figure 8.6 and let the filter node of stage B consist of two statements, S B 1 and S B 2 . For the sake of simplicity, we do not show the context in expansion node of stage A.

Following the cases of single statement, we schedule S B 1 and S B 2 as a macro statement S B , fusing the tile loop with those of stage A and C , and obtain a schedule tree regardless of the correctness similar to Figure 8.6. An expansion node is allowed to insert right underneath the filter node consisting of S B 1 and S B 2 , followed by a band node representing the map of their point loops and its child sequence node defining their execution order. We finally obtain a schedule tree as shown in Figure 8.10. Right trapezoid tiling can be handled in the same way. domain

Figure 8.10 -Schedule tree of multiple statements

Multi-dimensional Statements

Considering the cases of multi-dimensional statements, we are allowed to handle a single dimension each time without impacting on other dimensions, implying that handling the cases of multi-dimensional statements is a process of invoking the cases of one-dimensional statement recurrently, with both acute and right trapezoid tiling cases being taken into account.

Evaluations

We conduct experiments on both the PolyMage benchmarks and representative iterated stencils. The image processing pipelines covered by the PolyMage benchmarks are extracted from the Halide benchmarks, varying widely in structure and complexity.

Experimental Setup

We implement both trapezoid tiling techniques as well as the complementary transformations explained in Section 8.8 in PPCG, a polyhedral compiler that exploits parallelism and data locality of programs automatically. All transformations are applied automatically by default when passing --acute and --right flags to PPCG for switching between acute and right trapezoid tiling. The PPCG version we use is ppcg-0.07-26-g236d559. PPCG can take C programs as input, automatically generating OpenMP code on general-purpose multicores and CUDA code on heterogeneous accelerators.

The experiments are conducted on a 32-core, dual-socket workstation with an NVIDIA Quadro P6000 GPU. Each CPU is a 2.10GHz 16-core Intel Xeon(R) E5-2683 v4. We use the icc compiler (18.0.1) from Intel Parallel Studio XE 2018, with the flags -fast -qopenmp. CUDA code is compiled with the NVIDIA CUDA (9.1.95) toolkit with the -O3 optimization flag. Each benchmark is executed 11 times, of which the first run is discarded and the average of the remaining is recorded.

Image Processing Pipelines

The PolyMage framework takes a DSL as input and generates both naïve and optimized OpenMP codes. A naïve version is generated by PolyMage without schedule transformations and overlapped tiling, of which the sequential code is used as the baseline and also as the input PPCG. The image processing benchmarks used in our experiments are listed in Table 9.1, together with the number of stages and the actual execution times of each baseline. One Résumé De nos jours, l'optimisation des compilateurs est de plus en plus mise à l'épreuve par la diversité des langages de programmation et l'hétérogénéité des architectures. Le modèle polyédrique est un puissant cadre mathématique permettant aux programmes d'exploiter la parallélisation automatique et l'optimisation de la localité, jouant un rôle important dans le domaine de l'optimisation des compilateurs. Une limite de longue date du modèle réside dans sa restriction aux programmes affines à contrôle statique, ce qui a entraîné une demande émergente de prise en charge d'extensions non affines. Cela est particulièrement aigu dans le contexte d'architectures hétérogènes où une variété de noyaux de calcul doivent être analysés et transformés pour répondre aux contraintes des accélérateurs matériels et pour gérer les transferts de données à travers des espaces mémoire. Nous explorons plusieurs extensions non affines du modèle polyhédral, dans le contexte d'un langage intermédiaire bien défini combinant des éléments affines et syntaxiques. D'un côté, nous expliquons comment les transformations et la génération de code pour des boucles avec des limites de boucle dynamiques non dépendantes des données et dynamiques sont intégrées dans un cadre polyédrique, élargissant ainsi le domaine applicable de la compilation polyédrique dans le domaine des applications non affines. D'autre part, nous décrivons l'intégration du pavage en recouvrement pour les calculs de pochoir dans un cadre polyhédral général, en automatisant les transformations non affines dans la compilation polyhédrique. Nous évaluons nos techniques sur des architectures de CPU et de GPU, en validant l'efficacité des optimisations en effectuant une comparaison approfondie des performances avec des frameworks et des librairies écrites à la pointe de la technologie.

Mots Clés

Programmation parallèle, compilation polyédrique, parallélisation automatique Abstract Nowadays, optimizing compilers are increasingly challenged by the diversity of programming languages and heterogeneity of architectures. The polyhedral model is a powerful mathematical framework for programs to exploit automatic parallelization and locality optimization, playing an important role in the field of optimizing compilers. A long standing limitation of the model has been its restriction to static control affine programs, resulting in an emergent demand for the support of non-affine extensions. This is particularly acute in the context of heterogeneous architectures where a variety of computation kernels need to be analyzed and transformed to match the constraints of hardware accelerators and to manage data transfers across memory spaces. We explore multiple non-affine extensions of the polyhedral model, in the context of a welldefined intermediate language combining affine and syntactic elements. On the one hand, we explain how transformations and code generation for loops with non-affine, data-dependent and dynamic loop bounds are integrated into a polyhedral framework, extending the applicable domain of polyhedral compilation in the realm of non-affine applications. On the other hand, we describe the integration of overlapped tiling for stencil computations into a general polyhedral framework, automating non-affine transformations in polyhedral compilation. We evaluate our techniques on both CPU and GPU architectures, validating the effectiveness of the optimizations by conducting an in-depth performance comparison with state-of-the-art frameworks and manually-written libraries.