Efficient Search-Space Pruning for Integrated
Fusion and Tiling Transformations
Xiaoyang Gao1 , Sriram Krishnamoorthy1 , Swarup Kumar Sahoo1 , Chi-Chung
Lam1 , Gerald Baumgartner2 , J. Ramanujam3 , and P. Sadayappan1
1Department of Computer Science and Engineering
The Ohio State University, Columbus, OH 43210, USA{gaox,krishnsr,sahoo,clam,saday}@cse.ohio-state.edu
2
Department of Computer Science
Louisiana State University, Baton Rouge, LA 70803, USA
gb@csc.lsu.edu
3
Department of Electrical and Computer Engineering andCenter for Computation and Technology
Louisiana State University, Baton Rouge, LA 70803, USA
jxr@ece.lsu.edu

Abstract. Compile-time optimizations involve a number of transformations such as loop permutation, fusion, tiling, array contraction, etc.
Determination of the choice of these transformations that minimizes theexecution time is a challenging task. We address this problem in the context of tensor contraction expressions involving arrays too large to fit inmain memory. Domain-specific features of the computation are exploited
to develop an integrated framework that facilitates the exploration of theentire search space of optimizations. In this paper, we discuss the exploration of the space of loop fusion and tiling transformations in orderto minimize the disk I/O cost. These two transformations are integrated
and pruning strategies are presented that significantly reduce the numberof loop structures to be evaluated for subsequent transformations. The
evaluation of the framework using representative contraction expressions
from quantum chemistry shows a dramatic reduction in the size of thesearch space using the strategies presented.

1

Introduction

Optimizing compilers incorporate a number of loop transformations such as permutation, tiling, fusion, etc. Considerable work has been done on improving
locality and/or parallelism by loop fusion [8–11, 19]. Fusion often creates imperfectly nested loops, which are more complex to tile effectively than perfectly
nested loops. Several works have addressed the tiling of imperfectly nested loops[2, 20]. Although there has been much progress in developing unified frameworks
for modeling a variety of loop transformations [1, 2, 16], their use has so far been
restricted to optimization of indirect performance metrics such as reuse distance,
degree of parallelism, etc.

The development of model-driven optimization strategies that target direct
performance metrics remains a difficult task. In this paper, we address the problem in the specific domain of tensor contractions involving tensors too large to
fit into physical memory. We use certain properties of the computations in thisdomain to integrate various transformations and investigate pruning strategies
to reduce the search space to be explored.
The large sizes of the tensors involved require the development of out-ofcore implementations that orchestrate the movement of data between disk andmain memory. In this paper, we discuss the integration of loop fusion and tilingtransformations with the objective of minimizing disk I/O cost. We first divide
the input program into several independent loop nests, then enumerate the setof fusion structures of each loop nest. Then, a generalized tiling approach thatsignificantly reduces the number of loop structures to be explored is presented. It
also enables subsequent optimizations of I/O placements and loop permutations.
This approach enables an exploration of the entire search space using a realistic
performance model, without the need to resort to heuristics and search of alimited subspace of the search space to limit search time.
The rest of this paper is organized as follows. In the next section, we elaborate on the computational context of interest and introduce some preliminary
concepts; in addition, an overview of the program synthesis system and theoverall approach are given. Section 3 describes a tree partitioning algorithm.
In Section 4, we propose a loop structure enumeration algorithm and prove its
completeness. The reduction in the space of loop structures to be explored isshown for representative computations in Section 5. Conclusions are provided in
Section 6.

2

Background

The work presented in this paper is being developed in the context of the Tensor
Contraction Engine (TCE) program synthesis tool [3–6, 14]. The TCE takes as
input a high-level specification of a computation expressed as a set of tensor contraction expressions, and transforms it into efficient parallel code. The current
prototype of the TCE incorporates several compile-time optimizations which aretreated in a decoupled manner, with the transformations being performed in a
pre-determined sequence. In [12, 13], we presented an integrated approach to determine the tile sizes and I/O placements for a fixed loop structure. Techniques
to prune the search space of possible I/O placements, orderings, loop permutations and tiling for given a choice of fusion of tensor contractions were presented
in [18]. In this paper, we present a technique to enumerate the various fusion
structures and develop an algorithm to significantly reduce the number of loop
nests to be evaluated for each fusion structure.
2.1

