A

Scheduling and Tiling Reductions on Realistic machines

arXiv:1801.05909v1 [] 18 Jan 2018

Nirmal Prajapati
October 25, 2016

Computations, where the number of results is much smaller than the input data and are produced through
some sort of accumulation, are called Reductions. Reductions appear in many scientific applications. Usually,
reductions admit an associative and commutative binary operator over accumulation. Reductions are therefore highly parallel. Given unbounded fan-in, one can execute a reduction in constant/linear time provided
that the data is available. However, due to the fact that real machines have bounded fan-in, accumulations
cannot be performed in one time step and have to be broken into parts. Thus, a (partial) serialization of
reductions becomes necessary. This makes scheduling reductions a difficult and interesting problem.
There have been a number of research works in the context of scheduling reductions. We focus on the
scheduling techniques presented in [Gupta et al. 2002], identify a potential issue in their scheduling algorithm
and provide a solution. In addition, we demonstrate how these scheduling techniques can be extended to
“tile" reductions and briefly survey other studies that address the problem of scheduling reductions.

1. INTRODUCTION

Reductions are those computations in which an associative and commutative operator accumulate a set of points into a single value. Following example illustrates such a computation.
for i = 1 to N
A[i] += B[i-1];
The reduction operator is associative and commutative, which implies that accumulations
need not admit any order. Therefore, we should be able to exploit parallelism. With an
unbounded number of processors and unbounded fan-in operator, accumulations can be
done in a single time step. However, for a real machine, both the number of processors and
the fan-in are bounded. This necessitates ordering of accumulations.
Many scientific and engineering applications spend most of their execution time in nested
loops. The task of optimizing such nested loops involve dataflow analysis of the program [Feautrier 1991]. The dataflow analysis of the above program reflects loop-carried
data dependences that prevent parallelism. Therefore, scheduling reductions is a difficult
problem.
A multidimensional affine function that imposes a particular order of execution is called a
Schedule. A schedule must satisfy the precedence constraints imposed by the dependences.
Polyhedral model renders a powerful abstraction that enables precise reasoning for the legality of transformations. Iterations of the nested loops can be viewed as integer points in
a polyhedron. The computations impose dependence constraints which are represented as
affine functions of the indices. Reductions can also be described as Systems of Affine Recurrence Equations (SAREs) over polyhedral domains. Reduction dependences are implicit
in the SAREs. Such a representation allows us to focus on non-reduction dependences that
must be satisfied by the schedule.
[Redon and Feautrier 1994] present a scheduling technique that optimally schedules
reductions over CRCW PRAM model. They assumed that accumulations can happen in one
time step. This approach was extended by [Gupta et al. 2002] to work on realistic machines.
They invented a scheduling technique for machines with binary operators and exclusive
writes. They claim that on such machines, their scheduling method gives efficient solutions
with a constant factor slowdown compared to best possible schedules on a CRCW PRAM
model. We discover a flaw in their technique which violates their claim of “exclusive writes".
Using a counter example, we prove that their scheduling technique requires a machine with

concurrent writes. We solve this problem by introducing additional causality constraints and
show that the modified scheduling constraints guarantee a schedule with exclusive writes.
Furthermore, we extend the scheduling approach to tile reductions. Tiling [Wolfe 1987;
Irigoin and Triolet 1988] is a classic iteration space partitioning technique which combines a
set of points into tiles, where each tile can be executed atomically. Tiling comes in handy for
exploiting data locality [Wolf and Lam 1991], minimizing communication [Andonov et al.
2001; Xue 1997] and maximizing parallelism. Reductions are usually serialized before tiling.
Most of the tiling techniques such as [Bondhugula et al. 2008a; Doerfert et al. 2015] take
as input serialized reduction programs. Loop transformation techniques are used to find
tiling hyperplanes. In the majority of the cases, serializing imposes uniform dependences
which are tileable. Serializing reductions, however, negates the fact that accumulations can
be carried out in any order and imposes unnecessary intra-tile and inter-tile dependences.
The tiling techniques developed in this paper maximize parallelism as well as improve data
locality in a work efficient manner.
We also discuss other limitations of Gupta’s scheduling algorithm and suggest possible
improvements. We highlight the unexplored areas and address other related work in the
context of scheduling reductions.
The rest of the paper is organized as follows. Section 2 describes necessary background.
In section 3, we describe [Gupta et al. 2002] scheduling technique using examples. Section
4 exposes the flaw in their technique using a counter example which shows that “exclusive
writes" condition is violated and proposes a solution that is guaranteed to give schedules
for machines with “exclusive writes". Section 5 extends this scheduling technique to tile
reductions. Section 6 suggests other possible improvements, section 7 discusses related work
and section 8 concludes the paper.
2. BACKGROUND

