A Domain-Extensible Compiler with Controllable Automation of Optimisations
Thomas Kœhler
Submitted in fulfilment of the requirements for the
Degree of Doctor of Philosophy
School of Computing Science
College of Science and Engineering
University of Glasgow

2022

## Abstract

In high performance domains like image processing, physics simulation or
machine learning, program performance is critical. Programmers called
_performance engineers_ are responsible for the challenging task of optimising
programs. Two major challenges prevent modern compilers targeting
heterogeneous architectures from reliably automating optimisation. First,
domain-specific compilers such as Halide for image processing and TVM for
machine learning are difficult to extend with the new optimisations required
by new algorithms and hardware. Second, automatic optimisation is often unable
to achieve the required performance, and performance engineers often fall back
to painstaking manual optimisation.

This thesis shows the potential of the Shine compiler to achieve domain-
extensibility, controllable automation, and generate high performance code.
_Domain-extensibility_ facilitates adapting compilers to new algorithms and
hardware. _Controllable automation_ enables performance engineers to gradually
take control of the optimisation process.

The first research contribution is to add 3 code generation features to Shine,
namely: synchronisation barrier insertion, kernel execution, and storage
folding. Adding these features requires making novel design choices in terms
of compiler extensibility and controllability. The rest of this thesis builds
on these features to generate code with competitive runtime compared to
established domain-specific compilers.

The second research contribution is to demonstrate how extensibility and
controllability are exploited to optimise a standard image processing pipeline
for corner detection. Shine achieves 6 well-known image processing
optimisations, 2 of them not being supported by Halide. Our results on 4 ARM
multi-core CPUs show that the code generated by Shine for corner detection
runs up to 1.4$\times$ faster than the Halide code. However, we observe that
controlling rewriting is tedious, motivating the need for more automation.

The final research contribution is to introduce _sketch-guided equality
saturation_ , a semi-automated technique that allows performance engineers to
guide program rewriting by specifying rewrite goals as sketches: program
patterns that leave details unspecified. We evaluate this approach by applying
7 realistic optimisations of matrix multiplication. Without guidance, the
compiler fails to apply the 5 most complex optimisations even given an hour
and 60GB of RAM. With the guidance of at most 3 sketch guides, each 10 times
smaller than the complete program, the compiler applies the optimisations in
seconds using less than 1GB.

###### Contents

  1. Abstract
  2. Acknowledgements
  3. Declaration
  4. 1 Introduction
  5. 2 Motivating Background
    1. 2.1 High Performance Programming
    2. 2.2 The Compiler Extensibility Challenge
    3. 2.3 The Compiler Controllability Challenge
    4. 2.4 Conclusion
  6. 3 Code Generation in a Domain-Extensible Compiler
    1. 3.1 The Rise Language & Shine Compiler
    2. 3.2 Implicit Barrier Insertion for Shine
    3. 3.3 Explicit Kernel Execution for Shine
    4. 3.4 Explicit Storage Folding for Shine
    5. 3.5 Conclusion
  7. 4 Beyond Halide Scheduling with Controlled Rewriting
    1. 4.1 Domain-Extensibility with Control of Optimisations
    2. 4.2 The Harris Corner Detection Case Study
    3. 4.3 Optimising the Harris Corner Detection with Elevate
    4. 4.4 Evaluation of Runtime Performance
    5. 4.5 Conclusion
  8. 5 Sketch-Guided Equality Saturation
    1. 5.1 Background on Equality Saturation
    2. 5.2 Motivation for Semi-Automatic Optimisation
    3. 5.3 Sketch-Guided Equality Saturation
    4. 5.4 Efficient Equality Saturation for the Lambda Calculus
    5. 5.5 Evaluation of Lambda Calculus Encoding
    6. 5.6 Evaluation of Sketch Guidance
    7. 5.7 Conclusion
  9. 6 Discussion
    1. 6.1 Summary
    2. 6.2 Limitations
    3. 6.3 Ongoing & Future Work
  10. A Optimised Harris Corner Detection
    1. A.1 Comparable to Halide Reference
    2. A.2 Beyond Halide Reference
  11. B Optimised Matrix Multiplication
  12. C Handwritten Sketches and Selection of Discovered Programs
    1. C.1 Matrix Multiplication Sketches
    2. C.2 Rise Programs for the parallel Matrix Multiplication

## Acknowledgements

This thesis is the result of 3 short years of research, and 6 long months of
write-up. I am grateful to those who helped me throughout this journey.

Thanks to my supervisors, Michel Steuwer and Phil Trinder, for their
encouragements, guidance, and constructive feedback. Michel, you always cared
about my personal growth, I owe you for attending many interesting events and
learning to confidently promote my work. Phil, your helicopter view helped me
to focus on the bigger picture, strengthened my research and writing. I am
still using the mug from your last supply drop.

Thanks to everybody in the School of Computing Science in Glasgow for creating
such a welcoming environment. The pandemic made me realise how much I missed
the social lunches, coffee breaks and pub nights. I had the pleasure of many
refreshing interactions with Adrian Ramsingh, Dejice Jacob, Cristian Urlea and
Kyle Simpson in our shared F101 office.

Thanks to the Halide developers and egg developers for maintaining open-source
code, to Max Willsey for sharing his insights on equality saturation.

Thanks to the entire Lift, Rise and Elevate teams for their work and our
shared moments. Bastian Köpcke and Federico Pizzuti, you developed the Shine
compiler at my side, and shared many pre-pandemic whiteboard sessions.
Christof Schlaak, you unexpectedly shared the post-pandemic write-up
experience with me (or was it a table tennis competition?).

Above all, thanks to my loved ones, who gave me the strength to withstand this
journey. Gözel Shakeri, my partner, you gave meaning to work-life balance,
celebrated my successes, and relativized my failures. Madeleine Nœuveglise and
François Kœhler, my parents, you supported me unconditionally as I moved far
away from home. My friends, stay who you are. My levain, may you live a long
and sour life.

## Declaration

The work reported in this thesis is primarily my own unless otherwise
explicitly stated, and has not been submitted for any other degree or
qualification. This thesis unifies and expands work reported in the following
publications:

  * [1]

maxnames[1]

  * [2]

maxnames[2]

During the work on this thesis, I also co-authored the following publications:

  * [3]

maxnames[3]

  * [4]

maxnames[4]

  * [5]

maxnames[5]

##  Chapter 1 Introduction

Breakthroughs in artificial intelligence keep making our machines smarter:
they become better at diagnosing diseases, driving vehicles, predicting the
weather, understanding natural languages, and more. Part of this success is
due to machine learning models like artificial neural networks that have
increasingly more neural connections (aka parameters). Networks went from
millions of parameters (e.g. ResNet) to billions of parameters (e.g. Microsoft
Turing), and are anticipated to reach trillions of parameters. Such models
require massive computing power, which is a major issue [6] since physical
resources are limited (left of Figure 1.1). In the cloud, huge data centres
are taking an increasing toll on our planet: a recent study estimates that
training a single natural language processing neural network has a carbon
footprint equivalent to driving 4,311 km in a European car and would take 823
tree-months to offset [7]. At the edge, running advanced neural networks on
small devices with low energy supply is not yet possible. To sustain progress,
we need to use resources more efficiently.

Figure 1.1: (Left-to-Right) Physical resources are limited at all scales: from
the impact of the cloud on our planet, to battery limitations at the edge.
Different actors optimise resource usage: while domain scientists design new
algorithms, hardware architects design new hardware, and performance engineers
design new program optimisations.

In high performance domains like image processing, physics simulation or
machine learning, resources are used more efficiently through continuous
algorithm, hardware and program optimisation (right of Figure 1.1). Here
domain scientists are responsible for optimising algorithms, hardware
architects optimise hardware, and _performance engineers_ are programmers
responsible for optimising programs. The three optimisation aspects are inter-
dependent: in particular, for each new algorithm and hardware architecture,
new specialised program optimisations are required.

Optimising programs is challenging. On one side, skilled performance engineers
perform optimisations manually on low-level code (e.g. C, OpenCL), which takes
months and risks introducing bugs, slowing down algorithm and hardware
innovation. On the other side, modern compilers targeting heterogeneous and
parallel hardware do not reliably automate program optimisation, facing the
two following major challenges.

_Challenge 1) Extending compilers with new program optimisations is hard._

Traditional _domain-agnostic_ compilers such as LLVM [8] or GCC only support a
fixed, generic set of program abstractions and optimisations. As a result,
they do not automate optimisations for specific domains and performance
engineers must perform them manually (Section 2.1).

Established _domain-specific_ compilers such as Halide [9] for image
processing or TVM [10] for machine learning only support a fixed, specialised
set of program abstractions and optimisations. Although they successfully
automate many optimisations, extending them with new optimisations as the
domain evolves is difficult [11]. Domain evolution includes changes in
algorithms, changes in program optimisations, but also changes in specialised
hardware (e.g. machine learning accelerators [12]). Performance engineers are
not necessarily compiler writers, and falling back to manual optimisation is
often easier for them than modifying a compiler to achieve their goals. In
consequence, many state-of-the-art high-performance kernel libraries (e.g.,
BLAS, cuDNN, MKL) are still predominantly written by hand [13].

Recently, _domain-extensible_ compilers such as Delite [14], Lift [15, 16], or
AnyDSL [17] provide an extensible set of program abstractions and
optimisations, showing potential to mitigate _Challenge 1_ (Section 2.2).
However, domain-extensible compilers remain relatively immature and lack
adoption compared to established domain-specific compilers. For example,
previous work does not investigate whether Lift is able to generate code with
competitive runtime performance compared to established domain-specific
compilers.

_Challenge 2) Manual optimisation is tedious, automatic optimisation is
unsatisfactory._

Compilers such as Lift [15] or PolyMage [18] aim for _full automation_ of
optimisations, removing performance engineers from the optimisation process.
This is highly desirable in scenarios where performance engineers are not
available. However, full automation of optimisations is not always feasible or
even desirable as it may result in poor performance or may be too time-
consuming [19, 20, 21]. When compiler optimisations are unsatisfactory,
performance engineers often fall back to manual optimisation in order to
achieve their performance goals [22, 23, 24].

Compilers such as Halide offer _control_ of optimisations through schedules
[9, 10, 25]. However, schedules are challenging to write [26]. Successful
research goes into automatically generating schedules [27, 28, 29], but
sacrifices control for full automation of optimisations.

This thesis explores _controllable automation_ (Section 2.3.3) of
optimisations that embraces trade-offs between full automation and precise
control of optimisations. Controllable automation enables performance
engineers to gradually take control, instead of abruptly falling back to
manual optimisation when compiler optimisations are unsatisfactory, and has
the potential to mitigate _Challenge 2_ (Section 2.3).

##### Thesis Statement

This thesis shows the potential of a novel compiler design to achieve domain-
extensibility, controllable automation, and generate high performance code.
Domain-extensibility is combined with controlled rewriting in an image
processing case study to generate faster code than the established domain-
specific compiler Halide. Controllable automation is exploited in a linear
algebra case study to automatically explore an optimisation space while
providing the performance engineer with control over the optimisation outcome.

This thesis makes the following research contributions:

  1. 1.

Enhancing Code Generation in a Domain-Extensible Compiler (Chapter 3).
Three important code generation features are added to a domain-extensible
compiler with controllable automation of optimisations called Shine, making it
possible to:

     * •

generate correct and efficient synchronisation barriers

     * •

generate multiple computation kernels and the host code launching them

     * •

generate folded storage for temporary arrays

Adding these features requires making novel design choices in terms of
compiler extensibility and controllability. Crucially the following chapters
require these features to generate imperative code with competitive runtime
performance compared to established domain-specific compilers such as Halide.

  2. 2.

Going Beyond Halide Scheduling with Controlled Rewriting [1] (Chapter 4).
Domain-extensibility is combined with controlled rewriting to optimise a
standard image processing pipeline: the Harris corner detection [30]. We show
how rewriting is controlled in Shine to achieve 6 well-known image processing
pipeline optimisations, including 2 optimisations that are not supported by
Halide schedules. Our results on four mobile ARM multi-core CPUs, often used
for image processing tasks, show that the code generated using Shine for the
Harris corner detection is up to 1.4$\times$ (geomean of 1.27$\times$) faster
than Halide. However, we also observe that controlling rewriting is tedious,
motivating the following chapter that aims to lower performance engineer
effort with semi-automation.

  3. 3.

Proposing a Novel Semi-Automatic Optimisation Technique [2] (Chapter 5).
A new semi-automatic optimisation technique called _sketch-guided equality
saturation_ is developed, based on a fully automated technique called equality
saturation. Sketch-guiding allows performance engineers to guide program
rewriting by specifying rewrite goals as _sketches_ : program patterns that
leave details unspecified. We evaluate sketch-guided equality saturation by
applying 7 realistic optimisations of matrix multiplication. Unguided equality
saturation alone does not scale to the 5 most complex optimisations, even
given an hour and 60GB of RAM. With the guidance of at most 3 sketch guides,
each 10 times smaller than the complete program, the compiler applies the
optimisations in seconds using less than 1GB of RAM. We also explore how to
efficiently encode a polymorphically typed lambda calculus for equality
saturation. The runtime and memory consumption of unguided equality saturation
over lambda terms is reduced by orders of magnitude, which is also beneficial
for sketch-guided equality saturation.

##  Chapter 2 Motivating Background

This chapter motivates the contributions of this thesis by presenting
technical background and related work. The remaining background will be
introduced separately in the technical chapters. Similarly, more specific
related work will be presented throughout the thesis and discussed in Chapter
6.

Section 2.1 provides background on how performance engineers traditionally
program and optimise performance-demanding applications in low-level languages
like C and OpenCL.Performance engineers typically achieve orders of magnitude
performance improvements by manually exploring potential optimisations using
their expert knowledge. Unfortunately this optimisation process is time
consuming, risks introducing bugs, and needs to be repeated for each new
algorithm and hardware target.

To improve over this traditional optimisation process, we would like to
combine higher-level programming models with compilers that automate program
optimisations, but this is a challenging task.

Section 2.2 motivates research on _domain-extensible_ compilers. The compiler
_extensibility challenge_ is introduced: extending compilers with new
optimisations is hard. Existing compilers are categorised into traditional
domain-agnostic compilers; established domain-specific compilers; and emerging
domain-extensible compilers. The advantages and disadvantages of each category
of compilers are highlighted.

Section 2.3 motivates research on _controllable automation_ of optimisations.
The compiler _controllability challenge_ is introduced: even though manual
optimisation is tedious for performance engineers, automatic optimisation is
not always satisfactory. The advantages and disadvantages of automated
optimisation and controlled optimisation are discussed before introducing the
concept of controllable automation of optimisations.

###  2.1 High Performance Programming

####  2.1.1 Hardware Architectures are Evolving

