1

arXiv:1708.03183v1 [] 10 Aug 2017

Automated Tiling of Unstructured Mesh Computations with
Application to Seismological Modelling
FABIO LUPORINI, Imperial College London
MICHAEL LANGE, Imperial College London
CHRISTIAN T. JACOBS, University of Southampton
GERARD J. GORMAN, Imperial College London
J. RAMANUJAM, Louisiana State University
PAUL H. J. KELLY, Imperial College London
Sparse tiling is a technique to fuse loops that access common data, thus increasing data locality. Unlike
traditional loop fusion or blocking, the loops may have different iteration spaces and access shared datasets
through indirect memory accesses, such as A[map[i]] – hence the name “sparse”. One notable example of
such loops arises in discontinuous-Galerkin finite element methods, because of the computation of numerical
integrals over different domains (e.g., cells, facets). The major challenge with sparse tiling is implementation –
not only is it cumbersome to understand and synthesize, but it is also onerous to maintain and generalize, as it
requires a complete rewrite of the bulk of the numerical computation. In this article, we propose an approach
to extend the applicability of sparse tiling based on raising the level of abstraction. Through a sequence
of compiler passes, the mathematical specification of a problem is progressively lowered, and eventually
sparse-tiled C for-loops are generated. Besides automation, we advance the state-of-the-art by introducing:
a revisited, more efficient sparse tiling algorithm; support for distributed-memory parallelism; a range of
fine-grained optimizations for increased run-time performance; implementation in a publicly-available library,
SLOPE; and an in-depth study of the performance impact in Seigen, a real-world elastic wave equation solver
for seismological problems, which shows speed-ups up to 1.28× on a platform consisting of 896 Intel Broadwell
cores.
CCS Concepts: • Computing methodologies → Parallel algorithms; • Mathematics of computing →
Mathematical software performance; • Software and its engineering → Compilers; Domain specific
languages;
Additional Key Words and Phrases: Finite element method, unstructured mesh, compiler, performance optimization, loop fusion, loop tiling, sparse tiling
ACM Reference format:
Fabio Luporini, Michael Lange, Christian T. Jacobs, Gerard J. Gorman, J. Ramanujam, and Paul H. J. Kelly. 2017.
Automated Tiling of Unstructured Mesh Computations with Application to Seismological Modelling. 1, 1,
Article 1 (August 2017), 29 pages.
DOI: 10.1145/nnnnnnn.nnnnnnn
This work was supported by the Engineering and Physical Sciences Research Council through grants EP/I00677X/1,
EP/L000407/1, EP/I012036/1], by the Imperial College London Department of Computing, and by the Imperial College
London Intel Parallel Computing Centre (IPCC). The authors would like to thank the HPC Service Support team at Imperial
College London for their help with the Helen cluster. The authors would also like to thank Gheorghe-Teodor Bercea,
Lawrence Mitchell, and David Ham for their suggestions during the development of this project.
Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee
provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the
full citation on the first page. Copyrights for components of this work owned by others than the author(s) must be honored.
Abstracting with credit is permitted. To copy otherwise, or republish, to post on servers or to redistribute to lists, requires
prior specific permission and/or a fee. Request permissions from permissions@acm.org.
© 2017 Copyright held by the owner/author(s). Publication rights licensed to ACM. XXXX-XXXX/2017/8-ART1 $15.00
DOI: 10.1145/nnnnnnn.nnnnnnn
, Vol. 1, No. 1, Article 1. Publication date: August 2017.

1:2

1

F. Luporini et al.

INTRODUCTION

In many unstructured mesh applications, for example those approximating the solution of partial
differential equations (PDEs) using the finite volume or the finite element method, sequences of
numerical operators accessing common fields need to be evaluated. Usually, these operators are
implemented by iterating over sets of mesh elements and computing a kernel in each element.
In languages such as C or Fortran, the resulting sequence of loops is typically characterized
by heterogeneous iteration spaces and accesses to shared datasets (reads, writes, increments)
through indirect pointers, like A[map[i]]. One notable example of such operators/loops arises
in discontinuous-Galerkin finite element methods, in which numerical integrals over different
domains (e.g., cells, facets) are evaluated; here, A could represent a discrete function, whereas map
could store connectivity information (e.g., from mesh elements to degrees of freedom). In this
article, we devise compiler theory and technology to automate a sophisticated version of sparse
tiling, a technique to maximize data locality when accessing shared fields (like the A and map arrays
in the earlier example), which consists of fusing a sequence of loops by grouping iterations such
that all data dependencies are honored. The goal is to improve the overall application performance
with minimal disruption (none, if possible) to the source code.
Three motivating real-world applications for this work are Hydra, Volna and Seigen. Hydra [30]
is a finite-volume computational fluid dynamics application used at Rolls Royce for the simulation
of next-generation components of jet engines. Volna [11] is a finite-volume computational fluid
dynamics application for the modelling of tsunami waves. Seigen [17] aims to solve the elastic wave
equation using the discontinuous Galerkin finite element method for seismic exploration purposes.
All these applications are characterized by the presence of a time-stepping loop, in which several
loops over the mesh (thirty-three in Hydra, ten in Volna, twenty-five in Seigen) are repeatedly
executed. These loops are characterized by the irregular dependence structure mentioned earlier,
with for example indirect increments in one loop (e.g., A[m[i]] += f(...)) followed by indirect
reads in one of the subsequent loops (e.g., b = g(A[n[j]])). The performance achievable by Seigen
through sparse tiling will extensively be studied in Section 7.
Although our work is general in nature, we are particularly interested in supporting increasingly
sophisticated seismological problems that will be developed on top of Seigen. This has led to the
following strategic decisions:
Automation, but no interest in legacy codes Sparse tiling is an “extreme optimization”.
An implementation in a low level language requires a great deal of effort, as a thoughtful restructuring of the application is necessary. In common with many other low level
transformations, it also makes the source code impenetrable, affecting maintenance and
extensibility. We therefore aim for a fully automated system based on domain-specific
languages (DSLs), which abstracts sparse tiling through a simple interface (i.e., a single
construct to define a scope of fusible loops) and a tiny set of parameters for performance
tuning (e.g., the tile size). We are not interested in automating sparse tiling in legacy codes,
in which the key computational aspects (e.g., mesh iteration, distributed-memory parallelism) are usually hidden for software modularity, thus making such a transformation
almost impossible.
Unstructured meshes require mixed static/dynamic analysis Unstructured meshes are
often used to discretize the computational domain, since they allow for an accurate representation of complex geometries. Their connectivity is stored by means of adjacency lists
(or equivalent data structure), which leads to indirect memory accesses within the loop
, Vol. 1, No. 1, Article 1. Publication date: August 2017.

Automated Tiling of Unstructured Mesh Computations

1:3

nests. Indirections break static analysis, thus making purely compiler-based approaches
insufficient. Runtime data dependence analysis is essential for sparse tiling, so integration
of compiler and run-time tracking algorithms becomes necessary.
Realistic datasets not fitting in a single node Real-world simulations often operate on
terabytes of data, hence execution on multi-node systems is often required. We have
extended the original sparse tiling algorithm to enable distributed-memory parallelism.
Sparse tiling does not change the semantics of a numerical method – only the order in which
some iterations are executed. Therefore, if most sections of a PDE solver suffer from computational
boundedness and standard optimizations such as vectorization have already been applied, then
sparse tiling, which targets memory-boundedness, will only provide marginal benefits (if any).
Likewise, if a global reduction is present in between two loops, then there is no way for sparse
tiling to be applied, unless the numerical method itself is rethought. This is regardless of whether
the reduction is explicit (e.g., the first loop updates a global variable that is read by the second loop)
or implicit (i.e., within an external function, as occurs for example in most implicit finite element
solvers). These are probably the two greatest limitations of the technique; otherwise, sparse tiling
may provide substantial performance benefits.
The rest of the article is structured as follows: in Section 2 we present the abstraction on
which sparse tiling relies. We then show, in Section 3, examples of how the algorithm works on
shared- and distributed-memory systems. This is followed by the formalization of the algorithms
(Sections 4, 5) and the implementation of the compiler that automates sparse tiling (Section 6). The
experimentation is described in Section 7. A discussion on the limitations of the algorithms and the
future work that we expect to carry out in the years to come conclude the article.
2

THE LOOP CHAIN ABSTRACTION FOR UNSTRUCTURED MESH APPLICATIONS

The loop chain is an abstraction introduced in [21]. Informally, a loop chain is a sequence of loops
with no global synchronization points, with attached some extra information to enable run-time
data dependence analysis – necessary since indirect memory accesses inhibit common static approaches to loop optimization. The idea is to replace static with dynamic analysis, exploiting the
information carried by a loop chain. Loop chains must somehow be added to or automatically
derived (e.g., exploiting a DSL) from the input code. A loop chain will then be used by an inspector/executor scheme [32]. The inspector is an algorithm performing data dependence analysis using
the information carried by the loop chain, which eventually produces a sparse tiling schedule. This
schedule is used by the executor, a piece of code semantically equivalent to the original sequence of
loops (i.e., computing the same result) executing the various loop iterations in a different order.
Before diving into the description of the loop chain abstraction, it is worth observing two aspects.
• The inspection phase introduces an overhead. In many scientific computations, the data
dependence pattern is static – or, more informally, “the topology does not change over
time”. This means that the inspection cost may be amortized over multiple iterations of the
executor. If instead the mesh changes over time (e.g., in case of adaptive mesh refinement),
a new inspection must be performed.
• To adopt sparse tiling in a code there are two options. One possibility is to provide a library
and leave the application specialists with the burden of carrying out the implementation
(re-implementation in case of legacy code). A more promising alternative consists of raising
the level of abstraction: programs can be written using a DSL; loop chain, inspector, and
executor can then be automatically derived at the level of the intermediate representation.
As we shall see in Section 6, the tools developed in this article enable both approaches,
though our primary interest is in the automated approach (i.e., via DSLs).
, Vol. 1, No. 1, Article 1. Publication date: August 2017.

Process 0
core

non-exec
exec
owned

F. Luporini et al.
owned
exec
non-exec

1:4
Process 1
core

Fig. 1. Two partitions of a mesh distributed to two neighboring processes. The core region includes all
iterations that can be processed without reading halo data. The owned iterations can be processed only by
reading halo data. Exec halo is the set of iterations that are redundantly executed because they indirectly
access owned iterations. Finally, the non-exec halo region includes off-process iterations that only need to
be read to compute the exec halo iterations. Matching colors across the two processes represent identical
subsets of iterations in the non-partitioned mesh. The image was inspired by an example in [27].

These points will be further elaborated in later sections.
The loop chain abstraction was originally defined as follows:
• A loop chain L consists of n loops, L 0 , L 1 , ..., Ln−1 . There are no global synchronization
points in between the loops. Although there may be dependencies between successive
loops in the chain, the execution order of a loop’s iterations does not influence the result.
• D is a set of disjoint m data spaces, D 0 , D 1 , ..., Dm−1 . Each loop accesses (reads from, writes
to) a subset of these data spaces. An access can be either direct (e.g., A[i]) or indirect (e.g.,
A[map(i)]).
• R Ll →Dd (i) and WLl →Dd (i) are access relations for a loop Ll over a data space Dd ∈ D. They
indicate which locations in the data space Dd an iteration i ∈ Ll reads from and writes to,
respectively. A loop chain must provide all access relations for all loops. For example, if Ll
writes to the array A as A[B(i)] = f(...), then the loop chain will have to provide an
access relation B Ll →A (i).
We here refine this definition, and specialize it for unstructured mesh applications. This allows
the introduction of new concepts, which will be used to extend the sparse tiling algorithm presented
in [34]. Some terminology and ideas are inspired by the programming model of OP2, a library for
unstructured mesh applications [14] used to implement the already mentioned Hydra code.
• A loop chain L consists of n loops, L 0 , L 1 , ..., Ln−1 . There are no global synchronization
points in between the loops. Although there may be dependencies between successive
loops in the chain, the execution order of a loop’s iterations does not influence the result.
• S is a set of m disjoint sets, S 0 , S 1 , ..., Sm−1 . Sets are used to represent iteration and data
spaces. Possible sets are the topological entities in the mesh or the degrees of freedom
associated with a function.
A set S is logically split into three contiguous regions: core (S c ), boundary (S b ), and
non-exec (S ne ). Given a process P and a set S:
S c : the iterations of S that are owned by P.
S b : the boundary region can be seen as the union of two sub-regions, owned (S owned )
and exec (S ex ec ). As shown in Figure 1, S owned are iterations that belong to P which
, Vol. 1, No. 1, Article 1. Publication date: August 2017.

