MATHEMATICAL ENGINEERINGTECHNICAL REPORTS

Domain-Specific Optimization for Skeleton
Programs Involving Neighbor Elements

Kento EMOTO, Kiminori MATSUZAKI,
Zhenjiang HU and Masato TAKEICHIMETR 2007–05

February 2007

DEPARTMENT OF MATHEMATICAL INFORMATICSGRADUATE SCHOOL OF INFORMATION SCIENCE AND TECHNOLOGYTHE UNIVERSITY OF TOKYOBUNKYO-KU, TOKYO 113-8656, JAPANWWW page: http://www.i.u-tokyo.ac.jp/edu/course/mi/index e.shtml

The METR technical reports are published as a means to ensure timely dissemination ofscholarly and technical work on a non-commercial basis. Copyright and all rights therein
are maintained by the authors or by other copyright holders, notwithstanding that they
have offered their works here electronically. It is understood that all persons copying thisinformation will adhere to the terms and constraints invoked by each author’s copyright.
These works may not be reposted without the explicit permission of the copyright holder.

Domain-Specific Optimization for Skeleton Programs
Involving Neighbor Elements
Kento EMOTO, Kiminori MATSUZAKI, Zhenjiang HU and Masato TAKEICHI
Graduate School of Information Science and Technology
The University of Tokyo{emoto,kmatsu}@ipl.t.u-tokyo.ac.jp{hu,takeichi}@mist.i.u-tokyo.ac.jp

Abstract. Skeletal parallel programming enables us to develop parallel programs easily
by composing ready-made components called skeletons. However, a simply-composedskeleton program often lacks efficiency due to overheads of intermediate data structuresand communications. Many studies have focused on optimizations by fusing successive
skeletons to eliminate the overheads. Existing fusion transformations, however, are too
general to achieve adequate efficiency for some classes of problems. Thus, a specificfusion optimization is needed for a specific class. In this paper, we propose a specificoptimization of skeleton programs involving neighbor elements, which is often seen inscientific computations. We start with a normal form that abstracts the programs ofinterest. Then, we develop fusion rules that transform a skeleton program into thenormal form. Finally, we make efficient parallel implementation of the normal form.

1

Introduction

Recently, the increasing popularity of parallel machines like PC clusters and multi-coreCPUs attracts more and more users. However, development of efficient parallel programs
is difficult due to synchronization, interprocessor communications, and data distributionthat complicate the parallel programs. This difficulty calls for a methodology of parallel
programming with ease, and many studies have addressed themselves to it. As one promisingsolution, skeletal parallel programming [2, 10] has been proposed.
In skeletal parallel programming users develop parallel programs by composing skeletons,
which are abstracted basic patterns in parallel programs. Each skeleton is given as a higher
order function that takes concrete computations as its parameters, and conceals low-level
parallelism from users. Therefore, users can develop parallel programs with the skeletons in
a similar way to developing sequential programs.
Efficiency is one of the most important topics in the research of skeletal parallel programming. Since skeleton programs are developed in a compositional style, they often have
overheads of many loops and intermediate data. To make skeleton programs efficient, notonly each skeleton is implemented efficiently in parallel, but also optimizations over multiple
skeletons are necessary.
There have been several studies on the optimizations over multiple skeletons based on
fusion transformations [5–7, 9, 13], which were studied in depth in the field of functional
programming [4, 12]. In particular, general fusion optimizations [5, 7, 9, 13] have achieved
good results both in theory and in practice. For example, Hu et al. [7] proposed a set of
fusion rules based on a general form of skeletons named accumulate.
Although the general fusion optimizations so far are reasonably powerful, there is stilllarge room for further optimizations. Due to the generality of their fusion transformations,some overheads in skeleton programs are left through the general fusion optimizations. In
many cases such overheads can be removed if we provide a program-specific implementation.
Thus, some specific optimizations are important for efficiency of skeleton programs.
1

In this paper, we propose a specific optimization for skeleton programs that involve
neighbor elements, which is often seen in scientific computations. The proposed optimization
is based on the following strategy.
First, we formalize a normal form that captures the domain-specific computations.
Then, we develop fusion rules that transform a skeleton program into the normal
form. Finally, we provide an efficient parallel implementation of the normal form.
We formalized a normal form and fusion rules for the class of skeleton programs and developed a small system for fusing skeleton programs into the normal form implemented
efficiently in parallel. The experiment results show effectiveness of the domain-specific optimization.
The rest of this paper is organized as follows. Section 2 defines skeletons and our target
skeleton programs in this paper. It also explains a strategy for optimization of our target
skeleton programs. Section 3 designs a normal form for the target skeleton programs. Section 4 gives fusion rules for transformation of a skeleton program into the normal form.
Section 5 explains parallel implementation of the normal form. Section 6 expands our target
programs. Section 7 shows experiment results. Section 8 discusses the applicability of ourstrategy and related work. Section 9 concludes this paper.

2

Preliminaries

2.1

NotationsNotations in this paper follow that of Haskell [1], a pure functional language.
Function application is denoted by a space and the argument may be written withoutbrackets. Thus f a means f (a) in ordinary notation. Functions are curried, i.e. f a b means(f a) b. The function application binds more strongly than any other operator, so f a ⊗ bmeans (f a)⊗b, but not f (a⊗b). Function composition is denoted by ◦, so (f ◦g) x = f (g x)from its definition. The identity function, which is the identity of function composition, isdenoted by id. Binary operators are denoted by ⊕, ⊗ and so on. Binary operators can beused as functions by sectioning as follows: a ⊕ b = (a⊕) b = (⊕b) a = (⊕) a b. The identity ofa binary operator ⊕ is denoted by ı⊕ . Pairs are Cartesian products of plural data, written
like (x, y). It can be extended to the case of arbitrary number of elements like (x, y, z).
2.2

Lists and Basic Functions

Lists are constructed by three constructors: [ ] makes an empty list, [·] makes a singleton,and ++ concatenates two lists.
data List α = List α ++ List α| [·] α| []Here, the concatenation constructor ++ is associative. For simplicity, we denote List α by[α]. With the above constructors, a list of three elements a0 , a1 and a2 is constructed as[·] a1 ++ [·] a2 ++ [·] a3 . For simplicity, we denote a list by lining elements up between ‘[’
and ‘]’ separated by ‘,’ . So, the above list is denoted by [a1 , a2 , a3 ]. We also use a notationa : x = [a] ++ x to show the head of a list.
2

We define basic functions on lists that used in the rest of this paper.
length (a : x) = 1 + length x
length [ ]=0init (x ++ [a]) = x
init [ ]= []tail (a : x) = x
tail [ ]= []last (x ++ [a]) = a
head (a : x) = a
rev (a : x) = rev x ++ [a]rev [ ]
= []repeat a = a : repeat aat i (a : x) = if i = 0 then a else at (i − 1) x
Function length calculates the length of a list. Function init removes the last element ofa list, while tail removes the head element. Function last returns the last element of a list,
while head returns the head element. Function rev reverses a list. Function repeat generates
an infinite list of the given element. Function at returns the element at the given index.
2.3

Skeletons on Lists

We introduce skeletons on lists. Formal definitions of those skeletons are shown in Figure 1.
In the following, we explain skeletons using intuitive definitions and examples.
Skeleton map applies a given function f to each element of the list.
map f [a1 , . . . , an ] = [f a1 , . . . , f an ]For example, computation that calculates squares of each element is described with map asfollows.
map sqr [1, 2, 3, 4] = [1, 4, 9, 16]
Skeleton zip zips two lists of the same length.
zip [a1 , . . . , an ] [b1 , . . . , bn ] = [(a1 , b1 ), . . . , (an , bn )]Similarly, zip3 zips three lists. A composition of map and zip is often denoted by zipWith f =map (λ(x, y) → f x y) ◦ zip.
Skeleton shift shifts elements of a list to the right by one, and inserts the given value
e as the leftmost element. The rightmost element is thrown away. Similarly, skeleton shift
shifts elements to the left and inserts the given value as the rightmost element.
shift e [a1 , . . . , an ] = [a2 , . . . , an , e]shift e [a1 , . . . , an ] = [e, a1 , . . . , an−1 ]Skeleton reduce calculates a reduction of a list by the given associative binary operator⊕. Application of reduce to an empty list results in the identity ı⊕ of operator ⊕.
reduce (⊕) [a1 , . . . , an ] = a1 ⊕ · · · ⊕ an
For example, a summation of the elements are calculated by reduce as follows.
reduce (+) [1, 2, 3, 4] = 103(= 1 + 2 + 3 + 4)

map :: (α → β) → [α] → [β]map f [ ]= []map f (a : x) = f a : map f x
zip :: [α] → [β] → [(α, β)]zip [ ] [ ]= []zip (a : x) (b : y) = (a, b) : zip x yshift :: α → [α] → [α]shift e (a : x) = x ++ [e]shift :: α → [α] → [α]shift e (x ++ [a]) = e : x
reduce :: (α → α → α) → [α] → α
reduce (⊕) [ ]= ı⊕
reduce (⊕) (a : x) = a ⊕ reduce (⊕) xscan :: (α → α → α) → α → [α] → [α]scan (⊕) e [ ]= []scan (⊕) e (a : x) = let e0 = (e ⊕ a) in e0 : (scan (⊕) e0 x)
scanr :: (α → α → α) → α → [α] → [α]scanr (⊕) e [ ]= []scanr (⊕) e [a]= [a ⊕ e]scanr (⊕) e (a : x) = let (c : y) = scanr (⊕) e x in (a ⊕ c) : (c : y)Fig. 1. Definitions of Skeletons on Lists