In the quest to maximise runtime performance and energy efficiency, hardware
architects continuously come up with significant architectural changes. First,
hardware architectures became highly parallel [31, 32]. Nowadays, hardware
architectures are increasingly heterogeneous [33, 34]. Such hardware evolution
seriously impacts high performance programming practices. Performance
engineers must consider hardware parallelism and heterogeneity to achieve high
performance when optimising software [35].

In particular, various forms of parallelism are exploited:

  * •

_instruction-level parallelism_ , where multiple instructions are executed
simultaneously. The hardware commonly exploits this parallelism implicitly at
runtime (e.g. instruction pipelining, superscalar execution), but some
hardware explicitly exposes this parallelism to software instead (e.g. Very
Long Instruction Word).

  * •

_vector-level parallelism_ , where multiple vector elements are processed
simultaneously. This is also referred to as Single Instruction, Multiple Data
(SIMD). Many processors explicitly expose this parallelism. As a result,
software must use explicit instructions to exploit this parallelism.

  * •

_thread-level parallelism_ , where multiple software threads are executed
simultaneously. A thread is a programmed sequence of instructions that can be
scheduled independently to other threads on the available hardware cores. Care
must be taken when developing explicitly multi-threaded software, as it is a
difficult and error-prone task [36].

Two common hardware architectures are often combined in heterogeneous systems:
multi-core Central Processing Units (CPUs) and many-core Graphics Processing
Units (GPUs).

##### Multi-Core CPUs

Optimised for individual thread performance and low latency, multi-core CPUs
have a few sophisticated cores able to run relatively heavyweight threads.
CPUs dynamically adapt to various usage patterns through complex execution,
control flow and memory management. To exploit locality of memory accesses,
CPUs have a hierarchy of caches that automatically stores recently accessed
data closer to the cores. Three levels of caches named L1, L2 and L3 are
commonly used; where L1 is the closest to the cores. The closer to a core a
cache is, the smaller and faster it is. Even though data is moved implicitly
between caches, program optimisations can still improve cache usage by
improving the locality of accesses.

While CPUs focus on individual thread performance, they still expose
significant parallelism. It is standard for CPUs to execute 8 to 16 threads
simultaneously and to provide 128 to 256 bits of vector parallelism (i.e. 4 to
8 32-bit values) on top of instruction-level parallelism.

##### Many-Core GPUs

Optimised for high throughput, many-core GPUs hide high latency operations by
overlapping the execution of many relatively lightweight threads. Originally
specifically designed to accelerate 3D graphics rendering, GPUs have been
generalised to support more use cases and are becoming increasingly flexible
[37]. GPUs are now popular for massively parallel and regular computations in
many performance demanding domains.

It is standard for GPUs to contain hundreds or thousands of cores, grouped
hierarchically to share hardware resources during execution. For example,
while GPUs have a cache hierarchy, they also typically provide scratchpad
memory which is shared by a group of threads. A scratchpad is similar to a
cache, but is less sophisticated because data transfers are explicitly
programmed instead of implicit. In general, GPUs tend to avoid sophisticated
hardware techniques with high hardware cost and instead favor simpler
techniques with low hardware cost. This shift from implicit hardware logic to
explicit software programming makes program optimisation both more important
and difficult.

##### Other Processors

While this thesis focuses on CPUs and GPUs, heterogeneous systems also include
reconfigurable hardware such as Field-Programmable Gate Arrays (FPGAs) or
Coarse-Grained Reconfigurable Architectures (CGRAs); as well as highly
specialised hardware such as Digital Signal Processors (DSPs) or Tensor
Processing Units (TPUs). The contributions of this thesis may be combined with
related work that tackles such hardware [38, 39, 40]. For example, the D2A
methodology [41] addresses compilation challenges that are specific to
accelerators, such as the mismatch between fine-grained compiler IR
(Intermediate Representation) operations and coarse-grained accelerator
operations, and the need for a formal hardware specification. In the future,
processors are likely to become even more diverse, complicating program
optimisation tasks: performance engineers may need to exploit quantum [42],
neuromorphic [43] or optical processors [44].

####  2.1.2 The OpenCL Programming Model is too Low-Level

To achieve high performance on parallel and heterogeneous hardware, the
industry standard is to use relatively low-level programming models. In this
thesis, _low-level code_ is code that is specialised to a given hardware for
performance reasons. By extension, a _low-level programming model_ is a
programming model where non-specialised code performs poorly, i.e. the
programming model fails to provide what is also called performance portability
[45]. We specifically provide more background on the OpenCL programming model
[46, 47], an open standard targeting diverse processors that we will use in
this thesis.

The OpenCL standard consists of a C API to orchestrate computation from the
_host_ and a programming language to express computation on _OpenCL devices_.

##### Compute Model

Programs executed on a device are called _kernels_. A kernel is a program that
each thread executes in a _Single Program, Multiple Data_ (SPMD) fashion. The
host program may submit a kernel for execution on a device by specifying an
N-dimensional index space (where N goes from 1 to 3). For each point in the
index space, an instance of the kernel program called _work-item_ is executed.
Each point in the index space is a global identifier for the corresponding
work-item and can be queried by the kernel program to inform its execution.

The index space is decomposed into _work-groups_ of multiple work-items that
may share hardware resources for cooperation. Each work-group has a
corresponding identifier, and each work-item has a local identifier within its
work-group. This thread hierarchy closely corresponds to the compute hierarchy
found in most GPU architectures, allowing efficient runtime scheduling, and is
compatible with execution on other architectures such as CPUs.

##### Memory Model

The OpenCL memory model also closely corresponds to the memory hierarchy of
common GPU architectures, exposing different address spaces:

  * •

_Global memory_ is accessible to all work-items in all work-groups, and also
to the host program. On CPUs and GPUs, it typically corresponds to device RAM.

  * •

_Local memory_ is only accessible to work-items of the same work-group. It
usually corresponds to RAM on CPUs, and faster scratchpad memory on GPUs.

  * •

_Private memory_ is only accessible to a single work-item. On CPUs and GPUs,
it typically corresponds to registers, the fastest possible memory.

OpenCL also provides _constant memory_ , read-only global memory that remains
constant during kernel execution, however we will not exploit this type of
memory in this thesis.

##### Kernel Programming

The OpenCL C programming language, used to define kernels, is derived from the
C99 specification. On one hand, OpenCL C disallows pointers to functions,
recursion, and dynamic memory allocation. On the other hand, OpenCL C
introduces address space qualifiers, custom data types and built-in functions.

LABEL:fig:opencl-add shows a simple implementation of point-wise addition of
two arrays using an OpenCL kernel. In line 1, the kernel signature consists of
pointers referring to arrays in global memory (a, b, output) and an integer
representing the size of these arrays (n). The kernel program is executed by
multiple threads, and starts by querying the global identifier of the current
thread in line 2 (get_global_id(0)). Depending on the number of available
threads (get_global_size(0)), the program may loop depending on how many
elements the current thread needs to compute. An element is computed inside
the loop using indexing in line 3.

⬇

1kernel void add_v1(global float* a, global float* b, global float* output,
int n) {

2 for (int gid = get_global_id(0); gid < n; gid += get_global_size(0)) {

3 output[gid] = a[gid] + b[gid];

4 }

5}

Listing 1: Point-wise addition of two arrays in OpenCL

In practice, optimised kernel implementations are significantly more complex.
LABEL:fig:opencl-add-opt shows another implementation of the same computation.
The program leverages vector-level parallelism using vload4 and vstore4 to
transfer vectors of 4 single-precision floating-point values, as well as a
vectorised addition. The program loops are tiled to increase single-thread
workload and potentially improve data locality: chunks of 4 vectors are
computed together.

⬇

1kernel void add_v2(global float* a, global float* b, global float* output,
int n) {

2 for (int gid = get_global_id(0); gid < n/16; gid += get_global_size(0)) {

3 for (int i = 0; i < 4; i += 1) {

4 float4 av = vload4(gid * 4 + i, a);

5 float4 bv = vload4(gid * 4 + i, b);

6 vstore4(av + bv, gid * 4 + i, output);

7 }

8 }

9}

Listing 2: Point-wise addition of two arrays in OpenCL: tiled, vectorized.

Optimised code is typically orders of magnitude faster than naive code, but is
also more complex and specialised to a given hardware. For performance
engineers, manually optimising low-level code using techniques such as tiling
or vectorisation is time-consuming and risks introducing bugs. This is
aggravated by the fact that many optimisation decisions require global
thinking: naive compositions of locally optimised code often perform poorly.

##### Other Industrial Programming Standards

The _C_ and _C++_ languages are commonly used to program CPUs, in combination
with primitives such as _POSIX threads_ for multi-threading and _SIMD
instructions_ (e.g. SSE, AVX, NEON) for vectorisation. _CUDA_ is a popular
closed framework developed by NVIDIA for NVIDIA GPUs. It is also common to
build convenient abstractions on top of these standards. _SYCL_ allows single-
source programming in C++ for OpenCL, avoiding the need to write separate
OpenCL C code, and providing implicit memory transfers. _OpenMP_ and _OpenACC_
are used for quick code acceleration with compiler directives, such as
parallelising a for loop by preceding it with #pragma omp parallel for. All of
these programming models suffer from similar performance engineering
productivity and portability issues as OpenCL: non-specialised code performs
poorly compared to specialised code.

####  2.1.3 Higher-Level Programming Models are Challenging to Compile

Higher-level programming models and languages are a promising way to increase
productivity, performance and portability. We give concrete examples with
domain-specific languages and array programming languages before mentioning
common compilation challenges.

##### Domain-Specific Languages

Domain-Specific Languages (DSLs) enable combining convenient, hardware-
agnostic programming with high-performance for specific domains. This is a
proven success for domains such as image processing [9, 48, 18], signal
processing [49], linear algebra [50], machine learning [10] or partial
differential equations [51, 52].

For example, the point-wise addition of two arrays can be defined in the
Halide image processing DSL by writing a high-level _algorithm_ as in
LABEL:halide-add-algorithm. Crucially, DSLs separate concerns: domain
scientists can focus on writing high-level programs (what to compute) while
compilers and performance engineers optimise programs for performance (how to
compute). For example, LABEL:halide-add-algorithm could be compiled by Halide
into an OpenCL implementation (e.g. LABEL:fig:opencl-add, LABEL:fig:opencl-
add-opt) or into some other optimised implementation, depending on the target
hardware.

⬇

Var i; Func output;

output(i) = a(i) + b(i);

Listing 3: Point-wise addition of two arrays in Halide

##### Array Programming Languages

Array programming languages such as Accelerate [53], SaC [54], NOVA [55], Lift
[15], Futhark [56], or Dex [57] provide a middle-ground between highly
specialised languages and general-purpose languages. Multi-dimensional arrays
are a key abstraction for many performance-demanding application domains (e.g.
physics, chemistry, mechanics, image processing, data processing, linear
algebra, machine learning), as they can be used to represent widely used
tensors [58, 59].

Some array programming languages focus on index-based notations [54, 57],
similar to how the i index is used in the Halide algorithm of LABEL:halide-
add-algorithm. Other array programming languages focus on collective array
operations such as map or reduce [53, 55, 15], that can also be seen as
algorithmic skeletons [60]. For example, the point-wise addition of two arrays
can be defined in the Lift language as in LABEL:fig:lift-add. zip a b combines
two arrays a and b whose elements are added pairwise using map, and no
explicit indexing is required.

⬇

map ($\lambda$x. (fst x) + (snd x)) (zip a b)

Listing 4: Point-wise addition of two arrays in Lift

Array programming languages offer the potential to share common infrastructure
across multiple application domains while still being convenient and hardware-
agnostic.

##### Developing Compilers for High-Level Programs

Whether high-level programs are written in domain-specific or array languages,
building and maintaining compilers that generate high-performance code for a
variety of algorithmic domains and hardware architectures is challenging.

Industrial-strength compilers such as Halide for image processing and TVM for
machine learning are the result of years of engineering effort from multiple
contributors. The Halide repository contains 278k lines of code written by 167
contributors.111https://github.com/halide/Halide/commit/820ec1f963f06f53a1808eb9b4631f2031be7468
The TVM repository contains 794k lines of code written by 752
contributors.222https://github.com/apache/tvm/commit/178f82dc481bf31961206412c22dd5519a245b49
Further, both use other compiler technologies such as LLVM as backends. LLVM
8.0.1 comprises 6,887k lines of code written by 1,210 contributors, with an
estimated development cost of 529M$ [61].

The following sections motivate the need to address two specific compiler
challenges: the _extensibility challenge_ and the _controllability challenge_.

###  2.2 The Compiler Extensibility Challenge

####  2.2.1 Domain-Agnostic Compilers

Traditional, general-purpose compilers such as LLVM [8] or GCC are _domain-
agnostic_ : they only support a fixed, generic set of program abstractions and
optimisations.

LLVM is a widely adopted compiler framework consisting of production-quality
libraries for modular compiler construction. The LLVM IR (Intermediate
Representation) provides many generic abstractions that are useful across
multiple programming languages and hardware targets, such as types, functions
and exceptions. LLVM also provides many powerful optimisation passes, such as
dead code elimination, common subexpression elimination, function inlining,
strength reduction or target-specific instruction selection. LLVM is
successfully used as part of compilers for the C/C++ (clang compiler), Rust
and Swift languages, to name a few.

While the LLVM framework is modular and extensible with custom compiler
analyses and compiler passes, the LLVM IR is generic and even lower-level than
languages like C or OpenCL. Higher-level program constructs must be lowered by
compilers targeting LLVM, erasing domain-specific information which can be
hard or impossible to recover. As a result, domain-agnostic compilers
typically do not automate optimisations for specific domains. Instead,
performance engineers must perform them manually as seen in Section 2.1.2.

##### Optimised Libraries

Numerous _libraries of optimised functions_ (e.g. BLAS, cuDNN, MKL, NVIDIA
Performance Primitives, ARM Compute Library, Intel IPP, OpenCV) provide
optimised implementations of common computations. However, they are developed
by performance engineers at high cost. Limited engineering budget leads to
limited functionality and limited support for certain use cases or hardware.
Moreover, the composition of optimised functions through library calls is
often far from optimal because many optimisations cannot be applied across
library calls. For example, Weld [62] achieves order of magnitude performance
improvements by optimising across libraries and functions using a common
intermediate representation.

####  2.2.2 Domain-Specific Compilers

Established _domain-specific_ compilers only support a fixed, specialised set
of program abstractions and optimisations. Examples include Halide [9] and
PolyMage [18] for image processing; TVM [10] for machine learning; SPIRAL for
signal processing [49]; Firedrake for partial differential equations [52].