Automated Tiling of Unstructured Mesh Computations

1:5

are redundantly executed by some other processes; S ex ec are iterations from other processes which are redundantly executed by P. We will see that redundant computation
preserves atomic execution – a property that enables executing tiles without the need
for synchronization.
S ne : these are iterations of other processes that are communicated to P because they
need be read to correctly compute S b .
A set is uniquely identified by a name and the sizes of its three regions. For example, the
notation S = (vertices, C, B, N ) defines the vertices set, which comprises C elements
in the core region (iterations 0 to C − 1), B elements in the boundary region (iterations C to
C + B − 1), and N elements in the non-exec region (iterations C + B to C + B + N − 1).
• The depth is an integer indicating the extent of the boundary region of a set. This constant
is the same for all sets in S.
• M is a set of k maps, M 0 , M 1 , ..., Mk −1 . A map of arity a is a vector-valued function M :
Si → S j0 × S j1 × ... × S ja−1 that connects each element of Si to one or more elements in
S j . For example, if a triangular cell c is connected to three vertices v 0 , v 1 , v 2 , we have
M(c) = [v 0 , v 1 , v 2 ].
• A loop Li over the iteration space S is associated with d descriptors, D 0 , D 1 , ..., Dd −1 . A
descriptor D is a 2-tuple D = <M, mode>. M is either a map from S to some other sets
or the special placeholder ⊥, which indicates that memory accesses are direct to some
data associated with S itself. mode is one of [r , w, i], meaning that a memory access is
respectively of type read, write or increment.
There are a few crucial differences in this refined definition for the unstructured mesh case.
One of them is the presence of sets in place of data spaces. In unstructured mesh applications, a
loop tends to access multiple data spaces associated with the same set. The idea is to rely on sets,
rather than data spaces, to perform data dependence analysis. This can significantly improve the
inspection cost, because typically |S| << |D|. Another crucial difference is the characterization
of sets into the three regions core, boundary and non-exec. As we shall see, this separation is
essential to enable distributed-memory parallelism. The extent of the boundary regions is captured
by the depth of the loop chain. Informally, the depth tells how many extra “strips” of elements
are provided by the neighboring processes. This allows some redundant computation along the
partition boundary and also limits the depth of the loop chain (i.e., how many loops can be fused).
The role of the parameter depth will be clear by the end of Section 5.
3

LOOP CHAIN, INSPECTION AND EXECUTION EXAMPLES

Using the example in Figure 2a, we describe the actions performed by our sparse tiling inspector.
The inspector takes as input the loop chain illustrated in Figure 2b. We show two variants, for
shared-memory and distributed-memory parallelism. The value of the variable mode in Figure 2b
determines the variant to be executed.
3.1

Overview

The inspector starts with partitioning the iteration space of a seed loop, for example L 0 . Partitions
are used to initialize tiles: the iterations of L 0 falling in Pi – or, in other words, the edges in partition
Pi – are assigned to the tile Ti . Figure 3 displays the situation after the initial partitioning of L 0 for
a given input mesh. There are four partitions, two of which (P0 and P3 ) are not connected through
any edge or cell. These four partitions correspond to four tiles, [T0 , T1 , T2 , T3 ], with Pi = Ti .
As detailed in the next two sections, the inspection proceeds by populating Ti with iterations
from L 1 and L 2 . The challenge of this task is guaranteeing that all data dependencies – read after
, Vol. 1, No. 1, Article 1. Publication date: August 2017.

1:6

F. Luporini et al.

for t = 0 to T {
inspector = init_inspector (...) ;
// L 0 : loop over edges
for e = 0 to E {
x = X + e;
tmp_0 = tmp + edges2vertices [e + 0];
tmp_1 = tmp + edges2vertices [e + 1];
kernel1 (x , tmp_0, tmp_1) ;
}

// Three sets , edges , cells , and vertices
E = set ( inspector , " edges " , core_edges ,
boundary_edges , nonexec_edges , ...) ;
C = set ( inspector , " cells " , core_cells ,
boundary_cells , nonexec_cells , ...) ;
V = set ( inspector , " verts " , core_verts ,
boundary_verts , nonexec_verts , ...) ;

// L 1 : loop over cells
for c = 0 to C {
res = R + c;
tmp_0 = tmp + cells2vertices [c + 0];
tmp_1 = tmp + cells2vertices [c + 1];
tmp_2 = tmp + cells2vertices [c + 2];
kernel2 ( res , tmp_0 , tmp_1 , tmp_2 );
}

// Two maps , from edges to vertices and from cells
to vertices
e2vMap = map ( inspector , E , V , edges2verts , ...) ;
c2vMap = map ( inspector , C , V , cells2verts , ...) ;
// The loop chain comprises three loops
// Each loop has some descriptors
loop ( inspector , E , { ⊥ , "r"}, { e2vMap , "i" }) ;
loop ( inspector , C , { ⊥ , "r"}, { c2vMap , "i" }) ;
loop ( inspector , E , { e2vMap , "i" }) ;

// L 2 : loop over edges
for e = 0 to E {
tmp_0 = tmp + edges2vertices [e + 0];
tmp_1 = tmp + edges2vertices [e + 1];
kernel3 ( tmp_0, tmp_1) ;
}

// Now can run the inspector
return inspection ( mode , inspector , tile_size , ...)
;

}

(a) Example sequence of sparse-tilable loops.

(b) Loop chain for the example program.

Fig. 2. On the left, a “toy” program used as running example in Section 3 to illustrate the loop chain abstraction
and show how the sparse tiling algorithms (inspection, execution) work. Note that all parameters passed to the
kernels are pointers. On the right, a code snippet showing the loop chain corresponding to the program on the
left. The syntax is very close to the actual API of SLOPE, the sparse tiling library that we have implemented,
described in Section 6.

P1
P3

P0

P2

Fig. 3. Partitioning of the seed loop. The vertices are illustrated to make the connectivity of the mesh clear,
although they do not belong to any partition yet.

write, write after read, write after write – are honored. The output of the inspector is eventually
passed to the executor. The inspection carries sufficient information for computing sets of tiles in
parallel. Ti is always executed by a single thread/process and the execution is atomic; that is, it does
not require communication with other threads/processes. When executing Ti , first all iterations
from L 0 are executed, then all iterations from L 1 and finally those from L 2 .
, Vol. 1, No. 1, Article 1. Publication date: August 2017.

Automated Tiling of Unstructured Mesh Computations

1:7

C0
V1

V0
C1

C2
V2

(a) A snapshot of the mesh after tiling L 0 .
B

B

B

B

B

B
B
B

B
G

R

B
R

R

R

R

R

R

R
R

R

R
R

G

R

R
R

B

G

R

R

B

B

G

B

B

B

B

R
R

B

B
B

R

R

B

R

R
R

(b) The vertices are written by L 0 , so a projection must
be computed before tiling L 1 . Here, the projection is
represented as colored vertices.

R

R
R
R

(c) A snapshot of the mesh after tiling L 1 .

(d) A snapshot of the mesh after tiling L 2 .

Fig. 4. Four passes of the inspection algorithm for shared-memory parallelism.

3.2

Inspection for Shared-Memory Parallelism

Similarly to OP2, to achieve shared-memory parallelism we use coloring. Two tiles that are given
the same color can be executed in parallel by different threads. Two tiles can have the same color
if they are not connected, because this ensures the absence of race conditions through indirect
memory accesses during parallel execution. In the example we can use three colors: red (R), green
(G), and blue (B). T0 and T3 are not connected, so they are assigned the same color. The colored
tiles are shown in Figure 4a. In the following, with the notation Tic we indicate that the i-th tile has
color c.
To populate [T0G , T1B , T2R , T3G ] with iterations from L 1 and L 2 , we first have to establish a total
ordering for the execution of partitions with different colors. Here, we assume the following order:
green (G), blue (B), and red (R). This implies, for instance, that all iterations assigned to T1B must
be executed before all iterations assigned to T2R . By “all iterations” we mean the iterations from L 0
(determined by the seed partitioning) as well as the iterations that will later be assigned from tiling
L 1 and L 2 . We assign integer positive numbers to colors to reflect their ordering, where a smaller
number means higher execution priority. We can assign, for example, 0 to green, 1 to blue, and 2 to
red.
To schedule the iterations of L 1 to [T0G , T1B , T2R , T3G ], we need to compute a projection for any
write or local reduction performed by L 0 . The projection required by L 0 is a function ϕ : V → T
mapping the vertices in V – as indirectly incremented during the execution of L 0 , see Figure 2a –
to a tile Tic ∈ T. Consider the vertex v 0 in Figure 4b. v 0 has 7 incident edges, 2 of which belong to
T0G , while the remaining 5 to T1B . Since we established that G ≺ B, v 0 can only be read after T1B has
finished executing the iterations from L 0 (i.e., the 5 incident blue edges). We express this condition
by setting ϕ(v 0 ) = T1B . Observe that we can compute ϕ by iterating over V and, for each vertex,
applying the maximum function (MAX) to the color of the adjacent edges.
We now use ϕ to schedule L 1 , a loop over cells, to the tiles. Consider again v 0 and the adjacent cells [c 0 , c 1 , c 2 ] in Figure 4b. These three cells have in common the fact that they are
adjacent to both green and blue vertices. For c 1 , and similarly for the other cells, we compute
MAX(ϕ(v 0 ), ϕ(v 1 ), ϕ(v 2 )) = MAX(B, G, G) = B = 1. This establishes that c 1 must be assigned
, Vol. 1, No. 1, Article 1. Publication date: August 2017.

1:8

F. Luporini et al.

to T1B , because otherwise (c 1 assigned instead to T0G ) a read to v 0 would occur before the last
increment from T1B took place. Indeed, we recall that the execution order, for correctness, must
be “all iterations from [L 0 , L 1 , L 2 ] in the green tiles before all iterations from [L 0 , L 1 , L 2 ] in the blue
tiles”. The scheduling of L 1 to tiles is displayed in Figure 4c.
To schedule L 2 to [T0G , T1B , T2R , T3G ] we employ a similar process. Vertices are again written by
L 1 , so a new projection over V will be necessary. Figure 4d shows the output of this last phase.

B

B

B

B

B

B

B
B

G
G

G

G

G
G

(a) After tiling L 0

R
B

G

G

R
R

G
G

G

G

G

G

R

G
G

R

R
R

G

G

G
R

G

G

G
G

G

G

G

G

R

R

R

R
R

R

(b) After tiling L 1

G
R
R

(c) After tiling L 2

Fig. 5. Tiling the program in Figure 2a for shared-memory parallelism can lead to conflicts. Here, the two
green tiles eventually become adjacent, creating race conditions.

Conflicting Colors. It is worth noting how T2R “consumed” the frontier elements of all other tiles
every time a new loop was scheduled. Tiling a loop chain consisting of k loops has the effect of
expanding the frontier of a tile of at most k vertices. With this in mind, we re-inspect the loop
chain of the running example, although this time employing a different execution order – blue (B),
red (R), and green (G) – and a different seed partitioning. Figure 5 shows that, by applying the same
procedure described in this section, T0G and T3G will eventually become adjacent. This violates the
precondition that tiles can be given the same color, and thus run in parallel, as long as they are not
adjacent. Race conditions during the execution of iterations belonging to L 2 are now possible. This
problem will be solved in Section 5.
3.3

Inspection for Distributed-Memory Parallelism

In the case of distributed-memory parallelism, the mesh is partitioned and distributed to a set of
processes. Neighboring processes typically exchange (MPI) messages before executing a loop L j .
A message includes all “dirty” dataset values required by L j modified by any Lk , with Lk ≺ L j . In
the running example, L 0 writes to vertices, so a subset of values associated with border vertices
must be communicated prior to the execution of L 1 . To apply sparse tiling, the idea is to push all
communications to the beginning of the loop chain: as we shall see, this increases the amount
of data to be communicated at a time, but also reduces the number of synchronizations (only 1
synchronization between each pair of neighboring processes per loop chain execution).
From Section 2 it is known that, in a loop chain, a set is logically split into three regions, core,
boundary, and non-exec. The boundary tiles, which originate from the seed partitioning of the
boundary region, will include all iterations that cannot be executed until the communications have
terminated. The procedure described for shared-memory parallelism – now performed individually
by each process on a partition of the input mesh – is modified as follows:
(1) The core region of the seed loop L 0 is partitioned into tiles. Unless aiming for a mixed
distributed/shared-memory scheme, there is no need to assign identical colors to unconnected tiles, as a process will execute its own tiles sequentially. Colors are assigned
increasingly, with Ti given color i. As long as tiles with contiguous ID are also physically
, Vol. 1, No. 1, Article 1. Publication date: August 2017.