A reduction of variable X can be represented as
X = reduce(⊕, f, R)

(1)

where ⊕ is an associative and commutative accumulation operator, the body of the reduction
is some variable R, and f is a projection function that maps a subset of the points in the
domain of R, represented as Dom(R), to zX ∈ Dom(X). For every zX ∈ Dom(X), the
Parametrized Reduction Domain of zX , referred
  as P (zX ), is the set of points in R that are
z
mapped to zX by the function f = z → AP
+ cP where AP is a constant matrix, cP is
p
a constant vector and p is a vector of the size parameters.
2.1. Schedule

Schedule of a variable X is a vector that represents the time instant at which zX ∈ Dom(X)
is computed and is given by λX (zX ), a multidimensional affine function. For any two variables X and Y , if X(z) depends on Y (z 0 ) then Y (z 0 ) must be computed before X(z). The
dependence imposes the following causality constraints on the schedule
λX (z)  λY (z 0 ) + TeqX (z, z 0 )

(2)

where TeqX (z, z 0 ) is the time to compute the RHS of the equation X(z) after Y (z 0 )
becomes available.
[Redon and Feautrier 1994] scheduling algorithm assumes a CRCW PRAM with unbounded number of processors and unbounded fan-in operators. The reductions can, therefore, be performed in a single time step. The causality constraints for equation (1) on such
a machine are given by
2

zX ∈ Dom(X), zR ∈ P (zX ),
λX (zX )  λR (zR ) + 1

(3)

However, this technique is not applicable for scheduling reductions on real machines.
[Gupta et al. 2002] developed a technique that schedules reductions on realistic machines
with bounded fan-in (binary) operators and exclusive writes. They claim that their scheduling algorithm generates efficient schedules with a constant fold slow down compared to
optimal schedules obtained on a CRCW PRAM model.
The following section explains in detail the scheduling algorithm as presented in [Gupta
et al. 2002].
3. GUPTA’S SCHEDULING ALGORITHM

Let λR (zR ) be the schedule of zR ∈ Dom(R). λR (zR ) = t are the equitemporal hyperplanes,
or “slices", defined as the set of points in zR ∈ P (zX ) that become available for accumulation
at time t. A temporary variable T empX(zX , t) is defined as follows




      
AP
z
c
T empX = reduce ⊕, z →
+ P
,R
(4)
ΛR
p
αR
 
z
where λR (zR ) = ΛR R + αR is the schedule for R.
p
T empX(zX , t) are the partial accumulations of equitemporal hyperplanes in P (zX ). These
intermediate results are accumulated to get final answer zX ∈ Dom(X). Equation (1) is
modified as
X = reduce(⊕, (zX , t → zX ), T empX)

(5)

Consider, X as the following reduction:
X=

i=N,j=i
X

R(i, j)

(6)

i=0,j=0

Dom(R) = {i, j|0 ≤ j ≤ i ≤ N }
Assume, λR (zR ) = i are the equitemporal hyperplanes. With this information, we can
deduce that there are N equitemporal hyperplanes. i.e. Dom(T empX) has N elements where
the t-th element is a partial accumulation of the t-th equitemporal hyperplane in P (zX ).
Figure 1(b) shows the new set of equations obtained after decomposing the reduction as
shown in equations (4) and (5).
Correspondingly, the causality constraints in equation (2) are updated to accommodate
T empX.
f (zX ) ≤ t ≤ l(zX ),
λT empX (zX , t)  t + TeqT empX (zX , t)
λX (zX )  λT empX (zX , t) + TeqX (zX , t)

(7)
(8)

where f (zX ) and l(zX ) are the first and last time steps at which values in P(zX ) are
available. TeqT empX (zX , t) is the time to compute the reduction of all the values in the
equitemporal hyperplane (zX , t). For a machine with binary operators, size(zX , t) gives
the time to accumulate the equitemporal hyperplane at (zX , t). The number of time steps
3