Domain-specific compilers are successfully established in industry where they
automate many domain-specific optimisations and achieve impressive performance
results. Halide is used for some processing tasks in the Google Pixel camera,
Adobe Photoshop, and YouTube. TVM is used and developed by companies like AMD,
ARM, Microsoft, NVIDIA, and Samsung.

However, extending domain-specific compilers with new optimisations as the
domain evolves is difficult [11]. Domain evolution includes changes in
algorithms, changes in program optimisations, but also changes in specialised
hardware.

Compiler extensions are particularly difficult when they have an impact on the
entire compilation stack, requiring simultaneous expertise of high-level
algorithms, multiple compilation aspects, and low-level hardware targets.
Separating concerns facilitates extension, for example by separating the
definition of correct transformations from the search for best transformations
as shown later in this thesis.

##### Extending Halide is Hard

The Halide Development Roadmap from
2020333https://github.com/halide/Halide/issues/5055 highlights that extending
Halide is hard, listing unsolved questions such as:

  * •

"How do we make Halide easier to use for researchers wanting to […] extend it
[…]?"

  * •

"How do we make Halide more useful on current and upcoming hardware?"

  * •

"How do we make Halide more useful for new types of application?"

In fact, solutions to the extensibility challenge have been independently
researched in the course of this thesis by Halide authors [63, 13], showing
that our research direction is valuable and our motivations shared by the
broader community.

##### Extending TVM is Hard

The TVM Unity vision for 2022444https://tvm.apache.org/2021/12/15/tvm-unity
demonstrates that, currently, extending TVM is hard. Two problematic
boundaries are identified in the compilation stack. First, vertical boundaries
between successive abstraction layers. Each local decision made within a layer
has a global impact on the following layers. As a result, introducing support
for new hardware in the last layer requires re-thinking how decisions are made
in all prior layers, requiring global instead of local engineering efforts.
Second, horizontal boundaries between lowering strategies. For example, TVM
may decide to target optimised libraries or to generate code from scratch, but
typically cannot combine both approaches.

> “[…] boundaries are slowing down the pace of innovation in machine learning.
> New hardware accelerators are emerging with new levels of capability and
> performance, but harnessing them will require fluid collaboration between ML
> scientists, ML engineers, hardware vendors that these boundaries prevent. To
> cope with the rapid pace of change in ML systems, frameworks need to support
> incremental evolution: Incorporating new capabilities should require effort
> proportional to the change, not wholesale re-engineering at each level. ”
> (TVM Unity vision)

####  2.2.3 Domain-Extensible Compilers

Building DSLs and their compilers has long been recognized as a highly complex
task, which is why many projects aim to simplify the development of DSLs. For
example, Spoofax [64] is a language workbench simplyfing the development of
DSLs. The idea is to produce parsers, type checkers, compilers, interpreters
and even IDE support from declarative language specifications. Gammars are
used to declare syntax, rewrite rules to declare semantics. With Spoofax,
multiple domain-specific compilers are still constructed even if the task is
simplified. The constructed compilers may still be hard to adapt as their
domain evolves, as seen for established domain-specific compilers in Section
2.2.2. Another approach is to build _domain-extensible_ compilers providing an
extensible set of program abstractions and optimisations.

Delite [14] is a DSL framework providing a fixed set of generic parallel
patterns that DSLs can target. Domain-specific optimisations can be defined
using staging, a technique used to eliminate the cost of DSL abstractions
while lowering them to generic parallel patterns. Delite demonstrated that
multiple DSLs can be composed together while sharing common infrastructure
[65]. The Delite projet also identified that its staging mechanism is not
sufficient, and showed that using term rewriting techniques is beneficial for
certain use cases [66].

AnyDSL [17] is a more recent approach leveraging partial evaluation to replace
DSL compilers with self-specialising library code written in a language called
Impala. The idea of implementing languages as libraries has a long strand of
research predating AnyDSL, such as the Racket language that leverages macros
instead of partial evaluation [67].

Lift [15] combines a high-level functional language with an extensible rewrite
system to define an optimisation space and generate high-performance code.
Lift has demonstrated its ability to generate high-performance code by
deriving multiple low-level implementations from a single portable high-level
program, obtaining performance on par with highly tuned platform-specific
libraries on various GPUs [68]. Lift is extensible with rewrite rules and
language primitives by design, and has been successfully extended for stencil
computations [16] which are common in multiple domains including image
processing.

Overall, domain-extensible compilers show a lot of potential, but they remain
relatively immature and lack adoption compared to established domain-specific
compilers. This thesis builds on the foundations laid by Lift and addresses
some of its shortcomings. For example, previous work does not investigate
whether Lift is able to generate code with competitive runtime performance
compared to Halide for image processing as we do in Chapter 4. Lift also does
not address the compiler controllability challenge that we discuss in the next
section.

While we follow a rewrite-based approach, multiple extensible techniques can
complement each other, as demonstrated by Delite for staging and rewriting. We
believe that research on each technique is useful on its own.

##### Extending Compilers with Rewrite Rules

Many compilers before Lift already allowed programmers to express domain-
specific optimisations as rewrite rules. The Glasgow Haskell Compiler (GHC)
allows this, but only applies rewrite rules according to a simple strategy,
assuming that the right-hand side of a rule is always preferable to the left
[69]. The following pragma would instruct GHC to always rewrite map f (map g
xs) into map (f . g) xs:

⬇