Automated Tiling of Unstructured Mesh Computations

1:9

contiguous in the mesh, this assignment retains memory access locality when “jumping”
from executing Ti to Ti+1 .
(2) The same process is applied to the boundary region. Thus, a situation in which a tile includes
iterations from both the core and the boundary regions is prevented by construction. Further,
all tiles within the boundary region are assigned colors higher than those used for the core
tiles. This constrains the execution order: no boundary tiles will be executed until all core
tiles are computed.
(3) We map the whole non-exec region of L 0 to a single special tile, Tne . This tile has the highest
color and will actually never be executed.
boundary
Process 0

core

owned

boundary

non-exec
exec

exec

non-exec

owned

core

Process 1

Fig. 6. A snapshot of the two mesh partitions on Process 0 and Process 1 after inspecting the seed loop
L 0 for distributed-memory parallelism. On each process, there are five tiles in total: two in the core region
(green and violet), two in the boundary region (red and light blue), and Tne . The boundary tiles can safely
cross the owned and exec sub-regions (i.e., the private local iterations and the iterations to be redundantly
computed, respectively). However, no tile can include iterations from both the core and the boundary regions.

boundary
Process 0

core

owned

non-exec
exec

non-exec

boundary
exec

owned

core

Process 1

Fig. 7. A snapshot of the two mesh partitions on Process 0 and Process 1 at the end of the inspection
for distributed-memory parallelism. Tne expands over the boundary region, which minimizes the amount of
redundant computation to be performed. At the end of the execution phase, the orange edges will contain
“dirty values”, but correctness is not affected as the exec region only includes off-process data. The boundary
tiles expand over the core region: this is essential for correctness since none of the red and blue entities from
[L 0 , L 1 , L 2 ] can be executed until the MPI communications have terminated.

From this point on, the inspection proceeds as in the case of shared-memory parallelism. The
application of the MAX function when scheduling L 1 and L 2 makes higher color tiles (i.e., those
having lower priority) “expand over” lower color ones.
, Vol. 1, No. 1, Article 1. Publication date: August 2017.

1:10

F. Luporini et al.

In Figure 6, a mesh is partitioned over two processes and a possible seed partitioning and tiling
of L 0 illustrated. We observe that the two boundary tiles (the red and light blue ones) will expand
over the core tiles as L 1 and L 2 are tiled, which eventually results in the scheduling illustrated in
Figure 7. Roughly speaking, if a loop chain consists of n loops and, on each process, n − 1 extra
layers of iterations are provided (the exec regions in Figure 6), then all boundary tiles are correctly
computed.
The schedule produced by the inspector is subsequently used by the executor. On each process, the
executor starts with triggering the MPI communications required for the computation of boundary
tiles. All core tiles are then computed, since no data from the boundary region is necessary. Hence,
computation is overlapped with communication. As all core tiles are computed and the MPI
communications terminated, the boundary tiles can finally be computed.
Efficiency considerations. The underlying hypothesis is that the increase in data locality will
outweigh the overhead induced by the redundant computation and by the bigger volume of data
exchanged. This is motivated by several facts: (i) the loops being memory-bound; (ii) the core region
being much larger than the boundary region; (iii) the amount of redundant computation being
minimized through the special tile Tne , which progressively expands over the boundary region,
thus avoiding unnecessary calculations.
4

DATA DEPENDENCY ANALYSIS

The loop chain abstraction, described in Section 2, provides the information to construct an
inspector capable of analyzing data dependencies and thus build a legal sparse tiling schedule.
The dependencies between two different loops may be of type flow (read after write), anti (write
after read), or output (write after write). Further, as elaborated below, there may be “reduction
dependencies” between iterations of the same loop. In the following, assume that loop L x , having
iteration space S x , always comes before loop Ly , having iteration space Sy , in the loop chain. Let us
identify a descriptor of a loop L with mmode
S i →S j : this simply indicates that the loop L i has iteration
space Si and uses a map m to write/read/increment elements (respectively, mode ∈ {w, r , i}) in the
space S j .
The flow dependencies can then be enumerated by considering pairs of points (i® and j®) in the
iteration spaces of the two loops L x and Ly :
r
®
®
{i® → j® | i® ∈ S x ∧ j® ∈ Sy ∧ mw
S →S (i ) ∩ m S →S (j ) , ∅}.
x

z

y

z

Anti and output dependencies are defined in a similar way. The anti dependencies for all pairs of
loops L x and Ly are:
®
{i® → j® | i® ∈ S x ∧ j® ∈ Sy ∧ mrS →S (i®) ∩ mw
S →S (j ) , ∅}.
x

z

y

z

While the output dependencies between loops L x and Ly are:
w
®
{i® → j® | i® ∈ S x ∧ j® ∈ Sy ∧ mw
S →S (i ) ∩ m S
x

z

y →S z

(j®) , ∅}.

In essence, there is a storage-related data dependence between two iterations from different loops
(and therefore between the tiles they are placed in) when one of those iterations writes to a data
element and the other iteration reads from or writes to the same data element.
There are local reductions, or “reduction dependencies” between two or more iterations of the
same loop when those iterations “increment” the same location(s); that is, when they read, modify
with a commutative and associative operator, and write to the same location(s). The reduction
dependencies in L x are:
{i® → j® | i® ∈ S x ∧ j® ∈ S x ∧ miS →S (i®) ∩ miS →S (j®) , ∅}.
x

, Vol. 1, No. 1, Article 1. Publication date: August 2017.

z

x

z

Automated Tiling of Unstructured Mesh Computations

1:11

A reduction dependency between two iterations within the same loop indicates that those two
iterations must be executed atomically with respect to each other.
As seen in the example in Section 3, our inspector algorithm handles data dependencies, including
those between non-adjacent loops, by tracking projections. In the next section we formalize how
projections are constructed and used.
5

ALGORITHMS

The pseudo-code for the sparse tiling inspector is shown in Algorithm 1. Given a loop chain and a
“seed” tile size, the algorithm produces a schedule suitable for mixed distributed/shared-memory
parallelism. This schedule – in essence, a set of populated tiles – is used by the executor to perform
the sparse-tiled computation. The executor pseudo-code is displayed in Algorithm 2. In the next
two sections, we elaborate on the main steps of these algorithms. The notation is summarized in
Table 1.
ALGORITHM 1: The inspection algorithm

ALGORITHM 2: The executor algorithm

Input: The loop chain L = [L 0, L 1, ..., L n−1 ], a tile size t s
Output: A set of tiles T, populated with iterations from L

Input: A set of tiles T
Result: Execute the loop chain

seed ← 0;
Φ ← ∅;
C ←⊥;
σs e ed , T ← partition(S s e ed , t s );
seed_map ← find_map(S s e ed , L);
conflicts ← false;

Tc , Tb ← group_tiles_by_region(T);

do
color(T, seed_map);
for j = 1 to n − 1 do
project(L j−1 , σ j−1 , Φ, C );
σ j ← tile(L j , Φ);
assign(σ j , T);
end
if has_conflicts(C) then
conflicts ← true;
add_fake_connection(seed_map, C);
end
while conflicts;
compute_local_maps(T);
return T

start_MPI_comm();
foreach color do
foreach T ∈ Tc s.t. T.color == color do
execute_tile(T);
end
end
end_MPI_comm();
foreach color do
foreach T ∈ Tb s.t. T.color == color do
execute_tile(T);
end
end

ALGORITHM 3: Projection of a tiled loop

ALGORITHM 4: Building a tiling function

Input: A loop L j , a tiling function σ j , the projections set Φ, the conflicts
matrix C
Result: Update Φ and C

Input: A loop L j , the projections set Φ
Output: The tiling function σ j

foreach D ∈ L j .descriptors do
if D.map == ⊥ then
Φ = Φ ∪ σi ;
else
inverse_map ← map_invert(D.map);
S t , S j , values, offset ← inverse_map;
ϕ S t ←⊥;
for e = 0 to S t .size do
for k = offset[e] to offset[e + 1] do
Tl as t = T[values[k ]];
max_color ← MAX(ϕ S t [e].color, Tl as t .color);
if max_color , ϕ S t [e].color then
ϕ S t [e] ← Tl as t ;
end
end
end
update(C, T, ϕ S t );
Φ = Φ ∪ ϕSt ;
end
end

foreach D ∈ L j .descriptors do
if D.map == ⊥ then
σ j ← Φ[S j ];
else
arity ← D.map.arity;
ϕ S ← Φ[D . map .S j ];
for e = 0 to S j .size do
σ j [e] ← T⊥ ;
for k = 0 to arity do
adjacent_tile ← ϕ S [D . map . values[e ∗ arity + k ]];
max_color ← MAX(σ j [e].color,
adjacent_tile.color);
if max_color , σ j [e].color then
σ j [e] ← adjacent_tile;
end
end
end
end
end
return σ j

σ j ←⊥;

, Vol. 1, No. 1, Article 1. Publication date: August 2017.

1:12

F. Luporini et al.
Symbol

L
Lj
Sj
S jc , S jb , S jne
S
D
r, w, i
T
T[i]
ϕS
Φ
σj
ts

Meaning
The loop chain
The j -th loop in L
The iteration space of L j
The core, boundary, and non-exec regions of S j
A generic set in L
A descriptor of a loop
Possible values for D .mode
The set of all tiles
Accessing the i -th tile
A projection ϕ S : S → T
The set of all available projections
A tiling function σ j : S j → T for L j
seed tile size

Table 1. Summary of the notation used throughout the section.

5.1

Inspector

Choice of the seed loop. The seed loop Lseed is used to initialize the tiles. Theoretically, any
loop in the chain can be chosen as seed. Supporting distributed-memory parallelism, however, is
cumbersome if Lseed , L 0 . This is because more general schemes for partitioning and coloring
would be needed to ensure that no iterations in any S bj are assigned to a core tile. A limitation of
our inspector algorithm in the case of distributed-memory parallelism is that it must be Lseed = L 0 .
In the special case in which there is no need to distinguish between core and boundary tiles
(because a program is executed on a single shared-memory system), Lseed can be chosen arbitrarily.
If we however pick Lseed in the middle of the loop chain (L 0 ≺ ... ≺ Lseed ≺ ...), a mechanism
for constructing tiles in the reverse direction (“backwards”), from Lseed towards L 0 , is necessary.
In [34], we propose two “symmetric” algorithms to solve this problem, forward tiling and backward
tiling, with the latter using the MIN function in place of MAX when computing projections. For
ease of exposition, and since in the fundamental case of distributed-memory parallelism we are
imposing Lseed = L 0 , we here neglect this distinction1 .
Tile initialization. Let ts be the user-specified average tile size. The algorithm starts with partitionc
ing Sseed
into m subsets {P0 , P1 , ..., Pm−1 } such that |Pi | = ts (except possibly for Pm−1 ), Pi ∩ P j = ∅,
c
c
and ∪m−1
i=0 Pi = S seed . Among all possible legal partitionings, we choose the one that splits S seed
into blocks of ts contiguous iterations, with P 0 = {0, ..., ts − 1}, P 1 = {ts, ..., 2ts − 1}, and so on. We
b
analogously partition Sseed
into k subsets. We create m + k + 1 tiles, one for each of these partitions
ne
c ,T b , ...,T b
and one extra tile for Sseed . We therefore have T = {T0c , ...,Tm−1
,T ne }.
m
m+k −1 m+k
A tile Ti has four fields, as summarized in Table 2.
• The region is used by the executor to schedule tiles in a given order. This field is set right
c
b
after the partitioning of Lseed , as a tile (by construction) exclusively belongs to Sseed
, Sseed
,
ne
or Sseed .
• The iteration lists contain the iterations in L thatTi will have to execute. There is one iteration
list [Ti ]j for each L j ∈ L. At this stage of the inspection we have [Ti ]seed = [Ti ]0 = Pi ,
whereas still [Ti ]j = ∅ for j = 1, ..., n.
• Local maps may be used for performance optimization by the executor in place of the global
maps provided through the loop chain; this will be discussed in more detail in Section 7
and in the Supplementary Materials.
• The color gives a tile a scheduling priority. If shared-memory parallelism is requested,
adjacent tiles are given different colors (the adjacency relation is determined through the
maps available in L). Otherwise, colors are assigned in increasing order (i.e., Ti is given
1 The algorithm implemented in the library presented in Section 6.1 supports backwards tiling for shared-memory parallelism.

, Vol. 1, No. 1, Article 1. Publication date: August 2017.

Automated Tiling of Unstructured Mesh Computations

1:13

color i). The boundary tiles are always given colors higher than that of core tiles; the
non-exec tile has the highest color. The assignment of colors is carried by the function
color in Algorithm 1.
Field
region
iterations lists
local maps
color