Fig. 1: Gupta’s scheduling technique for the equation in (a).
required for linear accumulation of all values in any box B is defined as its size, which is
also equal to the 1-norm of its principal diagonal.
For the example in Figure 1(a), the size(zX , t) of a hyperplane is given by t = i. Figure 1(c) shows the iteration space and the table in Figure 1(d) give the sizes of the equitemporal hyperplanes at (zX , t).
The scheduling constraints in equations (7) and (8) get reduced to
size(zX , t)  λX (zX ) − t

(9)

The slack of an equitemporal hyperplane is defined as sl(zX , t) = λX (zX ) − t. The
scheduling constraints in (9) are further modified to
sl(zX , t)  size(zX , t)

(10)

Note that by definition size(zX , t) is always one-dimensional scalar. If slack is truly
multidimensional then the constraint in (10) is trivially satisfied. However, if the slack is
not multidimensional, then the value of the innermost dimension of slack should be greater
than the size(zX , t).
These causality constraints can be formulated as an integer linear program and solved
using a PIP [Feautrier 1988] solver. Using the above constraints to formulate causality for
the example in Figure 1, we obtain schedule λX (zX ) as shown in Figure 1(d).
For the same example, assume that λR (zR ) = i − j are the equitemporal hyperplanes.
Figure 2(a) shows the hyperplanes at which values in P (zX ) become available for accumu4

Fig. 2: Gupta’ scheduling technique for equation (6), given λR (zR ) = i − j.

lation. For this example, P (zX ) = Dom(R). Again, there are N equitemporal hyperplanes.
i.e. Dom(T empX) has N elements where the t − th element is a partial accumulation of
the t − th equitemporal hyperplane in zX .
The size(zX , t) of each hyperplane is shown in Figure 2(b). We obtain the schedule for
zX such that it requires a machine with concurrent writes. As shown in Figure 2(b), the
accumulations are scheduled at the same time step. In such cases, [Gupta et al. 2002]
suggests that we slow down the schedule by a factor of 2 in the innermost dimension of time
and obtain enough time for the accumulations in equation (5) on a machine with binary
operators. Figures 1(c) and (d) show how this modification easily solves the problem.
5

In the next section, we provide a counter example to show that this modification does
not solve the problem of concurrent writes in all cases.
4. COUNTER EXAMPLE

Consider, variable X as the following reduction equation:
X(i) =

j=i
X

R(i, j)

(11)

j=0

Dom(R) = {i, j|0 ≤ j ≤ i ≤ N }
Here, X is a one-dimensional variable. Dom(R) is again two-dimensional. Assume,
λR (zR ) = j are the equitemporal hyperplanes. With this information, we can see that
there are j equitemporal hyperplanes. i.e. the slice at (zX , t) has a single element. The
t − th element is a partial accumulation of the t − th equitemporal hyperplane in zX which
is a single point, therefore, there are no partial accumulations.
The size of all the equitemporal hyperplanes is equal to 1. Figure 3(c) shows the iteration
space of equation (11) and Figure 3(d) shows the desird schedule.
Using these constraints, we will now solve for t. Assume, we get λR (zR ) = t = 0. The
schedule λX (zX ) is calculated using the constraint in equation (10). This constraint is
trivially satisfied with a size of 1 for all equitemporal hyperplanes. Figures 3(e) and (f)
show the issue. Here the schedule is zero-dimensional and we can no longer slow it down by
a constant factor to accumulate the intermediate results. This violates the “exclusive write"
condition and requires a machine with Concurrent Writes!.
In the previous example, note that the number of values to be accumulated into zX
is more than then number of time steps between f (zX ) and l(zX ). This prevents linear
accumulation. We define size0 (zX ) as the total number of time steps required for linear
accumulation of all values in the reduction body T empX(zX , t) of X in equation (5). Let
the total number of time steps between f (zX ) and l(zX ) be TzX .
If the condition
∀zX ∈ Dom(X), TzX + 1 > size0 (zX )

(12)

is satisfied, then concurrent writes in zX can be avoided. Suppose, f (zX ) = 1 and l(zX ) =
N . Therefore, TzX = N . If size0 (zX ) < N , then the condition is trivially satisfied. However,
this constraint cannot be satisfied for the example (11) where TzX = 1 and size0 (zX ) = N .
Let the total number of equitemporal hyperplanes in P (zX ) be EzX . Again, if EzX <
size0 (zX ) < TzX , then exclusive writes cannot be guaranteed by equation (12).
Now let’s see what happens if the condition
∀zX ∈ Dom(X), EzX + 1 > size0 (zX )