{-# RULES

"map/map" forall f g xs.

map f (map g xs) = map (f . g) xs

#-}

While such a simple strategy can be effective in optimising programs, it falls
short for use cases where deciding which rewrite rule is beneficial when is
hard.

Stratego [70] is a language for defining customised rewriting strategies,
which is included in the Spoofax language workbench. Phobos is a compiler
frontend enabling the development of domain-specific languages by combining an
open term language with term rewriting [71].

On one hand, tools like GHC, Stratego and Phobos allow domain-extensbility via
rewrite rules but do not focus on high-performance code generation. On the
other hand, tools like Spiral [72, 49] successfully generate high-performance
code for heterogeneous hardware using rewrite rules, but are domain-specific
and not extensible. This thesis follows Lift’s ambition of combining domain-
extensibility with high-performance code generation for heterogeneous hardware
using rewrite rules.

###  2.3 The Compiler Controllability Challenge

####  2.3.1 Automatic Optimisation

Many compilers aim for fully automated optimisation, some of them lean towards
_greedy optimisation_ , while others lean towards _explorative optimisation_.

##### Greedy Optimisation

Greedy optimisation is often achieved using a fixed pipeline of optimisation
_passes_ , relying on expert heuristics to make local decisions that are
assumed to be beneficial for final performance.

For example, the LLVM framework [8] provides many transformation and
optimisation passes, and allows defining custom passes. In GHC, user-defined
rewrite rules are applied greedily [69]. Many other projects rely on greedy
optimisation passes, such as Accelerate [73] and the Shir compiler that adopts
ideas from Lift to target FPGAs [40].

Typically, greedy optimisation is fast, but may result in poor performance as
it gets stuck in local minima and may poorly predict performance benefits. The
problem of deciding when to apply which optimisation pass is known as the
_phase ordering problem_ , and has a huge impact on the performance of the
rewritten program.

##### Explorative Optimisation

Explorative optimisation is a more holistic approach, where a space of
possible implementations arising from different combinations of optimisations
is explored, in order to make a decision.

Random sampling, or Monte Carlo methods, can be used to navigate the space of
possible implementations. Using this approach, finding an implementation with
high enough performance can be very time consuming. Lift [15] uses random
sampling, and optimising a single convolution takes hours to reach peak
performance on a GPU [16, 74].

Performance models can be used to choose from the space of possible
implementations. Some performance models are analytical, describing
performance as an equation, as done in TyTra [38, 75], SPIRAL [72, 49] and
Telamon [76]. Analytical performance models are typically fast to evaluate,
estimating a program’s performance without executing it. Creating an accurate
performance model for complex systems is a challenging task, which is why
empirical performance models are also used, predicting future performance by
learning from previous performance data [77]. Although performance modelling
is making great progress, it remains a standing challenge that is the focus of
multiple papers each year [78, 79, 80, 81, 82].

Automatic, explorative optimisation is an active research field and novel
techniques are constantly developed to explore various optimisation spaces,
such as iterative compilation [83], adaptive compilers [84], or equality
saturation [85, 86].

##### Advantages of Automatic Optimisation

Fully automated optimisation is invaluable if performance engineers are not
available, or are allocated to other tasks. State-of-the-art automation of
tedious, low-level implementation decisions such as register allocation or
instruction selection delivers satisfying performance for most applications,
and frees up development time for higher-level optimisations where a
performance engineer can have higher impact. Valuable human time is saved,
trading it for machine time. Both inside and outside of computer science,
automation has the potential to increase task feasibility, productivity and
quality [87].

##### Disadvantages of Automatic Optimisation

Fully automated optimisation may result in poor performance or require too
much compilation time [19, 20, 21]. Crucially, control is usually sacrificed.
The compilation flags offered by compilers like clang, that uses LLVM for
C-like languages, only provide very limited control. Therefore, when compiler
optimisations are unsatisfactory, performance engineers often fall back to
manual optimisation in order to achieve their performance goals [22, 23, 24].

####  2.3.2 Controlled Optimisation

Instead of manually optimising low-level code, or entrusting optimisation to a
black box automatic compiler, performance engineers may control semi-automatic
compilers that effectively act as optimisation assistants. Techniques for
controlled optimisation include transformation scripts, scheduling APIs, and
rewriting strategies.

##### Transformation scripts

The polyhedral framework [88] is a powerful technique modelling loop nests as
polyhedra, enabling sophisticated dependency analysis, and providing tools to
explore the space of valid transformations. URUK _transformation scripts_ were
proposed to define and apply sequences of polyhedral transformations [89, 90].
The CHiLL framework extended this idea to more complex loop transformations
[91]. An example of CHiLL script that applies loop transformations such as
tile(statement_id, loop_level, tile_size) is given in LABEL:fig:chill-trans-
script. Loo.py [92] is a more recent programming system embedded in Python and
inspired by the polyhedral model, that provides a library of loop
transformations.

⬇

permute([3,1,2])

tile(0,2,TJ)

tile(0,2,TI)

tile(0,5,TK)

datacopy(0,3,2,[1])

datacopy(0,4,3)

unroll(0,4,UI)

unroll(0,5,UJ)

Listing 5: CHiLL transformation script used to optimise matrix multiplication
in [91].

##### Scheduling APIs

Halide [9] has popularised the development of compilers that offer scheduling
APIs. In this setting, a _schedule_ defines how to optimise an _algorithm_
that defines what to compute at a high level of abstraction, separating the
two concerns. The Halide scheduling API focuses on exposing important trade-
offs between parallelism, redundancy, and locality; by defining when and where
functions are computed and stored. LABEL:halide-add-schedule gives a simple
example of Halide schedule. Following the success of Halide, the TVM [10],
Fireiron [25] and Tiramisu [93] compilers apply the same principle to
different domains.

⬇

output.split(i, gid, i, 16)

.parallel(gid)

.vectorize(i, 4);

Listing 6: Schedule for the point-wise addition Halide algorithm from
LABEL:halide-add-algorithm, loosely comparable to the OpenCL implementation
from LABEL:fig:opencl-add-opt.

##### Rewriting Strategies

As discussed in the previous section, languages such as Stratego [70] enable
the definition of custom _rewriting strategies_ [94]. As this thesis follows a
rewrite-based approach to tackle the extensibility challenge, rewriting
strategies are of particular interest to tackle the controllability challenge.
Chapter 4 will demonstrate how rewriting strategies, written in a language
called Elevate [3], can be used to apply transformations that are beyond what
is possible with Halide schedules on a case study.

##### Advantages of Controlled Optimisation

Transformation scripts, schedules and rewriting strategies all enable
performance engineers to take control over the optimisation process. They
expose a structured optimisation space to the performance engineer, who
explores it in order to achieve his performance goals. Compared to manual
optimisation of low-level code, controlled optimisation saves development
time, and given a correct compiler implementation, avoids introducing bugs
during optimisation.

##### Disadvantages of Controlled Optimisation

Controlled optimisation remains challenging, as argued for schedules in [26].
Further, systems built around these techniques typically do not allow smooth
trade-offs between automation and control, as this trade-off is built-in.
Therefore, controlled optimisation may be tedious and take too much
performance engineering effort compared to automatic optimisation.

####  2.3.3 Controllable Automation of Optimisations

This thesis explores _controllable automation_ of optimisations that promotes
trade-offs between automation and control of optimisations. If a compiler
supports both extremes, as well as a spectrum in-between, performance
engineers should be able to gradually take control of optimisations depending
on their performance requirements and time budget.

Properties like domain-extensibility and controllable automation are hard to
quantify on a scale, as they represent informal design principles. However, it
is possible to identify whether different approaches provide controllable
automation of specific choices or not, as in Table 2.1.

Approach |  Parameter Choice (tile size) |  Optimisation Choice (tiling)
---|---|---
Halide Schedule | ✗ | ✗
TVM Schedule Template | ✓ | ✗
This Thesis | ✓ | ✓
Table 2.1: This thesis allows controllable automation of optimisation choices
like tiling, as opposed to current Halide or TVM approaches.

Halide schedules neither provide controllable automation of parameters (e.g.
tile size), nor of optimisations (e.g. tiling). Successful research goes into
_autoscheduling_ , the challenge of automatically generating schedules, both
for Halide [27, 28, 29] and TVM [95, 96]. However to the best of our
knowledge, autoscheduling techniques offer no controllable automation:
performance engineers cannot constrain autoscheduling to apply a specific
tiling optimisation, or not use specific tile sizes. This contrasts with the
conclusion of the first Halide paper: _"the ultimate solution must allow a
smooth trade off between inference when it is sufficient, and sparse
programmer control when it is necessary"_ [9].

TVM _schedule templates_ provide controllable automation of parameters, but
not of optimisations. For example, performance engineers may delegate the
tuning of some numerical parameters, such as tile sizes, while manually
specifying other parameters in a schedule template [10]. However, this
approach does not allow delegating optimisations that significantly change
code structure, such as deciding whether to use tiling, or how to fuse
operators [95].

By contrast, the techniques presented in this thesis enable controllable
automation of both parameters and optimisations (Chapter 5).

In the polyhedral community, there exists multiple trade-offs between
automation and control. Transformation scripts can be used to structure
automatic empirical optimisation [97, 98]. Powerful heuristics have also been
developed to automate optimisation, without using transformation scripts, in
PLuTo [99], PPCG [100] and PolyMage [18, 101]. Ways to increase interactivity
have been explored, _"enabling users to examine, refine, replay and even
design complex optimizations semi-automatically in partnership with the
compiler"_ [102]. Program optimisation tactics, loosely comparable to
rewriting strategies, have been developed independently during this thesis for
the polyhedral model [103, 104]. While the polyhedral framework is amenable to
controllable automation, we focus on term rewriting techniques instead in this
thesis. Term rewriting is more flexible than the polyhedral framework [105],
which gets its strength from a restricted but carefully structured
transformation space.

In the previous section, we already identified term rewriting as a compelling
technique to tackle the extensibility challenge. Term rewriting is also a
compelling technique to tackle the controllability challenge. Indeed, the
definition of rewrite rules can be clearly separated from their application
during rewriting. Given rewrite rules, we can decide to fully automate their
application as done in Lift with random sampling, to precisely control their
application as done in Stratego with rewriting strategies as we will leverage
in Chapter 4, or trade-off between the two as we will explore in Chapter 5.

###  2.4 Conclusion

This chapter motivated the contributions of this thesis by presenting
technical background and related work.

Section 2.1 discussed how, in a context of evolving hardware architectures,
the programming models commonly used in industry, such as OpenCL, are too low-
level. Higher-level programming models such as domain-specific languages or
array programming languages have been successful at increasing productivity
and portability. However, building and maintaining compilers that generate
high-performance code for a variety of algorithmic domains and hardware
architectures, while minimising engineering efforts, is challenging.

Section 2.2 discussed the compiler _extensibility challenge_. Domain-agnostic
compilers such as LLVM are typically not able to automate important domain-
specific optimisations. Extending domain-specific compilers such as Halide or
TVM with new optimisations is hard. Domain-extensible compilers show potential
to mitigate the extensibility challenge, but remain relatively immature
compared to domain-specific compilers. Research is needed to further develop
domain-extensible compiler technology.

Finally, Section 2.3 discussed the compiler _controllability challenge_.
Automatic optimisation is invaluable if performance engineers are not
available, but may result in poor performance or be too time consuming for
some applications. Alternatively, performance engineers may take control of
the optimisation process. However, controlling the optimisation process is
challenging as there is typically no smooth transition between automation and
control of optimisations. Research is needed to enable performance engineers
to gradually take control instead of abruptly falling back to manual
optimisation.

> Imagination is the beginning of creation. You imagine what you desire, you
> will what you imagine and at last you create what you will.

— George Bernard Shaw

##  Chapter 3 Code Generation in a Domain-Extensible Compiler

The previous chapter identified term rewriting as a compelling technique to
tackle the compiler extensibility and controllability challenges. Lift [15]
addresses the _extensibility challenge_ by combining a high-level functional
language with an extensible term rewriting system [68]. However, Lift does not
address the _controllability challenge_ , and optimising a single convolution
takes hours to reach peak performance on a GPU [16, 74]. Shine is a novel
compiler inspired by Lift, with the additional goal to provide controllability
by exploring trade-offs between automation and control of rewrite rule
applications. A key concern for Shine is how to generate imperative code from
rewritten functional programs, with competitive runtime with established
domain-specific compilers such as Halide and TVM. Code generation in Shine
does not aim to be extensible and controllable, but rather to predictably
preserve the optimisation choices encoded during rewriting.

This chapter presents the design and implementation of three important code
generation features. These features are crucial for Shine to generate faster
code than Halide on an image processing case study via controlled rewriting
(Chapter 4), and similarly fast code as TVM on a linear algebra case study via
semi-automated rewriting (Chapter 5).

Section 3.1 introduces the Shine compiler and its Rise language, both the
result of collaboration between multiple researchers [4, 106]. The rest of
this chapter presents work that is solely my own. Section 3.2 contributes a
synchronisation _barrier insertion_ algorithm that does not need to be
modified when extending Rise patterns, contrasting with the barrier
elimination algorithm of Lift [107]. While barrier insertion is implicit and
automated, the next two features add new Rise patterns explicitly exposing
implementation choices to rewriting, allowing design space exploration and
external control. Section 3.3 demonstrates how a Rise pattern is added for
explicit _kernel execution_ , and how Shine is modified to generate
corresponding imperative code. Section 3.4 discusses how Rise patterns are
added for explicit _storage folding_ , and how Shine is modified to generate
corresponding imperative code. Section 3.5 concludes.

###  3.1 The Rise Language & Shine Compiler

Shine is a domain-extensible compiler inspired by the Lift compiler, its
implementation is open-source.111https://github.com/rise-lang/shine Figure 3.1
gives an overview of the Shine compiler, which is meant to be a bridge between
domain-specific languages (left of Figure 3.1) and hardware targets (right of
Figure 3.1). A program written in a domain-specific language is translated
into a _high-level Rise program_ (this step is left for future work and not
the topic of this thesis). Then, the high-level program is rewritten into a
_low-level Rise program_ where implementation decisions are explicitly encoded
(Section 3.1.1). Finally, imperative code is generated and can be executed on
a given hardware target (Section 3.1.3). The focus of this chapter (Chapter 3)
is to enhance imperative code generation from low-level Rise programs (red
rectangle in Figure 3.1).

Figure 3.1: Rise & Shine: A Domain-Extensible Compiler Design

Rise is a functional data-parallel language inspired by the Lift language.
More precisely, Rise is a lambda calculus extended with higher-order functions
called _patterns_ (e.g. map or reduce) and a restricted form of dependent
typing (Section 3.1.2). Functional languages are referentially transparent,
meaning that an expression is equivalent to its value. There are no side
effects, and semantics-preserving rewrite rules can be easily defined to
encode program optimisations.

The rest of this section gives a complete overview of the compilation stack,
to contextualise the following sections.

##### Relationship to Lift

Rise & Shine started as an effort to re-implement Lift in a more principled
way by following the typed lambda calculi formalism [108]. In Shine, patterns
are implemented as higher-order functions, and functions as expressions.
Adding a pattern only requires providing a name and a type. This contrasts
from the Lift implementation, where adding a pattern also requires defining
how to infer its type (Table 3.1). Both implementations then evolved
separately as the original research group forked to explore two different
directions. There are no plans to merge the implementations, although it would
be possible.

language | functions are expressions | function type | type inference
---|---|---|---
Lift | no | no | pattern-specific code
Rise | yes | yes | generic for all functions
Table 3.1: Rise follows the typed lambda calculus formalism more closely than
Lift.

####  3.1.1 Optimising Programs via Rewriting

The Shine compiler takes as input a high-level Rise program describing _what_
to compute, rather than _how_ to compute. For example, a vector dot product is
represented as the high-level Rise program dot:

⬇

def dot a b = zip a b $\vartriangleright$ map ($\lambda$x. (fst x) $\times$
(snd x)) $\vartriangleright$ reduce + 0

Listing 7: dot program in Rise

The zip a b pattern combines two vectors a and b whose elements are multiplied
pairwise using map before they are summed using reduce + 0. The triangle
symbol (x $\vartriangleright$ y) indicates the chaining of operations via
function application (y x) or composition ($\lambda$z. y (x z)).

The high-level dot program does not encode how it is executed. We could
parallelise the map, store the intermediate result, and perform a sequential
reduce. Alternatively, we could avoid intermediate storage by fusing map and
reduce into a single sequential reduction. Many more options are possible, and
such choices are encoded explicitly by applying rewrite rules.

To achieve a fused version avoiding the intermediate results, a compiler
writer or performance engineer writes the reduceMapFusion rewrite rule.
Currently, the rule is written in an informal syntax, future work may design a
formal language for Rise rewrite rules. This rule states that mapping a
function f over an array before reducing the array is equivalent to reducing
the array while applying f on the go. Note that the reduction must be
performed sequentially because the reduction operator is not commutative
anymore.

⬇

rule reduceMapFusion = map f $\vartriangleright$ reduce g init

${\color[rgb]{0.19921875,0.1328125,0.53515625}\definecolor[named]{pgfstrokecolor}{rgb}{0.19921875,0.1328125,0.53515625}\bm{\mapsto}}$
reduceSeq ($\lambda$acc x. g acc (f x)) init

Listing 8: reduceMapFusion rewrite rule.

Deciding which rewrite rule to apply where is challenging, and will be
discussed in Chapters 4 and 5. For now, consider that reduceMapFusion is
applied to map $\times$ $\vartriangleright$ reduce + 0 in the high-level dot
product program, replacing it with the right hand side of the rewrite and
producing the following low-level program:

⬇

def dotSeq a b = zip a b $\vartriangleright$ reduceSeq ($\lambda$acc x. acc +
(fst x) $\times$ (snd x)) 0

Listing 9: dotSeq program in Rise, equivalent to dot from LABEL:dot-rise.

Generating the equivalent C or OpenCL code from this representation is
conceptually straightforward, but has some technical challenges that have been
explored in prior work [107, 106] and are discussed in Section 3.1.3. Some
patterns generate loops (e.g. reduceSeq), others generate scalar expressions
(e.g. +, $\times$, 0), and others affect indexing (e.g. zip, fst, snd). The C
function dotSeqC is generated from dotSeq, implementing the dot product with a
sequential reduction loop as expected:

⬇

1void dotSeqC(float* output, int n, float* a, float* b) {

2 float acc;

3 acc = 0.0f;

4 for (int i = 0; i < n; i++) {

5 acc = acc + (a[i] * b[i]);

6 }

7 output[0] = acc;

8}

Listing 10: dotSeqC prorgam in C, generated from LABEL:dot-seq-rise

####  3.1.2 High-Level and Low-Level Rise Programs

Rise is a functional language with anonymous functions (written as $\lambda$x.
e where x is a variable name and e a Rise expression), familiar function
application (written as e e), identifiers and literals. The language is
embedded in Scala, allowing meta-programming and the definition of macros,
i.e. Scala code that will generate Rise programs. Rise also defines a set of
high-level patterns to describe computations as shown in LABEL:fig:high-level-
patterns together with their types. Formally, patterns are higher-order
functions that can depend on natural numbers (nat) and data types (data). Rise
has no support for general recursion or iteration, and instead relies on
extensible and higher-level array patterns that explicitly expose data
parallelism [109, 110]: map applies a function to each element of an array;
reduce reduces all elements of an array to a single value given a binary
reduction function. Multi-dimensional array reshaping is common in array
languages [54], Rise has patterns such as split, join, or transpose for this
purpose.

⬇

1$+$ $\mid$ $\times$ : t $\to$ t $\to$ t

2map : (s $\to$ t) $\to$ n.s $\to$ n.t

3reduce : (t $\to$ t $\to$ t) $\to$ t $\to$ n.t $\to$ t

4split : (n: nat) $\to$ (n $\times$ m).t $\to$ m.n.t

5join : n.m.t $\to$ (n $\times$ m).t

6transpose : n.m.t $\to$ m.n.t

7slide : (sz sp: nat) $\to$ (sp $\times$ n + sz - sp).t $\to$ n.sz.t

8zip : n.s $\to$ n.t $\to$ n.(s $\times$ t)

9fst : (s $\times$ t) $\to$ s

10snd : (s $\times$ t) $\to$ t

Listing 11: Selection of Rise high-level patterns and their type

We write s $\to$ t for a function type with input of type $s$ and output of
type $t$, n.t for an array type with $n$ elements of type $t$, s $\times$ t
for a pair type with component types $s$ and $t$. To avoid the challenge of
storing functions (or closures) in heterogeneous hardware memory, the type
system ensures that function types cannot be stored in memory, only data types
(e.g. it is challenging for GPUs).

The Shine compiler rewrites a high-level Rise program into a low-level Rise
program that describes _how_ the result is computed, encoding implementation
decisions explicitly. Rise’s low-level patterns (LABEL:fig:low-level-patterns)
indicate specific implementation decisions.

⬇

1mapSeq : (s $\to$ t) $\to$ n.s $\to$ n.t

2reduceSeq : (s $\to$ t $\to$ s) $\to$ s $\to$ n.t $\to$ s

3mapGlobal(dim) : (s $\to$ t) $\to$ n.s $\to$ n.t

4toMem : (a: addr) $\to$ t $\to$ t

5asVector : (m: nat) $\to$ (n $\times$ m).t $\to$ n.<m>t

6asScalar : n.<m>t $\to$ (n $\times$ m).t

7vectorFromScalar : t $\to$ <m>t

Listing 12: Selection of Rise low-level patterns and their type

For example, mapSeq and reduceSeq respectively implement map and reduce with
sequential loops. Some low-level patterns are specific to the target
programming model (such as OpenCL) or hardware architecture (such as SIMD
vector support). For OpenCL, mapGlobal introduces parallelism by parallelising
across the dimension dim of global threads. The toMem pattern is used to
explicitly encode storing an expression in the given address space in memory
(a: addr). Other patterns enable SIMD vectorisation (e.g. asVector). A vector
type with $m$ elements of type $t$ is written as <m>t.

####  3.1.3 Code Generation through DPIA

To generate imperative code (e.g. C or OpenCL) from a low-level Rise program,
a hybrid functional-imperative language called DPIA is used as an intermediate
language. The DPIA language [106] is a variation of idealised ALGOL [111] that
separates program _phrases_ into three categories: functional _expressions_ ,
imperative _commands_ , and imperative _acceptors_. One can think of commands
as C statements, and of acceptors as C pointers that can be manipulated and
ultimately used for writing to a memory location. DPIA basically uses the same
functional constructs and data types as Rise. Imperative constructs are added,
such as assignment (p = q), sequential composition (;), memory allocation
(new), and loops (for). The type system reflects the separation into
expressions (exp[dt, rw]), acceptors (acc[dt]), and commands (comm). Functions
and pairs can freely combine phrases of all categories.

A selection of DPIA patterns is shown in LABEL:fig:dpia-functional and
LABEL:fig:dpia-imperative. While Rise introduces patterns as higher-order
functions with implicit type-level arguments more suitable for rewriting, DPIA
introduces patterns as fully applied values with explicit type-level arguments
more suitable for translation and code generation.

⬇

zip(n: nat, s: data, t: data, w: access,

a: exp[n.s, w], b: exp[n.t, w]): exp[n.(s x t), w]

fst(s: data, t: data, x: exp[(s x t), read]): exp[s, read]

snd(s: data, t: data, x: exp[(s x t), read]): exp[t, read]

mapSeq(n: nat, s: data, t: data, f: exp[s, read] $\to$ exp[t, write],

x: exp[n.s, read]): exp[n.t, write]

reduceSeq(n: nat, a: addr, s: data, t: data,

f: exp[s, read] $\to$ exp[t, read] $\to$ exp[s, write],

init: exp[s, write], x: exp[n.t, read]): exp[s, read]

toMem(a: addr, t: data, x: exp[t, write]): exp[t, read]

idx(n: nat, t: data, i: exp[idx[n], read], x: exp[n.t, read]): exp[t, read]

// ^ written as x[i]

Listing 13: Selection of functional DPIA patterns

⬇

assign(t: data, p: acc[t], q: exp[t, read]): comm // written as p = q

seq(c1: comm, c2: comm): comm // written as c1 ; c2

new(a: addr, t: data, body: (exp[t, read] x acc[t]) $\to$ comm): comm

for(n: nat, body: exp[idx[n], read] $\to$ comm): comm

zipAcc1(n: nat, s: data, t: data, x: acc[n.(s x t)]): acc[n.s]

zipAcc2(n: nat, s: data, t: data, x: acc[n.(s x t)]): acc[n.t]

idxAcc(n: nat, t: data, i: exp[idx[n], read], a: acc[n.t]): acc[t]

Listing 14: Selection of imperative DPIA patterns

The type system of DPIA enforces important invariants and assumptions.
Expression types include an access annotation, which can either be read or
write and enforces that DPIA programs must explicit encode how to write in
what memory. A read annotation signifies that the value can be read from
memory while a write annotation indicates that the value must first be written
to memory before being readable. To convert a write expression into a read
expression, the toMem pattern can be used and specifies a choice of address
space.

Loop patterns like for iterate over bounded indices of type idx[n], meaning
that the index value belongs to the interval $[0;n[$. This is important to
statically enforce safe array accesses and to remove the need for runtime
checks.

##### From Rise to DPIA

Going from low-level Rise to functional DPIA requires two steps which we will
not detail: (1) inferring access annotations to create DPIA types from Rise
types, the work of Bastian Köpcke; (2) translating patterns to be fully
applied values with explicit type-level arguments, a personal engineering
contribution.

To illustrate the second step, $\lambda$f x. mapSeq(n, s, t, f, x) would
result from translating mapSeq, extracting n, s, t from the instantiated Rise
type.

##### Translation to Imperative

Translating functional DPIA into imperative DPIA is performed by two mutually
recursive translation functions (LABEL:fig:types-acc-con-trans).

⬇

$\mathcal{A}$(e: exp[dt, write], out: acc[dt]): comm

$\mathcal{C}$(e: exp[dt, read], k: exp[dt, read] $\to$ comm): comm

Listing 15: Type of the acceptor and continuation translations

The _acceptor translation_ $\mathcal{A}$ produces a command which writes a
functional write expression to the memory represented by an acceptor
(LABEL:fig:acceptor-trans). The _continuation translation_ $\mathcal{C}$
produces a command which reads from a functional read expression and calls a
_continuation_ function to continue the translation as required
(LABEL:fig:continuation-trans). To start the overall translation process, we
need to know where to write the program result. For this, we generate an
output acceptor according to the data type computed by the functional
expression before invoking the $\mathcal{A}$ translation: imperative =
$\mathcal{A}$(functional, output). The output will later correspond to a
runtime parameter in the generated C function (as in Section 3.1.1), or OpenCL
kernel.

⬇

$\mathcal{A}$(zip(n, s, t, write, a, b), out) =

$\mathcal{A}$(a, zipAcc1(n, s, t, out)) ;

$\mathcal{A}$(b, zipAcc2(n, s, t, out))

$\mathcal{A}$(mapSeq(n, s, t, f, x), out) = $\mathcal{C}$(x, $\lambda$xT.

for(n $\lambda$i. $\mathcal{A}$(f(idx(n, s, i, xT), idxAcc(n, t, i, out)))))

$\mathcal{A}$(reduceSeq(n, a, s, t, f, init, x), out) =

$\mathcal{C}$(reduceSeq(n, a, s, t, f, init, x))($\lambda$r.
$\mathcal{A}$(r)(out))

Listing 16: Selection of acceptor translations

⬇

$\mathcal{C}$(zip(n, s, t, read, a, b), k) =

$\mathcal{C}$(a, $\lambda$aT. $\mathcal{C}$(b, $\lambda$bT. k(zip(n, s, t,
read, aT, bT))))

$\mathcal{C}$(fst(s, t, read, x), k) = $\mathcal{C}$(x, $\lambda$xT. k(fst(s,
t, read, xT)))

$\mathcal{C}$(snd(s, t, read, x), k) = $\mathcal{C}$(x, $\lambda$xT. k(snd(s,
t, read, xT)))

$\mathcal{C}$(reduceSeq(n, a, s, t, f, init, x), k) =

$\mathcal{C}$(x, $\lambda$xT. new(a, s, $\lambda$(accE, accA).

$\mathcal{A}$(init, accA) ;

for(n, $\lambda$i. $\mathcal{A}$(f(accE, idx(n, t, i, xT)), accA)) ;

k(accE)))

$\mathcal{C}$(toMem(a, t, x), k) = new(a, t, $\lambda$(tmpE, tmpA).
$\mathcal{A}$(e, tmpA) ; k(tmpE))

Listing 17: Selection of continuation translations

For example, translating the dotSeq program (LABEL:dot-seq-rise) from Rise to
functional DPIA and further to imperative DPIA leads to the following
translation steps (with syntactic sugar):

⬇

$\mathcal{A}$(reduceSeq ($\lambda$a x. a + (fst x) $\times$ (snd x)) 0 (zip a
b), output)

=

$\mathcal{C}$(reduceSeq ($\lambda$a x. a + (fst x) $\times$ (snd x)) 0 (zip a
b))($\lambda$r. $\mathcal{A}$(r)(output))

=

$\mathcal{C}$(zip a b, $\lambda$z. new (exp[f32, read] x acc[f32]) in
$\lambda$(accE, accA). (

$\mathcal{A}$(0, accA);

for n $\lambda$i.

$\mathcal{A}$(accE + (fst z[i]) $\times$ (snd z[i]), accA);

output = accE))

=

new (exp[f32, read] x acc[f32]) in $\lambda$(accE, accA). (

$\mathcal{A}$(0, accA);

for n $\lambda$i.

$\mathcal{A}$(accE + (fst (zip a b)[i]) $\times$ (snd (zip a b)[i]), accA);

output = accE)

=

new (exp[f32, read] x acc[f32]) in $\lambda$(accE, accA). (

accA = 0.0f;

for n $\lambda$i.

(accA = accE + ((fst (zip a b)[i]) * (snd (zip a b)[i])));

output = accE)

##### Imperative Passes

The use of DPIA in the Shine compiler enables transformation passes to be
applied at the imperative DPIA level, as we will show in Section 3.2. This is
valuable because imperative DPIA has more convenient abstractions and richer
semantics compared to C or OpenCL: for example, multi-dimensional arrays can
be used and precise type information has not been erased yet.

##### Code Generation

Finally, code such as C or OpenCL is generated from the imperative DPIA
program. One important transformation which happens at this stage is the
translation of indexing patterns (e.g. fst, snd, zip) into array indexing, and
the flattening [112] (or linearisation [113]) of multi-dimensional arrays into
one-dimensional arrays. The C code generated for the dotSeq program
(LABEL:dot-seq-rise) was shown in LABEL:dot-seq-c-code.

###  3.2 Implicit Barrier Insertion for Shine

When multiple threads access the same memory concurrently and at least one
thread is writing, synchronisation is required. Without synchronisation,
operation ordering is non-deterministic, and the computation may produce
incorrect results: this type of bug is called a _data race_ [114].

This section contributes DPIABI, a synchronisation barrier insertion algorithm
for Shine. DPIABI transforms imperative DPIA programs and does not need to be
modified when extending functional Rise patterns, as opposed to the barrier
elimination algorithm of Lift [107, 115]. Following Lift’s design, barrier
insertion is implicit and not controllable by rewriting, allowing Rise
programs to ignore low-level synchronisation requirements.

####  3.2.1 Synchronisation in OpenCL

In OpenCL kernels, calling barrier(flag) synchronises work-items from the same
work-group. The flag argument specifies which address spaces require a memory
fence to ensure correct ordering of memory operations (CLK_LOCAL_MEM_FENCE for
local memory and CLK_GLOBAL_MEM_FENCE for global memory). We must also ensure
that all work-items execute the barrier, to avoid undefined behaviour and
potential _deadlocks_ [116] where some work-items are indefinitely waiting for
others to encounter the barrier:

> “All the work-items of a work-group must execute the barrier before any are
> allowed to continue execution beyond the barrier. Note that the work-group
> barrier must be encountered by all work-items of a work-group executing the
> kernel or by none at all. ” (OpenCL 1.2 Specification)

OpenCL does not support synchronisation across work-groups inside a kernel,
instead multiple kernels must be launched. We will discuss how to add support
for launching multiple kernels in Section 3.3 and focus on work-item barriers
in this section.

####  3.2.2 Barrier Insertion as an Imperative DPIABI Transformation

Because functional Rise programs ignore low-level synchronisation
requirements, code generation must take care of inserting work-item barriers.
Barriers are relatively expensive [117, 118, 119], to maximise performance of
the generated code, we seek to minimise the number of barriers executed at
runtime while ensuring correct synchronisation.

Previous work on Lift follows a pessimistic approach, first inserting barriers
at the end of all parallel map patterns to avoid any data-race (e.g. after
mapGlobal, mapWorkGroup and mapLocal patterns), then eliminating barriers as
an optimisation [107, 115]. There are two problems with this approach. First,
the pessimistic approach is not sufficient to guarantee correctness because it
does not ensure that all work-items execute the barriers. Following the
pessimistic approach without any optimisation generates programs with
undefined behaviour, as expanded on in Section 3.2.5. Second, Lift’s barrier
elimination is analysing how functional patterns such as mapLocal are
composed, and their interaction with other patterns that might result in data
sharing between work-items (split, join, etc). We argue that this requires the
definition of too many special cases (for example, special treatment is
required when encountering several mapGlobal inside the arguments of a zip, or
nested mapLocal patterns). Additionally, such an approach may require
modifying the barrier elimination optimisation when new functional patterns
are added, reducing the ease of compiler extension.

Therefore, we follow instead an optimistic approach, inserting barriers as
required in imperative DPIA programs. Before generating OpenCL kernel code
from imperative DPIA, the following imperative passes are applied in sequence:

  1. 1.

Inject work item sizes. The OpenCL number of groups, local size and global
size values are inlined if they are statically known.

  2. 2.

Flag private array loops. Loops iterating over private arrays are flagged to
be unrolled, which helps OpenCL perform register allocation.

  3. 3.

Unroll loops. The unrolled loops have been flagged in the previous pass, or
explicitly during rewriting.

  4. 4.

Simplify nats. Tries to simplify certain arithmetic expressions (e.g. $x\times
0=0$), such simplifications are not only applied here but also throughout the
compilation process.

  5. 5.

Insert memory barriers. The topic of this section realised in the DPIABI
algorithm.

  6. 6.

Hoist memory allocations. In OpenCL, global and local memory cannot be
allocated while the kernel is running, it must be allocated upfront. This pass
hoists such memory allocations as required (e.g. outside of loops).

  7. 7.

Adapt kernel parameters. Kernel parameters in global or local memory which
have a scalar type in DPIA are represented as arrays of size 1 for OpenCL.

####  3.2.3 Definition of DPIABI

DPIABI is defined by analysing reads and writes to memory, making conservative
barrier insertions to ensure synchronisation according to the observed data
dependencies.

The DPIABI analysis records memory reads and writes in two mutable
dictionaries from memory identifier to address space, one representing reads
from outer scope allocations, and one representing work-group parallel writes
to outer scope allocations. This corresponds to the following record structure
(a case class is used in Scala):

DPIABI is defined recursively in LABEL:lst:barrier-insertion. The DPIABI
function is the entry point and calls
$\text{DPIA}_{\text{BI}}^{\text{fresh\\_rec}}$, which creates empty data of
type D before calling $\text{DPIA}_{\text{BI}}^{\text{rec}}$. Barriers are
represented by the barrier DPIA command and may be inserted in-between
sequences (inserted for seq in Line 30), or at the end of loop bodies
(inserted for for or parFor in Line 41).

⬇

1def DPIABI(p: DPIAcomm): DPIAcomm =
$\text{DPIA}_{\text{BI}}^{\text{fresh\\_rec}}$(p, Map())1

2

3def $\text{DPIA}_{\text{BI}}^{\text{fresh\\_rec}}$(p: DPIAcomm,

4 allocs: Map[Ident, Address]): (DPIAcomm, D) =

5 def d = D(MutMap(), MutMap())

6 ($\text{DPIA}_{\text{BI}}^{\text{rec}}$(p, allocs, d), d)

7

8def $\text{DPIA}_{\text{BI}}^{\text{rec}}$(p: DPIAcomm, allocs: Map[Ident,
Address], d: D): DPIAcomm =

9 match p

10 for(n, $\lambda$x. body) =>

11 for(n, $\lambda$x. $\text{DPIA}_{\text{BI}}^{\text{loop\\_body}}$(body,
allocs, d, MutMap()))

12 parFor(level, ..., out, $\lambda$x o. body) =>

13 parFor(level, ..., out, $\lambda$x o.
$\text{DPIA}_{\text{BI}}^{\text{loop\\_body}}$(body, allocs, d,

14 if level == local then collect writes by inspecting out

15 else MutMap()))

16 new(adr, dt, $\lambda$x. body) =>

17 def allocs2 = if adr != private then allocs + (x -> adr) else allocs

18 new(adr, dt, $\lambda$x. $\text{DPIA}_{\text{BI}}^{\text{rec}}$(body,
allocs2, d))

19 assign(dt, lhs, rhs) =>

20 collect reads into d.reads by inspecting rhs

21 assign(dt, lhs, rhs)

22 seq(a, b) =>

23 def (a2, ad) = $\text{DPIA}_{\text{BI}}^{\text{fresh\\_rec}}$(a, allocs)

24 def (b2, bd) = $\text{DPIA}_{\text{BI}}^{\text{fresh\\_rec}}$(b, allocs)

25 def dependencies = (ad.reads.keys $\cap$ bd.wg_writes.keys)

26 $\cup$ (bd.reads.keys $\cap$ ad.wg_writes.keys)

27 extend d.reads with ad.reads and bd.reads

28 extend d.wg_writes with ad.wg_writes and bd.wg_writes

29 if dependencies non empty then

30 seq(a2, seq(make barrier according to dependencies, b2))

31 else

32 seq(a2, b2)

33 [...]

34

35def $\text{DPIA}_{\text{BI}}^{\text{loop\\_body}}$(p: DPIAcomm, allocs:
Map[Ident, Address], d: D,

36 outer_wg_writes: MutMap[Ident, Address]): DPIAcomm =

37 def p2 = $\text{DPIA}_{\text{BI}}^{\text{rec}}$(p, allocs, d)

38 def dependencies = d.wg_writes.keys $\cap$ d.reads.keys

39 if dependencies non empty then

40 clear d.reads and d.wg_writes

41 seq(p2, make barrier according to dependencies)

42 else

43 extend d.wg_writes with outer_wg_writes

44 p2

Listing 18: Pseudo-code for the DPIABI barrier insertion algorithm

To analyse dependencies, reads and work-group parallel writes are collected in
lines 20 and 14. The parFor pattern corresponds to parallel loops (e.g.
parFor(local, ...) corresponds to mapLocal), and differently from sequential
loops, explicitly mentions its global output (out) and its local output (o),
making analysis easier. Information about memory allocations and their address
space is also collected in line 17, since this information is critical to
select the right barrier flags (CLK_LOCAL_MEM_FENCE and CLK_GLOBAL_MEM_FENCE).

####  3.2.4 Evaluating the Correctness and Efficiency of DPIABI

##### Experimental Setup

The correctness and efficiency of barrier insertion is evaluated on 48
programs (38 unit tests and 10 benchmarks), by observing and comparing the
code generated by Lift and Shine. The generated code is considered correct if
it respects the OpenCL 1.2 specification. The generated code is considered
more efficient if it executes less barriers at runtime, or uses less flags.

The 38 unit tests mainly come from the Lift
repository222https://github.com/lift-
project/lift/blob/5e8a18df48ab791ae66016007e18826047a015f7/src/test/opencl/generator/TestBarrier.scala,
with a few additions. The unit tests compose various patterns in various ways:
parallel and sequential map patterns to introduce threads and loops, toMem to
introduce temporary memory, and indexing patterns to introduce data sharing
between threads.

The 10 benchmarks come from the experimental evaluation of Lift presented in
[107]. The benchmarks target GPU hardware and represent different domains
including physics simulations, statistics, and linear algebra. Two benchmarks
were ported but elided from the results, as they are not fairly comparable:
"NBody NVIDIA" and "MM NVIDIA". Shine allocates additional memory for these
programs due to orthogonal issues.

##### Correctness and Efficiency Results

Table 3.2 summarises the results, identifying 6 differences between the code
generated by Shine and the one generated by Lift:

Impact | Difference | Unit tests | Benchmarks | Total
---|---|---|---|---
Lift incorrect | Incorrect flags |  13 (34%) |  1 (10%) |  14 (29%)
Lift more efficient | Better analysis |  13 (34%) |  0 (0%) |  13 (27%)
Single iteration loops |  6 (16%) |  1 (10%) |  7 (15%)
Rise more efficient | Better analysis |  0 (0%) |  2 (20%) |  2 (4%)
Barrier position |  1 (2%) |  0 (0%) |  1 (2%)
Trade-off | Sync. vs alloc. |  7 (18%) |  0 (0%) |  7 (15%)
Total |  38 (100%) |  10 (100%) |  48 (100%)
Table 3.2: Number and percentage of unit tests and benchmarks that exhibit
each difference identified in the code generated by Shine and Lift. Shine
generates correct barriers in all cases, which is not the case of Lift.

  1. 1.

_Incorrect flags ( Lift incorrect)._ While Shine generates correct barriers
for all programs, Lift generates barriers with incorrect flags for 14 programs
(29%). For the test shown in Figure 3.2, a barrier on global memory is
generated instead of a barrier on local memory. In practice, the program still
computes the expected result, however the OpenCL 1.2 specification is not
respected and this could lead to undefined behaviour with a different OpenCL
implementation. With DPIABI, Shine inserts barriers with correct flags: one in
line 5 to protect the read-after-write dependency on local memory, and one in
line 10 to protect the write-after-read dependency between outer loop
iterations.

⬇

1in $\vartriangleright$ mapWorkGroup(0) (

2 mapLocal(0) ($\lambda$x. x) $\vartriangleright$

3 toMem local $\vartriangleright$

4 slide 3 1 $\vartriangleright$

5 mapLocal(0) sum

6))

⬇

1for (int wg0 = get_group_id(0); wg0 < n; wg0 += get_num_groups(0)) {

2 for (int l0 = get_local_id(0); l0 < m; l0 += get_local_size(0)) {

3 [...] // read from global input; write to local memory

4 }

5 barrier(CLK_LOCAL_MEM_FENCE);

6 for (int l1 = get_local_id(0); l1 < m-2; l1 += get_local_size(0)) {

7 [...] // read from local memory; write to global output

8 }

9- barrier(CLK_GLOBAL_MEM_FENCE);

10+ barrier(CLK_LOCAL_MEM_FENCE);

11}

Figure 3.2: Example low-level Rise program and its generated OpenCL code. In
red (-), barrier that would be generated by Lift. In green (+), barrier
generated by Shine. The second barrier flag is incorrect with Lift, but
correct with Shine.

  2. 2.

_Better analysis ( Lift more efficient)._ DPIABI only tracks memory accesses
for entire allocations. Lift generates fewer barriers for 13 programs (27%) by
reasoning about data sharing in the functional program. This difference is not
visible in the benchmarks that exclusively use local memory to share data
between work-items. However, integrating a more precise dependency analysis
into DPIABI, for example by tracking accessed intervals using symbolic natural
numbers, would enable Shine to generate as many or fewer barriers than Lift in
all observed cases.

  3. 3.

_Single iteration loops ( Lift more efficient)._ Lift avoids generating
barriers for loops whose body is only executed once in 7 programs (15%), while
Shine generates superfluous barriers. Shine could avoid generating such
barriers by eliminating single iteration loops before barrier insertion,
detecting such loops during barrier insertion, or by differentiating between
two scenarios which are currently conflated both in Lift and Rise: (1) the
loop will be executed again; or (2) the loop will be exited.

  4. 4.

_Better analysis ( Rise more efficient)_. Rise avoids generating barriers that
Lift generates, or uses more precise flags, in 2 programs (4%). This is
because Lift fails to reason about specific combinations of functional
patterns, while Shine handles these combinations naturally at the imperative
level, without any special treatment.

  5. 5.

_Barrier position ( Rise more efficient)_. Lift only inserts barriers at the
end of parallel map patterns. For 1 program (2%) this strategy is not flexible
enough and leads to inefficient barriers compared to Shine. Figure 3.3 shows
an example where Shine is able to generate barriers at the end of sequential
loops, which results in a program executing $2\times(m-1)\times n$ times fewer
barriers at runtime ($m-1$ less per sequential loop execution).

  6. 6.

_Synchronisation vs allocation (Trade-off)_. For 7 (15%) programs, Lift
allocates more memory than Rise. As a result, Lift generates less barriers
because the memory is not re-used while Rise generates more barriers because
the memory is re-used. It is unclear which trade-off is best without
benchmarking a specific program on a specific hardware target. It may even be
that a fine-grain combination of different options for different allocations
is optimal. Future work may explore this trade-off.

⬇

1 in $\vartriangleright$ mapWorkGroup(0) (

2 mapSeq (mapLocal(0) ($\lambda$x. x)) $\vartriangleright$

3 toMem local $\vartriangleright$

4 map (slide 3 1) $\vartriangleright$

5 mapSeq (mapLocal(0) sum)

6 ))

⬇

1for (int wg0 = get_group_id(0); wg0 < n; wg0 += get_num_groups(0)) {

2 for (int i = 0; < m; i += 1) {

3 for (int l0 = get_local_id(0); l0 < o; l0 += get_local_size(0)) {

4 [...] // read from global input; write to local memory

5 }

6- barrier(CLK_LOCAL_MEM_FENCE);

7 }

8+ barrier(CLK_LOCAL_MEM_FENCE);

9 for (int i = 0; < m; i += 1) {

10 for (int l1 = get_local_id(0); l1 < o-2; l1 += get_local_size(0)) {

11 [...] // read from local memory; write to global output

12 }

13- barrier(CLK_GLOBAL_MEM_FENCE);

14 }

15+ barrier(CLK_LOCAL_MEM_FENCE);

16}

Figure 3.3: Example low-level Rise program and its generated OpenCL code. In
red (-), barrier that would be generated by Lift. In green (+), barrier
generated by Shine. The second barrier flag is only correct with Shine. The
code generated by Shine is also more efficient: $2\times(m-1)\times n$ times
less barriers will be executed at runtime.

####  3.2.5 Limitations of DPIABI

In addition to the "_Better analysis_ " and "_Single iteration loops_ " issues
(Items 2 and 3 above), there are four other limitations worth mentioning.

##### Barrier Reachability

In some cases, the generated barriers might not be reached by all work-items,
which would result in undefined behaviour. By optimising barrier placement on
the 48 programs from Section 3.2.4, this issue is avoided in 11 unit tests and
3 benchmarks compared to following the naive pessimistic approach of Lift with
barrier elimination optimisations turned off. However, even when optimising
barrier placement, the issue remains in 2 unit tests. An example is given in
Figure 3.4, where part of the work-items might not enter the loop in line 3.
Previous work on Lift suffered from similar limitations and implemented a mix
of compilation time and runtime checks to report the issue. For Shine, an
additional imperative DPIA pass could be implemented to fix the code as
illustrated in Figure 3.4.

⬇

1 in $\vartriangleright$ mapWorkGroup(1) (mapWorkGroup(0) (

2 mapLocal(1) (

3 mapLocal(0) ($\lambda$x. x) $\vartriangleright$

4 toMem(local) $\vartriangleright$

5 slide 3 1 $\vartriangleright$

6 mapLocal(0) sum

7 ))))