Possible values
core, boundary, non-exec
one list of iterations [Ti ] j for each L j ∈ L
one list of local maps for each L j ∈ L; one local map for each map used in L j
an integer representing the execution priority

Table 2. The tile data structure.

Populating tiles by tracking data dependencies. To schedule a loop to tiles we use projections. A
projection is a function ϕ S : S → T. Initially, the projections set Φ is empty. Each time a loop is tiled,
Φ may be added some new projections or old projections may be updated. Φ, and consequently the
tiling functions for all loops in L, are derived incrementally (within the loop at line 1 in Algorithm 1)
starting from σseed : Sseed → T, the tiling function of Lseed . In the following, we discuss in detail
how projections and tiling functions are constructed.
Deriving a projection from a tiling function. Algorithm 3 takes as input (the descriptors of) L j
and its tiling function σ j : S j → T to update Φ. The algorithm also updates the conflicts matrix
C ∈ Nm×m , which indicates whether two tiles having the same color will become adjacent once
L j+1 is tiled.
A projection tells what tile a set element logically belongs to at a given point of the inspection. A
new projection ϕ S is needed if the elements of S are written by a loop. Let us consider the non-trivial
case in which writes or increments occur indirectly through a map M : S j → S t0 × S t1 × ... × S ta−1 .
To compute ϕ St , we first determine the inverse map (an example is shown in Figure 8). Then, we
iterate over all elements of S t and, for each e ∈ S t , we determine the last tile that writes to e, say
Tl ast . This is accomplished by applying the MAX function over the color of the tiles accessing e.
We finally simply set ϕ St [e] = Tl ast .
9
map:

3

7

9

1

1

3

inverse
map:
7

1

1

1

oﬀset [3]

oﬀset [7]

oﬀset [9]

Fig. 8. Representation of an inverse map. The original map shows that the triangular cell 1 is adjacent to
three vertices, namely 3, 7, and 9. The inverse map associates vertices to cells. Since the mesh is unstructured,
different vertices can be incident to a different number of cells. The array offset determines the distance
between two consecutive vertices in the inverse map. For instance, all entries in the inverse map between
offset[3] and offset[4] are cells incident to vertex 3 (in the illustrated case, the first of these is cell 1).

Deriving a tiling function from the available projections. Using Φ, we compute σ j as described in
Algorithm 4. The algorithm is similar to the projection of a tiled loop, with the main difference
being that now we use Φ to schedule iterations correctly. Finally, σ j is inverted and the iterations
are added to the corresponding iteration lists [Ti ]j , for all Ti ∈ T.
, Vol. 1, No. 1, Article 1. Publication date: August 2017.

1:14

F. Luporini et al.

Detection of conflicts. If C indicates the presence of at least one conflict, say between Ti 1 and Ti 2 ,
we add a “fake connection” between these two tiles and loop back to the coloring stage. Ti 1 and Ti 2
are now connected, so they will be assigned different colors.
On the history of the algorithm. The first algorithm for generalized sparse tiling inspection was
introduced in [34]. In this section, a new, enhanced version of that algorithm has been presented.
In essence, the major differences are: (i) support for distributed-memory parallelism; (ii) use of
coloring instead of a task graph for tile scheduling; (iii) speculative inspection with backtracking if
a coloring conflict is detected; (iv) use of sets, instead of datasets, for data dependency analysis; (v)
use of inverse maps for parallelization of the projection and tiling routines; (vi) computation of
local maps. Most of these changes contributed to reduce the inspection cost, as discussed later.
5.2

Executor

The sparse tiling executor is illustrated in Algorithm 2. It consists of four main phases: (i) exchange of
halo regions amongst neighboring processes through non-blocking communications; (ii) execution
of core tiles (in overlap with communication); (iii) wait for the termination of the communications;
(iv) execution of boundary tiles.
As explained in Sections 3.3, a sufficiently deep halo region enables correct computation of
the boundary tiles. Further, tiles are executed atomically, meaning that all iterations in a tile are
computed without ever synchronizing with other processes. The depth of the boundary region,
which affects the amount of off-process data to be redundantly computed, increases with the number
n of loops to be fused. In the example in Figure 6, there are n = 3 loops, and three “strips” of extra
vertices are necessary for correctly computing the fused loops without tile-to-tile synchronizations.
We recall from Section 2 that the depth of the loop chain indicates the extent of the boundary
region. This parameter imposes a limit to the number of fusible loops. If L includes more loops
than the available boundary region – that is, if n > depth – then L will have to be split into
shorter loop chains, to be fused individually. As we shall see (Section 6), in our inspector/executor
implementation the depth is controlled by the Firedrake’s DMPlex module.
6

IMPLEMENTATION: SLOPE, PYOP2, AND FIREDRAKE

The implementation of automated sparse tiling is distributed over three open-source software
modules.
Firedrake An established framework for the automated solution of PDEs through the finite
element method [28].
PyOP2 A module used by Firedrake to apply numerical kernels over sets of mesh components.
Parallelism is handled at this level.
SLOPE A library for writing inspector/executor schemes, with primary focus on sparse tiling.
PyOP2 uses SLOPE to apply sparse tiling to loop chains.
The SLOPE library is a side-product of the present research. The interplay amongst Firedrake,
PyOP2 and SLOPE is outlined in Figure 9 and discussed in more detail in the following sections.
6.1

SLOPE: a Library for Sparse Tiling Irregular Computations

SLOPE is an open source software that provides an interface to build loop chains and to express
inspector/executor schemes for sparse tiling2 .
The loop chain abstraction implemented by SLOPE has been formalized in Section 2. In essence, a
loop chain comprises some sets (including the separation into core, boundary, and non-exec regions),
2 SLOPE

is available at https://github.com/coneoproject/SLOPE

, Vol. 1, No. 1, Article 1. Publication date: August 2017.

Automated Tiling of Unstructured Mesh Computations
1A

UFL
specification

1B

PyOP2
program

loop _chain

loop _chain

Firedrake

PyOP2

2

1:15
1C
3

p
y
t
h
o
n

C program

Slope

4

Fig. 9. Sparse tiling in the Firedrake-PyOP2-SLOPE framework. There are three ways of sparse tiling a loop
chain: decorating a Firedrake program (1A), decorating a sequence of loops in a PyOP2 program (1B), writing
both the loop chain and the inspector/executor codes explicitly in C through calls to SLOPE (1C). Both (1A)
and (1B) use the loop_chain interface (details in Section 6.2). The kernels generated within a loop_chain
are pre-processed in PyOP2 (2) and forwarded to SLOPE through its Python interface (3). SLOPE now has
access to the loop chain, so it can generate an inspector/executor scheme and return it to PyOP2 (4). The
inspector is compiled and executed. The result, a schedule (i.e., the output of Algorithm 1), is cached and
used as input to the executor. Each time the same loop_chain is encountered in a Firedrake/PyOP2 program,
the corresponding schedule is reused.

maps between sets, and a sequence of loops. Each loop has one or more descriptors specifying
what and how different sets are accessed. The example in Figure 2b illustrates the interface exposed
by SLOPE. It implements the algorithms in Section 5. Further, it provides additional features to
estimate the effectiveness and to verify the correctness of sparse tiling, including generation of
VTK files (suitable for Paraview [2]), to visualize the partitioning of the mesh into colored tiles, as
well as insightful inspection summaries (showing, for example, number and average size of tiles,
total number of colors used, time spent in critical inspection phases).
In the case of shared-memory parallelism, the following sections of code are parallelized through
OpenMP:
Inspection The projection and tiling algorithms; in particular, the loops over set elements in
Algorithm 3 and Algorithm 4).
Execution The computation of same colored tiles; that is, the two loops over tiles in Algorithm 2.
6.2

PyOP2: Lazy Evaluation and Interfaces

PyOP2 is a Python library offering abstractions to model an unstructured mesh – in terms of sets
(e.g. vertices, edges), maps between sets (e.g., a map from edges to vertices to express the mesh
topology), and datasets associating data to each set element (e.g. 3D coordinates to each vertex)
– and applying numerical kernels to sets of entities [28]. In this section, we focus on the three
relevant contributions to PyOP2 made through our work: (i) the interface to identify loop chains;
(ii) the lazy evaluation mechanism that allows loop chains to be built; (iii) the interaction with
SLOPE to automatically build and execute inspector/executor schemes.
To apply sparse tiling to a sequence of loops, the loop_chain interface was added to PyOP2. This
interface, exemplified in Figure 10, is also exposed to the higher layers, for example in Firedrake. In
the listing, the name uniquely identifies a loop chain. Other parameters (most of them optional) are
useful for performance evaluation and performance tuning. Amongst them, the most important are
the tile_size and the fusion_scheme. The tile_size specifies the initial average size for the
seed partitions. The fusion_scheme allows to specify how to break a long sequence of loops into
smaller loop chains, which makes it possible to experiment with a full set of sparse tiling strategies
without having to modify the source code.
, Vol. 1, No. 1, Article 1. Publication date: August 2017.

1:16

F. Luporini et al.
with loop_chain ( name , tile_size , fusion_scheme , ...) :
< some PyOP2 parallel loops are expressed / generated here >

Fig. 10. The loop_chain interface in PyOP2.

PyOP2 exploits lazy evaluation of parallel loops to generate an inspector/executor scheme.
The parallel loops encountered during the program execution – or, analogously, those generated
through Firedrake – are pushed into a queue, instead of being executed immediately. The sequence
of parallel loops in the queue is called the trace. If a dataset f needs to be read, for example because
a user wants to inspect its values or a global linear algebra operation is performed, then the trace is
traversed – from the most recent parallel loop to the oldest one – and a new sub-trace produced.
The sub-trace includes all parallel loops that must be executed to evaluate f correctly. The sub-trace
can then be executed or further pre-processed.
All loops in a trace that were created within a loop_chain scope are sparse tiling candidates. In
detail, the interaction between PyOP2 and SLOPE is as follows:
(1) Figure 10 shows that a loop_chain defines a new scope. As this scope is entered, a stamp
s 1 of the trace is generated. This happens “behind the scenes”, because the loop_chain is
a Python context manager, which can execute pre-specified routines prior and after the
execution of the body. As the loop_chain’s scope is exited, a new stamp s 2 of the trace is
computed. All parallel loops in the trace generated between s 1 and s 2 are placed into a
sub-trace for pre-processing.
(2) The pre-processing consists of two steps: (i) “simple” fusion – consecutive parallel loops
iterating over the same iteration space that do not present indirect data dependencies are
merged; (ii) generation of a loop chain representation for SLOPE.
(3) In (ii), PyOP2 inspects the sequence of loops and translates all relevant data structures (sets,
maps, loops) into a format suitable for the SLOPE’s Python interface. C code implementing
an inspector for the loops in the loop_chain is returned by SLOPE. PyOP2 compiles and
executes this code, which results in an inspection for the loop chain.
(4) A “software cache” mapping loop_chains to inspections is used. This whole process needs
therefore be executed only once for a given loop_chain.
(5) The executor is built in an analogous way to the inspector.
6.3

Firedrake/DMPlex: the S-depth Mechanism for Extended Halo Regions