(13)

is satisfied. There would be enough time steps for the accumulation of the intermediate
values into the final answer zX . However, if there is only one equitemporal hyperplane like
our example (11), then the constraint in (13) will not be satisfied. Neither the number of
equitemporal hyperplanes EzX nor the total number of time steps TzX can be guaranteed
to be more than size0 (zX ).
Therefore, we suggest the following scheduling constraints
λX (zX )  size0 (zX )

(14)

in addition to the constraints in (10). Linear accumulations are now guaranteed on a machine
with bounded fan-in.
The above scheduling technique can be further optimized to get better schedules. In the
next section, we will see how this scheduling technique can be extended to tile reductions
in order to maximize parallelism and improve data locality.
6

Fig. 3: Gupta’s scheduling technique violates Exclusive Write condition for equation (11).

7

Fig. 4: [Gupta et al. 2002] scheduling technique for equation in (a), given λR (zR ) = i.
5. TILING REDUCTIONS

Tiling or blocking computations is a strategy of dividing the iteration space into tiles where
each tile is a set of points [Wolfe 1987; Irigoin and Triolet 1988]. A tiling is considered to be
legal if there are no dependence based cycles between tiles and if all tiles can be executed
atomically.
→
−
Let θX ( i ) define a set of tiling hyperplanes that tile the iteration space of a variable
X. For any two variables X and Y , if X depends on Y then the following tiling legality
constraint must be satisfied for all the dependencies between X and Y
→
−
→
−
θX ( i ) − θY ( i ) ≥ 0

(15)

The above condition ensures the legality of tiling as shown in [Bondhugula et al. 2008b].
With the knowledge that accumulations can be carried out in any order, we can eliminate
the reduction dependences from the dependence set. However, (15) must hold for all other
dependencies.
8

Consider the equation shown in Figure 4(a). Figure 4(b) shows the decomposition of X
as per equation (4). Figures 4(c) and (d) show the iteration space and schedule obtained
using the formulation discussed in Section 3.
We provide an incremental approach to finding tiling hyperplanes for reductions. We first
show how equitemporal hyperplanes can be tiled, succeeded by tiling P (zX ) and finally
suggest possible tilings for the reduction body R.
5.1. Tiling Equitemporal Hyperplanes

Tiling an equitemporal hyperplane (zX , t) using any tiling hyperplane is a legal tiling. This is
due to the fact that all values in an equitemporal hyperplane are available for accumulation
at the same time and that all of them contribute to a single value in T empX. Therefore, the
tiling legality condition (15) holds for any tiling hyperplane. We are left with many possible
choices. We choose orthogonal tiling hyperplanes with tiles of size s in every dimension.
We introduce a new variable T ileT empX such that (zX , t, b) ∈ Dom(T ileT empX ) maps to
the bth tile in T empX, where (zX , t) ∈ Dom(T empX).

T ileT empX = reduce ⊕, z →

AP
ΛR
γ

!

z
p
s

!



c
+ P
αR

!

!
,R

(16)

where γ is a function that divides every dimension of the equitemporal hyperplane at
(zX , t) into tiles of size s such that 1 ≤ b ≤ τ (zX , t), where τ (zX , t) is the total number
of tiles in the slice at (zX , t). With this definition of variable T ileT empX , equation (4) is
modified as
T empX = reduce(⊕, (zX , t, b → zX , t), T ileT empX )

(17)

T empX(zX , t) is the accumulation of τ (zX , t) tiles. Equation (5) remains unchanged.
Consider the equation shown in Figure 5(a). Figure 5(b) shows the decomposition of X
as per equations (16) and (17). We now obtain causality constraints for equations (16) and
(17). The precedence constraints on T ileT empX state that
1 ≤ b ≤ τ (zX , t), f (zX )  t  l(zX ), zx ∈ Dom(X)
λT ileT empX (zX , t, b)  t + TeqT ileT empX (zX , t, b)

(18)

λT empX (zX , t)  λT ileT empX (zX , t, b) + TeqT empX (zX , t)

(19)

In an equitemporal hyperplane, all the tiles can be executed simultaneously.
TeqT ileT empX (zX , t, b) is given by the size of the tile. We assume that tile size is s in every dimension. Let the number of dimension of the equitemporal hyperplane be d(zX , t).
Hence, size of a tile (zX , t, b) will be given by d(zX , t) × s. We get
λT ileT empX (zX , t, b)  t + [d(zX , t) × s]