Skeleton scan accumulate the intermediate results of a reduction.
scan (⊕) e [a1 , . . . , an ] = [b1 , . . . , bn ]where bi = e ⊕ a1 ⊕ · · · ⊕ ai
For example, an accumulation of a summation is described with scan as follows.
scan (+) 2 [1, 2, 3, 4] = [3, 5, 8, 12]Skeleton scanr is the reverse of scan.
scanr (⊕) e [a1 , . . . , an ] = [c1 , . . . , cn ]where ci = ai ⊕ · · · ⊕ an ⊕ e
For example, a reverse accumulation of a summation is described with scanr as follows.
scanr (+) 2 [1, 2, 3, 4] = [12, 11, 9, 6]

2.4

Target Skeleton Programs

We formalize skeleton programs that are targets of this paper. Our targets are skeleton
programs that involve neighbor elements using combination of shift, shift, zip and map.
We also deal with skeleton programs that perform one accumulation or one reduction byscan (scanr) or reduce to the result of the above programs. We define the former program4

as Program, a program with accumulation as ProgramS , and a program with reduction as
ProgramR . A program of ProgramS can have arbitrary number of map after scan or scanr.
data Program α = map (β → α) (Program β)| shift α (Program α)| shift α (Program α)| zip (Program β) (Program γ)| [α]data ProgramS α = scan (α → α → α) α (Program α)| scanr (α → α → α) α (Program α)| map (β → α) (ProgramS β)
data ProgramR α = reduce (α → α → α) (Program α)Here, Program α and ProgramS α are programs that generate lists of which elements havetype α, while ProgramR α is a program that generates a value of type α. Type β in map isbound locally. Types β and γ in zip are bound by the relation α = (β, γ). For simplicity,
we do not use binding of variables in the above skeleton programs. Since the above skeleton programs are inputs for optimization algorithms, we distinguish skeletons in the above
skeleton programs from skeletons used in algorithms by attaching underlines.
Evaluation of the above defined program is given by evalP , evalPS and evalPR . Thisevaluation is the same as that removes the underlines of skeletons.
evalP :: Program α → [α]evalP (map f x) = map f (evalP x)
evalP (shift e x) = shift e (evalP x)
evalP (shift e x) = shift e (evalP x)
evalP (zip y z)
= zip (evalP y) (evalP z)
evalP (x)=xevalPS :: ProgramS α → [α]evalPS (scan (⊕) e x) = scan (⊕) e (evalP x)
evalPS (scanr (⊕) e x) = scanr (⊕) e (evalP x)
evalPS (map f x)
= map f (evalPS x)
evalPR :: ProgramR α → αevalPR (reduce (⊕) x) = reduce (⊕) (evalP x)
In the rest of this paper, we first focus on Program α to explain our idea. Then, we discuss
ProgramS α and ProgramR α in Section 6 as extensions of the optimization of Program α.
2.5

A Running Example

An example of our target programs, which involves neighbor elements, is computation of anumerical solution of differential equations by difference methods. We use a simple programfor difference method as our running example.
We consider a wave-equation as a concrete differential equation.
∂u
∂u
= −C
∂t
∂x
This equation describes propagation of waves. To calculate the propagation of waves, we
replace differential terms by differences. The value of u at time n and at location i is denotedby uni .

un+1
− uni
−C
i
a−2 uni−2 + a−1 uni−1 + a0 uni + a1 uni+1=∆t∆xwhere (a−2 , a−1 , a0 , a1 ) = (1/6, −1, 1/2, 1/3)
5