Firedrake uses DMPlex [23] to handle meshes. DMPlex is responsible for partitioning, distributing
over multiple processes, and locally reordering a mesh. The MPI parallelization is therefore managed
through Firedrake/DMPlex.
During the start-up phase, each MPI process receives a contiguous partition of the original
mesh from DMPlex. The required PyOP2 sets, which can represent either topological components
(e.g., cells, vertices) or function spaces, are created. As intuitively shown in Figure 1, these sets
distinguish between multiple regions: core, owned, exec, and non-exec. Firedrake initializes the
four regions exploiting the information provided by DMPlex.
To support the loop chain abstraction, Firedrake must be able to allocate arbitrarily deep halo
regions. Both Firedrake and DMPlex have been extended to support this feature [19]. A parameter
called S-depth (the name has historical origins, see for instance [7]) regulates the extent of the halo
regions. A value S-depth = n indicates the presence of n strips of off-process data elements in each
set. The default value is S-depth = 1, which enables computation-communication overlap when
, Vol. 1, No. 1, Article 1. Publication date: August 2017.

Automated Tiling of Unstructured Mesh Computations

1:17

executing a single loop at the price of a small amount of redundant computation along partition
boundaries. This is the default execution model in Firedrake.
7
7.1

PERFORMANCE EVALUATION
The Seigen Computational Framework

Seigen is a seismological modelling framework capable of solving the elastic wave equation on
unstructured meshes. Exploiting the well-known velocity-stress formulation [39], the seismic model
is expressible as two first-order linear PDEs, which we refer to as velocity and stress. These
governing equations are discretized in space through the discontinuous-Galerkin finite element
method. The evolution over time is obtained by using a fourth-order explicit leapfrog scheme based
on a truncated Taylor series expansion of the velocity and stress fields. The particular choice of
spatial and temporal discretizations has been shown to be non-dissipative [8]. More details can
be found in [17]. Seigen, which is built on top of Firedrake, is part of OPESCI, an ecosystem of
software for seismic imaging based on automated code generation [15].
Seigen has a set of test cases, which differ in various aspects, such as the initial conditions of the
system and the propagation of waves. However, they are all based upon the same seismological
model; from a computational viewpoint, this means that, in a time step, the same sequence of loops
is executed. In the following, we focus on the explosive_source test case (see the work by [12]
for background details).
7.2

Implementation and Validation

In a time loop iteration, eight linear systems need to be solved, four from velocity and four from
stress. Each solve consists of three macro-steps: assembling a global matrix A; assembling a global
vector b; computing x in the system Ax = b. There are two global “mass” matrices, one for velocity
and one for stress. Both matrices are time invariant, so they are assembled before entering the
time loop, and block-diagonal, as a consequence of the spatial discretization employed (a block
belongs to an element in the mesh). The inverse of a block-diagonal matrix is again block-diagonal
and is determined by computing the inverse of each block. The solution of the linear system Ax = b,
expressible as x = bA−1 , can therefore be evaluated by looping over the mesh and computing a
“small” matrix-vector product in each element, where the matrix is a block in A−1 . Assembling
the global vectors boils down to executing a set of loops over mesh entities, particularly over
cells, interior facets, and exterior facets. Overall, twenty-five loops are executed in a time loop
iteration. Thanks to the hierarchy of “software caches” employed by Firedrake, the translation from
mathematical syntax into loops is only performed once.
Introducing sparse tiling into Seigen was relatively straightforward. Three steps were required:
(i) embedding the time stepping loop in a loop_chain context (see Section 6.2), (ii) propagating
user input relevant for performance tuning, (iii) creating a set of fusion schemes. A fusion scheme
establishes which sub-sequences of loops within a loop_chain will be fused and the respective seed
tile sizes. If no fusion schemes were specified, all of the twenty-five loops would be fused using a
default tile size. As we shall see, operating with a set of small loop chains and heterogeneous tile
sizes is often more effective than fusing long sequences of loops.
Seigen has several mechanisms to validate the correctness of the seismological model and the
test cases. The numerical results of all code versions (with and without tiling) were checked and
compared. Paraview was also used to verify the simulation output. Snapshots of the simulation
output are displayed in Figure 11.
, Vol. 1, No. 1, Article 1. Publication date: August 2017.

1:18

F. Luporini et al.

(a) Right after the explosion.

(b) Wave propagation.

(c) Final snapshot.

Fig. 11. Three phases of the explosive_source test case in Seigen, following an explosion at a point source.

7.3

Computational Analysis of the Loops

We here discuss computational aspects of the twenty-five fusible loops. The following considerations
derive from an analytical study of the data movement in the loop chain, extensive profiling through
the Intel VTune Amplifier tool [16], and roofline models (available in [17]).
Our initial hypothesis was that Seigen would have benefited from sparse tiling. Not only does
data reuse arise within single loops (e.g., by accessing vertex coordinates from adjacent cells), but
also across consecutive loops, through indirect data dependencies. This seemed to make Seigen a
natural fit for sparse tiling. The eight “solver” loops perform matrix-vector multiplications in each
mesh element. It is well established that linear algebra operations of this kind are memory-bound.
Four of these loops arise from velocity, the others from stress. There is significant data reuse
amongst the four velocity loops and amongst the four stress loops, since the same blocks in the
global inverse matrices are accessed. We therefore hypothesized performance gains if these loops
were fused through sparse tiling.
We also observed that because of the particular mesh employed, the exterior facet loops, which
implement the boundary conditions of the variational problems, have negligible cost. However, the
cells and facets loops do have significant cost, and data reuse across them arises. Six loops iterate
over the interior mesh facets to evaluate facet integrals, which ensure the propagation of information
between adjacent cells in discontinuous-Galerkin methods. The operational intensity of these loops
is much lower than that of cell loops, and memory-boundedness is generally expected. Consecutive
facet and cell integral loops share fields, which creates cross-loop data reuse opportunities, thus
strengthening the hypothesis about the potential of sparse tiling in Seigen.
All computational kernels generated in Seigen are optimized through COFFEE [24], a system
used in Firedrake that, in essence, (i) minimizes the operation count by restructuring expressions
and loop nests, and (ii) maximizes auto-vectorization opportunities by applying transformations
such as array padding and by enforcing data alignment.
7.4

Setup and Reproducibility

There are two parameters that we can vary in explosive_source: the polynomial order of the
method, q, and the input mesh. We test the spectrum q ∈ {1, 2, 3, 4}. To test higher polynomial
orders, changes to both the spatial and temporal discretizations would be necessary. For the spatial
discretization, the most obvious choice would be tensor product function spaces, which at the
moment of writing is still unavailable in Firedrake. We use as input a two-dimensional rectangular
domain of fixed size 300×150 tessellated with unstructured triangles (explosive_source only
supports two-dimension domains); to increase the number of elements in the domain, thus allowing
to weak scale, we vary the mesh spacing h.
The generality of the sparse tiling algorithms, the flexibility of the loop_chain interface, and
the S-depth mechanism made it possible to experiment with a variety of fusion schemes without
, Vol. 1, No. 1, Article 1. Publication date: August 2017.

Automated Tiling of Unstructured Mesh Computations

1:19

changes to the source code. Five fusion schemes were devised, based on the following criteria: (i)
amount of data reuse, (ii) amount of redundant computation over the boundary region, (iii) memory
footprint (the larger, the smaller the tile size to fit in cache). The fusion schemes are summarized
in Table 4. The full specification, along with the seed tile size for each sub loop chain, is available
at [46].
The experimentation was conducted on two platforms, whose specification is reported in Table 3.
The two platforms, Erebus (the “development machine”) and Helen (a Broadwell-based architecture
in the Helen cluster [33]) were idle and exclusive access had been obtained when the runs were
performed. Support for shared-memory parallelism is discontinued in Firedrake, so only distributedmemory parallelism with 1 MPI process per physical core was tested. The MPI processes were
pinned to cores. The hyperthreading technology was tried, but found to be generally non-profitable.
The code used for running the experiments was archived with the Zenodo data repository service:
Firedrake [42], PETSc [43], PETSc4py [44], FIAT [41], UFL [49], TSFC [48], PyOP2 [45], COFFEE [40],
SLOPE [47], and Seigen [46].
System
Node
DRAM

Erebus
1x4-core
Intel I7-2600 3.4GHz

Helen
2x14-core
Intel Xeon E5-2680 v4 2.40GHz

16 GB

128 GB (node)

Cache hierarchy

L1=32KB, L2=256KB, L3=8MB

L1=32KB, L2=256KB,
L3/socket=35MB

Compilers

Intel icc 16.0.2

Intel icc 16.0.3

Compiler flags

-O3 -xHost -ip

-O3 -xHost -ip

MPI version

Open MPI 1.6.5

SGI MPT 2.14

Table 3. Systems specification.

Fusion
scheme

Number of
loop chains

fs1

3

fs2

8

fs3

6

fs4

3

fs5

2

Criterion
Fuse costly cells
and facets loops
More aggressive
than fs1
fs2, include all
solver loops
More aggressive
than fs3
velocity, stress

S-depth
2
2
2
3
4

Table 4. Fusion schemes summary.

In the following, for each experiment, we collect three measurements.
Overall completion time – OT Used to compute the maximum application speed-up when
sparse tiling is applied.
Average compute time – ACT Sparse tiling impacts the kernel execution time by increasing
data locality. Communication is also influenced, especially in aggressive fusion schemes:
the rounds of communication decrease, while the data volume exchanged may increase.
ACT isolates the gain due to increased data locality from (i) the communication cost and (ii)
any other action performed in Firedrake (executing Python code) between the invocation
of kernels. This value is averaged across the processes.
Average compute and communication time – ACCT As opposed to ACT, the communication cost is also included in ACCT. By comparing ACCT and ACT, the communication
overhead can be derived.
As we shall see, all of these metrics will be essential for a complete understanding of the sparse
tiling performance.
To collect OT, ACT and ACCT, the following configuration was adopted. All experiments were
executed with “warm cache”; that is, with all kernels retrieved directly from the Firedrake’s software
cache of compiled kernels, so code generation and compilation times are not counted. All of the
non-tiled explosive_source tests were repeated three times. The minimum times are reported
(negligible variance). The cost of global matrix assembly – an operation that takes place before
entering the time loop – is not included in OT. Firedrake needs to be extended to assemble blockdiagonal matrices directly into vectors (an entry in the vector would represent a matrix block).
Currently, this is instead obtained in two steps: first, by assembling into a CSR matrix; then, by
, Vol. 1, No. 1, Article 1. Publication date: August 2017.

1:20

F. Luporini et al.
h
1.0

1.2

q

OT original (s)

OT tiled (s)

ΩOT

Ω ACT

Ω ACCT

1
2
3
4
1
2
3
4

687
1453
3570
7057
419
870
1937
4110

596
1200
2847
5827
365
715
1549
3390

1.15
1.21
1.25
1.21
1.15
1.22
1.25
1.21

1.17
1.25
1.28
1.22
1.17
1.26
1.28
1.23

1.16
1.25
1.27
1.22
1.16
1.26
1.27
1.22

Table 5. OT, ACT and ACCT on Erebus, with 1 MPI process per core.

explicitly copying the diagonal into a vector (a Python operation). The assembly per se never
takes more than 3 seconds, so it was reasonable to exclude this temporary overhead from our
timing. The inspection cost due to sparse tiling is included in OT, and its overhead will be discussed
appropriately in a later section. Extra costs were minimized: no check-pointing, only two I/O
sessions (at the beginning and at the end of the computation), and minimal logging. The time loop
has a fixed duration, while the time step depends on the mesh spacing h to satisfy the CourantFriedrichs-Lewy necessary condition (i.e., CFL limit) for the numerical convergence. In essence,
finer meshes require proportionately smaller time steps to ensure convergence.
We now proceed with discussing the achieved results. Below, a generic instance of explosive_source
optimized with sparse tiling will be referred to as a “tiled version”, otherwise the term “original
version” will be used.
7.5

Single-node Experimentation