(20)

as the constraints for λT ileT empX (zX , t, b).
TeqT empX (zX , t) of T empX, defined as the size(zX , t), is the time it takes to accumulate
all the partial answers produced by each tile, which is also equal to τ (zX , t). Therefore,
causality constraints on the schedule for equation (17) can be formulated as:
λT empX (zX , t)  t + τ (zX , t) + [d(zX , t) × s]
9

(21)

Fig. 5: Extending the scheduling techniques for tiling equitemporal hyperplanes for the
equation in (a), given λR (zR ) = i.
Similarly, we deduce the following as the causality constraints for equation (5).
λX (zX )  t + τ (zX , t) + [d(zX , t) × s]
0

λX (zX )  size (zX )

(22)
(23)

Figures 5(c) and (d) show the iteration space and schedule obtained using the formulation
discussed above. Using tile size s = 3 and N = 9, we get the number of tiles τ = 3 and
hence λX (zX ) can be scheduled as early as 15.
Note, when s does not equally divide every dimension of an equitemporal hyperplane, we
get partial tiles whose size is smaller than full tiles. Therefore, above causality constraints
are satisfied trivially for partial tiles.
With the above formulation, there are many possible choices for tile size s. If we choose
s = 1, then there will be only one point in each tile. We provide a cost function that leads
to good tile size and maximizes parallelism.
10

5.1.1. Towards finding good solutions. To minimize the total execution time, reductions of an
equitemporal hyperplane can be implemented using a binary-tree-like algorithm. However,
this is not work-efficient. To accumulate N elements, binary-tree algorithms take logN steps.
Brent’s theorem suggests N/logN parallel instances where each instance performs logN
work [Brent 1999]. Later, all N/logN parallel instances contribute to the final accumulation
of the partial results. The cost is now given by (N/logN )∗logN = N , which is work efficient.
Our proposed cost function provides such work efficient solutions.
We seek to minimize the following cost function Let,

ω = [d(zX , t) × s] − τ (zX , t)
minimize(ω)

(24)

The tile size s used in Figures 5(c) and (d) reflect the result of the cost function (24).
Note that our tile size optimization function assumes equal tile size in every dimension.
This restriction can be lifted to enable rectangular tiling.
This technique of tiling equitemporal hyperplanes applies to those cases where the equitemporal hyperplane is at least one dimensional. If the hyperplanes are zero-dimensional
then above formulation does not apply. This motivates tiling across equitemporal hyperplanes.
5.2. Tiling Parametrized Reduction Domain of zX

As shown above, tiling equitemporal hyperplanes is straightforward. Let us consider tiling
the Parameterized reduction domain of zX . If there are no dependencies between equitemporal hyperplanes, then orthogonal tiling hyperplane can be chosen and all tiles can be
launched independently. The partial answers can be accumulated to get the final answer
zX .
However, if there exists dependence in t, then orthogonal tiling hyperplanes might not
be legal. Notice, the only dependences that affect tiling legality are between equitemporal
hyperplanes (as per the definition of slices, all values become available for accumulation
in an equitemporal hyperplane at the same time). The problem of tiling P (zX ) is thus
reduced to time tiling. A well-known approach is to enforce forward communication only
constraint as presented in [Griebl et al. 2005]. If we want to maximize parallelism and
minimize synchronization, then we can use the time-partitioning technique of [Lim and Lam
1998]. If minimizing communication is also desired optimization, then the cost optimization
techniques of [Bondhugula et al. 2008a] can be used. If the dependences are uniform, then
time tiling techniques for stencil computations such as [Tang et al. 2011; Grosser et al.
2014; Bondhugula et al. 2016] can be used for maximizing parallelism together with the
concurrent start.
The dependences in t impose additional constraints on tile sizes. The tile size along
each tiling hyperplane must be greater than the length of the longest dependence in the
hyperplane. Using these tiling hyperplanes and additional dependence based constraints,
we can now optimize for tile size with our cost function.
We can eliminate the variable T empX and use only one variable T ileX to tile P (zX ).
−
T ileX is defined using an affine function →
γ which represents the tiling hyperplanes and the
tile size s. X is now an accumulation of the partial answers produced by each tile.
Additional analysis is needed to tile the reduction body. Let R be the reduction body of
X. We want to tile all the d dimensions of the reduction body using the same approach
as mentioned above. In order to do so, the schedule of tiles must admit the schedules of
both variables R and X. Tiling hyperplanes such that the tiling legality constraint (15) is
satisfied can be found. Reduction body can be partitioned using these tiling hyperplanes to
decompose the reduction.
11