Here, we use a difference based on two elements on the left and an element on the right.
Rearranging the above equation, we get the following recurrence equation.
un+1= c−2 uni−2 + c−1 uni−1 + c0 uni + c1 uni+1i
A skeleton program next that computes the values at the next time form the current valuesof u in parallel is given as follows. Here, we consider simple boundary conditions: u n−1 = bl0 ,un0 = bl1 and unN +1 = br for a fixed N .
0 = map (c ×) (shift b (shift b u))
next u = let v−2
−2
l0
l1
0
v−1 = map (c−1 ×) (shift bl1 u)
v00 = map (c0 ×) u
v10 = map (c1 ×) (shift br u)00 )v− = map add (zip v−2v−1v+ = map add (zip v00 v10 )in map add (zip v− v+ )0 , v 0 , v 0 , and v 0 for readability, although the definition of
We use intermediate variables v−21−10target programs Program does not have variables. The correspondence of the program next
and the above recurrence equation is as follows. First, to generate a list corresponding to
the first term c−2 uni−2 , we apply two shift to shift the elements, and apply map to multiply
the coefficient c−2 . Similarly, lists corresponding to the second through the fourth terms are
generated by using shift, shift and map. Then, zipping these four lists by zip and adding
elements by map add , we obtain the final result. Since each of four lists are shifted by shiftand shift, ith element of the resulting list is c−2 ui−2 + c−1 ui−1 + c0 ui + c1 ui+1 .
In the rest of this paper, we explain our idea by using this example next.

2.6

A Strategy for Optimization

We propose specific optimization for our target skeleton programs according to the following
strategy. This strategy is based on the observation that domain-specific skeleton programs
are often developed with a fixed set of skeletons composed in a specific manner.
1. Design a normal form that abstracts target computations.
2. Develop fusion rules that transform a skeleton program into the normal form.
3. Implement the normal form efficiently in parallel.
In designing a normal form, we should have the following regards in mind. A normal form is
specified to describe any computation of target programs but should not be too general. A
normal form specific to the target programs enables us to develop efficient implementation
for it. In addition, a normal form should be closed under the fusion rules to maintain the
result of optimization in the form. Once we formalize a normal form with fusion rules
and efficient implementation, we can perform the optimization easily: we first transform askeleton program into the normal form with the fusion rules, and then we translate theprogram in the normal form to an efficient program.

3

Normal Form of Target Skeleton Programs

In this section, we formalize a normal form that describes the computation of Program αthat involves neighbor elements by using shift, shift, zip and map. The objective of thisnormal form is to perform main part of such computation in one loop.
6shift b0 (shift b1 u)

b0

b1

shift b1 u

b1

u

−
→
u[0]
−
→
u[1]

−
→
u[0]
−
→
u[1]
−
→
u[2]

shift b2 u

[[

[

l1

,

−−−→
u[i − 2]−−−→u[i − 1]
→
−
u [i]
−−−→
u[i + 1]l2 ]cecomputation withintroduced constantscommon computation
except for the index i

←
−
u[2]
←
−
u[1]
←
−
u[0]
b2

r1 ] ]][

computation withintroduced constants

Fig. 2. Regions of elements calculated by triples for next. A list of computational trees for elementson the left edge [l1 , l2 ]. A common computational tree for elements in the center part ce. A list ofcomputational trees for elements on the right edge [r1 ]. Elements involved in the computation are
boxed by dashed lines.

3.1

Example of Normal Form

Consider calculating the running example next in a single loop. By a simple observation,
its computation can be divided into three parts according to forms of computations ofelements: the computation of the center part elements where only elements of the given
list are involved, the computation of the left part where constants introduced by shift areinvolved, and the computation of the right part where constants introduced by shift areinvolved (shown in Fig. 2).
For the example next, each element in the center part is calculated by the following→−
computation. Here, u [i] denotes ith element form the left of u.
−−−→
−−−→
→
−
−−−→
ce = add (add (c−2 ×u[i − 2], c−1 ×u[i − 1]), add (c0 ×u [i], c1 ×u[i + 1]))This computation is common among the elements in the center part.
On the other hand, the leftmost element is calculated by the following.
−→−→
l1 = add (add (c−2 ×bl0 , c−1 ×bl1 ), add (c0 ×u[0], c1 ×u[1]))−→−→This computation involves variables (u[0] and u[1]) and constants (bl0 and bl1 ) introducedby shift. Similarly, the second element and the rightmost element involve constants. Here,←−u [i] denotes ith element from the right of u.
−
→
−
→
−
→
l2 = add (add (c−2 ×bl1 , c−1 ×u[0]), add (c0 ×u[1], c1 ×u[2]))
←
−
←
−
←
−
r1 = add (add (c−2 ×u[2], c−1 ×u[1]), add (c0 ×u[0], c1 ×br ))Summarizing these observations, whole computation of next can be denoted by thefollowing triple of computations.
– a list of computational trees for elements on the left edge [l1 , l2 ]– common computational tree for elements in the center part ce– a list of computational trees for elements on the right edge [r1 ]−−−→Figure 3 shows these computational trees. Here, u1 denotes the element u[i − 1] on theleft of the ith element. Similarly, u1 denotes the element on the right, and u2 denotes
the second element on the left. Figure 2 shows a region that each member of the triple
7




+
+

 c × c × c × c ×01
 −2 −1
b l0

b l1



+

+

−
→
u[0]

−
→
u[1]

+

+

c−2 × c−1 × c0 ×
b l1

−
→
u[0]

−
→
u[1]




c1 × 

+++c−2 × c−1 × c0 ×

c1 ×

u2 u1

u1

−
→
u[2]

u



+



+
+


 c × c × c × c × 01
 −2 −1

←
−
u[2]

←
−
u[1]

←
−
u[0]

br

Fig. 3. Triple describing computation of next.

calculates. To obtain the result from this triple, we calculate each element on the edges by its
computational tree, and calculate elements on the center part by the common computational
tree against indices in a single loop.
Generally, such a triple denotes a computation of a target skeletal program. Thus, inthe next section, we formalize this triple as a normal form of our target skeletal programs.
Transformation of a skeleton program into a normal form is shown in Section 4.
3.2

Definition of Normal Form

A normal form is defined as a triple of two lists of computational trees for elements on theedges and a common computational tree for the center part.
type NForm α = ([Tree α], Tree α, [Tree α])
We denote this triple by using special brackets like [[ ls, zms, rs ]]. Here, ls is the list ofcomputational trees for the left edge, zms is the common computational tree for the center
part, and rs is the list of computational trees for the right edge. Computational trees innormal forms are defined as follows.
data Tree α = Node ((β, γ) → α) (Tree β) (Tree γ)| Leafv (β → α) (Var β)| Leafc αdata Var α = Var [α] Int| Fix [α] Int Direction| Holedata Direction = FromL | FromR
Here, types β and γ are bound locally. A node of the tree holds left and right children
zipped by zip and a function applied by successive map. There are two kinds of leaves: Leaf cdenotes a constant introduced by shift or shift, and Leafv denotes a list variable and holds
a function applied by successive map. Access to a list is formalized by Var α: Var is used
is index access of the center part, and Fix is used in computations edge elements where
the indices are fixed. Thus, Var holds the list and the amount of shifting, while Fix holds
the list, the index and the origin of the index. The origin of the fixed index is specified by
Direction. A special value Hole is used in parallel implementation.
−→
For example, add (c−2 ×bl1 , c−1 ×u[0]) is described as follows.
Node add (Leafc (c−2 ×bl1 )) (Leafv (c−1 ×) (Fix u 0 FromL))
A part of common tree c−2 ×u2 is described as follows.
Leafv (c−2 ×) (Var u (−2))
Note that negative value of the amount of shifting means shifting to the left.
83.3

Semantics of Normal Form

We give a semantics of the normal form by sequential evaluation of the normal form. The
center part is calculated by a single loop with the common computational tree. Each element
on both edges is calculated by its own computational tree.
First, we define the evaluation of computational trees by evaluation function eval T0 and
evalT as follows. This evaluation uses an auxiliary function eval V defined below.
evalT0 :: Tree α → α
evalT0 x = evalT x 0
evalT :: Tree α → Int → α
evalT (Node f l r) i = f (evalT l i, evalT r i)
evalT (Leafv f v) i = f (evalV v i)
evalT (Leafc c) i
=c
The evaluation function evalT0 for elements on the left and right edges is defined by evalTignoring the index. The evaluation function evalT performs computation according to thedefinition of computational trees. The auxiliary function eval V processes index accessing of
input lists. The result of evalV for Hole is not defined.
evalV :: Var α → Int → α
evalV (Var u s) i
= at (i − s) u
evalV (Fix u s FromL) = at s u
evalV (Fix u s FromR) = at s (rev u)
Using the above evaluation functions, we define a sequential program eval that evaluates
the normal form. Here, the lengths of involved lists are supposed to be n.
eval :: NForm α → [α]eval [[ ls, zms, rs ]] = map evalT0 ls ++ map (evalT zms) idces ++ map evalT0 rs
where l = length ls ; r = length rs
idces = [l..(n − r − 1)]
Each element on both edges is calculated by its own computational tree using eval T0 defined
above. The center part is calculated by a single loop (map (eval T zms)) with the common
computational tree.

4

Fusion Rules for Transformation to Normal FormIn this section, we define rules to transform a skeleton program Program α into a normal
form.
4.1

Transformation with Fusion Rules

The transformation of a skeleton program into a normal form is performed by the following
function compile and fusion rules one by one.
compile :: Program α → NForm α
compile (map f x) = fuseMap f (compile x)compile (shift e x) = fuseShift  e (compile x)compile (shift e x) = fuseShift  e (compile x)compile (zip x y)
= fuseZip (compile x) (compile y)compile (x)= [[ [ ], Leafv id (Var x 0), [ ] ]]Fusion rules for each skeletons are defined as follows. Figure 4 illustrates the fusion rules.
9map fshift e⇓⇓(a) Fusion of map composes the given function to eachtree of the normal form. The structure of the normal
form is not changed.

(b) Fusion of shift discards the rightmost computational tree and introduces a new constant tree tothe leftmost. The amount of shifting in the tree forcenter part is updated by one.

zip
shift e⇓⇓(c) Fusion of shift discards the leftmost computational
tree and introduces a new constant tree to the rightmost. The amount of shifting in the tree for center
part is updated by −1.

(d) Fusion of zip instantiates the shaded parts so thatthe sizes of lists of trees for edges become uniform,and zips the corresponding trees.

Fig. 4. An image of fusion rules. Rectangles show the resulting lists. The three parts separated byvertical lines correspond to the triple of a normal form. Changed parts are shaded.

Fusion Rule for map Fusion of skeleton map is performed by composing the given function
to roots of computational trees.
fuseMap :: (α → β) → NForm α → NForm β
fuseMap f [[ ls, zms, rs ]] = [[ map (comp f ) ls, comp f zms, map (comp f ) rs ]]Composition of a function is defined by the following comp.
comp :: (α → β) → Tree α → Tree βcomp f (Node g l r) = Node (f ◦ g) l r
comp f (Leafv g v) = Leafv (f ◦ g) v
comp f (Leafc c)= Leafc (f c)
For non-constant roots, comp composes the given function to the function held in the root.
For constant roots (i.e. constant leaves), comp applies the given function to the constant togenerate a new constant root.
Fusion Rules for shiftfi and shiftfl Fusion of skeletons shift and shift is performed byinsertion and deletion of the leftmost and the rightmost trees, and update of the amount ofshifting.
fuseShift  :: α → NForm α → NForm α
fuseShift  e [[ ls, zms, rs ]] = [[ tail ls, slide (−1) zms, rs ++ [Leafc e] ]]fuseShift  :: α → NForm α → NForm α
fuseShift  e [[ ls, zms, rs ]] = [[ [Leafc e] ++ ls, slide 1 zms, init rs ]]For the left-shift, the fusion discards the leftmost computational tree and introduces aconstant computational tree to the rightmost. Then, it updates the amount of shifting in thecommon computational tree for the center part. This update is performed by the following
slide. The right-shift is similar to the left-shift.
slide :: Int → Tree α → Tree αslide d (Node f l r) = Node f (slide d l) (slide d r)
slide d (Leafv f (Var x s)) = Leafv f (Var x (s + d))
slide d x = x
10Fusion Rule for zip Fusion of skeleton zip needs unification of the lengths of lists tobe zipped. Thus, the common trees for the center parts are instantiated to expand thelists of computational trees for edges. Then, each pair of corresponding trees is zipped byintroducing a new node with id function.
fuseZip :: NForm α → NForm β → NForm (α, β)fuseZip [[ ls 1 , zms 1 , rs 1 ]] [[ ls 2 , zms 2 , rs2 ]]= let zms = Node id zms 1 zms 2
ls = trim FromL ls 1 ls 2 zms 1 zms 2
rs = trim FromR (rev rs 1 ) (rev rs 2 ) zms 1 zms 2
in [[ ls, zms, rev rs ]]The function trim defined below unifies the lengths of lists by instantiating common tree
zms k by a function insts defined below.
trim :: Direction → [Tree α] → [Tree β] → Tree α → Tree β → [Tree (α, β)]trim d ts 1 ts 2 zms 1 zms 2= let n1 = length ts 1 ; n2 = length ts 2
(ts 01 , ts 02 ) = (ts 1 ++ insts d zms 1 n1 n2 , ts 2 ++ insts d zms 2 n2 n1 )in zipWith (Node id) ts 01 ts 02
The instantiation is performed by the following insts and inst.
insts :: Direction → Tree α → Int → Int → [Tree α]
insts d zms s e = map (inst d zms) [s..(e − 1)]inst d (Node f l r) i = Node f (inst d l i) (inst d r i)
inst d (Leafv f (Var x s)) i = let s0 = case d of FromL → s; FromR → −sin Leafv f (Fix x (−s0 + i) d)
inst d x i = x
These four fusion rules and the base case rule can transform any skeleton program definedby Program into the normal form. We conclude this fact as a theorem.
Theorem 1. Any skeleton program defined by Program can be transformed into the normal
form by using the four fusion rules and the base case rule.
Proof. This is proved by induction on the structure of Program. The base case is shown by
the transformation of an input list. Induction cases are shown by the four fusion rules.
Complete proof is shown in Appendix Atu
4.2

Example of Transformation

As a brief explanation of the rules, we transform the example next into the normal form.
For readability, we use a brief notation used in the previous examples instead of the data
structures.
Most simplest case is the transformation of the argument list u. A list u needs only the→−common computational tree u [i] that is just the element of u.
u ⇒ [[ [ ], id u, [ ] ]]Here, id u is the brief notation of Leafv id (Var u 0).
Next, we transform shift bl1 u. This shift introduces the constant bl1 to the leftmost
element, a new computational tree of the constant bl1 is introduced to the normal form.
shift bl1 [[ [ ], id u, [ ] ]] ⇒ [[ [bl1 ], id u1 , [ ] ]]11

Also, the amount of shifting in the common tree is updated by 1.
Then, we fuse map (c−1 ×) to the above result.
map (c−1 ×) [[ [bl1 ], id u1 , [ ] ]] ⇒ [[ [c−1 ×bl1 ], c−1 ×u1 , [ ] ]]The constant bl1 is replaced by c−1 × bl1 , and the function c−1 × is composed to id held in
the root of the common tree. Since id is the identity of function composition, id is removed
in the result.
Similarly, other applications of shift, shift and map result in the following normal
forms.
map (c0 ×) u ⇒ [[ [ ], c0 × u, [ ] ]]map (c1 ×) (shift br u) ⇒ [[ [ ], c1 × u1 , [c1 × br ] ]]map (c−2 ×) (shift bl0 (shift bl1 u)) ⇒ [[ [c−2 × bl0 , c−2 × bl1 ], c−2 × u2 , [ ] ]]In the second transformation, a new tree is introduced by shift to the rightmost, andthe amount of shifting in the common tree is updated by 1 to the left. After the lasttransformation, there are two computational trees for elements on the left edges introducedby two shift, and the amount of shifting in the common tree becomes 2.
Next, we perform fusion of zip to transform zip v00 v10 , i.e.
zip [[ [ ], c0 × u, [ ] ]] [[ [ ], c1 × u1 , [c1 × br ] ]].
Since the lengths of lists of two normal forms to be zipped are not the same (i.e. [ ] and[c1 × br ] for the right edges), we have to unify the lengths by instantiating the common trees.
Instantiation means to fix the indices in the common trees for elements on the edges. The
result of instantiation and zip is as follows.
←−[[ [ ], (c0 × u, c1 × u1 ), [(c0 × u[0], c1 × br )] ]]←−Here, instantiation of the common tree of the first normal form c0 ×u results in c0 ×u[0], and
it is zipped with the rightmost tree of the second normal form to make the new rightmost
tree.
0
0 .
Similarly, we obtain the following normal form from zip v−2
v−1−→[[ [(c−2 × bl0 , c−1 × bl1 ), (c−2 × bl1 , c−1 × u[0])], (c−2 × u2 , c−1 × u1 ), [ ] ]]Continuing these fusions, we finally obtain the following normal form for the example
next.
−
→
−
→
[[[add (add (c−2 × bl0 , c−1 × bl1 ), add (c0 × u[0], c1 × u[1])),
−
→
−
→
−
→
add (add (c−2 × bl1 , c−1 × u[0]), add (c0 × u[1], c1 × u[2]))],add (add (c−2 × u2 , c−1 × u1 ), add (c0 × u, c1 × u1 )),
←
−
←
−
←
−
[add (add (c−2 × u[2], c−1 × u[1]), add (c0 × u[0], c1 × br ))] ]]5

Parallel Implementation of Normal FormIn this section, we give a parallel implementation of the normal form.
12

(1)
(2)

.

&

dist

⇓ eval 0⇓ eval 0

(3)

&

(4)

. &

.
combine(1) Lists in a normal form are divided into two parts to generate two normal forms. Then, these normal
forms are distributed to processors.
(2) Each processor calculates elements that can be calculated by local data. Elements on the edges
cannot be calculated because they need data on other processors.
(3) To calculate the elements on the edges, communication is performed to complete the data.
(4) Perform calculation of the elements on the edges. Combine the edge elements to the other elements.

Fig. 5. An image of parallel implementation of the normal form (two processors).

5.1

Parallel Implementation of Normal Form

Based on parallel implementation of existing skeletons [8], we consider parallel implementation of the normal form by four steps: (1) distribution of a normal form (input lists), (2)the first local computation, (3) global communication, (4) the second local computation.
We explain the idea of parallel implementation using the example next. Figure 5 shows animage of parallel implementation of the normal form using two processors.
First, we distribute the normal form among processors. The input list u is divided into
two parts: u = u1 ++ u2 . There are two processors, and each processor has a part of thedivided list. Letting the original normal form be [[ ls, zms, rs ]], each processor has one of thefollowing distributed normal forms.
nf 1 = [[ ls, zms 1 , [ ] ]] ; nf 2 = [[ [ ], zms 2 , rs ]]Here, the list of computational trees for the left edge is held in the first normal form nf 1 ,while the list of computational trees for the right edge is held in the last normal form
nf 2 . The common computational tree zms k of the normal form nf k is created from zms
by replacing the input list u with the part of the list uk . Thus, for the example next, thedistributed normal form nf k is as follows.
add (add (c−2 ×uk2 , c−1 ×uk1 ), add (c0 ×uk , c1 ×uk1 ))Next, in the first local step, processors calculate own partial results in parallel. Since
elements on the edges of the distributed normal forms needs elements of both u1 and u2 , theseelements cannot be calculated in this phase. For example, the rightmost element calculated
by nf 1 needs an element of u2 , which is underlined in the following expression.
←
−
←
−
←
−
−
→
add (add (c−2 ×u1 [2], c−1 ×u1 [1]), add (c0 ×u1 [0], c1 ×u2 [0]))Since these elements cannot be calculated in this local phase, computational trees for theseelements are held until the global computation phase. Here, a hole • is substituted to anelement that is not available in this phase.
←
−
←
−
←
−
add (add (c−2 ×u1 [2], c−1 ×u1 [1]), add (c0 ×u1 [0], c1 ×•))Similarly, the left edge of nf 2 generates the following computational tree with holes.
−→
add (add (c−2 ×•, c−1 ×•), add (c0 ×•, c1 ×u2 [0]))13Note that this tree complement the tree generated from nf 1 . Generally, letting the maximum
amounts of shifting to the left and the right be l and r, the number of computational trees
with holes generated on each edge is l + r. Trees on both sides of an edge are complementary
to each other.
Third, in the global communication step, neighboring processors communicate incomplete trees with holes to each other to complete those trees. Of course, the first local computation can hide the time of this communication phase.
Fourth, in the second local step, each processor calculates the elements on the edges
with the completed trees to complete the resulting distributed list.
For the example next, nf 1 generates the following trees with holes,
←
−
←
−
←
−
add (add (c−2 ×u1 [2], c−1 ×u1 [1]), add (c0 ×u1 [0], c1 ×•)),←−←−add (add (c−2 ×u1 [1], c−1 ×u1 [0]), add (c0 ×•, c1 ×•)),←−add (add (c−2 ×u1 [0], c−1 ×•), add (c0 ×•, c1 ×•)),while nf 2 generates the following trees.
−→add (add (c−2 ×•, c−1 ×•), add (c0 ×•, c1 ×u2 [0]))−→−→add (add (c−2 ×•, c−1 ×•), add (c0 ×u2 [0], c1 ×u2 [1]))
−
→
−
→
−
→
add (add (c−2 ×•, c−1 ×u2 [0]), add (c0 ×u2 [1], c1 ×u2 [2]))Zipping these trees, we obtain the following complete computational trees for elements on
the edge.
←
−
←
−
←
−
−
→
add (add (c−2 ×u1 [2], c−1 ×u1 [1]), add (c0 ×u1 [0], c1 ×u2 [0]))
←
−
←
−
−
→
−
→
add (add (c−2 ×u1 [1], c−1 ×u1 [0]), add (c0 ×u2 [0], c1 ×u2 [1]))
←
−
−
→
−
→
−
→
add (add (c−2 ×u1 [0], c−1 ×u2 [0]), add (c0 ×u2 [1], c1 ×u2 [2]))After these four steps, these completed results are gathered to the root processor, or becomea new input to another normal form. Distribution of input will be skipped in the latter case.
Summarizing the above ideas, we get the following parallel implementation of the normal
form.
eval = globalReduction ◦ map (localEval ) ◦ distFirst, a normal form is distributed among processors by dist (step (1)). Next, each processor
evaluates a part of distributed normal forms by localEval (step (2)). Then, these results are
combined globally by globalReduction (step (3) and (4)). Local computation localEval takes
a distributed normal form, and generates a triple of its partial result (pls, cs, prs), where cs
is the calculated elements on the center part, and pls and prs are computational trees withholes for elements on left and right edges. Global computation globalReduction communicatespartial results (pls 1 , cs 1 , prs 1 ) and (pls 2 , cs 2 , prs 2 ) of both sides of an edge, then generates a
new partial result (pls 1 , cs 1 ++glue prs 1 pls 2 ++cs 2 , prs 2 ). Here, glue denotes the complement
process. Then, the final result is obtained by extracting the center result of the final partial
result.
5.2

Formalization of Parallel Implementation of Normal Form

We formalize the parallel implementation explained in the previous section. In the rest of
this section, p means the number of processors and n means the length of the lists involvedin the computation.
14

First, we define the distribution dist that divides a normal form into p parts.
dist :: Int → NForm α → [NForm α]
dist p [[ ls, zms, rs ]] = let divs = division p nzmss = distribute zms divs
lss = ls : take (p − 1) (repeat [ ])rss = take (p − 1) (repeat [ ]) ++ [rs]in zip3 lss zmss rss
Here, division p n calculates the division of input lists, and distribute zms divs distributes
the input lists according to the division. This process generates distributes a computational
tree zms i that is generated by replacing a list u (= u1 ++ · · · ++ ui ++ · · · ++ up ) in the originaltree zms with ui , and zms i is distributed to ith processor and held in the normal form nf i
(this is generated by zip3 lss zmss rss). The list of computational trees for the left edge isheld in the first normal form nf 1 , while the list of computational trees for the right edge isheld in the last normal form nf p .
Next, we define the triple of the partial result PResult.
type PResult α β = ([Tree α], β, [Tree α])
We use a special brackets hh pls, cs, prs ii to denote a partial result. Here, cs is the calculated
elements on the center part, and pls and prs are computational trees with holes for elements
on left and right edges. We abstract the type of calculated elements on the center part as β
for generality.
Next, we define an auxiliary function maxShift that calculates the maximum amount ofshifting.
maxShift :: Tree α → (Int, Int)
maxShift (Node l r) = let (l1 , r1 ) = maxShift l(l2 , r2 ) = maxShift rin (max l1 l2 , max r1 r2 )maxShift (Leafv (Var s)) = if s < 0 then (−s, 0) else (0, s)= (0, 0)maxShift
We define a general function gEval 0 to perform the local computation.
gEval 0 :: (NForm α → [Int] → β) → NForm α → PResult α β
gEval 0 fc [[ ls, zms, rs ]] = hh pls, cs, prs ii
whereidces = [r..(n − l − 1)](l, r) = maxShift zms
cs = fc (ls, zms, rs) idces
pls = map (instP FromL zms) [(−l)..(r − 1)]
prs = map (instP FromR zms) (rev [(−r)..(l − 1)])This gEval 0 takes a function to calculate the center part of the partial result, and generates
a partial result consisting of the result computed by the given function and two lists ofcomputational trees with holes generated by the following function instP .
instP :: Direction → Tree α → Int → Tree α
instP d (Node f l r) i = Node f (instP d l i) (instP d r i)
instP d (Leafv f (Var x s)) i = let s0 = case d of FromL → s; FromR → −si0 = (−s0 + i)
in if (i0 ≤ (length x))||(i0 < 0) then Leafv f Holeelse Leafv f (Fix x i0 d)
instP d x i = x
15This instP is the same of inst except that it substitutes a hole for an element that is notavailable in this phase.
Using gEval 0 , we define the function eval 0 that performs local computation of the normal
form as follows.
eval 0 :: NFormα → PResultα[α]eval 0 = gEval 0 fcMap
where fcMap [[ ls, zms, rs ]] idces= map evalT0 ls ++ map (evalT zms) idces ++ map evalT0 rs
The defined function fcMap evaluates the computational tree zms against indices given by
gEval 0 , and returns a list of the resulting elements.
Next, we define the process of completion of a computational tree from two trees complementary to each other.
combine :: Tree α → Tree α → Tree α
combine (Node f1 l1 r1 ) (Node f2 l2 r2 ) = Node f1 (combine l1 l2 ) (combine r1 r2 )combine (Leafv f1 Hole) (Leafv f2 v)= Leafv f2 v
combine (Leafv f1 v) (Leafv f2 Hole)) = Leafv f1 v
combine (Leafc c) (Leafc c0 )= Leafc cThis function completes the tree by filling holes with values held in the other tree. Using
this auxiliary function, we define an operator (·,·) for global reduction.
((·,·) ) :: (α → α → α) → (β → α) → PResult β α → PResult β α → PResult β αhh pls 1 , cs 1 , prs 1 ii (⊕,f ) hh pls 2 , cs 2 , prs 2 ii = hh pls 1 , cs, prs 2 ii
where es = zipWith ((evalT0 ◦) ◦ combine) prs 1 pls 2
cs = cs 1 ⊕ cata (⊕) f es ⊕ cs 2This operator takes an associative operator and a function, and calculates the reduction of
the elements on edges. This reduction is done by the following general function.
cata :: (α → α → α) → (β → α) → [β] → αcata (⊕) f [ ]= ı⊕cata (⊕) f (a : x) = f a ⊕ cata (⊕) f x
Since the computation of the normal form generates a list, we use (++,[·]) that performs
reduction by ++ and [·].
Using the functions defined above, we define the parallel implementation of normal form
parEval as follows.
parEval :: Int → NForm α → [α]parEval p=extract ◦ reduce ((++,[·]) ) ◦ map eval 0 ◦ dist p
The last function extract extracts the center value from the final partial result.
extract :: PResult α β → β
extract hh pls, cs, prs ii = cs