Hundreds of single-node experiments were executed on Erebus, which was easily accessible in
exclusive mode and much quicker to use for VTune profiling. The rationale of these experiments
was to assess how sparse tiling impacts the application performance by improving data locality.
We only show parallel runs at maximum capacity (1 MPI process per physical core); the benefits
of sparse tiling in sequential runs tend to be negligible (if any), because (i) the memory access
latency is only marginally affected when a large proportion of bandwidth is available to a single
process, (ii) hardware prefetching is impaired by the small iteration space of tiles, (iii) translation
lookaside buffer (TLB) misses are more frequent due to tile expansion. Points (ii) and (iii) will be
elaborated upon.
Table 5 reports execution times and speed-ups, indicated with the symbol Ω, of the tiled version
over the original version for the three metrics OT, ACT and ACCT. We report the best speed-up
obtained after varying a number of parameters (tile size, fusion scheme and other optimizations
discussed below).
The parameters that we empirically varied to obtain these results were: (i) the fusion scheme,
fsX , X ∈ {1, 2, 3, 4, 5}, see Table 4; (ii) the seed tile size – for each fs and q, up to four different
values chosen to maximize the likelihood of fitting in L2 or L3 cache, were tried3 .
Further, a smaller experimentation varying (i) type of indirection maps (local or global, see
Section 5) and (ii) tile shape (through different mesh partitioning algorithms), as well as introducing
(iii) software prefetching and (iv) extended boundary region (to minimize redundant computation)
led to the conclusion that these optimizations may either improve or worsen the execution time, in
ways that are too difficult to predict beforehand. We therefore decided to exclude these parameters
from the full search space. This simplifies the analysis that will follow and also allowed the
3A

less extensive experimentation with “more adverse” tile sizes showed that: (i) a very small value causes dramatic
slow-downs (up to 8× slower than the original versions); (ii) larger values cause proportionately greater performance drops.
, Vol. 1, No. 1, Article 1. Publication date: August 2017.

fs4

fs5

140

250 320
Tile size factor
[q = 1, h = 1.2]
fs1

fs2

fs3

400

fs4

fs5

140

250 320
Tile size factor

400

1.30
1.25
1.20
1.15
1.10
1.05
1.00

1.30
1.25
1.20
1.15
1.10
1.05
1.00

[q = 2, h = 1.0]
fs1

fs2

fs3

fs4

fs5

70

140 200
Tile size factor
[q = 2, h = 1.2]
fs1

fs2

fs3

300

fs4

fs5

70

140 200
Tile size factor

300

1.30
1.25
1.20
1.15
1.10
1.05
1.00

1.30
1.25
1.20
1.15
1.10
1.05
1.00

[q = 3, h = 1.0]
fs1

fs2

fs3

fs4

fs5

ACT speed-up

fs3

1:21

45

60
Tile size factor
[q = 3, h = 1.2]
fs1

fs2

fs3

75

fs4

fs5

ACT speed-up

fs2

ACT speed-up

fs1

ACT speed-up

1.30
1.25
1.20
1.15
1.10
1.05
1.00

[q = 1, h = 1.0]

ACT speed-up

1.30
1.25
1.20
1.15
1.10
1.05
1.00

ACT speed-up

ACT speed-up

ACT speed-up

Automated Tiling of Unstructured Mesh Computations

45

60
Tile size factor

75

1.30
1.25
1.20
1.15
1.10
1.05
1.00

1.30
1.25
1.20
1.15
1.10
1.05
1.00

[q = 4, h = 1.0]
fs1

20

fs2

fs3

45
Tile size factor
[q = 4, h = 1.2]
fs1

20

fs4

fs2

fs3

45
Tile size factor

fs5

70

fs4

fs5

70

Fig. 12. Search space exploration on Erebus, with h ∈ {1.0, 1.2} and q ∈ {1, 2, 3, 4}. Each plot shows the
average compute time (ACT) speed-up achieved by multiple tiled versions over the original (non-tiled) version.
The seed tile size of a loop chain in an fs, in terms of seed loop iterations, is the product of the “tile size
factor” (x-axis) and a pre-established multiplier (an integer in the range [1, 4]; full list available at [46]).

execution of the whole test suite in less than five days. The interested reader is invited to refer to
the Supplementary Materials attached to this article for more information.
In Figure 12, we summarize the results of the search space exploration. A plot shows the ACT
speed-ups achieved by multiple tiled versions over the non-tiled version for given q and h. In these
single-node experiments, the ACCT trend was basically identical (as one can infer from Table 5),
with variations in speed-up smaller than 1%.
PyOP2 was enhanced with a loop chain analyzer (LCA) capable of estimating the best- and
worst-case tile memory footprint, as well as the percentage of data reuse ideally achievable4 . We use
this tool, as well as Intel VTune, to explain the results shown in Figure 12. We make the following
observations.
• fs, unsurprisingly, is the parameter having largest impact on the ACT. By influencing the
fraction of data reuse convertible into data locality, the amount of redundant computation
and the data volume exchanged, fusion schemes play a fundamental role in sparse tiling.
This makes automation much more than a desired feature: without any changes to the
source code, multiple sparse tiling strategies could be studied and tuned. Automation is one
of our major contributions, and this performance exploration justifies the implementation
effort.
• There is a non-trivial relationship between ACT, q and fs. The aggressive fusion schemes
are more effective with high q – that is, with larger memory footprints – while they tend
to be less efficient, or even deleterious, when q is low. The extreme case is fs5, which fuses
two long sequences of loops (twelve and thirteen loops each). In Figure 12 (Erebus), fs5 is
never a winning choice, although the difference between fs3/fs4 and fs5 decreases as q
grows. If this trend continued with q > 4, then the gain from sparse tiling could become
increasingly larger.
4 It

is part of our future plans to integrate this tool with SLOPE, where the effects of tile expansion can be taken into account
to provide better estimates.
, Vol. 1, No. 1, Article 1. Publication date: August 2017.

1:22

F. Luporini et al.
• A non-aggressive scheme fuses only a few small subsets of loops. As discussed in later
sections, these fusion schemes, despite affording larger tile sizes than the more aggressive
ones (due to the smaller memory footprint), suffer from limited cross-loop data reuse. For
fs1, LCA determines that the percentage of reusable data in the three fused loop chains
decreases from 25% (q = 1) to 13% (q = 4). The drop is exacerbated by the fact that no
reuse can be exploited for the maps. Not only are these ideal values, but also a significant
number of loops are left outside of loop chains. The combination of these factors motivate
the lack of substantial speed-ups. With fs2, a larger proportion of loops are fused and the
amount of shared data increases. The peak ideal reuse in a loop chain reaches 54%, which
translates into better ACTs. A similar growth in data reuse can be appreciated in more
aggressive fusion schemes, with a peak of 61% in one of the fs5’s loop chains. Nevertheless,
the performance of fs5 is usually worse than fs4. As we shall clarify in Section 7.8, this
is mainly due to the excessive memory footprint, which in turn leads to very small tiles.
Speculatively, we tried running a sixth fusion scheme: a single loop chain including all of
the 25 loops in a time iteration. In spite of an ideal data reuse of about 70%, the ACT was
always significantly higher than all other schemes.
• Figure 12 shows the ACT for a “good selection” of tile size candidates. Our approach was as
follows. We took a very small set of problem instances and we tried a large range of seed
tile sizes. Very small tile sizes caused dramatic slow-downs, mostly because of ineffective
hardware prefetching and TLB misses. Tile sizes larger than a fraction of L3 cache also led
to increasingly higher ACTs. If we consider q = 4 in Figure 12, we observe that the ACT of
fs2, fs3, and fs4 grows when the initial number of iterations in a tile is as big as 70. Here,
LCA shows that the tile memory footprint is, in multiple loop chains, higher than 3MB,
with a peak of 6MB in fs3. This exceeds the proportion of L3 cache that a process owns
(on average), which explains the performance drop.

7.6

Multi-node Performance

Weak-scaling experiments were carried out on thirty-two nodes split over two racks in the Helen
cluster at Imperial College London [33, 38]. The node architecture as well as the employed software
stack are summarized in Table 3. The rationale of these experiments was to assess whether the
change in communication pattern and amount of redundant computation caused by sparse tiling
could affect the run-time performance.
For each q in the usual range [1−4], fs3 generally resulted in the best performance improvements,
due to its trade-off between gained data locality and restrained redundant computation (whose
effect obviously worsen as q grows). Figure 13 summarizes the ACCT speed-ups achieved by the
best tiled version (i.e., the one found by empirically varying the tile size, with same tile size factor
as in Figure 12) over the original version. The weak scaling trend is remarkable, with only a small
drop in the case q = 1 when switching from one rack (448 cores) to two racks (896 cores), which
disappears as soon as the per-process workload becomes significant. For instance, with q = 2, each
process already computes over more than 150k degrees of freedom for the velocity and stress fields.
The peak speed-up over the original version, 1.28×, is obtained in the test case q = 3 when running
with 448 processes.
7.7 Negligible Inspection Overhead
As explained in Section 7.4, the inspection cost was included in OT. In this section, we quantify
this overhead in a representative problem instance. In all other problem instances, the overhead
was either significantly smaller than or essentially identical (for reasons discussed below) to that
, Vol. 1, No. 1, Article 1. Publication date: August 2017.

Automated Tiling of Unstructured Mesh Computations
Weak scaling on Helen
q=1 q=2 q=3

1.4

ACCT speed-up (original / tiling)

1:23
q=4

1.2
1.0
0.8
0.6
0.4
0.2
0.0
28
(119K)

56
(240K)

112
(477K)

224
(953K)

Number of processes
(Mesh elements)

448
(1.91M)

896
(3.80M)

Fig. 13. Weak scaling performance of Seigen’s explosive_source on the Helen cluster. The ACCT speed-up
is relative to a single process. Results for q ∈ {1, 2, 3, 4} are displayed. The simulations ran for 1000 timesteps.

reported here. We consider explosive_source on Erebus with h = 1.0, q = 4, and fs5. With this
configuration, the time step was ∆t = 481 · 10−6 (we recall from Section 7.4 that ∆t is a function
of h). Given the simulation duration T = 2.5, in this test case 5198 time steps were performed. A
time step took on average 1.15 seconds. In each time step, twenty-five loops, fused as specified by
fs5, are executed. We know that in fs5 there are two sub loop chains, which respectively consist
of thirteen and twelve loops. To inspect these two loop chains, 1.4 and 1.3 seconds were needed
(average across the four MPI ranks, with negligible variance). Roughly 98% of the inspection time
was due to the projection and tiling functions, while only 0.2% was spent on the tile initialization
phase (see Section 5). These proportions are consistent across other fusion schemes and test cases.
After 200 time steps (less than 4% of the total) the inspection overhead was already close to 1%.
Consequently, the inspection cost, in this test case, was eventually negligible. The inspection cost
increases with the number of fused loops, which motivates the choice of fs5 for this analysis.
7.8 On the Main Performance Limiting Factor
The tile structure, namely its shape and size, is the key factor affecting sparse tiling in Seigen.
The seed loop partitioning and the mesh numbering determine the tile structure. The simplest
way of creating tiles consists of “chunking” the seed iteration space every ts iterations, with ts
being the initial tile size (see Section 5). Hence, the chunk partitioning inherits the original mesh
numbering. In Firedrake, and therefore in Seigen, meshes are renumbered during initialization
applying the reverse Cuthill-McKee (RCM) algorithm. Using chunk partitioning on top of an RCMrenumbered mesh has the effect of producing thin, rectangular tiles, as displayed in Figure 14a.
This dramatically affects tile expansion, as a large proportion of elements will lie on the tile border.
There are potential solutions to this problem. The most promising would consist of using a Hilbert
curve, rather than RCM, to renumber the mesh. This would lead to more regular polygonal tiles
when applying chunk partitioning. Figures 14b and 14c show the tile structure that we would
ideally want in Seigen, from a Hilbert-renumbered mesh produced outside of Firedrake. As later
discussed, introducing support for Hilbert curves in Firedrake is part of our future work.
The memory footprint of a tile grows quite rapidly with the number of fused loops. In particular,
the matrices accessed in the velocity and stress solver loops have considerable size. The larger
the memory footprint, the smaller ts for the tile to fit in some level of cache. Allocating small tiles
, Vol. 1, No. 1, Article 1. Publication date: August 2017.

1:24

F. Luporini et al.

(a) chunk partitioning in Seigen.

(b) Ideal Hilbert partitioning.

(c) Zoomed in Hilbert partitioning.

Fig. 14. Representation of seed iteration space partitioning strategies via SLOPE and Paraview.

has unfortunately multiple implications. First, the proportion of iterations lying on the border grows,
which worsen the tile expansion phenomenon discussed, thus affecting data locality. Secondly, a
small ts impairs hardware prefetching, since the virtual address streams become more irregular.
Finally, using small tiles implies that a proportionately larger number of tiles are needed to “cover”
the iteration space; in the case of shared-memory parallelism, this increases the probability of
coloring conflicts, which result in higher inspection costs.
8