6. FURTHER IMPROVEMENTS

We identify the following potential areas of improvements to the scheduling techniques
presented in this paper.
(1) Reiterating over the scheduling approach discussed in section 4, the causality constraints
are formulated by making some assumptions regarding λR and then these constraints
are simultaneously resolved to find schedules for all the variables. Hence, while formulating the causality constraints it is not possible to recognize the exact equitemporal
hyperplanes without the knowledge of λR . The scheduling techniques do not show how
to make an optimal choice of equitemporal hyperplanes while formulating the causality
constraints. The analysis is, therefore, sub-optimal and can be improved.
(2) The suggested scheduling and tiling techniques do not consider the program size parameters. Reconsider the example in Figure 4 (a) with modified size parameters such
that 0 ≤ i ≤ M and 0 ≤ i ≤ N . Assume λR (zR ) = i as the equitemporal hyperplanes as
shown in Figure 4 (c). The number of equitemporal hyperplanes will be given by M and
the size of an equitemporal hyperplane will be given by N . The causality constraints
in the equation (10) will impose the condition that the slack must be greater than or
equal to the size. Either M ≥ N or M < N . Without the knowledge of the values
of N and M , it is not possible to find a schedule that satisfies both the inequalities.
Assuming that the values of N and M are known and that M > N , after solving if we
get λR (zR ) = t = j which used N as the size of equitemporal hyperplanes, we get an
incorrect schedule of X. Therefore, it becomes necessary to consider size parameters.
(3) In the situations where equitemporal hyperplanes have different slacks, it is suggested
that equitemporal hyperplanes be scheduled in the decreasing order of slack. However,
the method does not discuss ordering of equitemporal hyperplanes when they have the
same slack.
(4) Furthermore, if the reduction operator is not commutative then accumulations must
admit an order. The scheduling techniques presented in this paper can be extended
to consider non-commutative operators. For example, a reduction computation with
an associative but non-commutative operator can be tiled using the techniques presented in Section 5 with additional constraint such that accumulations within a tile are
lexicographically ordered. The accumulations of partial answers into the final answer
can also be ordered in the lexicographically increasing order of tiles. Such an ordering
imposed by non-commutative operator might slow down the schedule by a constant
factor. Parallelism can, however, be exploited irrespective of the commutativity of the
operator.
7. RELATED WORK

The scheduling technique of [Karp et al. 1967] solved the problem of scheduling Systems
of Uniform Recurrence Equations (SUREs). Using polyhedral [Rajopadhye et al. 1986]
presented a technique for synthesizing systolic architectures from recurrence rquations which
enable scheduling Affine Recurrence Equations. [Feautrier 1992a; 1992b] give closed form
schedules as affine functions of the indices of a nested loop program. The reader is referred
to the book [Darte et al. 2000] which details scheduling algorithms for recurrence equations.
The problem of finding schedules in the presence of reductions was initially tackled by [Redon and Feautrier 1994]. They assumed a CRCW PRAM model where accumulations can
be carried out in a single time step. They also show that their technique can be extended
to machines with a bounded number of processors by serializing reductions using a partialbinary-tree algorithm. However, they did not show how this can be done efficiently. Building
over their scheduling technique, [Gupta et al. 2002] developed an algorithm to determine
effective serialization of reductions to achieve the fastest possible linear schedules on an
exclusive writes machines with bounded fan-in. Scheduling SAREs do not need to consider
12

memory based dependence like [Doerfert et al. 2015; Sato and Iwasaki 2011; Bondhugula
et al. 2008a] and hence provide more flexibility.
While scheduling SAREs with reductions, the reduction dependences are implicit which
also allows for maximal parallelism. Other techniques such as [Pugh and Wonnacott 1994;
Stock et al. 2014] make the dependences explicit to improve parallelization. Automatic
Parallelization technique such as Polly’s polyhedral optimizer [Doerfert et al. 2015] tries to
achieve parallelism by introducing privatization.
The problem of finding optimal schedules are directed towards optimizing some cost
function that miniminze latency or delay, or maximize fine-grained parallelism [Feautrier
1992a; 1992b; Redon and Feautrier 1994]. Tiling improves data locality and works such as
[Bondhugula et al. 2008a] process loops where reductions are serialized.
In certain cases, it becomes necessary to find piecewice linear schedules. It is, however,
difficult to determine the pieces automatically. In the paper [Wonnacott et al. 2015], the
authors show how to find the optimal piece-wise schedule for Optimal String Parenthesization problem and use the Mostly-Tileable technique for tiling. The schedule was, however,
found by hand.
8. CONCLUSION