6

Expansion of Target Programs

Based on the results for Program, we expand our target programs to Program S , which
includes accumulation by scan or scanr, and ProgramR , which includes reduction by reduce(see Section 2).
16

6.1

Target Programs with Accumulation

First, we expand our target programs to ProgramS . A target skeleton program of ProgramS
has one accumulation by scan or scanr after the computation (i.e. Program) that involves
neighbor elements using combination of shift, shift, zip and map. The target skeleton
program also has arbitrary number of map after the accumulation.
An Example of Target Programs with Accumulation As an example of the target
programs, consider a program to solve a tridiagonal linear system of equations. The outputof the program is xs = [x1 , . . . , xn ] that satisfies the following linear equations for the given
coefficients ds = [d1 , . . . , dn ], es = [e1 , . . . , en ], fs = [f1 , . . . , fn ], bs = [b1 , . . . , bn ].

 


b1
x1
d1 f1
 e2 d2 f2
 x2   b2 

 


 e3 d3 f3
 x2   b2 

 



 ..  =  .. 
.. ..
..






.
.
.

 .   . en−1 dn−1 fn−1 xn−1  bn−1 
bn
xn
en dn
|

{z

}|

A

{z
x

}

=

| {z }b

A skeleton program that solves this tridiagonal linear system of equations is given as follows.
This program is based on the LU decomposition of the coefficient matrix A = LU [11].
solveTS ds es fs bs= let us = map g2 (scan (⊗1 ) eye (map g1 (zip ds (map mul (zip es (shift 0 fs))))))
ms = map div (zip es (shift ∞ us))ys = map g5 (scan (⊗2 ) eye 0 (map g3 (zip ms bs)))xs = map g5 (scanr (⊗3 ) eye 0 (map g4 (zip us (zip ys fs))))in xs
where (a11 , a12 , a21 , a22 ) ⊗1 (b11 , b12 , b21 , b22 )= (a11 ∗ b11 + a12 ∗ b21 , a11 ∗ b12 + a12 ∗ b22 ,a21 ∗ b11 + a22 ∗ b21 , a21 ∗ b12 + a22 ∗ b22 )(a11 , a21 ) ⊗2 (b11 , b21 ) = (a11 ∗ b11 , a21 ∗ b11 + b21 )(a11 , a12 ) ⊗3 (b11 , b12 ) = (a11 ∗ b11 , a11 ∗ b12 + a12 )
g1 (d, ef ) = (d, 1, −ef , 0)
g2 (a11 , a12 , a21 , a22 ) = a11 /a12
g3 (m, b) = (−m, b)
g4 (u, (y, f )) = (−f /u, y/u)g5 (a11 , a12 ) = a12
mul a b = a ∗ b
div a b = a/beye = (1, 0, 0, 1)
eye 0 = (1, 0)
Here, us corresponds to the upper triangular matrix U , ms corresponds to the lower triangular matrix L, ys is the solution of the linear equation Lys = bs, and xs the solution of thelinear equation Axs = bs. The operators ⊗1 , ⊗2 and ⊗3 are multiplication of 2 × 2 matrices,
although ⊗2 and ⊗3 omit the half of the elements. This program consists of three parts:
the first part performs the LU decomposition with scan to obtain us and ms, the second
part performs the forward substitution by scan to obtain ys, and the third part performs
the backward substitution by scanr to obtain xs.
17Normal Form with Accumulation We extend the normal form of Program to hold theoperator, the direction and the element of the accumulation, and the functions applied by
the last map. Thus, the extended normal form NFormS is defined as follows.
type NFormS α = (NForm β, (Direction, β → β → β, β, β → α))
In the extended normal form ([[ ls, zms, rs ]], (d, ⊕, e, f )), [[ ls, zms, rs ]] specifies the computation before the accumulation, d is the direction of accumulation, ⊕ is the associative binary
operator used in the accumulation, e is the initial element of the accumulation, and f is afunction applied to each element after the accumulation. The direction d is FromL for theaccumulation by scan, and FromR for scanr.
For example, the example solveTS is described by the following three normal forms.
us ⇒ ([[ [g1 (ds , mul (es , 0))], g1 (ds , mul (es , fs 1 )), [ ] ]], (FromL, ⊗1 , eye, g2 ))ys ⇒ ([[ [g3 (div (es , ∞), bs )], g3 (div (es , us 1 ), bs ), [ ] ]], (FromL, ⊗2 , eye 0 , g5 ))xs ⇒ ([[ [ ], g4 (us , (ys , fs )), [ ] ]], (FromR, ⊗3 , eye 0 , g5 ))Note that the computation of ms is absorbed by the computation of ys.
Fusion Rules with Accumulation We extend the fusion rules to handle scan, scanr andmap after accumulation. The transformation of a skeleton program into an extended normal
form is performed by the following function compileS and fusion rules one by one.
compileS :: ProgramS α → NFormS α
compileS (scan (⊕) e x) = fuseScan (⊕) e (compile x)compileS (scanr (⊕) e x) = fuseScanr (⊕) e (compile x)compileS (map f x)
= fuseMapS f (compileS x)
Transformation of the computation before the accumulation is done by compile defined in
Section 4.
Fusion rules for accumulations are as follows.
fuseScanr :: (α → α → α) → α → NForm α → NFormS α
fuseScanr (⊕) e [[ ls, zms, rs ]] = ([[ ls, zms, rs ]], (FromR, (⊕), e, id))
fuseScan :: (α → α → α) → α → NForm α → NFormS α
fuseScan (⊕) e [[ ls, zms, rs ]] = ([[ ls, zms, rs ]], (FromL, (⊕), e, id))
Each rule appends the operator ⊕ and the initial element e to the normal form, and marks
the direction of the accumulation by FromL or FromR. The last element of the extended
normal form is initialized by the identity function id.
The fusion rule for map after accumulation is as follows.
fuseMapS :: (β → α) → NFormS β → NFormS α
fuseMapS f = (nf , (d, (⊕), e, g)) = (nf , (d, (⊕), e, f ◦ g))
This rule merely composes the given function f to the function g in the extended normalform.
It is obvious that any target program can be transformed into an extended normal formby the above fusion rules and the fusion rules defined in Section 4.
18