Computational Context

In the class of computations considered, the final result can be expressed using
a collection of multi-dimensional summations of the product of several input

arrays. For example, we consider a transformation used in quantum chemistry
to transform a set of two-electron integrals from an atomic orbital (AO) basis
to a molecular orbital (MO) basis:
X
B(a, b, c, d) =C1(d, s) × C2(c, r) × C3(b, q) × C4(a, p) × A(p, q, r, s)
p,q,r,sHere, all arrays would be initially stored on disk. The indices p, q, r and s havethe same range N. The indices a, b, c and d have the same range V. Typicalvalues for N range from 60 to 1300; the value for V is usually between 50 and
1000.
The calculation of B is done in four steps to reduce the number of floating
point operations.
X
T 1(a, q, r, s) =C4(a, p) × A(p, q, r, s)p

T 2(a, b, r, s) =X

C3(b, q) × T 1(a, q, r, s)X

C2(c, r) × T 2(a, b, r, s)X

C1(d, s) × T 3(a, b, c, s)q

T 3(a, b, c, s) =r

B(a, b, c, d) =s

The sequence of contractions in this form can be represented by an operation
tree as shown in Fig. 1(a). The leaves correspond to the input arrays and theroot corresponds to the output array. The intermediate arrays and output arrayare produced by the tensor contraction of their immediate children. The edges
in the operation tree represent the producer-consumer relationship between contractions.
Assuming that the available memory space is less than V 4 (which is 3TB for
V = 800), any one of the arrays A, T 1, T 2, T 3 and B is too large to entirely
fit in memory. Therefore, if the computation is implemented as a succession offour independent steps, the intermediates T 1, T 2 and T 3 have to be written to
disk after they are produced, and read from disk before they are used in thenext step. Furthermore, the amount of disk access volume could be much larger
than the total volume of the data on disk. Since none of these arrays can be
fully stored in memory, it is not possible to read each element only once fromdisk. Suitable fusion of the common loops between producing and consuming
contractions can reduce the size of the intermediate array, making it feasible toretain it in memory. Henceforth, the term intermediate node will be used to refer
to both the intermediate array produced in the corresponding interior node, andthe contraction that produces it. The reference shall be clear from the context.
Given a choice of fusion, an intermediate node not fused with its parent
divides the operation tree into two parts, both of which can be evaluated independently. Such an intermediate node is called a cut-point. A cut-point node isassumed to be resident on disk. A connected operation tree without any interiorB = SUM(T3*C1)

for a,r,q,s,pt1a,q,r,s += Ap,q,r,s ∗ C4a,p

T3 = SUM(T2*C2)
T3

C1

for a,b,r,s,qt2a,b,r,s += t1a,q,r,s ∗ C3b,q

T2 = SUM(T1*C3)
T2

C2

for a,b,c,r,s

t3a,b,c,s += t2a,b,r,s ∗ C2c,r

T1 = SUM(A*C4)

A

C3

for
 a,b,c,d,s
Ba,b,c,d += t3a,b,c,s ∗ C1d,sC4(a) Operation tree for the four-index
transform(b) Corresponding unfused code structure

Fig. 1. Operation tree and unfused code structure for the four-index transform

cut-points is called a fused sub-tree. The divided operation tree for the four-index
transform corresponding to T 1 being a cut-point is shown in Fig. 2(a). The loop
nesting tree (LNT) represents the loop structure of a fused sub-tree. Each node
in a LNT is labeled by the indices of a set of fully permutable loops appearing
together at the same level in the imperfectly nested loop structure. Loops inthe children nodes are surrounded by loops in the parent node. Fig. 2(b) shows
two possible LNT’s for the two fused subtrees in Fig. 2(a), respectively. The
corresponding code structure is shown in Fig. 2(c).
2.2Overall approach