We studied previous works that address the problem of scheduling SAREs in the presence of
reductions. We show that method of scheduling reductions developed in [Gupta et al. 2002]
has an error in the formulation of causality constraints which leads to concurrent writes.
We exposed this error with an example and provided a solution that guarantees exclusive
writes. The scheduling techniques presented in this paper gives optimal linear schedules.
Above all, reductions remain memory-bound computations. Therefore, exploiting data
locality using tiling techniques can improve the performance [Wolf and Lam 1991]. Using the knowledge of reduction operator being associative and commutative, we extended
Gupta’s scheduling technique to scheduling as well as tiling reductions. Tiling is also useful for coarse-grained parallelism [Lim and Lam 1998; Xue 2000]. We demonstrated that
tiling the equitemporal hyperplanes renders maximal parallelism. When the accumulations
are serialized, like most other techniques do, then similar parallelism can not be achieved
because serializing imposes an execution order on tiles.
The tile size optimization technique presented in this paper maximizes parallelism; employs data locality and provides work efficient solutions which also reduces the total number
of synchronizations at the same time. This is achieved by reducing the number of elements
that contribute to final accumulation.
REFERENCES
R. Andonov, S. Balev, S. Rajopadhye, and N. Yanev. 2001. Optimal Semi-Oblique Tiling. In Proceedings
of the thirteenth annual ACM symposium on Parallel algorithms and architectures (SPAA ’01). ACM,
New York, NY, USA.
U. Bondhugula, V. Bandishti, and I. Pananilath. 2016. Diamond Tiling: Tiling Techniques to Maximize
Parallelism for Stencil Computations. IEEE Transactions on Parallel and Distributed Systems PP, 99
(2016), 1–1. DOI:http://dx.doi.org/10.1109/TPDS.2016.2615094
U. Bondhugula, A. Hartono, J. Ramanujam, and P. Sadayappan. 2008a. PLuTo: A Practical and Fully
Automatic Polyhedral Program Optimization System. In ACM Conference on Programming Language
Design and Implementation. ACM SIGPLAN, Tuscon, AZ, 101–113.
Uday Bondhugula, Albert Hartono, J. Ramanujam, and P. Sadayappan. 2008b. A Practical Automatic
Polyhedral Parallelizer and Locality Optimizer. In Proceedings of the 29th ACM SIGPLAN Conference
on Programming Language Design and Implementation (PLDI ’08). ACM, New York, NY, USA, 101–
113. DOI:http://dx.doi.org/10.1145/1375581.1375595
Richard P. Brent. 1999. Some Parallel Algorithms for Integer Factorization.. In EuroPar’99 Parallel Processing (Lecture Notes in Computer Science, No. 1685), P. Amestoy, P. Berger, M. Daydé, I. Duff,
V. Frayssé, L. Giraud, and D. Ruiz (Eds.). Springer-Verlag, 1–22.
A. Darte, Y. Robert, and F. Vivien. 2000. Scheduling and Automatic Parallelization. Birkhäuser.

13