RELATED WORK

The data dependence analysis that we have developed in this article is based on the loop chain
abstraction, which was originally presented in [21]. This abstraction is sufficiently general to
capture data dependencies in programs structured as arbitrary sequences of loops, particularly to
create inspector/executor schemes for many unstructured mesh application. Inspector/executor
strategies were first formalized by [32]. They have been used to exploit data reuse and to expose
shared-memory parallelism in several studies [9, 10, 20, 36].
Sparse tiling is a technique based upon inspection/execution. The term was coined by [36, 37] in
the context of the Gauss-Seidel algorithm and also used in [35] in the Moldyn benchmark. However,
the technique was initially proposed by [10] to parallelize computations over unstructured meshes,
taking the name of unstructured cache blocking. In this work, the mesh was initially partitioned; the
partitioning represented the tiling in the first sweep over the mesh. Tiles would then shrink by one
layer of vertices for each iteration of the loop. This shrinking represented what parts of the mesh
could be accessed in later iterations of the outer loop without communicating with the processes
executing other tiles. The unstructured cache blocking technique also needed to execute a serial
clean-up tile at the end of the computation. [1] also developed an algorithm very similar to sparse
tiling, to parallelize Gauss-Seidel computations. The main difference between [36, 37] and [10] was
that in the former work the tiles fully covered the iteration space, so a sequential clean-up phase at
the end could be avoided. All these approaches were either specific to individual benchmarks or
not capable of scheduling across heterogeneous loops (e.g., one over cells and another over degrees
of freedom). These limitations had been addressed in [34].
The automated code generation technique presented in [29] examines the data affinity among
loops and performs partitioning with the goal of minimizing inter-process communication, while
maintaining load balancing. This technique supports unstructured mesh applications (being based
on an inspector/executor strategy) and targets distributed memory systems, although it does not
exploit the loop chain abstraction and does not introduce any sort of loop reordering transformation.
Automated code generation techniques based on polyhedral compilers have been applied to
structured grid benchmarks or proxy applications [5]. However, there has been very little effort in
providing evidence that these tools can be effective in real-world applications. Time-loop diamond
, Vol. 1, No. 1, Article 1. Publication date: August 2017.

Automated Tiling of Unstructured Mesh Computations

1:25

tiling was applied in [4] to a proxy application, but experimentation was limited to shared-memory
parallelism. In [31], instead, an approach based on raising the level of abstraction, similar to the
one presented in this paper, is described and evaluated. The experimentation is conducted using
realistic stencil codes ported to the OPS library. The main difference with respect to our work is
the focus on structured grids (i.e., different types of numerical methods are targeted).
In structured codes, multiple layers of halo, or “ghost” elements, are often used to reduce
communication [3]. Overlapped tiling exploits the very same idea: trading communication for
redundant computation along the boundary [50]. Several works tackle overlapped tiling within
single regular loop nests (mostly stencil-based computations), for example [6, 22, 25]. Techniques
known as “communication avoiding” [9, 26] also fall in this category. To the best of our knowledge,
overlapped tiling for unstructured mesh applications has only been studied analytically, by [13].
Further, we are not aware of any prior techniques for automation.

9

FUTURE WORK

Using a Hilbert curve numbering will lead to dramatically better tile shapes, thus mitigating the
performance penalties due to tile expansion, TLB misses and hardware prefetching described in
Section 7. This extension is at the top of our future work priorities.
Shared-memory parallelism was not as carefully tested as distributed-memory parallelism. First
of all, we would like to replace the current OpenMP implementation in SLOPE with the MPI Shared
Memory (SHM) model introduced in MPI-3. Not only does a unified programming model provide
significant benefits in terms of maintainability and complexity, but the performance may also be
greater as suggested by recent developments in the PETSc community. Secondly, some extra work
would be required for a fair comparison of this new hybrid MPI+MPI programming model with
and without sparse tiling.
The experimentation was carried out on a number of “conventional” Intel Xeon architectures;
we aim to repeat the experimentation on the Intel’s Knights Landing soon.
Finally, a cost model for automatic derivation of fusion schemes and tile sizes is still missing.

10 CONCLUSIONS
Sparse tiling aims to turn the data reuse in a sequence of loops into data locality. In this article,
three main problems have been addressed: the specialization of sparse tiling to unstructured mesh
applications via a revisited loop chain abstraction, automation via DSLs, effective support for
shared- and distributed-memory parallelism. The major strength of this work lies in the fact that
all algorithmic and technological contributions presented derive from an in-depth study of realistic
application needs. The performance issues we found through Seigen would never have been exposed
by a set of simplistic benchmarks, as often used in the literature. Further experimentation will be
necessary when 3D domains and high-order discretizations will be supported by Seigen. In essence,
the performance experimentation shows systematic speed-ups, in the range 1.10×-1.30×. This is
hopefully improvable by switching to Hilbert curve numberings and by exploiting shared memory
through a suitable paradigm. Finally, our opinion is that sparse tiling is an “extreme” optimization:
at least for unstructured mesh application, it is unlikely that it will lead to speed-ups in the order
of magnitudes. However, through automation via DSLs, and with suitable optimization and tuning,
it may still play a key role in improving the performance of real-world computations.

, Vol. 1, No. 1, Article 1. Publication date: August 2017.

1:26
11
11.1

F. Luporini et al.
SUPPLEMENTARY MATERIALS
Summary of the Optimizations Attempted in Seigen

A number of potential optimizations were attempted when applying sparse tiling to Seigen. The
experimentation with these optimizations was, however, inconclusive. Although performance
improvements were observed in various problem instances, in a significant number of other cases
either a detrimental effect or no benefits were noticed at all. Below we briefly discuss the impact of
four different execution strategies: (i) use of global maps, (ii) variation in tile shape, (iii) software
prefetching, (iv) use of an extended boundary region to minimize redundant computation.
Global and local maps Algorithm 1 computes so called local maps to avoid an extra level
of indirection in the executor. Although no data reuse is available for the local maps (each
fused loop has its own local maps), there might be benefits from improved hardware
prefetching and memory latency. We compared the use of global and local maps (i.e., the
former are normally constructed by Firedrake and provided in the loop chain specification,
the latter are computed by SLOPE), but no definitive conclusion could be drawn, as both
performance improvements and deteriorations within a 5% range were observed.
Tile shape The seed loop partitioning and the mesh numbering determine the tile shape.
The simplest way of creating tiles consists of “chunking” the seed iteration space every
ts iterations, with ts being the initial tile size (see Section 5.1). The chunk partitioning
inherits the original mesh numbering. In Firedrake, and therefore in Seigen, meshes are
renumbered during initialization applying the reverse Cuthill-McKee (RCM) algorithm.
Using chunk partitioning on top of an RCM-renumbered mesh has the effect of producing
thin, rectangular tiles. This affects tile expansion, as a large proportion of elements lie on
the tile border. There are potential solutions to this problem. The most promising would
consist of using a Hilbert curve, rather than RCM, to renumber the mesh. This would lead
to more regular polygonal tiles when applying chunk partitioning. Unfortunately, Hilbert
curves are at the moment not supported in Firedrake. An alternative consists of partitioning
the seed iteration space with a library like METIS [18] before applying RCM. However,
we discovered that this approach too was not exempt from side effects. The main cause
was increased translation lookaside buffer (TLB) misses, which occur whenever the CPU
cannot retrieve the mapping to the physical page corresponding to a given virtual page.
Since the page table has a hierarchical structure, handling a TLB miss usually requires
multiple accesses to memory. Hence, TLB misses are much more costly than cache misses.
Sparse tiling increases the TLB miss/hit ratio because of the fragmented streams of virtual
addresses. This is evident (and more pronounced) when the tile size is small, in which
case a TLB miss is quite likely to occur when jumping to executing a new loop. This
problem is exacerbated by the metis partitioning (in contrast to chunk), which leads to
irregular tile shapes. Here, tile expansion may eventually incorporate iterations living in
completely different virtual pages. VTune experimentation with q = 1 and q = 2 versions
of explosive_source showed that chunk- and metis-based sparse tiling suffer from an
increase in TLB misses of roughly 16% and 35%, respectively. To mitigate this issue, we
explored the possibility of using larger virtual pages through Linux’s Transparent Huge
Pages mechanism, which was enabled to automatically allocate memory in virtual pages
of 2MB (instead of the default 4KB) – as long as the base array addresses were properly
aligned. However, no significant differences were observed, and a deeper investigation is
still necessary.
, Vol. 1, No. 1, Article 1. Publication date: August 2017.

Automated Tiling of Unstructured Mesh Computations

1:27

Software prefetching In a loop, there is usually more than a single stream of memory accesses amenable to hardware prefetching (e.g., accesses to the indirection maps; direct
accesses to data values; indirect accesses to data values if the mesh has a good numbering).
Sparse tiling, unfortunately, impairs hardware prefetching for two reasons: (i) the virtual
addresses streams are considerably shorter; (ii) tile expansion introduces irregularities in
these streams. Software prefetching can be used together with hardware prefetching to
minimize memory stalls. PyOP2 and SLOPE have been extended to emit intrinsics instructions to prefetch the iteration i’s maps and data values while executing the iteration i − k
at distance k. No compelling evidence that this further transformation could systematically
improve the performance was found.
Extended boundary region The special non-exec tile Tne (see Sections 3 and 5) reduces the
amount of redundant computation in long loop chains by expanding over boundary tiles.
There are two ways of creating Tne : either an extra layer of data is added to the boundary
region (e.g., see Figure 6) or during inspection, by searching for mesh boundaries. The
current implementation only supports the first option. Manually deriving Tne would be not
only algorithmically complex, but also potentially very expensive.
REFERENCES
[1] M. F. Adams and J. Demmel. 1999. Parallel Multigrid Solver Algorithms and Implementations for 3D Unstructured
Finite Element Problem. In Proceedings of SC99. Portland, Oregon.
[2] Utkarsh Ayachit. 2015. The ParaView Guide (Full Color Version): A Parallel Visualization Application (paraview 4.3 ed.).
Kitware, Incorporated. http://www.amazon.com/exec/obidos/redirect?tag=citeulike07-20&path=ASIN/1930934300
[3] Federico Bassetti, Kei Davis, and Daniel J. Quinlan. 1998. Optimizing Transformations of Stencil Operations for Parallel
Object-Oriented Scientific Frameworks on Cache-Based Architectures. In Proceedings of the Second International
Symposium on Computing in Object-Oriented Parallel Environments (ISCOPE ’98). Springer-Verlag, London, UK, UK,
107–118. http://dl.acm.org/citation.cfm?id=646894.709707
[4] Uday Bondhugula, Vinayaka Bandishti, Albert Cohen, Guillain Potron, and Nicolas Vasilache. 2014. Tiling and
Optimizing Time-iterated Computations on Periodic Domains. In Proceedings of the 23rd International Conference on
Parallel Architectures and Compilation (PACT ’14). ACM, New York, NY, USA, 39–50. DOI:https://doi.org/10.1145/
2628071.2628106
[5] Uday Bondhugula, Albert Hartono, J. Ramanujam, and P. Sadayappan. 2008. A Practical Automatic Polyhedral
Parallelizer and Locality Optimizer. In Proceedings of the 2008 ACM SIGPLAN Conference on Programming Language
Design and Implementation (PLDI ’08). ACM, New York, NY, USA, 101–113. DOI:https://doi.org/10.1145/1375581.
1375595
[6] Li Chen, Zhao-Qing Zhang, and Xiao-Bing Feng. 2002. Redundant computation partition on distributed-memory
systems. In Algorithms and Architectures for Parallel Processing, 2002. Proceedings. Fifth International Conference on.
252–260. DOI:https://doi.org/10.1109/ICAPP.2002.1173583
[7] A. T. Chronopoulos. 1991. s-Step Iterative Methods for (Non)Symmetric (In)Definite Linear Systems. SIAM J. Numer.
Anal. 28, 6 (1991), 1776–1789. DOI:https://doi.org/10.1137/0728088 arXiv:http://dx.doi.org/10.1137/0728088
[8] Sarah Delcourte, Loula Fezoui, and Nathalie Glinsky-Olivier. 2009. A high-order Discontinuous Galerkin method for
the seismic wave propagation. ESAIM: Proceedings 27 (2009), 70–89. DOI:https://doi.org/10.1051/proc/2009020
[9] James Demmel, Mark Hoemmen, Marghoob Mohiyuddin, and Katherine Yelick. 2008. Avoiding Communication in
Sparse Matrix Computations. In Proceedings of International Parallel and Distributed Processing Symposium (IPDPS).
IEEE Computer Society.
[10] Craig C. Douglas, Jonathan Hu, Markus Kowarschik, Ulrich Rüde, and Christian Weiß. 2000. Cache Optimization for
Structured and Unstructured Grid Multigrid. Electronic Transaction on Numerical Analysis (February 2000), 21–40.
[11] Denys Dutykh, Raphaël Poncet, and Frédéric Dias. 2011. The VOLNA code for the numerical modeling of tsunami
waves: Generation, propagation and inundation. European Journal of Mechanics - B/Fluids 30, 6 (2011), 598 – 615.
Special Issue: Nearshore Hydrodynamics.
[12] W. W. Garvin. 1956. Exact Transient Solution of the Buried Line Source Problem. Proceedings of the Royal Society of
London, Series A 234, 1199 (1956). DOI:https://doi.org/10.1098/rspa.1956.0055
[13] M. B. Giles, G. R. Mudalige, C. Bertolli, P. H. J. Kelly, E. Laszlo, and I. Reguly. 2012. An Analytical Study of Loop
Tiling for a Large-Scale Unstructured Mesh Application. In Proceedings of the 2012 SC Companion: High Performance
, Vol. 1, No. 1, Article 1. Publication date: August 2017.