The program synthesis system takes an operation tree representing a set of tensor
contractions as input, and generates an efficient loop structure with explicit disk
I/O statements to implement the computation. The optimization process may
be viewed in terms of the following steps.
1. Operation Tree Partitioning: In this step, we divide the original operation
tree into several fused subtrees by identifying cut-points. The optimal loop
structures for the subtrees are independent of each other, and are determined
separately.
2. Loop Structures Enumeration: For each fused subtree, we enumerate candidate loop structures to be evaluated, as a set of LNT’s.
3. Intra-Tile Loop Placements: For a given LNT, we tile all loops at each nodeand propagate intra-tile loops to all the nodes below it.
4. Disk I/O Placements and Orderings: We then explore various possible placements and orderings of disk I/O statements for each disk array in a tiled loop
structure with a pruning strategy to determine the best placement and ordering.

T1 = SUM(A*C4)

Aa,r,q,s,t (T1)

C4

loopN est2 :B = SUM(T3*C1)a,b

T3 = SUM(T2*C2)
T3

C1

r,s

T2 = SUM(T1*C3)
T2

C2

loopN est1 :for a,r,q,s,pt1a,q,r,s += Ap,q,r,s ∗ C4a,pq (T2)c (T3)

for a,b

for r,s

q
 for
 
  t2r,s += t1a,q,r,s ∗ C3b,qc,d,s (B)   for ct3c,s += t2r,s ∗ C2c,r for c,d,s
Ba,b,c,d += t3c,s ∗ C1d,s(b) Loop nesting trees
T1

C3(a) Divided operation trees(c) Corresponding code
structure

Fig. 2. Representations involved in generation of a fused code structure.

5. Tile Size Selection: For each combination of loop transformations and I/Oplacements, the I/O cost is formulated as a non-linear optimization problem
in terms of the tile sizes. The tile sizes that minimize the disk I/O cost are
determined using a general-purpose non-linear optimization solver.
6. Code Generation: We calculate the disk access cost for each solution obtained, and generate code for the one with the minimal disk I/O cost.
The possible choices of fused subtrees are first enumerated. This is explainedin Section 3. Given a fused sub-tree, the optimal loop structure and the corresponding cost can be determined by the following steps: 1) enumerating all
candidate loop structures; 2) enumerating placements and ordering of disk I/O
statments; 3) determining the tile size to minimize the disk I/O cost for eachcombination; and 4) selecting the program structure with the minimal disk I/O
cost. The algorithm for enumerating candidate loop structures is discussed in
Section 4. The search space of disk I/O placements and orderings, loop permutations and tile sizes is pruned and modeled as a non-linear optimization problem
in [18], which is then solved to determine the disk I/O cost. In this paper, we
focus on determination of the fused sub-trees and the enumeration of candidate
loop nesting trees to be evaluated.

3

Tree Partitioning

In this section, we discuss the procedure to enumerate the set of all fused subtrees
to be evaluated. In general, fusing a loop between the producer of an intermediatearray and its consumer eliminates the corresponding dimension of the array andreduces the array size. If the array fits in memory after fusion, no disk I/O

is required for that array. On the other hand, if the array does not fit in thephysical memory after fusion, the disk I/O cost will remain the same and there
is no improvement in locality. Therefore, fusion of any loops corresponding toan intermediate node is assumed to cause the resulting intermediate to residein memory. Alternatively, an intermediate node not fused with its parent (cutpoint) is assumed to reside in disk.
An arbitrary operation tree with M intermediate nodes theoretically has
O(2M ) possible fused sub-trees, but not all of them are legal. If both the childrenof an intermediate node are fused with it, then the loops corresponding to thesummation indices in the given node must be the outermost loops; and it can
not be fused with its parent anymore. Thus, either the node itself or one of its
children must be a cut-point.
Based on this property, we can restrict the number of top sub-trees to O(M 2 ).
The algorithm to enumerate the fused sub-trees rooted at a given node is shownin Algorithm 1. It proceeds in a bottom-up fashion, constructing all fused subtrees rooted at a given node from those of its children. Given a node t with twochildren lef t and right, we can extend a fused sub-tree from either lef t or rightto include the given node. These sub-trees can further be extended to include thegiven node’s parent. Besides, the given node can be considered as a cut-point.
In this scenario, all possible pairs of left and right fused sub-trees may form avalid fused sub-tree for the given node. The field t.T reeSet represents the set of
fused sub-trees which can be extended to include the parent of t.