Parallel Implementation of Normal form with Accumulation Based on parallel
implementation of skeleton scan, we design parallel implementation of the normal form.
Parallel implementation of skeleton scan is given as follows.
scan (⊕) e =reduce (++) ◦ zipWithP (postScan (⊕))◦ ((prescan (⊕) e ◦ map last)4id)◦ map (scan (⊕) ı⊕ ) ◦ dist p
where postScan (⊕) a x = map (a⊕) xHere, zipWithP f (x, y) = zipWith f x y and (f 4g) x = (f x, g x). First, this implementation
performs local accumulation in parallel by scan. Then, it performs global accumulation byprescan defined below. Finally, postScan adds the accumulated value calculated by prescan
to each element of the result of the local accumulation.
prescan :: (α → α → α) → α → [α] → [α]prescan (⊕) e [ ]= []prescan (⊕) e (a : x) = e : prescan (⊕) (e ⊕ a) x
Based on the parallel implementation of scan, we consider the following parallel implementation of a normal form with accumulation. We only show the implementation of thenormal form with accumulation by scan. The implementation for accumulation by scanr issimilar to that by scan.
parEvalScan :: Int→NFormS α→[α]parEvalScan p (nf , (FromL, ⊕, e, f )) = reduce (++) ◦ zipWithP (evalPostScan (⊕) f )◦ ((prescan ((⊕,id) ) ([ ], e, [ ]) ◦ map takeLast)4id)◦ map (evalScan 0 (⊕) ı⊕ ) ◦ dist p $ nf
Here, takeLast hh pls, cs, prs ii = hh pls, last cs, prs ii. Basic structure is the same as the implementation of scan. Main difference is that local and global computation deal with triples
of partial results defined. The local computation evalScan 0 that generates a partial result isdefined as follows. The generated partial result holds the result of local accumulation in thecenter of the triple.
evalScan 0 :: (α → α → α) → α → NForm α → PResult α [α]evalScan 0 (⊕) e = gEval 0 fcScan
where fcScan [[ ls, zms, rs ]] idces = sls ++ scs ++ srs
where sls = scata (⊕) e evalT0 ls
e0 = last (e : sls)
scs = scata (⊕) e0 (evalT zms) idces
e00 = last (e0 : scs)
srs = scata (⊕) e00 evalT0 rs
Here, scata defined below performs accumulation and evaluation of the computational tree
at the same time.
scata :: (α → α → α) → α → (β → α) → [β] → [α]scata (⊕) e f [ ]= []scata (⊕) e f (a : x) = let e0 = (e ⊕ f a)in e0 : scata (⊕) e0 f x
19