1:28

[14]

[15]

[16]
[17]
[18]
[19]

[20]

[21]

[22]

[23]

[24]

[25]

[26]

[27]
[28]

[29]

[30]

[31]
[32]
[33]

F. Luporini et al.
Computing, Networking Storage and Analysis (SCC ’12). IEEE Computer Society, Washington, DC, USA, 477–482. DOI:
https://doi.org/10.1109/SC.Companion.2012.68
M. B. Giles, G. R. Mudalige, Z. Sharif, G. Markall, and P. H. J. Kelly. 2011. Performance Analysis of the OP2 Framework
on Many-core Architectures. SIGMETRICS Perform. Eval. Rev. 38, 4 (March 2011), 9–15. DOI:https://doi.org/10.1145/
1964218.1964221
Gerard Gorman, Marcos de Aguiar, David Ham, Felix Herrmann, Christian Jacobs, Paul Kelly, Michael Lange, Chris
Pain, Matthew Piggott, Spencer Sherwin, Felippe Vieira Zacarias, Mike Warner, Fabio Luporini, and Navjot Kukreja.
2015. OPESCI project web page. http://www.opesci.org. (2015).
Intel Corporation. 2016. Intel VTune Performance Analyzer. software.intel.com/en-us/intel-vtune-amplifier-xe. (2016).
Christian T. Jacobs, Michael Lange, Fabio Luporini, and Gerard J. Gorman. 2017. Application of code generation to
high-order seismic modelling with the discontinuous Galerkin finite element method. (2017). In preparation.
George Karypis and Vipin Kumar. 2011. MeTis: Unstructured Graph Partitioning and Sparse Matrix Ordering System,
Version 5.0. http://www.cs.umn.edu/~metis. (2011).
Matthew G. Knepley, Michael Lange, and Gerard Gorman. 2015. Unstructured Overlapping Mesh Distribution in
Parallel. Submitted to ACM Transactions on Mathematical Software (TOMS) (2015). arXiv:1506.06194 http://arxiv.org/
abs/1506.06194
Christopher D. Krieger and Michelle Mills Strout. 2012. Executing Optimized Irregular Applications Using Task Graphs
Within Existing Parallel Models. In Proceedings of the Second Workshop on Irregular Applications: Architectures and
Algorithms (I A3 ) held in conjunction with SC12.
Christopher D. Krieger, Michelle Mills Strout, Catherine Olschanowsky, Andrew Stone, Stephen Guzik, Xinfeng
Gao, Carlo Bertolli, Paul Kelly, Gihan Mudalige, Brian Van Straalen, and Sam Williams. 2013. Loop Chaining: A
Programming Abstraction For Balancing Locality and Parallelism. In Proceedings of the 18th International Workshop on
High-Level Parallel Programming Models and Supportive Environments (HIPS). Boston, Massachusetts, USA.
Sriram Krishnamoorthy, Muthu Baskaran, Uday Bondhugula, J. Ramanujam, Atanas Rountev, and P Sadayappan.
2007. Effective Automatic Parallelization of Stencil Computations. SIGPLAN Not. 42, 6 (June 2007), 235–244. DOI:
https://doi.org/10.1145/1273442.1250761
Michael Lange, Lawrence Mitchell, Matthew G. Knepley, and Gerard J. Gorman. 2016. Efficient mesh management in
Firedrake using PETSc-DMPlex. SIAM Journal on Scientific Computing 38, 5 (2016), S143–S155. DOI:https://doi.org/10.
1137/15M1026092 arXiv:cs.MS/1506.07749
Fabio Luporini, Ana Lucia Varbanescu, Florian Rathgeber, Gheorghe-Teodor Bercea, J. Ramanujam, David A. Ham,
and Paul H. J. Kelly. 2015. Cross-Loop Optimization of Arithmetic Intensity for Finite Element Local Assembly. ACM
Trans. Archit. Code Optim. 11, 4, Article 57 (Jan. 2015), 25 pages. DOI:https://doi.org/10.1145/2687415
Jiayuan Meng and Kevin Skadron. 2009. Performance Modeling and Automatic Ghost Zone Optimization for Iterative
Stencil Loops on GPUs. In Proceedings of the 23rd International Conference on Supercomputing (ICS ’09). ACM, New
York, NY, USA, 256–265. DOI:https://doi.org/10.1145/1542275.1542313
Marghoob Mohiyuddin, Mark Hoemmen, James Demmel, and Katherine Yelick. 2009. Minimizing communication
in sparse matrix solvers. In Proceedings of the Conference on High Performance Computing Networking, Storage and
Analysis (SC ’09). ACM, Article 36, 12 pages. DOI:https://doi.org/10.1145/1654059.1654096
Florian Rathgeber. 2014. Productive and efficient computational science through domain-specific abstractions. Ph.D.
Dissertation. Imperial College London.
Florian Rathgeber, David A. Ham, Lawrence Mitchell, Michael Lange, Fabio Luporini, Andrew T. T. Mcrae, GheorgheTeodor Bercea, Graham R. Markall, and Paul H. J. Kelly. 2016. Firedrake: Automating the Finite Element Method by
Composing Abstractions. ACM Trans. Math. Softw. 43, 3, Article 24 (Dec. 2016), 27 pages. DOI:https://doi.org/10.1145/
2998441
Mahesh Ravishankar, John Eisenlohr, Louis-Noël Pouchet, J. Ramanujam, Atanas Rountev, and P. Sadayappan. 2012.
Code generation for parallel execution of a class of irregular loops on distributed memory systems. In Proc. Intl. Conf.
on High Perf. Comp., Net., Sto. & Anal. Article 72, 11 pages.
I. Z. Reguly, G. R. Mudalige, C. Bertolli, M. B. Giles, A. Betts, P. H. J. Kelly, and D. Radford. 2016. Acceleration of a
Full-Scale Industrial CFD Application with OP2. IEEE Transactions on Parallel and Distributed Systems 27, 5 (May 2016),
1265–1278. DOI:https://doi.org/10.1109/TPDS.2015.2453972
István Z Reguly, Gihan R Mudalige, and Mike B Giles. 2017. Loop Tiling in Large-Scale Stencil Codes at Run-time with
OPS. arXiv preprint arXiv:1704.00693 (2017).
Joel H. Salz, Ravi Mirchandaney, and Kay Crowley. 1991. Run-Time Parallelization and Scheduling of Loops. IEEE
Trans. Comput. 40, 5 (1991), 603–612.
Imperial College High Performance Computing Service. 2017. The cx2/Helen cluster. (2017). DOI:https://doi.org/10.
14469/hpc/2232

, Vol. 1, No. 1, Article 1. Publication date: August 2017.

Automated Tiling of Unstructured Mesh Computations

1:29

[34] M.M. Strout, F. Luporini, C.D. Krieger, C. Bertolli, G.-T. Bercea, C. Olschanowsky, J. Ramanujam, and P.H.J. Kelly. 2014.
Generalizing Run-Time Tiling with the Loop Chain Abstraction. In Parallel and Distributed Processing Symposium,
2014 IEEE 28th International. 1136–1145. DOI:https://doi.org/10.1109/IPDPS.2014.118
[35] Michelle Mills Strout, Larry Carter, and Jeanne Ferrante. 2003. Compile-time Composition of Run-time Data and
Iteration Reorderings. In Proc. ACM SIGPLAN Conf. Prog. Lang. Des. & Impl. (PLDI). ACM, New York, NY, USA.
[36] Michelle Mills Strout, Larry Carter, Jeanne Ferrante, Jonathan Freeman, and Barbara Kreaseck. 2002. Combining
Performance Aspects of Irregular Gauss-Seidel via Sparse Tiling. In Proceedings of the 15th Workshop on Languages
and Compilers for Parallel Computing (LCPC). Springer.
[37] Michelle Mills Strout, Larry Carter, Jeanne Ferrante, and Barbara Kreaseck. 2004. Sparse Tiling for Stationary Iterative
Methods. International Journal of High Performance Computing Applications 18, 1 (February 2004), 95–114.
[38] TOP500. 2016. cx2/Helen cluster specification in the TOP500 ranking. https://www.top500.org/system/178845. (2016).
[39] Jean Virieux. 1986.
P-SV wave propagation in heterogeneous media; velocity-stress finitedifference method.
Geophysics 51, 4 (1986), 889–901.
DOI:https://doi.org/10.1190/1.1442147
arXiv:http://geophysics.geoscienceworld.org/content/51/4/889.full.pdf
[40] Zenodo/COFFEE. 2017. coneoproject/COFFEE: A Compiler for Fast Expression Evaluation. (July 2017). DOI:https:
//doi.org/10.5281/zenodo.836678
[41] Zenodo/FIAT. 2017. firedrakeproject/fiat: The Finite Element Automated Tabulator. (July 2017). DOI:https://doi.org/
10.5281/zenodo.836679
[42] Zenodo/Firedrake. 2017. firedrakeproject/firedrake: an automated finite element system. (July 2017). DOI:https:
//doi.org/10.5281/zenodo.836680
[43] Zenodo/PETSc. 2017. firedrakeproject/petsc: Portable, Extensible Toolkit for Scientific Computation. (July 2017). DOI:
https://doi.org/10.5281/zenodo.836685
[44] Zenodo/PETSc4Py. 2017. firedrakeproject/petsc4py: The Python interface to PETSc. (July 2017). DOI:https://doi.org/
10.5281/zenodo.836684
[45] Zenodo/PyOP2. 2017. OP2/PyOP2: Framework for performance-portable parallel computations on unstructured
meshes. (July 2017). DOI:https://doi.org/10.5281/zenodo.836688
[46] Zenodo/Seigen. 2017. OPESCI/Seigen: An elastic wave equation solver for seismological problems based on the finite
element method . (Aug. 2017). DOI:https://doi.org/10.5281/zenodo.840000
[47] Zenodo/SLOPE. 2017. coneoproject/SLOPE: A run-time system to tile irregular loops. (July 2017). DOI:https:
//doi.org/10.5281/zenodo.836738
[48] Zenodo/TSFC. 2017. firedrakeproject/tsfc: The Two Stage Form Compiler firedrakeproject/tsfc: The Two Stage Form
Compiler. (July 2017). DOI:https://doi.org/10.5281/zenodo.836677
[49] Zenodo/UFL. 2017. firedrakeproject/ufl: The Unified Form Language. (July 2017). DOI:https://doi.org/10.5281/zenodo.
836683
[50] Xing Zhou, Jean-Pierre Giacalone, María Jesús Garzarán, Robert H. Kuhn, Yang Ni, and David Padua. 2012. Hierarchical
Overlapped Tiling. In Proceedings of the Tenth International Symposium on Code Generation and Optimization (CGO
’12). ACM, New York, NY, USA, 207–218. DOI:https://doi.org/10.1145/2259016.2259044

, Vol. 1, No. 1, Article 1. Publication date: August 2017.