4

Loop Structure Enumeration

In this section, we first present an algorithm that can generate a set of loop
structures for a fused sub-tree. Then, we present the result that for any loop
structure S of the fused sub-tree, we can find a corresponding loop structure S ′in the generated set, so that S ′ can be transformed to S by use of a multi-level
tiling strategy.
4.1

Enumeration Algorithm

In the previous section, we showed that a fused sub-tree must be in one of these
two forms:– All contractions form a chain called a contraction chain. For instance, Fig. 1
is such an operation tree in which the contraction chain is T 1, T 2, T 3, B.
– The contractions form two chains joining at the root node. In this case,the contraction chain is connected by these two chains. An example ofsuch an operation tree is shown in Fig. 3, in which the contraction chain
is T 1, T 2, B, T 3, T 4.
Given an operation tree that has n contraction nodes t1 , t2 , ..., tn , let ti .indicesdenote all loop indices surrounding the contraction node ti . First, we create acontraction chain of the operation tree. It corresponds to a sequence of perfectly

Algorithm 1 EnumerateFusedSubtrees(t: the root of a subtree) returns T reeSet
t1 = the left child of t; t2 = the right child of t; T reeSet = empty//Only one subtree
if both t1 and t2 are input nodes thenCreate a new Tree T r with T r.CutpointSet = ∅Insert T r into T reeSet
end if//Extending subtrees from the child not an input
if t1 is an input node and t2 is an intermediate node thenchildSet = t2 .T reeSet
Create a new Tree T r with T r.CutpointSet = {t2}Insert T r into T reeSet
end if
if t2 is an input node, and t1 is an intermediate node thenchildSet = t1 .T reeSet
Create a new Tree T r with T r.CutpointSet = {t1}Insert T r into T reeSet
end iffor each subtree st in childSet do
Create a new Tree T r with T r.CutpointSet = st.CutpointSet
Insert T r into T reeSet
end fort.T reeSet = T reeSet//Entending subtrees from either child, and cutting another child off
if both t and t2 are intermediate nodes thenchildSet1 = t1 .T reeSet
for each subtree st in childSet1 do
Create a new Tree T r with T r.CutpointSet = {st.CutpointSet, t2}
Insert T r into T reeSet
end forchildSet2 = t2 .T reeSet
for each subtree st in childSet2 do
Create a new Tree T r with T r.CutpointSet = {st.CutpointSet, t1}
Insert T r into T reeSet
end forCreate a new Tree T r with T r.CutpointSet = {t1, t2}
Insert T r into T reeSet
t.T reeSet = T reeSet//Merging subtrees from both children, and extending the result
for each pair of subtrees st1 in childSet1 and st2 in childSet2 do
Create a new Tree T r
T r.CutpointSet = {st1.CutpointSet, st2.CutpointSet}
Insert T r into T reeSet
end forend ifB = SUM(T2*T3)

T2 = SUM(T1*C3)
T3

T1 = SUM(A*C4)AC4T3 = SUM(C2*T4)

C3C2

T4 = SUM(D*C1)

DC1

Fig. 3. An operation tree with two chains

nested loops. Many different choices exist for the ordering of the fusions within
this sequence of perfectly nested loop nests. Each of the perfectly nested loops
corresponding to a contraction can be considered an independent loop nestingtree. The fusion of sub-trees producing and consuming an intermediate array
creates an imperfectly nested loop nest, in which some of the common loops are
merged. The process of construction of the loop nesting trees corresponding toa fused sub-tree can be modeled as a paranthesization problem. Consider thesequence of contraction nodes T1, T2, T3, and B in the operation tree shown in
Fig. 1. ((T 1(T 2 T 3))B) corresponds to a parenthesization in which the contractions producing T3 and consuming T3 are fused first and the resulting loop nest
is fused with the contractions producing T1 and B, in that order. Fig. 4 shows
one possible parenthesization for the four-index transform and the corresponding
loop nesting tree.
We enumerate all possible parenthesizations of the contraction chain. For
each parenthesization, a maximally fused loop structure is created by a recursive construction procedure. We call it maximally fused since, in the construction
procedure, each intermediate node will have its indices fused as much as possible
with its parent. The construction procedure is shown in Algorithm 2. It takes
a parenthesization P as input, and generate a corresponding LNT. A parenthesization of a contraction chain with n nodes has n − 1 pairs of parentheses. Eachpair of parentheses includes two elements, left and right element. Each element
is either a single contraction node or a parenthesization of a sub-chain within apair of parentheses.
Fig. 4 illustrates this proceduce for the ((T 1(T 2 T 3))B) parenthesization of
the four-index transform.