The final local computation evalPostScan (⊕) calculates accumulation of elements on theedges and adds the accumulated value to the result of the local computation.
evalPostScan :: (α → α → α) → (α → β)PResult α α → PResult α [α] → [β]evalPostScan (⊕) f hh epls, e, eprs ii hh pls, cs, prs ii = res
where fe (epr , pl ) = evalT0 (combine epr pl ))sls = scata (⊕) e fe (zip eprs pls)
e0 = last (e : sls)
scs = map (f ◦ e0 ⊕) cs
res = map f sls ++ scs
6.2

Target Programs with Reduction

Next, we expand our target programs to ProgramR .
Normal Form with Reduction We extend the normal form of Program to hold theoperator of the reduction. Thus, the extended normal form NFormR is defined as follows.
type NFormR α = (NForm β, β → β → β)In the extended normal form ([[ ls, zms, rs ]], ⊕), [[ ls, zms, rs ]] specifies the computation before
the accumulation, and ⊕ is the associative binary operator used in the reduction.
Fusion Rules with Reduction We extend the fusion rules to handle reduce. The transformation of a skeleton program into an extended normal form is performed by the following
function compileR and fusion rules one by one.
compileS :: ProgramS α → NFormS α
compileR (reduce (⊕) x) = fuseReduce (⊕) (compile x)
Transformation of the computation before the reduction is done by compile defined in Section 4.
The fusion rule for reduction is as follows.
fuseReduce :: (α → α → α) → NForm α → NFormR α
fuseReduce (⊕) [[ ls, zms, rs ]] = ([[ ls, zms, rs ]], (⊕))
The rule appends the operator ⊕ to the normal form.
It is obvious that any target program can be transformed into an extended normal formby the above fusion rule and the fusion rules defined in Section 4.
Parallel Implementation of Normal Form with Reduction In this section, we give
parallel implementation for ProgramR that perform reduction by reduce to the result of anormal form.
These programs can be executed in parallel by performing reduction, instead of generation of a list, in the implementation of a normal form. That is, parallel implementation is
obtained by giving the reduction operator to gEval 0 and (·,·) .
The function evalReduce 0 that performs local reduction is defined as follows.
evalReduce 0 :: (α → α → α) → NForm α → PResult α αevalReduce 0 (⊕) = gEval 0 fcReduce
where fcReduce [[ ls, zms, rs ]] idces= cata (⊕) evalT0 ls ⊕ cata (⊕) (evalT zms) idces ⊕ cata (⊕) evalT0 rs
20

Table 1. Running times and speedups of parallel programs against the number of processors. A
speedup is one with respect to a sequential program.
#processors
nexttime (s)speedupnext opt time (s)speedup

1
2
4
8
16
24
32
48
64
210.25 100.84 48.12 24.41 13.31 8.86 6.52 4.70 3.50
0.09
0.20 0.41 0.81 1.49 2.24 3.04 4.23 5.67
19.86
9.64 4.93 2.44 1.26 0.87 0.70 0.54 0.47
1.00
2.06 4.03 8.14 15.79 22.76 28.26 36.73 42.22

r[0] = c_2*b0 + c_1*b1 + c0*u[0] + c1*u[1];r[1] = c_2*b1 + c_1*u[0] + c0*u[1] + c1*u[2];for(int i = 2; i < n-1; i++) {r[i] = c_2*u[i-2] + c_1*u[i-1] + c0*u[i] + c1*u[i+1];}r[n-1] = c_2*u[n-3] + c_1*u[n-2] + c0*u[n-1] + c1*b2;
Fig. 6. A sequential program of next

This function performs reduction and evaluation of computational trees at the same timeby cata (⊕) evalT0 and cata (⊕) (evalT zms). Using this function, parallel implementationof a normal form with reduction is defined as follows.
parEvalReduce :: Int → NFormR α → α
parEvalReduce p (nf , ⊕) = extract ◦ reduce ((⊕,id) ) ◦ map (evalReduce 0 (⊕)) ◦ dist p $ nf

7

Experiment Result

We implemented a small domain-specific optimizer for the target programs. The system
reads a skeleton program written with our skeleton library SkeTo [8], and generates anoptimized C++ code of the program. For the examples next and solveTS , we measured
running times of skeleton programs and optimized programs. We used a PC cluster where
each of the nodes connected with Gigabit Ethernet has a CPU of Intel R Xeon R 2.80GHz
and 2GB memory, with Linux 2.4.21, GCC 4.1.1, and mpich 1.2.7.
7.1

Result of next

Table 1 shows measured running times and speedups. Running time is of applying next 100
times to an input list of 10,000,000 elements. A speedup is a ratio of running time of aparallel program to running time of a sequential program (shown in Figure 6).
The optimized program (next opt) achieves ten times faster running time than the original skeleton program, and the same running time as a sequential program on one processor.
This improvement was gained by elimination of redundant intermediate data and communications. Also, the optimized program achieves good speedups against the number ofprocessors. These results show effectiveness of the proposed optimization.
7.2

Result of solveTS

Table 2 shows measured running times and speedups. Running time is of applying solveTS
10 times to an input list of 1,000,000 elements. A speedup is a ratio of running time of aparallel program to running time of solveTS on one processor.
The optimized program (solveTS opt) achieves about 20% faster running time than theoriginal skeleton program. This improvement was gained by elimination of redundant intermediate data and communications. Also, the optimized program achieves good speedups
against the number of processors. These results show effectiveness of the proposed optimization.
21