⬇

1for (int wg0 = get_group_id(1); wg0 < n; wg0 += get_num_groups(1)) {

2 for (int wg1 = get_group_id(0); wg1 < m; wg1 += get_num_groups(0)) {

3- for (int l0 = get_local_id(1); l0 < o; l0 += get_local_size(1)) {

4+ for (int l0 = get_local_id(1); l0 < ctt(o); l0 += get_local_size(1)) {

5+ if (l0 < o) {

6 for (int l1 = get_local_id(0); l1 < p; l1 += get_local_size(0)) {

7 [...] // read from global input; write to local memory

8 }

9+ }

10 barrier(CLK_LOCAL_MEM_FENCE);

11+ if (l0 < o) {

12 for (int l2 = get_local_id(0); l2 < p-2; l2 += get_local_size(0)) {

13 [...] // read from local memory; write to global output

14 }

15+ }

16 barrier(CLK_LOCAL_MEM_FENCE);

17 }

18 }

19}

Figure 3.4: Example low-level Rise program and its generated OpenCL code where
barriers may not be encountered by all work-items. In red (-), buggy code that
would be generated by both Lift and Shine. In green (+), a potential fix where
the ctt function rounds up a number to a multiple of the involved work-items.

##### Extensibility of Indexing Patterns

To simplify DPIABI further and separate concerns, it would be valuable to
eliminate indexing patterns (e.g. split, join) beforehand. This would enable
extending indexing patterns without having to adjust the "collect reads"
operation, which would only need to deal with idx patterns (line 20 of
LABEL:lst:barrier-insertion).

##### Extensibility of Imperative Patterns

Adding a new imperative DPIA pattern requires adding a new case in
LABEL:lst:barrier-insertion, reducing the ease of extensibility. However, we
argue that this improves on previous Lift work where barrier elimination
happened at the functional level. In Shine, it is possible to add new
functional patterns without having to adapt DPIABI, as long as no new
imperative patterns are required.

##### Formal Verification

No proof of correctness is provided. Future work could provide an algorithm
that provably prevents any data race or undefined behaviour for any well-
formed low-level Rise program.

####  3.2.6 Related Work

##### Barriers for Collective Patterns

When using data-parallel collective patterns like map, a common approach is to
insert synchronisation barriers after collective operations. This is the
approach taken by the array-based languages Lift [15], Futhark [56], and SaC
[120].

##### Barriers for Imperative Programs

When transforming imperative code, using dependency analyses to reason about
synchronisation barriers is common [121]. Both barrier insertion [117, 122,
123, 124, 125] and barrier elimination [126, 127, 128, 129] have been studied,
often using advanced dependency analyses, such as polyhedral analyses. The
Tiramisu compiler explicitly exposes synchronisation decisions to the
performance engineer, checking validity using the polyhedral model [130].

##### Barriers in Domain-Specific Compilers

The problem of inserting correct and efficient barriers is barely discussed by
the Halide and TVM papers. For example, the TVM paper [10] simply mentions
that "memory synchronization barriers must be properly inserted to guarantee
that shared loaded data is visible to consumers". Following papers show that,
for some use cases, the TVM barrier insertion pass is unsatisfactory, leading
to poor performance and forcing users to modify the compiler [131].

The novelty of DPIABI is to insert barriers using imperative dependency
analyses in a domain-extensible compiler for a functional language based on
data-parallel collective patterns. We showed how, in the Shine domain-
extensible compiler, the design of DPIABI facilitates extending functional
Rise patterns and can even lead to generating more efficient barriers (Figure
3.3).

####  3.2.7 Summary

This section contributes DPIABI, a synchronisation barrier insertion algorithm
for Shine that transforms imperative DPIA programs. Crucially, DPIABI does not
need to be modified when extending Rise patterns, as opposed to the barrier
elimination algorithm of Lift [107, 115]. Following Lift’s design, barrier
insertion is implicit and not controllable by rewriting, allowing Rise
programs to ignore low-level synchronisation requirements.

The correctness and efficiency of DPIABI is evaluated on 38 unit tests and 10
benchmarks, mostly taken from prior Lift work. We identify 6 differences in
the code generated by Shine and Lift (Table 3.2), and observe that DPIABI
fixes Lift bugs in 13 unit tests and 1 benchmark. There is only 1 benchmark
where Shine inserts a barrier that Lift eliminates, and we provide a clear
pathway to improve DPIABI to generate more efficient barriers than Lift on all
48 unit tests and benchmarks.

###  3.3 Explicit Kernel Execution for Shine

Any high-level language that allows offloading computations using the OpenCL
or CUDA programming models must deal with executing kernels from the host.
Shine allows Rise programs to be compiled for CPUs by generating C code, or to
various hardware by generating OpenCL kernel code. However, so far Rise
programs cannot be compiled for a host CPU while offloading computations to
various hardware devices. Even if no computations are performed on the host,
host-side calls to the OpenCL C API are tightly coupled to the use of OpenCL
kernels. The host is responsible for dynamically allocating global/local
memory, launching kernels, and transferring data between host and device.
Launching multiple kernels is also the only way to synchronise across work-
groups, if strictly following the OpenCL 1.2 specification (Section 3.2.1).

This section enables explicit kernel execution in Rise by adding the oclRun
pattern. Crucially, this design makes kernel execution controllable by
rewriting in the domain-extensible Shine compiler, enabling to explore
different kernel decompositions of a high-level Rise program via rewriting.
The oclRun pattern allows a single low-level Rise program to express both
host-side and device-side computations. Shine is modified to generate
imperative code for multiple OpenCL kernels, as well as the necessary C host
code to launch them. To simplify host code generation, we implement a thin
runtime called LRA that abstracts over some of the OpenCL runtime details. The
focus is to support executing kernels on a single OpenCL device, support for
multiple devices is left for future work.

####  3.3.1 Adding a Rise Pattern for Kernel Execution

To make device execution explicit, while keeping the interface as simple as
possible, we introduce an oclRun pattern. This pattern is an identity function
which takes 3D local and global thread counts as additional input, it has the
following type:

⬇

oclRun : (ls0, ls1, ls2, gs0, gs1, gs2 : nat) $\to$ t $\to$ t

The input expression of oclRun will be computed on the device by generating
and calling a corresponding kernel. We also provide syntactic sugar to enable
providing 1D, 2D or 3D local and global thread counts, where the thread counts
of unspecified dimensions are set to 1:

⬇

oclRun(ls: LocalSize, gs: GlobalSize) = oclRun ls.x ls.y ls.z gs.x gs.y gs.z

This new pattern can be used to write Rise programs such as the one in
LABEL:fig:kernels-example, where store is similar to toMem without the OpenCL
address space parameter, and with a function parameter to provide name
binding. First, the value 3 is added to each element of in in parallel on the
device (line 4), storing the result into x1. Second, the value 1 is subtracted
to each element of x1 sequentially on the host (line 5), storing the result
into x2. Finally, the value 1 is added to each element of x2 in parallel on
the device (line 6), producing the final output.

⬇

1store : (s $\to$ t) $\to$ s $\to$ t

2

3$\Lambda$n : nat. $\lambda$in : n.i32. (

4 oclRun(LocalSize(2), GlobalSize(n/2)) (mapGlobal (add 3) in)
$\vartriangleright$ store ($\lambda$x1.

5 mapSeq ($\lambda$y. y-1) x1 $\vartriangleright$ store ($\lambda$x2.

6 oclRun(LocalSize(4), GlobalSize(n/4)) (mapGlobal (add 1) x2)))

7)

Listing 19: Example low-level Rise program mixing host computation (middle
mapSeq call) and device computation (both oclRun calls).

####  3.3.2 Creating a Lightweight Runtime Abstraction (LRA)

Low-level OpenCL host code needs to explicitly manage devices, kernels and
memory in great detail. To simplify code generation, a lightweight runtime
abstraction called LRA is designed, delegating some implementation decisions
to the runtime implementation. This runtime is meant to be simplistic and to
enable quick prototyping of our ideas, future work may implement or re-use
more sophisticated runtime libraries, if more advanced features are required.

##### Buffer Abstraction

The most important abstraction provided by LRA is the Buffer abstraction which
offers a unified and consistent view on memory that may reside on the host or
on the device, as required during program execution (LABEL:fig:buffer-abs).
This is not a new idea, as libraries such as SPOC [132], SkePU [133] and
SkelCL [134] provide a similar abstraction.

A Buffer is created with createBuffer and deleted with destroyBuffer. Creating
a buffer requires specifying its size in bytes and the type of accesses that
it should support. AccessFlags express a combination of reads/writes from the
host/device. To access a buffer, explicit synchronisation functions must be
called. hostBufferSync allows accessing a buffer from the host, returning a
void* pointer that may be used in regular C code. deviceBufferSync allows
accessing a buffer from the device, returning a DeviceBuffer structure that
may be used as a kernel parameter. The runtime user (e.g. the Shine compiler)
is responsible for avoiding data races while using this interface.333reminder:
multiple threads should not access the same memory concurrently if at least
one thread is writing

⬇

1typedef enum {

2 HOST_READ = 1 << 0,

3 HOST_WRITE = 1 << 1,

4 DEVICE_READ = 1 << 2,

5 DEVICE_WRITE = 1 << 3,

6} AccessFlags;

7

8Buffer createBuffer(Context ctx, size_t byte_size, AccessFlags access);

9void destroyBuffer(Context ctx, Buffer b);

10void* hostBufferSync(Context ctx, Buffer b,

11 size_t byte_size, AccessFlags access);

12DeviceBuffer deviceBufferSync(Context ctx, Buffer b,

13 size_t byte_size, AccessFlags access);

Listing 20: Interface for the Buffer abstraction.

We provide two different Buffer implementations: one with memory copies, and
one without. With memory copies, the buffer is allocated both in host memory
and in device memory. Memory copies are issued between both allocations as
required to maintain consistency. Without memory copies, the buffer is
allocated in memory that is accessible by both host and device. Memory
synchronizations are issued as required to maintain consistency. In both
cases, memory copies and synchronisations are performed lazily by keeping
track of dependencies between accesses at runtime in a very simple and coarse-
grain manner. Each implementation requires no more than 80 lines of C code,
illustrating how lightweight the abstractions are.

##### Context Abstraction

As already visible in LABEL:fig:buffer-abs, LRA uses a Context abstraction
(LABEL:fig:context-abs). The context is used to manage global state required
for OpenCL kernel execution. A Context may be created using
createDefaultContext or createContext, and deleted using destroyContext.
createContext enables specifying the OpenCL platform and OpenCL device desired
for execution, if multiple ones are available.

⬇

1Context createDefaultContext();

2Context createContext(const char* platform_subname,

3 const char* device_type_str);

4void destroyContext(Context ctx);

Listing 21: Interface for the Context abstraction.

##### Kernel Abstraction

Finally, LRA provides a Kernel abstraction (LABEL:fig:kernel-abs). A Kernel is
loaded from OpenCL kernel source code (loadKernelFrom...) encoded in a string,
or in a file. A Kernel is launched using launchKernel, which requires
specifying global and local thread counts as well as providing kernel
arguments (e.g. scalars, device buffers). It is deleted using destroyKernel.

⬇

1Kernel loadKernelFromSource(Context ctx, const char* name,

2 const char* source, size_t length);

3Kernel loadKernelFromFile(Context ctx, const char* name, const char* path);

4void launchKernel(Context ctx, Kernel k,

5 const size_t* global_size, const size_t* local_size,

6 size_t arg_count, const KernelArg* args);

7void destroyKernel(Context ctx, Kernel k);

Listing 22: Interface for the Kernel abstraction.

In addition to the two possible Buffer implementations, we provide a single
implementation for Context and Kernel, which consists of about 300 lines of
boilerplate C code.

##### Code Generation Objective

Given LRA, the code we want to generate from LABEL:fig:kernels-example is
depicted in LABEL:fig:kernels-code. The kernel source codes are provided as
string constants (e.g. k0_source), each one corresponding to an oclRun call. A
add3_t structure holds the execution environment for the computation (i.e. the
two Kernels), a add3_init function initializes this environment and a
add3_destroy function deletes it. Finally, a add3_run function runs the actual
computation: allocating Buffers, performing synchronisations, launching
kernels and performing host-side computations.

Crucially, if the Rise program is rewritten differently, the interface of the
generated code does not change (i.e. the signature of add3_init, add3_destroy
and add3_run). The Buffer abstraction is particularly important to enable
this: inputs and outputs may freely reside on either host or device without
having to generate different code, as the runtime will take care of both cases
seamlessly.

⬇

1const char k0_source[] =

2"__kernel␣__attribute__␣((reqd_work_group_size(2,␣1,␣1)))"

3"void␣k0(global␣int*␣restrict␣output,␣int␣n,␣const␣global␣int*␣restrict␣in){"

4"␣␣for␣(int␣gl_id␣=␣get_global_id(0);␣gl_id␣<␣n;␣gl_id␣=␣gl_id␣+␣(n␣/␣2))␣{"

5"␣␣␣␣output[gl_id]␣=␣((int)3)␣+␣in[gl_id];"

6"␣␣}"

7"}";

8

9const char k1_source[] =

10"__kernel␣__attribute__␣((reqd_work_group_size(4,␣1,␣1)))"

11"void␣k1(global␣int*␣restrict␣output,␣int␣n,␣const␣global␣int*␣restrict␣in){"

12"␣␣for␣(int␣gl_id␣=␣get_global_id(0);␣gl_id␣<␣n;␣gl_id␣=␣gl_id␣+␣(n␣/␣2))␣{"

13"␣␣␣␣output[gl_id]␣=␣((int)1)␣+␣in[gl_id];"

14"␣␣}"

15"}";

16

17typedef struct add3_t {

18 Kernel k0;

19 Kernel k1;

20} add3_t;

21

22void add3_init(Context ctx, add3_t* self){

23 (*self).k0 = loadKernelFromSource(ctx, "k0", k0_source, sizeof(k0_source) -
1);

24 (*self).k1 = loadKernelFromSource(ctx, "k1", k1_source, sizeof(k1_source) -
1);

25}

26

27void add3_destroy(Context ctx, add3_t* self){

28 destroyKernel(ctx, (*self).k0);

29 destroyKernel(ctx, (*self).k1);

30}

31

32void add3_run(Context ctx, add3_t* self, Buffer moutput, int n, Buffer min){

33 Buffer mx1 = createBuffer(ctx, ..., HOST_READ | DEVICE_WRITE);

34 {

35 DeviceBuffer b0 = deviceBufferSync(ctx, mx1, ..., DEVICE_WRITE);

36 DeviceBuffer b2 = deviceBufferSync(ctx, min, ..., DEVICE_READ);

37 launchKernel(ctx, (*self).k0, {n / 2, 1, 1}, {2, 1, 1}, 3, {b0, n, b2});

38 }

39 Buffer mx2 = createBuffer(ctx, ..., HOST_WRITE | DEVICE_READ);

40 {

41 int32_t* x2 = (int32_t*)hostBufferSync(ctx, mx2, ..., HOST_WRITE);

42 int32_t* x1 = (int32_t*)hostBufferSync(ctx, mx1, ..., HOST_READ);

43 for (int i = 0; i < n; i = 1 + i) {

44 x2[i] = x1[i] - ((int32_t)1);

45 }

46 }

47 {

48 DeviceBuffer b0 = deviceBufferSync(ctx, moutput, ..., DEVICE_WRITE);

49 DeviceBuffer b2 = deviceBufferSync(ctx, mx2, ..., DEVICE_READ);

50 launchKernel(ctx, (*self).k1, {n / 2, 1, 1}, {4, 1, 1}, 3, {b0, n, b2});

51 }

52 destroyBuffer(ctx, mx2);

53 destroyBuffer(ctx, mx1);

54}

Listing 23: Code using LRA that should be generated from LABEL:fig:kernels-
example

####  3.3.3 Compiling Kernel Executing Programs in Shine

To generate code targeting LRA from kernel executing Rise programs, we take
inspiration from partial compilers [135] to modify Shine. A _partial compiler_
delegates intermediate compilation tasks to child (partial) compilers before
combining their generated intermediate code into the final code. We implement
the partial C+OpenCL Shine compiler that delegates C host code generation and
OpenCL kernel code generation to dedicated child compilers (Figure 3.5).

Figure 3.5: The partial C+OpenCL Shine compiler delegates tasks to child
compilers.

##### Partial C+OpenCL Compiler

Programs with explicit kernel executions such as the one from
LABEL:fig:kernels-example can be partially compiled by decomposing them into
two simpler compilation tasks: host code generation and kernel code
generation. In the initial program, a kernel definition with name k$i$ is
generated for the $i$-th oclRun pattern through _lambda lifting_ [136]: all
free variables in the input expression of oclRun are eliminated by introducing
explicit parameters (LABEL:fig:kernels-k0 and LABEL:fig:kernels-k1). Each k$i$
can be compiled separately, while the associated oclRun pattern is replaced
with a kernelCall pattern that treats k$i$ as a regular callable function:
producing the host program which can also be compiled separately
(LABEL:fig:kernels-host).

Once separately compiled by child compilers, kernel code and host code are put
together in a single file as done in LABEL:fig:kernels-code.

⬇

1$\Lambda$n : nat. $\lambda$in : n.i32. mapGlobal (add 3) in

Listing 24: Rise k0 kernel program for LABEL:fig:kernels-example

⬇

1$\Lambda$n : nat. $\lambda$x2 : n.i32. mapGlobal (add 1) x2

Listing 25: Rise k1 kernel program for LABEL:fig:kernels-example

⬇

1$\Lambda$n : nat. $\lambda$in : n.i32. (

2 kernelCall("k0", LocalSize(2), GlobalSize(n/2), n, in) $\vartriangleright$
store (

3 mapSeq ($\lambda$y. y-1) x1 $\vartriangleright$ store ($\lambda$x2.

4 kernelCall("k1", LocalSize(4), GlobalSize(n/4), n, x2)))

5)

Listing 26: Rise host program for LABEL:fig:kernels-example

##### OpenCL Kernel Compiler

Generating kernel code does not require anything new: the existing Shine
OpenCL kernel compiler is simply re-used.

##### C Host Compiler

To generate host code, the existing Shine C compiler requires extensions to
target LRA and support compiling the kernelCall patterns (which become
kernelCallCmd imperative DPIA commands after a straightforward translation to
imperative).

First, we add an imperative DPIA compilation pass that, given a host program
which uses plain arrays, introduces LRA Buffers as needed, using the newBuffer
pattern:

⬇

k: (exp[buffer[dt], read] $\times$ acc[buffer[dt]]) $\to$ comm): comm