Algorithm 2 Construction(P )//Given a parenthesization, the algorithm map it to a maximally fused loop structure
in LNTl = P.lef t
r = P.right
if l is a parenthesization thenlt = Construction(lef t)
else if l is a contraction thenlt = Create a new LNT node
lt.indices = l.indiceslt.children = null
lt.contraction = l {lt is a leaf, which includes a contraction node in it}
end if
if r is a parenthesization thenrt = Construction(right)
else if r is a contraction thenrt = Create a new LNT node
rt.indices = r.indicesrt.children = null
rt.contraction = r {rt is a leaf, which includes a contraction node in it}
end ifcomindices = lt.indices ∩ rt.indiceslt.indices = lt.indices − comindicesrt.indices = rt.indices − comindiceslnt = Create a new LNT node
lnt.indices = comindiceslnt.children = {lt, rt}
return lnt

4.2

CompletenessIn this section, we state results that are useful in proving that the set of maximally fused loop structures generated by the enumeration algorithm above canrepresent all loop structures of a fused subtree.
Given an arbitrary loop nesting tree lnt, we can map it to a maximal fused
loop nesting tree lnt′ , which is generated by the enumeration algorithm above
and can be translated to lnt with proper multi-level tiling strategy. The mapping
algorithm consists of two steps:
1. Take lnt as input, and create a parenthesization P of the contraction chain
using the generation routine provided in Algorithm 3.
2. Apply the construction procedure in Algorithm 2 on P to generate a maximally fused loop structure lnt′ .
Obviously, lnt′ is the set of maximally fused loop structures generated by theenumeration algorithm. We note that lnt′ can be translated to lnt by sinking
indices at upper levels down.

Parenthesization(T2 T3)

LNT
a,b,r,sq (T2)c (T3)(T1 (T2 T3))a,r,sp,q (T1)bq (T2)c (T3)((T1 ( T2 T3)) B)a,sb,c,d (B)a,r,srp,q (T1)bq (T2)c (T3)

Fig. 4. Construction of a maximally fused loop structure for a particular parenthesization of the four-index transform.

Remark 1. For any pair of contraction nodes ti and tj , let common(lnt, ti , tj )be defined as the loops shared by ti and tj in lnt. We have common(lnt, ti , tj )⊆ common(lnt′ , ti , tj ).
Remark 2. If common(lnt, ti , tj ) ⊂ common(lnt′ , ti , tj ), then we can transform lnt′ to form lnt′′ by sinking indices down, so that common(lnt, ti , tj ) =common(lnt′′ , ti , tj )Applying the sinking operation in Remark 2 for each pair of contraction
nodes (ti , tj ), we can transform lnt′ to lnt′′ , which satisfies the condition:∀(ti , tj ), common(lnt, ti , tj ) = common(lnt′′ , ti , tj ). After that, if a node r has
no indices in r.indices, we remove r from lnt′′ , and put all children of r to its
parent. Then, lnt′′ is same as lnt.
Using a multi-level tiling strategy, a maximally fused loop strcuture can betransformed into an arbitrarily fused loop structure by appropriate choice of tile

Algorithm 3 Parenthesize(lnt)//Given an LNT, the algorithm map it to a corresponding parenthesization
if lnt.children 6= null thenP = nullfor each child c in lnt.children do
P ′ = Parenthesize(c)if P is null thenP = P′else
P = new Parenthesization(P , P ′ )end ifend forelse
P = c.contraction {c is a leaf and includes a contraction node}end ifreturn P