Table 2. Running times and speedups of parallel programs against the number of processors. A
speedup is one with respect to solveTS on one processor.
#processors
1
2
4
8
16
32
64
solveTStime (s) 117.03 69.92 40.54 20.21 10.14 4.97 2.79
speedupp 1.00 1.67 2.89 5.79 11.54 23.56 41.92
solveTS opt time (s) 62.20 60.81 30.47 15.26 7.66 3.91 2.18
speedupp 1.88 1.92 3.84 7.67 15.28 29.92 53.61

8

Discussion and Related WorkOne of the simplest fusion optimizations so far uses a general form called cataJ [9]. This cataJ
can describe the skeleton map and the reduction skeleton reduce, and consists of two things:
a function applied to each element of the input list, and an associative binary operator used
to reduce those elements. The general form cataJ can describe any computation written ascomposition of any number of map and at most one reduce at the last. Thus, cataJ is thought
to be a normal form of skeleton programs described with such compositions.
Hu et al. [7] proposed a more general fusion optimization using a general form called
accumulate with some fusion rules. This accumulate can describe skeleton scan, which calculates an accumulation of the input list with an associative binary operator, as well as
map and reduce. The general form accumulate essentially consists of two pairs of a function
and an associative binary operator, and it is thought to be a normal form of skeleton programs described with compositions of the skeletons. This accumulate is powerful so that it
can describe shift and shift too. However, describing shift and shift with accumulatecause some overheads due to lack of special consideration of elements on the edges. The
overheads are extensions of elements to treat all elements uniformly, and logarithmic steps
of interprocessor communications for general implementation of accumulation. Thus, weneed to consider a specific fusion optimization, i.e. a normal form, fusion rules and efficient
implementation.
Our normal form extends these fusion optimizations with special consideration of elements on the edges, which are introduced by shift and shift. The normal form separates
computation of elements on the edges from that of center elements, so that it does not
introduce the overheads accumulate needed.
Grelck et al. [6] proposed an optimization based on fusion of a general skeleton called
WITH-loop. Since their WITH-loop are based on index accessing, their method can deal
with shifting operation such as shift and shift. Also, their fusion can handle nesting useof skeletons. However, skeletons that perform accumulation or reduction, in which a regionof elements required in computation of an element varies, cannot be described with WITHloop.
Generalized scan (gen scan) proposed by Fischer et al. [3] can perform accumulation thatinvolves neighbor elements. We can describe shifting operations such as shift and shiftwith their gen scan. We think programs that perform accumulation after combinations ofmap, zip, shift and shift can be described by gen scan. However, since the condition forparallel implementation of gen scan is very complicated, formalization of fusion rules ofgen scan is very difficult.

9

Conclusion

In this paper, we proposed an optimization technique for skeletons programs that involve
neighbor elements, which cannot gain adequate efficiency with existing fusion techniques.
We give a normal form that abstracts computations that involve neighbor elements, fusion22

rules of skeletons into a normal form, and an efficient parallel implementation of a normal
form. The optimization is performed by transforming a skeleton program with fusion rules
into a normal form that has efficient implementation. Experimental results show the optimized program can be executed efficiently in parallel, and show effectiveness of proposed
optimization. We consider extension of this technique to multidimensional data structuresand nesting use of skeletons as our future work.

References
1. R. S. Bird. Introduction to Functional Programming using Haskell. Prentice Hall Series in
Computer Science. Prentice Hall, 2nd edition, April 1998.
2. M. Cole. Algorithmic Skeletons: Structural Management of Parallel Computation. Research
Monographs in Parallel and Distributed Computing. MIT Press, 1989.
3. P. F. Fischer, F. P. Preparata, and J. E. Savage. Generalized scans and tridiagonal systems.
Theoretical Computer Science, 255(1–2):423–436, 2001.
4. A. J. Gill, J. Launchbury, and S. L. P. Jones. A short cut to deforestation. In FPCA ’93
Conference on Functional Programming Languages and Computer Architecture. Copenhagen,Denmark, 9–11 June 1993, pages 223–232. ACM Press, 1993.
5. S. Gorlatch, C. Wedler, and C. Lengauer. Optimization rules for programming with collective
operations. In 13th International Parallel Processing Symposium / 10th Symposium on Paralleland Distributed Processing (IPPS / SPDP ’99), 12-16 April 1999, San Juan, Puerto Rico,
Proceedings, pages 492–499. IEEE Computer Society, 1999.
6. C. Grelck and S.-B. Scholz. Merging compositions of array skeletons in SaC. Parallel Computing,32(7–8):507–522, 2006.
7. Z. Hu, H. Iwasaki, and M. Takeichi. An accumulative parallel skeleton for all. In D. L. Métayer,
editor, Programming Languages and Systems, 11th European Symposium on Programming,
ESOP 2002, held as Part of the Joint European Conference on Theory and Practice of Software,
ETAPS 2002, Grenoble, France, April 8-12, 2002, Proceedings, volume 2305 of Lecture Notes in
Computer Science, pages 83–97. Springer, 2002.
8. K. Matsuzaki, H. Iwasaki, K. Emoto, and Z. Hu. A library of constructive skeletons for sequential
style of parallel programming. In InfoScale ’06: Proceedings of the 1st international conference on
Scalable information systems, volume 152 of ACM International Conference Proceeding Series,
page 13. ACM Press, 2006.
9. K. Matsuzaki, K. Kakehi, H. Iwasaki, Z. Hu, and Y. Akashi. A fusion-embedded skeleton
library. In M. Danelutto, M. Vanneschi, and D. Laforenza, editors, Euro-Par 2004 Parallel
Processing, 10th International Euro-Par Conference, Pisa, Italy, August 31 – September 3, 2004,
Proceedings, volume 3149 of Lecture Notes in Computer Science, pages 644–653. Springer, 2004.
10. F. A. Rabhi and S. Gorlatch, editors. Patterns and Skeletons for Parallel and Distributed
Computing. Springer-Verlag, 2002.
11. H. S. Stone. An efficient parallel algorithm for the solution of a tridiagonal linear system of
equations. Journal of the ACM, 20(1):27–38, 1973.
12. P. Wadler. Deforestation: Transforming programs to eliminate trees. In H. Ganzinger, editor,
ESOP ’88, 2nd European Symposium on Programming, Nancy, France, March 21-24, 1988,
Proceedings, volume 300 of Lecture Notes in Computer Science, pages 344–358. Springer, 1988.
13. C. Wedler and C. Lengauer. On linear list recursion in parallel. Acta Informatica, 35(10):875–909,
October 1998.

23

A

Proof of Theorem 1

In the following, we use an abbreviation of the index accessing for readability: hxi i = at i x.
We assume that the length of an input list is n. The goal of this proof is to show the equation
evalP prog = eval (compile prog) for any prog.
A.1

Base Case

What we have to prove is the following equation for a list x.
evalP x = eval [[ [], Leafc id (Var x 0), [] ]](1)

The ith element of the left-hand sides is as follows.
hevalP xii ={ definition of evalP }hevalP xii

The ith element of the right-hand sides is as follows.
heval [[ [], Leafc id (Var x 0), [] ]]ii ={ definition of eval }hmap evalT (Leafc id (Var x 0)) [0..n − 1]ii
= { ith element }evalT (Leafc id (Var x 0)) i= { definition of evalT and the identity function }evalV (Var x 0) i
= { definition of evalV and ith element }hxii

Thus, the following equation holds.
hevalP xii = heval [[ [], Leafc id (Var x 0), [] ]]ii
Since this equation holds for i ∈ [0..n − 1], the equation (1) holds.
A.2

Inductive Case for map

What we have to prove is the following equation for a function f and a program prog. Here,compile prog = [[ ls, zms, rs ]].
evalP (map f prog) = eval [[ map (comp f ) ls, comp f zms, map (comp f ) rs ]](2)

To prove this equation, we first show a lemma.
Lemma 1. Let f be a function, t be a computational tree, and i be an index. Then, thefollowing equation holds.
evalT (comp f t) i = f (evalT t i)
Proof. This is shown by induction on Tree.
The Node case:{ definition of comp }evalT (Node (f ◦ g) l r) i
= { definition of evalT }(f ◦ g) (evalT l i, evalT r i)
= { function composition }f (g (evalT l i, evalT r i))= { definition of evalT }f (evalT (Node g l r) i)

evalT (comp f (Node g l r)) i =24

The Leafv case:{ definition of comp }evalT (Leafv (f ◦ g) v) i
= { definition of evalT }(f ◦ g) (evalV v i)
= { function composition }f (g (evalV v i))= { definition of evalT }f (evalT (Leafv g v) i)

evalT (comp f (Leafv g v)) i =The Leafc case:{ definition of comp }evalT (Leafc (f c)) i= { definition of evalT }f c= { definition of evalT }f (evalT (Leafc c) i)

evalT (comp f (Leafc c)) i =tuNow, we show the equation (2). In the following, l = length ls, r = length rs, andidces = [l..(n − r − 1)].
evalP (map f prog){ definition of evalP }map f (evalP prog)= { induction hypothesis }map f (eval [[ ls, zms, rs ]])= { definition of eval }map f (map evalT0 ls ++ map (evalT zms) idces ++ map evalT0 rs)= { definition of map and its distributivity: map h ◦ map g = map (h ◦ g) }map (f ◦ evalT0 ) ls ++ map (f ◦ evalT zms) idces ++ map (f ◦ evalT0 ) rs
= { Lemma 1 and definition of evalT0 }map (evalT0 ◦comp f ) ls ++ map (evalT (comp f zms)) idces ++ map (evalT0 ◦comp f ) rs
= { distributivity of map and definition of eval }eval [[ map (comp f ) ls, comp f zms, map (comp f ) rs ]]=Thus, the equation (2) holds.
A.3