Sub-programs that only perform host-side computation are kept as-is, they can
still use plain arrays and the pre-existing C compilation, we call them _host
execution sections_ and introduce a corresponding hostExecution imperative
DPIA pattern:

To identify host execution sections, the algorithm used is fairly similar to
the barrier insertion algorithm from LABEL:lst:barrier-insertion. Host and
device side reads/writes are tracked to insert hostExecution patterns as
required.

Finally, C code generation from imperative DPIA is extended to support the
three new patterns. For each kernelCallCmd, deviceBufferSync and launchKernel
LRA calls are generated in a C block for scoping ({}). For each newBuffer,
createBuffer and destroyBuffer calls are generated in a block. For each
hostExecution, hostBufferSync calls are generated in a block also containing
the code generated for the host execution body. As desired, the code generated
for the Rise program in LABEL:fig:kernels-example corresponds to the C+OpenCL
program in LABEL:fig:kernels-code, modulo small syntactic details.

####  3.3.4 Evaluating the Boilerplate Avoided through Kernel Execution

##### Experimental Setup

The impact of using the oclRun Rise pattern for kernel execution is evaluated
by measuring how many lines of host code are automatically generated to launch
the kernels, rather than written by hand. As a case study, we consider
optimising the Harris corner (and edge) detector [30], a well established
image processing pipeline.