sizes. Multi-level tiling can transform the LNT of a loop structure as follows.
Each loop present in the root is split into two components, an inter-tile loop andan intra-tile loop. The intra-tile loop is placed on child nodes of the root. Thenthe loops present at each of the child nodes, including the intra-tile loops from theroot, are again split and intra-tile loops are placed on their respective child nodes.
This process is porformed recursively until leaf nodes are encountered. The loop
structure corresponding to the LNT can also be transformed accordingly. Fig. 5
shows the tiling of loop a in the LNT in Fig. 4 and the relationship between
different tiles, where a.range represents the range of loop a.
The sinking operation in an LNT can be modeled as a multi-level tiling in
the loop structure. If we tile a fused loop with a tile size equal to its loop range,
it leads to the same result as sinking the loop index from the original node to its
children. Let S and S ′ be loop structures representd by lnt and lnt′ respectively.
Since we can transform lnt′ to lnt by sinking operations, we can also transformS ′ to S by suitable multi-level tiling. We use an example to show the details of
the transformation procedure below.
An arbitrary fully fused loop structure S for the four-index transform isshown in Fig. 6(a), and the corresponding maximally fused loop structure S ′may be seen in Fig. 6(b). After we apply multi-level tiling, S ′ is translated to
the form shown in Fig. 7(a). In addition, if we set ranges of inter-tile loops asshown below, and remove all loops with range = 1, S ′ can be rewritten in theform shown in Fig. 7(b), which is exactly the same as S. The indexing of theintermediate arrays has been shown in a generic fashion.
aT2 = aT3 = sT1 = sT2 = sT3 = rT2 = qT1 = 1; aT1 = a.range; rI1 = r.range

aT1,saT1 .range × aI1 .range = a.range
aI1, b,c,d (B)a,r,saT2, raI2, p,q (T1)

aT3 .range × aI3 .range = aI2 .range

aT3,baI3,q (T2)

aT2 .range × aI2 .range = aI1 .range

aI3,c (T3)(a) Multi-level tiling loop a(b) Range of different level tiles

Fig. 5. An example of multi-level tiling in LNT
for
 a
for r

 for q,s,p  t1s,q += Ap,q,r,s ∗ C4a,p

b,s,q
 for
  t2b,r,s += t1s,q ∗ C3b,q for b,c,r,s t3b,c,s += t2b,r,s ∗ C2c,r for b,c,d,s
Ba,b,c,d += t3b,c,s ∗ C1d,s(a) Arbitrary fused loop structure: S

for a,sfor r
 for q p
  for    t1 += Ap,q,r,s ∗ C4a,p   for b 
t2b += t1 ∗ C3b,q  for b,c 
t3b,c += t2b ∗ C2c,r for b,c,d
Ba,b,c,d += t3b,c ∗ C1d,s(b) Maximally fused loop structure: S’Fig. 6. An arbitrary loop structure and the corresponding maximally fused structure

4.3

Complexity

The total number of loop structures generated by the enumeration algorithm
is the same as the number of parenthesizations of the contraction chain. For acontraction chain with n nodes, the number of all possible parenthesizations iscalled the nth Catalan Number. It is exponential in n, and the upper bound isO(4n /n3/2 ). In contrast, the number of possible loop structures is potentiallyexponential in the total number of distinct loop indices in the n intermediate
nodes, a considerably larger number. The fused operation tree is not very longfor most representative computations. In most practical applications, a fused
subtree usually has no more than 5 contractions in a single chain. Note that the
nth Catalan Number is not very large when n is small. The first six Catalan
Numbers are listed here: 1, 1, 2, 5, 14, 42, ....

for aT1 ,sT1

for
 aT1
for rT1 ,aT2 ,sT2

for rT1
 for qT1 ,rT2 ,aT3 ,sT3

 for
 
p,qI1 ,sI3   for p,qI1 ,rI2 ,aI3 ,sI3
  t1aI,qI,rI,sI + = Ap,q,r,s ∗ C4a,p   t1aI,qI,rI,sI + = Ap,q,r,s ∗ C4a,pb,qI1 ,sI3
 for
   for b, qI ,rI ,aI ,sI1