Inductive Cases for shiftfi and shiftfl

What we have to prove is the following equations for an element e and a program prog.
Here, compile prog = [[ ls, zms, rs ]].
evalP (shift e prog) = eval [[ tail ls, slide (−1) zms, rs ++ [Leafc e] ]](3)evalP (shift e prog) = eval [[ [Leafc e] ++ ls, slide 1 zms, init rs ]](4)

To prove this equation, we first show a lemma.
Lemma 2. Let d be an integer, i be an index, and t be a computational tree. Then, thefollowing equation holds.
evalT (slide d t) i = evalT t (i − d)25

(5)Proof. This is shown by induction on Tree and Var .
The Node case:{ definition of slide }evalT (Node f (slide d l) (slide d r))) i= { definition of evalT }f (evalT (slide d l) i, evalT (slide d r) i)= { induction hypothesis }f (evalT l (i − d), evalT r (i − d))= { definition of evalT }evalT (Node f l r) (i − d)evalT (slide d (Node f l r)) i =The Leafv with Var case:{ definition of slide }evalT (Leafv f (Var x (s + d))) i= { definition of evalT and evalV }f (at (i − (s + d)) x)= { arithmetic }f (at ((i − d) − s) x)= { definition of evalT and evalV }evalT (Leafv f (Var x s)) (i − d)evalT (slide d (Leafv f (Var x s))) i =Since other cases of Leafv and the case of Leafc ignore the index i, the equation (5) holds inthese cases.
Thus, the equation (5) holds.
tuNow, we show the equation (3). In the following, l = length ls and r = length rs. First,
we assume ls is not empty, i.e. l > 0.
evalP (shift e prog)= { definition of evalP }shift e (evalP prog)= { induction hypothesis }shift e (eval [[ ls, zms, rs ]])= { definition of eval }shift e (map evalT0 ls ++ map (evalT zms) [l..(n − r − 1)] ++ map evalT0 rs)= { definition of shift, and ls being not empty }tail (map evalT0 ls) ++ map (evalT zms) [l..(n − r − 1)] ++ map evalT0 rs ++ [e]= { definition of map and evalT0 }map evalT0 (tail ls) ++ map (evalT zms) (map (+1) [(l − 1)..(n − (r + 1) − 1)])++map evalT0 (rs ++ [Leafc e])= { distributivity of map, and Lemma 2 }map evalT0 (tail ls) ++ map (evalT (slide (−1) zms)) [l0 ..(n − r 0 − 1)]++map evalT0 (rs ++ [Leafc e])
wherel0 = length tail ls, r 0 = length (rs ++ [Leafc e])= { definition of eval }eval [[ tail ls, slide (−1) zms, rs ++ [Leafc e] ]]26

Next, we assume ls = [ ].
evalP (shift e prog){ the same as the previous calculation }shift e (map evalT0 ls ++ map (evalT zms) [l..(n − r − 1)] ++ map evalT0 rs)= { definition of shift, and ls being empty }tail (map (evalT zms) [0..(n − r − 1)]) ++ map evalT0 rs ++ [e]= { definition of tail }map (evalT zms) ([1..(n − r − 1)]) ++ map evalT0 rs ++ [e]= { definition of map and evalT0 }map (evalT zms) (map (+1) [0..(n − (r + 1) − 1)]) ++ map evalT0 (rs ++ [Leafc e])= { distributivity of map, Lemma 2, and tail ls = [ ] }map evalT0 (tail ls) ++ map (evalT (slide (−1) zms)) [l0 ..(n − r 0 − 1)]++map evalT0 (rs ++ [Leafc e])00wherel = tail ls, r = length (rs ++ [Leafc e])= { definition of eval }eval [[ tail ls, slide (−1) zms, rs ++ [Leafc e] ]]=Thus, the equation (3) holds.
The equation (4) is shown similarly.
A.4

Inductive Case for zip

What we have to prove is the following equation for two programs prog 1 and prog 2 . Here,compile prog 1 = [[ ls 1 , zms 1 , rs 1 ]] and compile prog 2 = [[ ls 2 , zms 2 , rs 2 ]].
evalP (zip prog 1 prog 2 ) = eval [[ ls, zms, rev rs ]]where zms = Node id zms 1 zms 2
ls = trim FromL ls 1 ls 2 zms 1 zms 2
rs = trim FromR (rev rs 1 ) (rev rs 2 ) zms 1 zms 2(6)

To prove this equation, we first show a lemma.
Lemma 3. Let i be an index, and t be a computational tree. Then, the following equation
holds.
evalT0 (inst FromL i) = evalT t i(7)

Proof. This is shown by induction on Tree and Var .
The Node case:{ definition of inst }evalT0 (Node f (inst FromL l i) (inst FromL r i))= { definition of evalT0 }f (evalT0 (inst FromL l i), evalT0 (inst FromL r i))= { induction hypothesis }f (evalT l i, evalT r i)
= { definition of evalT }evalT (Node f l r) i

evalT0 (inst FromL (Node f l r) i) =27

The Leafv with Var case:{ definition of inst }evalT0 (Leafv f (Fix x (−s + i) FromL))= { definition of evalT0 }f (at (−s + i) x)= { definition of evalT }evalT (Var x s) ievalT0 (inst FromL (Leafv f (Var x s)) i) =Since other cases of Leafv and the case of Leafc merely return the given tree and this tree
ignores the index, the equation (7) holds in these cases.
Thus, the equation (7) holds.
tu
We also use the following lemma.
Lemma 4. Let i be an index, and t be a computational tree. Then, the following equation
holds.
evalT0 (inst FromR i) = evalT t (n − 1 − i)tuProof. Similar to the proof of Lemma 3.

In the following, l1 = length ls 1 , r1 = length rs 1 , l2 = length ls 2 and r2 = length rs 2 .
First, we assume l1 ≥ l2 and r1 ≥ r2 .
evalP (zip prog 1 prog 2 )= { definition of evalP }zip (evalP prog 1 , evalP prog 2 )= { induction hypothesis }zip (eval [[ ls 1 , zms 1 , rs 1 ]]) (eval [[ ls 2 , zms 2 , rs 2 ]])= { definition of eval and zip, letting ls 1 = ls 11 ++ ls 12 and rs 1 = rs 11 ++ rs 12 }zip (map evalT0 ls 11 ) (map evalT0 ls 2 )++ zip (map evalT0 ls 11 ) (map (evalT zms 2 ) [l2 ..(l1 − 1)])++ zip (map (evalT zms 1 ) [l1 ..(n − r1 − 1)]) (map (evalT zms 2 ) [l1 ..(n − r1 − 1)])++ zip (map evalT0 rs 11 ) (map (evalT zms 2 ) [(n − r1 )..(n − r2 − 1)])++ zip (map evalT0 rs 12 ) (map evalT0 rs 2 )= { Lemma 3 and Lemma 4 }zip (map evalT0 ls 11 ) (map evalT0 ls 2 )++ zip (map evalT0 ls 11 ) (map evalT0 (map (inst FromL zms 2 ) [l2 ..(l1 − 1)])++ zip (map (evalT zms 1 ) [l1 ..(n − r1 − 1)]) (map (evalT zms 2 ) [l1 ..(n − r1 − 1)])++ zip (map evalT0 rs 11 ) (map evalT0 (map (inst FromR zms 2 ) [(r1 − 1)..r2 ])++ zip (map evalT0 rs 12 ) (map evalT0 rs 2 )= { introducing (Node id), and definition of eval and evalT0 }map evalT0 (zipWith (Node id) ls 11 ls 2 )++ map evalT0 (zipWith (Node id) ls 11 (map (inst FromL zms 2 ) [l2 ..(l1 − 1)]))++ map (evalT (Node id zms 1 zms 2 )) [l1 ..(n − r1 − 1)]++ map evalT0 (zipWith (Node id) rs 11 (map (inst FromR zms 2 ) [(r1 − 1)..r2 ]))++ map evalT0 (zipWith (Node id) rs 12 rs 2 )28

{ combining edge elements }map evalT0 (zipWith (Node id) ls 1 (ls 2 ++ map (inst FromL zms 2 ) [l2 ..(l1 − 1)]))++ map (evalT (Node id zms 1 zms 2 )) [l1 ..(n − r1 − 1)]++ map evalT0 (zipWith (Node id) rs 1 (map (inst FromR zms 2 ) [(r1 − 1)..r2 ] ++ rs 2 ))= { definition of eval and trim }eval [[ ls, zms, rev rs ]]where zms = Node id zms 1 zms 2
ls = trim FromL ls 1 ls 2 zms 1 zms 2
rs = trim FromR (rev rs 1 ) (rev rs 2 ) zms 1 zms 2=

The other cases (l1 ≥ l2 ∧ r1 < r2 , l1 < l2 ∧ r1 ≥ r2 and l1 < l2 ∧ r1 < r2 ) are similarly
shown.
Thus, the equation (6) holds.

29