Figure 3.6 depicts the Harris algorithm for a grayscale input image. Given the
input on the left, point-to-point operators (multiplications $\times$,
coarsity) and $3\times 3$ convolutions (sobel operators $S_{x}$ and $S_{y}$,
sums $+$) are combined to detect corners and edges highlighted in the output
on the right. As a composition of point wise and stencil operators, the Harris
detector is more complex than its individual parts, and exposes more
optimisation opportunities.

Figure 3.6: Harris corner detection computation flow, example image taken from
Halide.

We consider well-known optimisations presented in [23]. Table 3.3 shows how
many different OpenCL kernels need to be implemented to explore a space of
important optimisations. In particular, different decisions for an
optimisation called operator fusion lead to different kernel groupings. 4
different ways to fuse operators are considered, each variant offering a
different trade-off between memory accesses and arithmetic complexity. 5
different ways to optimise each resulting kernel are also considered (e.g.
tiling, vectorisation).

In total, benchmarking this relatively small optimisation space requires
implementing 50 different kernels, each associated with slightly different
host code to launch them.

##### Avoided Boilerplate Results

Table 3.4 shows how many lines of code are handwritten and generated with or
without the kernel execution feature (i.e. the oclRun pattern). Without
oclRun, 1K lines of host code are handwritten to launch the generated kernels.
This corresponds to about 20 lines per kernel launch, and does not even count
the code required to compose the kernels into the entire Harris processing
pipeline. With oclRun, 1.2K lines of host code are automatically generated.
Therefore, the kernel execution feature enables exploring the same
optimisation space without having to write 1K lines of boilerplate code.