2
33
  t2
 
aI,b,rI,sI + = t1aI,qI,rI,sI ∗ C3b,q
t2aI,b,rI,sI + = t1aI,qI,rI,sI ∗ C3b,q 
b,c,rI1 ,aI2 ,sI2
 for
  for b, c, rI1 , aI2 ,sI2
  t3aI,b,c,sI + = t2aI,b,rI,sI ∗ C2c,r 
t3aI,b,c,sI + = t2aI,b,rI,sI ∗ C2c,r for b,c,d,sI1
 for aI1 , b, c, d, sI1

Ba,b,c,d + = t3aI,b,c,sI ∗ C1d,s
Ba,b,c,d + = t3aI,b,c,sI ∗ C1d,s(a) After inserting intra-tile loops(b) After selecting proper tile counts

Fig. 7. Translate S’ to S by multi-level tiling strategy

5

Experimental Results

The enumeration algorithm discussed in Section 4.1 generates a set of candidates
loop structures to be considered for data locality optimization. Without thisalgorithm, and generalized tiling, the set of loop structures to be evaluated might
be too large, precluding their complete evaluation and necessitating the use ofheuristics.
We evaluate the effectiveness of our approach using the following tensor contractions from representative computations from the quantum chemistry domain.
1. Four-index transform (4index): Introduced in Section 2.
2. CCSD: The second and the third computations are from the class of Coupled Cluster (CC) equations [7, 15, 17] for ab initio electronic structure modeling. The sequence of tensor contraction expressions extracted from thiscomputation is shown as follows:P
S(j, i, b, a) = l,k (A(l, k, b, a)
P P× ( d ( c (B(d, c, l, k) × C(i, c)) × D(j, d)))3. CCSDT: This is a more accurate CC model. A sub-expression from the
CCSDT theory is:P
S(h3, h4, p1, p2)= p9,h6,h8 (y ooovvv(h8, h6, h4, p9, p1, p2)×PP
h10 t vo(p9, h10) ×
p7 (t vo(p7, h8)×P(tvo(p5,
h6)
×v
oovv(h10,
h3,
p7,
p5))
p5
We evaluated the fused subtree corresponding to the entire operation tree
without any cut-points. The number of all possible loop structures and thenumber of candidate loop structures enumerated by our approach are shownTable 1. Effectiveness of pruning of loop structures.
#Contractions #Loop structures Reduction
Total Pruned
4index
CCSD
CCSDT

4
3
4

241
69
182

5
2
5

98%
97%
98%

in Table 1. It can be seen that a very large fraction of the set of possible loop
structures, up to 98%, is pruned away using the approach developed in thispaper.

6

Conclusions

In this paper we addressed the problem of optimizing the disk access cost of
tensor contraction expressions by applying loop transformations. We discussed
approaches to partitioning of the operation tree into fused sub-trees and generating a small set of maximally-fused loop structures that cover all possible
imperfectly nested fused loop structures. The approach was evaluated on a setof computations representative of the targeted quantum chemistry domain anda significant reduction was demonstrated in the number of loop structures to beevaluated.
Acknowledgments This work is supported in part by the National Science Foundation through awards 0121676, 0121706, 0403342, 0508245, 0509442, 0509467,
and 0541409.

References
1. N. Ahmed, N. Mateev, and K. Pingali. Synthesizing transformations for locality
enhancement of imperfectly nested loops. In Proc. of ACM Intl. Conf. on Supercomputing, 2000.
2. N. Ahmed, N. Mateev, and K. Pingali. Tiling imperfectly-nested loops nests. In
Proc. of SC 2000, 2000.
3. G. Baumgartner, D.E. Bernholdt, D. Cociorva, R. Harrison, S. Hirata, C. Lam,
M. Nooijen, R. Pitzer, J. Ramanujam, and P. Sadayappan. A High-Level Approach
to Synthesis of High-Performance Codes for Quantum Chemistry. In Proc. of SC
2002, November 2002.
4. D. Cociorva, G. Baumgartner, C. Lam, J. Ramanujam P. Sadayappan, M. Nooijen,
D. Bernholdt, and R. Harrison. Space-Time Trade-Off Optimization for a Class ofElectronic Structure Calculations. In Proc. of ACM SIGPLAN PLDI 2002, pages
177–186, 2002.
5. D. Cociorva, X. Gao, S. Krishnan, G. Baumgartner, C. Lam, P. Sadayappan, and
J. Ramanujam. Global Communication Optimization for Tensor Contraction Expressions under Memory Constraints. In Proc. of IPDPS, 2003.

6. D. Cociorva, J. Wilkins, G. Baumgartner, P. Sadayappan, J. Ramanujam, M. Nooijen, D. E. Bernholdt, and R. Harrison. Towards Automatic Synthesis of HighPerformance Codes for Electronic Structure Calculations: Data Locality Optimization. In Proc. of the Intl. Conf. on High Performance Computing, volume 2228,
pages 237–248. Springer-Verlag, 2001.
7. T. Crawford and H. F. Schaefer III. An Introduction to Coupled Cluster Theory
for Computational Chemists. In K. Lipkowitz and D. Boyd, editor, Reviews in
Computational Chemistry, volume 14, pages 33–136. John Wiley, 2000.
8. C. Ding and K. Kennedy. Improving effective bandwidth through compiler enhancement of global cache reuse. J. Parallel Distrib. Comput., 64(1):108–134,
2004.
9. G. Gao, R. Olsen, V. Sarkar, and R. Thekkath. Collective Loop Fusion for ArrayContraction. In Proc. of the Fifth LCPC Workshop, 1992.
10. K. Kennedy. Fast greedy weighted fusion. In Proc. of ACM Intl. Conf. on Supercomputing, 2000.
11. K. Kennedy and K. S. McKinley. Maximizing loop parallelism and improving data
locality via loop fusion and distribution. In Proc. of Languages and Compilers for
Parallel Computing, pages 301–320. Springer-Verlag, 1993.
12. S. Krishnan, S. Krishnamoorthy, G. Baumgartner, C. Lam, J. Ramanujam, P. Sadayappan, and V. Choppella. Efficient synthesis of out-of-core algorithms using anonlinear optimization solver. In Proc. of IPDPS, page 34b, 2004.
13. S. Krishnan, S. Krishnamoorthy, G. Baumgartner, C. Lam, J. Ramanujam, P. Sadayappan, and V. Choppella. Efficient synthesis of out-of-core algorithms using
a nonlinear optimization solver. Journal of Parallel and Distributed Computing,66(5):659–673, May 2006.
14. C. Lam. Performance Optimization of a Class of Loops Implementing MultiDimensional Integrals. PhD thesis, The Ohio State University, Columbus, OH,
August 1999.
15. T. J. Lee and G. E. Scuseria. Achieving chemical accuracy with coupled cluster theory. In S. R. Langhoff, editor, Quantum Mechanical Electronic Structure
Calculations with Chemical Accuracy, pages 47–109. Kluwer Academic, 1997.
16. A. W. Lim and M. S. Lam. Maximizing Parallelism and Minimizing Synchronization with Affine Partitions. Parallel Computing, 24(3-4):445–475, May 1998.
17. J. M. L. Martin. Benchmark Studies on Small Molecules. In P. v. R. Schleyer, P. R.
Schreiner, N. L. Allinger, T. Clark, J. Gasteiger, P. Kollman, and H. F. Schaefer
III, editors, Encyclopedia of Computational Chemistry, volume 4, pages 115–128.
John Wiley, 1998.
18. S. K. Sahoo, S. Krishnamoorthy, R. Panuganti, and P. Sadayappan. Integrated loop
optimizations for data locality enhancement of tensor contraction expressions. In
Proc. of Supercomputing (SC 2005), 2005.
19. S. Singhai and K. S. McKinley. Loop Fusion for Parallelism and Locality. In
Proc. of Mid-Atlantic States Student Workshop on Programming Languages andSystems, 1996.
20. Y. Song and Z. Li. New Tiling Techniques to Improve Cache Temporal Locality.
In Proc. of ACM SIGPLAN PLDI, 1999.