Johannes Doerfert, Kevin Streit, Sebastian Hack, and Zino Benaissa. 2015. Polly’s Polyhedral Scheduling
in the Presence of Reductions, In Fifth International Workshop on Polyhedral Compilation Techniques
in conjunction with HiPEAC 2015. CoRR abs/1505.07716 (2015). http://arxiv.org/abs/1505.07716
Paul Feautrier. 1988. Parametric integer programming. In RAIRO Recherche Operationnelle. 22:243– 268.
P. Feautrier. 1991. Dataflow analysis of array and scalar references. International Journal of Parallel Programming 20, 1 (Feb 1991), 23–53.
P. Feautrier. 1992a. Some efficient solutions to the affine scheduling problem, Part I, One-dimensional
Time. Technical Report 28. Labaratoire MASI, Institut Blaise Pascal.
P. Feautrier. 1992b. Some efficient solutions to the affine scheduling problem, Part II, Multidimensional
Time. Technical Report 78. Labaratoire MASI, Institut Blaise Pascal.
Martin Griebl, Paul Feautrier, and Armin Größlinger. 2005. Forward Communication Only Placements and
Their Use for Parallel Program Construction. Springer Berlin Heidelberg, Berlin, Heidelberg, 16–30.
DOI:http://dx.doi.org/10.1007/11596110_2
Tobias Grosser, Albert Cohen, Justin Holewinski, P. Sadayappan, and Sven Verdoolaege. 2014. Hybrid
Hexagonal/Classical Tiling for GPUs. In CGO. Orlando, FL, 66.
Gautam Gupta, Sanjay Rajopadhye, and Patrice Quinton. 2002. Scheduling Reductions on Realistic Machines. In Proceedings of the Fourteenth Annual ACM Symposium on Parallel Algorithms and Architectures (SPAA ’02). ACM, New York, NY, USA, 117–126. DOI:http://dx.doi.org/10.1145/564870.564888
F. Irigoin and R. Triolet. 1988. Supernode Partitioning. In the 15th ACM SIGPLAN-SIGACT symposium
on Principles of programming languages (POPL ’88). ACM, New York, NY, USA, 319–329.
R.M. Karp, R.E. Miller, and S. Winograd. 1967. The organization of computations for uniform recurrence
equations. Journal of the ACM (JACM) 14, 3 (1967), 563–590.
Amy W. Lim and Monica S. Lam. 1998. Maximizing parallelism and minimizing synchronization with affine
partitions.. In Parallel Computing. DOI:http://dx.doi.org/10.1016/S0167-8191(98)00021-0
W. Pugh and D. Wonnacott. 1994. Static Analysis of Upper and Lower Bounds on Dependences and
Parallelism. ACM Trans. Program. Lang. Syst. 16, 4 (1994), 1248–1278.
Sanjay V. Rajopadhye, S. Purushothaman, and Richard Fujimoto. 1986. On Synthesizing Systolic Arrays
from Recurrence Equations with Linear Dependencies. In Proceedings of the Sixth Conference on Foundations of Software Technology and Theoretical Computer Science. Springer-Verlag, London, UK, UK,
488–503. http://dl.acm.org/citation.cfm?id=646824.706926
Xavier Redon and Paul Feautrier. 1994. Scheduling Reductions. In Proceedings of the 8th International Conference on Supercomputing (ICS ’94). ACM, New York, NY, USA, 117–125.
DOI:http://dx.doi.org/10.1145/181181.181319
Shigeyuki Sato and Hideya Iwasaki. 2011. Automatic Parallelization via Matrix Multiplication. SIGPLAN
Not. 46, 6 (June 2011), 470–479. DOI:http://dx.doi.org/10.1145/1993316.1993554
Kevin Stock, Martin Kong, Tobias Grosser, Louis-Noël Pouchet, Fabrice Rastello, J. Ramanujam, and P.
Sadayappan. 2014. A Framework for Enhancing Data Reuse via Associative Reordering. In Proceedings
of the 35th ACM SIGPLAN Conference on Programming Language Design and Implementation (PLDI
’14). ACM, New York, NY, USA, 65–76. DOI:http://dx.doi.org/10.1145/2594291.2594342
Yuan Tang, Rezaul Alam Chowdhury, Bradley C. Kuszmaul, Chi-Keung Luk, and Charles E. Leiserson.
2011. The Pochoir Stencil Compiler. In Proceedings of the 23rd ACM symposium on Parallelism in
algorithms and architectures (SPAA ’11). ACM, New York, NY, USA, 12.
M. E. Wolf and M. Lam. 1991. A Data Locality Optimizing Algorithm. In ACM SIGPLAN Conference on
Programming Language Design and Implementation (PLDI). Totonto, CA.
M. J. Wolfe. 1987. Iteration space tiling for memory hierarchies. Parallel Processing for Scientific Computing
(SIAM) (1987), 357–361.
David Wonnacott, Tian Jin, and Allison Lake. 2015. Automatic Tiling of “Mostly-Tileable" Loop Nests.
In Fifth International Workshop on Polyhedral Compilation Techniques in conjunction with HiPEAC
2015 (IMPACT ’15).
J. Xue. 1997. Communication-Minimal Tiling of Uniform Dependence Loops. J. Parallel and Distrib. Comput. 42, 1 (January 1997), 42–59.
Jingling Xue. 2000. Loop Tiling for Parallelism. Kluwer Academic Publishers, Norwell, MA, USA.

14