kernel grouping | variants per group | kernel variants
---|---|---
$[S_{x}],[S_{y}],[\times],[+],[coarsity]$ | 5 | 25
$[S_{x},S_{y},\times],[+,coarsity]$ | 5 | 10
$[S_{x},S_{y}],[\times,+,coarsity]$ | 5 | 10
$[S_{x},S_{y},\times,+,coarsity]$ | 5 | 5
_total_ | _total_ | 50
Table 3.3: Four different kernel groupings (a group is represented as $[]$), and the number of differently optimised kernel variants that are considered (e.g. vectorised, tiled). oclRun? | handwritten lines | generated lines
---|---|---
✗ | >1000 | 0
✓ | 0 | >1200
Table 3.4: Without oclRun, 1K lines of host code are handwritten to launch
each kernel for the Harris corner detection design space exploration case
study.

####  3.3.5 Related work

##### Kernel Execution in Domain-Specific Compilers

Halide and PolyMage support kernel executions, including exploring different
kernel groupings [137, 101], but are domain-specific compilers for image
processing.

##### Kernel Execution in General-Purpose Languages

Many general-purpose languages support executing kernels to offload
computations on GPUs. OCaml has the SPOC library [132], Haskell has Accelerate
[53], C/C++ has SYCL, OpenMP and OpenACC.

Typically, the host code is directly written in such languages, leveraging
runtime libraries without necessarily generating C host code. The Rise
language is not general-purpose, which is why we generate C host code.
Additionally, Rise and the oclRun pattern are designed for rewrite-based
design space exploration, which is not a design goal for the acceleration
libraries mentioned above.

##### Implicit Kernel Execution in High-Level Languages

Many high-level languages support executing kernels to accelerate
computations. The common strategy is to make kernel execution implicit and
automated by heuristics, as in MapCG [138], NOVA [55] SaC [139, 140], X10
[141] and Lime [142]. MATLAB is an exception, as it allows the manual
insertion of compiler directives in imperative code, similar to OpenACC [143,
144]. Both approaches differ from our explicit and functional oclRun pattern,
designed for rewrite-based design-space exploration.

##### Lift

Host code generation support has been independently added to Lift in [145].
However, the approach taken introduces two challenges for rewriting without
discussing them, as the paper focuses on manually written low-level programs.
First, the introduced equivalent to our oclRun pattern (called OclKernel) does
not take a single input expression, but rather a function as well as $n$
explicit arguments. It is unclear how a rewrite system would introduce such a
function and its arguments. By contrast, our design abstracts over the kernel
function, inferring it via lambda lifting, and allowing a simple rewrite rule
definition:

⬇

rule introOclRun(ls: LocalSize, gs: GlobalSize) =

(x : dt)
${\color[rgb]{0.19921875,0.1328125,0.53515625}\definecolor[named]{pgfstrokecolor}{rgb}{0.19921875,0.1328125,0.53515625}\bm{\mapsto}}$
oclRun(ls, gs) x

where dt : data

Second, explicit data transfers (as well as in-place computations) are
introduced. By contrast, our design keeps data transfers implicit at the
functional level for simplicity, and rewriting is not burdened with
introducing efficient and correct transfers.

####  3.3.6 Summary

In this section, we demonstrate how the oclRun pattern is added to Rise to
enable explicit kernel execution. To the best of our knowledge, Shine is the
first domain-extensible compiler that makes kernel execution controllable by
rewriting high-level functional programs. Different kernel decompositions of a
high-level Rise program can be explored via rewriting, producing a single low-
level Rise program representing both host-side and device-side computations.
Shine is modified to generate imperative code for multiple OpenCL kernels, as
well as the necessary C host code to launch them. Using oclRun in a relatively
simple Harris corner detection case study enables replacing 1K lines of
previously handwritten host code with 1.2K lines of automatically generated
code (Table 3.4). The kernel execution feature saves performance engineers
from writing boilerplate code and provides a stable interface for the
generated code, which is not the case if Shine generates a single OpenCL
kernel without host code.

###  3.4 Explicit Storage Folding for Shine

Stencil computations are used in a wide range of domains. A stencil
computation processes each element of an N-dimensional input array by
accessing a fixed neighborhood pattern (the stencil) to produce an output
element. When optimising stencils, it is important to leverage spatio-temporal
locality [23]. For example, to compute the next output of a 1D stencil during
iteration $i$ of sequential execution, only the last $m$ stencil inputs are
required. Instead of naively storing all stencil inputs in a temporary storage
$T$, storage folding optimisations only store the last $m$ inputs, improving
memory usage. We consider two storage folding optimisations:

  * •

_Circular buffering_ stores the last $m$ temporary results in memory, using
modulo indexing: $T[j]$ is stored in $M[j\mod m]$ with $j\in[i;i+m-1]$.

  * •

_Register rotation_ stores the last $m$ temporary results in registers,
rotating them between computation iterations: $T[i],\ldots,T[i+m-1]$ is stored
in registers $r_{0},\ldots,r_{m-1}$.

Circular buffering is supported by most image processing specific compilers,
including Halide. However, register rotation is not supported by Halide
schedules.444https://github.com/halide/Halide/issues/2905

This section enables explicit storage folding in Rise by adding the
circularBuffer and rotateValues patterns, along with 2 additional ones. As
with kernel execution in the last section, introducing these patterns makes
storage folding controllable and explorable by rewriting in the domain-
extensible Shine compiler. Chapter 4 will introduce storage folding patterns
during rewriting to generate faster code than Halide on a case study.

####  3.4.1 Adding Rise Patterns for Storage Folding

To explicitly encode storage folding decisions in functional programs, we
introduce two new functional patterns, called circularBuffer and rotateValues,
and with following types:

⬇

circularBuffer : (a: addr) $\to$ (alloc m: nat) $\to$

(s $\to$ t) $\to$ (n + m - 1).s $\to$ n.m.t

rotateValues : (a: addr) $\to$ (m: nat) $\to$

(t $\to$ t) $\to$ (n + m - 1).t $\to$ n.m.t

Both patterns have similar types, as they derive from slide m 1 : (n + m -
1).t $\to$ n.m.t, a sliding window pattern that was introduced to support
stencils in Lift [16]. Both patterns take an address space a as parameter,
controlling where temporary memory will be allocated.

circularBuffer takes an additional alloc parameter, allowing to allocate
memory for more than $m$ values (e.g. using a power of two makes modulo
indexing cheaper). It takes an s $\to$ t function specifying how to load
values into the circular buffer, potentially performing computations.

rotateValues takes a t $\to$ t function specifying how to write values into
memory. This function must be an identity function as it will be used to
rotate values between iterations. We talk about value rotation instead of
register rotation because the use of registers is not guaranteed. Regular
variables are used in the generated OpenCL code, and the OpenCL compiler may
decide to allocate registers for them if the address space is private.

Both patterns implement a finite _stream_ , where a sequence of values is made
available one at a time, sequentially. We add two other stream based
functional patterns, with following types:

⬇

iterateStream : (s $\to$ t) $\to$ n.s $\to$ n.t

mapStream : (s $\to$ t) $\to$ n.s $\to$ n.t

Like mapSeq, both patterns represent concrete implementations of map and have
exactly the same type. However, while mapSeq operates on arrays, mapStream
operates on streams, and iterateStream transforms streams into arrays. The
lack of stream type in Rise is discussed in Section 3.4.3.

Two simple examples using the new storage folding patterns in Rise are shown
in Figure 3.7. Both examples compute a simple stencil summing windows of 3
input values, corresponding to the high-level program slide 3 1
$\vartriangleright$ map (reduce + 0). The code that we want to generate from
the two example Rise programs is depicted in Figure 3.8.

⬇

1iterateStream (reduceSeq + 0)

⬇

1iterateStream (reduceSeq + 0)

Figure 3.7: Example Rise programs using storage folding patterns.

⬇

1int tmp[3];

2for (int i = 0; i < 2; i += 1) {

3 tmp[i] = input[i];

4}

5

6for (int i = 0; i < n-2; i += 1) {

7 tmp[2] = input[2 + i];

8

9 int acc = 0;

10 for (int j = 0; j < 3; j += 1) {

11 acc += tmp[j];

12 }

13

14 output[i] = acc;

15 tmp[0] = tmp[1];

16 tmp[1] = tmp[2];

17}

⬇

1int tmp[3];

2for (int i = 0; i < 2; i += 1) {

3 tmp[i] = input[i];

4}

5

6for (int i = 0; i < n-2; i += 1) {

7 tmp[(2 + i) % 3] = input[2 + i];

8

9 int acc = 0;

10 for (int j = 0; j < 3; j += 1) {

11 acc += tmp[(i + j) % 3];

12 }

13

14 output[i] = acc;

15}

Figure 3.8: Example OpenCL programs corresponding to the Rise programs of
Figure 3.7.

####  3.4.2 Adding a New Stream Translation to DPIA

To generate code for storage folding patterns, we extend DPIA with a new
translation from functional to imperative. The _stream translation_
$\mathcal{S}$ produces a command which reads from a functional read expression
e and calls a stream continuation function sk to continue the translation as
required (LABEL:fig:stream-trans-ty). The stream continuation function is
called by providing a function of type (i: nat) $\to$ (exp[dt, read] $\to$
comm) $\to$ comm, that should stream the $i$-th value by calling the inner
continuation function of type exp[dt, read] $\to$ comm.

⬇

$\mathcal{A}$(e: exp[dt, write], out: acc[dt]): comm

$\mathcal{C}$(e: exp[dt, read], k: exp[dt, read] $\to$ comm): comm

$\mathcal{S}$(e: exp[n.dt, read],

sk: ((i: nat) $\to$ (exp[dt, read] $\to$ comm) $\to$ comm) $\to$ comm): comm

Listing 27: Types of translations, including the new stream translation

The acceptor translation of iterateStream triggers a stream translation
(LABEL:fig:iterate-stream-trans). This makes sense because the pattern writes
to an array output by reading from a stream input. A selection of stream
translations is shown in LABEL:fig:stream-trans-sel. It excludes the stream
translation for circularBuffer which is similar to the one for rotateValues.

⬇

$\mathcal{A}$(iterateStream(n, dt1, dt2, f, array), out) =

$\mathcal{S}$($\lambda$(next : (i: nat) $\to$ (exp[dt1, read] $\to$ comm)
$\to$ comm).

for(n, $\Lambda$i. next i ($\lambda$(x : exp[dt1, read]).

$\mathcal{A}$(f x, out @ i))))

Listing 28: The acceptor translation of iterateStream calls $\mathcal{S}$

⬇

1$\mathcal{S}$(mapStream(n, dt1, dt2, f, input), sk) =

2 $\mathcal{S}$(input, $\lambda$(nextInput : (i: nat) $\to$ (exp[dt1, read]
$\to$ comm) $\to$ comm).

3 sk($\Lambda$i. $\lambda$(k : exp[dt2, read] $\to$ comm).

4 nextInput i ($\lambda$(x : exp[dt1, read]). $\mathcal{C}$(f x, k))

5 ))

6$\mathcal{S}$(rotateValues(a, n, m, dt, wrdt, input), sk) =

7 $\mathcal{S}$(input, $\lambda$(nextInput : (i: nat) $\to$ (exp[dt1, read]
$\to$ comm) $\to$ comm).

8 new(a, m.dt, $\lambda$rs.

9 for(m - 1, $\Lambda$i.

10 nextInput i ($\lambda$(x : exp[dt, read]). $\mathcal{A}$(wrdt x, rs.wr @
i)));

11 sk($\Lambda$i. $\lambda$(k : exp[m.dt, read] $\to$ comm).

12 nextInput (i + m - 1) ($\lambda$(x : exp[dt, read]).

13 $\mathcal{A}$(wrdt x, rs.wr @ (m - 1)));

14 k(rs.rd);

15 for(m - 1, $\Lambda$i. $\mathcal{A}$(wrdt (rs.rd @ (i + 1)), rs.wr @ i))

16 )))

Listing 29: Selection of stream translations

For mapStream, the stream translation of the input is first called in line 2,
giving access to nextInput. In line 3, the stream continuation sk is called to
define how to stream the $i$-the output value. The continuation k should be
called to return this output value. Finally in line 4, nextInput is called to
consume the next input value x, which is used to produce the next output value
f x. $\mathcal{C}$(f x, k) is used to translate f x to imperative.

For rotateValues, the stream translation is a bit more involved. A temporary
buffer is allocated to hold window values in line 8. The first $m-1$ window
values are stored in the buffer in lines 9-10, before sk allows streaming
values. To stream an output value, which is a window of $m$ values:

  * •

The $m$-th window value, or $i+m$-th input value, is loaded in lines 12-13.

  * •

The window of $m$ values is returned in line 14, by calling the continuation
k.

  * •

The values are rotated in the buffer to prepare for the next iteration in line
15.

To write values to the buffer, acceptor translations such as
$\mathcal{A}$(wrdt x, rs.wr @ i) are called, since wrdt x is a functional
program that needs to be translated to imperative.

It is also possible to transform any array into a stream. If there is no
specialised stream translation available, we fallback on a generic one, shown
in LABEL:fig:stream-trans-generic.

⬇

$\mathcal{S}$(e : exp[n.dt, read], sk) =

sk($\Lambda$i. $\lambda$(k : exp[dt, read] $\to$ comm). $\mathcal{C}$(e @ i,
k))

Listing 30: Generic stream translation used as a fallback

Using this newly defined stream translation, the code generated for Figure 3.7
corresponds to the desired code from Figure 3.8, modulo small syntactic
details. It is interesting to note that no new imperative patterns are
required to support storage folding.

####  3.4.3 Limitations of Storage Folding in Shine

The proposed storage folding implementation has two notable limitations.

##### Stream Type

The concept of finite streams is introduced, but is not reflected into the
type system. As a result, it is possible to construct invalid programs that