Kento Emoto, Kiminori Matsuzaki, Zhenjiang Hu, and Masato Takeichi Department of Mathematical Informatics, University of Tokyo {emoto,kmatsu,hu,takeichi}@ipl.t.u-tokyo.ac.jp Abstract. Computations on two-dimensional arrays such as matrices and images are one of the most fundamental and ubiquitous things in computational science and its vast application areas, but development of efficient parallel programs on two-dimensional arrays is known to be hard. To solve this problem, we have proposed a skeletal framework on two-dimensional arrays based on the theory of constructive algorithmics. It supports users, even with little knowledge about parallel machines, to develop systematically both correct and efficient parallel programs on two-dimensional arrays. In this paper, we apply our framework to the matrix-convolutions often used in image filters and difference methods. We show the efficacy of the framework by giving a general parallel program for the matrix-convolutions described with the skeletons, and a theorem that optimizes the general program into an application-specific one.
1
Introduction
Computations on two-dimensional arrays, such as matrix computations, image process- ing, and difference methods, are both fundamental and ubiquitous in scientific computa- tions and other application areas [7, 11, 15]. However, development of efficient parallel programs on two-dimensional arrays is known to be a hard task due to the necessity of considering data allocation, synchronization and communication between processors. Skeletal parallel programming is one promising solution to the situation [5, 16]. In this model, users build parallel programs by composing ready-made components (called skeletons ) implemented efficiently in parallel for various parallel architectures. Since low-level parallelism is concealed in the skeletons, users can obtain a comparatively ef- ficient parallel program without needing technical details of parallel computers or being conscious of parallelism explicitly. We have proposed a skeletal framework on two-dimensional arrays [9], based on the theory of constructive algorithmics (also known as Bird-Meertens Formalism ) [2,4]. Our framework provides users, even with little knowledge about parallel machines, with a concise way to describe safe and efficient parallel computations over two-dimensional arrays, and theorems for deriving and optimizing parallel programs. The main features of our framework are: (1) a novel use of the abide-tree representation [2] in develop- ing parallel programs for manipulating two-dimensional arrays; (2) a strong support for systematic development of both efficient and correct parallel programs in a highly abstract way; (3) an efficient implementation of basic skeletons in C++ and MPI on PC
K. Emoto et al.
clusters, guaranteeing that programs composed with these parallel skeletons can run ef- ficiently in parallel. To develop parallel programs in our framework, users construct a simple and general program that covers a class of problems, derive its efficient version using general techniques such as fusion, tupling and generalization, and then instantiate the general program to solve concrete problems. Usually, this derivation is summarized as a theorem (tool). In this paper, we give a domain-specific tool and show the efficacy of the frame- work. We focus on computations known as matrix-convolutions [12], in which each element in the resulting array depends on its surrounding elements. This set of compu- tations includes important and fundamental problems such as image filters, difference methods and the N -body problem (although this last problem seems more difficult than the others, it merely refers to not only the nearest neighbors but all the surrounding elements). The most general form mconv is described with three components: mconv f shrink = map f o map shrink o surrounds . Here, surrounds gathers all the surrounding elements for each element, shrink picks the necessary parts up from those gathered elements, and f calculates the resulting el- ement from them. This general form is parameterized by the two functions shrink and f , and users can solve many problems by specifying suitable ones. For example, users can develop a sharpen-filter by choosing the function shrink that reduces the surround- ings into a 3 × 3 matrix, and the function f that calculates the weighted sum of the nine values. We can further optimize instances of the general program to application-specific ones with the surrounding theorem . The main contributions of this paper are as follows. – We show the general parallel program for the matrix-convolutions described with parallel skeletons. Users can solve their problems as its instance. - We give the surrounding theorem that enables users to get an efficient program easily. The experimental results show that the derived program can be executed efficiently in parallel. Technical details of this paper are available in the master's thesis [8].
Notations
2
Notation in this paper follows that of Haskell [3], a pure functional language that can describe both algorithms and algorithmic transformation concisely. Function application is denoted by a space and the argument may be written without brackets. Thus, f a means f ( a ) in ordinary notation. Functions are curried, i.e. func- tions take one argument and return a function or a value, and the function application associates to the left. Thus, f a b means ( f a ) b . The function application binds more strongly than any other operator, so f a ⊗ b means ( f a ) ⊗ b , but not f ( a ⊗ b ). Function composition is denoted by ০, so ( f o g ) x = f ( g x ) from its definition. Binary operators can be used as functions by sectioning as follows: a ⊕ b = (a⊕) b = (⊕ b ) a = (⊕) a b . Two binary operators ≪ and ≫ are defined by a ≪ b = a, a ≫ b = b. Pairs are Cartesian products of plural data, written like ( x, y ). A function that applies functions f and g respectively to the elements of a pair ( x, y ) is denoted by ( f × g ). Thus, (f × g) (x,y) = (f x,g y).
2
Skeletal Framework on Two-Dimensional Arrays
3
In this section, we introduce our parallel skeletal framework on two-dimensional ar- rays [9] based on the theory of constructive algorithmics [2,4].
Abide-trees for Two-Dimensional Arrays
3.1
To represent two-dimensional arrays, we define the abide-trees, which are built up by three constructors | · | (singleton), ← (above) and $ (beside) following the idea in [2]. data  AbideTree α = | | α (AbideTree α) ⊕ (AbideTree α) (AbideTree α) φ (AbideTree α) Here, | · | a , or abbreviated as | a | , means a singleton array of a , i.e. a two-dimensional array of a single element a . For two-dimensional arrays x and y of the same width, x �y means that x is located above y . Similarly, for arrays x and y of the same height, x φ y means that x is located on the left of y . Moreover, ← and Φ are associative operators and satisfy the following abide (a coined term from above and beside) property. Definition 1  (Abide Property). Two binary operators ⊕ and ➤ are said to satisfy the abide property or to be abiding, if the following equation is satisfied: (x & u) ⊕ (y & v) = (x ⊕ y) ⊕ (u ⊕ v) . In the rest of the paper, we will assume that x has the same width as y when x ⊕ y appears, and that u has the same height as v for u φ v . Note that one two-dimensional array may be represented by many abide-trees, but these abide-trees are equivalent because of the abide property of Φ and Φ . For example, we can express the following 2×2 two-dimensional array by two equivalent abide-trees.
This is in sharp contrast to the quadtree representation of matrices [10], which does not allow such freedom. From the theory of constructive algorithmics [4], it follows that each constructively built-up data structure (i.e., algebraic data structure) is equipped with a powerful com- putation pattern called homomorphism. Definition 2  ((Abide-tree) Homomorphism). A function h is said to be an abide-tree homomorphism, if it is defined as follows for a function f  and binary operators ⊕, ⑧.
For notational convenience, we write (| f, ⊕, ⊗|) to denote h. When it is clear from the context, we just call ( f , ⊕ , ⊗ l ) homomorphism. Note that ⊕ and ⊗ in ( f , ⊕ , ⊗ l ) should be associative and satisfy the abide property, inheriting the properties of → and $. Intuitively, a homomorphism ( f, ⊕, ⊗ ) is a function to replace the constructors | · |, ⊕ and φ in an input abide-tree by f , ⊕ and ⊗ respectively.

Parallel Skeletons on Two-Dimensional Arrays
3.2
We introduce the parallel skeletons map, reduce, zipwith, scan and scanr for manip- ulating two-dimensional arrays. In the theory of constructive algorithmics [2,4], these functions are known to be the most fundamental computation components for manip- ulating algebraic data structures and for being glued together to express complicated computations. Intuitive definitions of the skeletons are shown in Fig. 1. All the skele- tons are implemented efficiently in parallel and their costs are shown in Table 1. The skeletons map and reduce are two special cases of homomorphism. The skele- ton map applies a function f to each element of a two-dimensional array while keeping the shape of the structure. The skeleton reduce collapses a two-dimensional array to a value using two abiding binary operators ⊕ and ⊗ . They are defined formally as map f = (|| o f, e, e (), and reduce( ⊕, 8) = (| id, ⊕, 8). The skeleton zipwith, an extension of map, takes two arrays of the same shape, applies a function f to corresponding elements of the arrays and returns a new array of the same shape. The skeletons scan and scanr, extensions of reduce, hold all values generated in reducing an array by reduce. The scan generates the result of reducing
4
Surrounding Theorem
imt sharpen_filter(imt **b, imt **a, imt m, imt for(imt i = 0; i < m; i++) b[i][j] = f(a[i][j], a[i-1][j], a[i+1][j], a[i-1][j+1] , a[i-1][j-1] , a[i+1][j+1] , a[i+1][j-1]); imt f(imt c, imt m, int s, int e. int v, int ne, int av, (-1) $ 5 return 5�c (-1) $n ◆ (-1) $c ✦ 0 $ 5 V : 1 [−1] $ v Fig. 2. C++ Code of the Sharpen Filter (Sequential Program) the upper-left subarray, while the scanr generates that of the lower-right subarray. We omit the formal definition of zipwith, scan and scanr for the space limitation.
4
Developing Parallel Programs for Matrix-Convolutions
In this section, focusing on the matrix-convolutions such as image filters and difference methods, we give the general form described with parallel skeletons, and then give the theorem to get optimized program from the general form. The matrix-convolution is computation in which each element of the resulting array depends on the surrounding elements. For example, the sharpen-filter that sharpens the input image is one instance of the matrix-convolution. A pixel of the resulting image is the weighted sum of the surrounding pixels of the input image. Similarly, the difference method is another instance of matrix-convolution since it calculates the new value of each point from the old values of the surrounding points. We show a code in C++ for the sharpen-filter in Fig. 2, to give a concrete image of the problems dealt with here.
A General Form Described with Parallel Skeletons
4.1
As argued in the introduction, the most general form of this kind of computation is thought to consist of three components: gathering all the surrounding elements of each element to it, shrinking those to the necessary amount, and applying a function to get a S new element from them. Thus, the program is described as follows: mconv f shrink = map f o map shrink o surrounds . The idea of our general from is illustrated in Fig. 3 that shows an image of execution of the sharpen-filter: (1) surrounds gathers all the surrounding elements for each element, (2) shrink picks the necessary parts up from those gathered elements, and (3) f cal- culates the resulting element from them. This general form has clear correspondences to the code in Fig. 2. The function f corresponds to f of the code, shrink corresponds to which elements are the arguments passed to f , and surrounds corresponds to for- loops. Thus, users can easily write their programs using the general form. This general form is parameterized by the two functions shrink and f , and users can solve many problems by specifying application-specific ones, as shown below. The function surrounds , which is commonly used in those problems, has two-phase calcu- lation as follows: (1) calculation of the parts of the northwest (i.e. c, n, w and nw ) by scan, and (2) that of the other parts by scanr. Its definition is as follows. surrounds = scanr(⊕, ⊕,) o map f, o scan(⊕, ⊕, ⊕, o map ff where ff a = (a, Nil, Nil, Nil) , nw, �w, �nw, (ca, na, wa, nwa) ⊕ f (cь, nь, wь, nwь) = ( cь , na ⊕ ca | ⊕ nь, (ca, na, wa, nwa) ➤ f (cь, nь, wь, nwь) = ( waφ Ca φ wb, nwaφ naφ nwb fr (c,n,w,nw) = (c,n,Nil,Nil,w,Nil,nw,Nil,Nil) (ca, na, sa, ca, wa, nea, nwa, sea, swa) ⊕ r (cb, nb, sb, eb, wb, neb, nwb, seb, swb) na , saやC LやS L wa , nea , nwa , seaやebやseb, swaやwhやswe (ca, na, sa, ca, wa, nea, nwa, sea, swa) Sr (cb, nb, sb, eb, wb, neb, nwb, seb, swb) na , sa , Ea Φ | Cb | Φ Eb, Wa , nea φ nb φ neb, nwa, sea φ sb φ seb, swa Here, Nil is a special value to indicate that there is no value, and we treat it as an identity of ➪ and Φ for simplification of the notation. Thus, Nilつx = x, x ∞ Nil = x, Nil φ x  =  x , and x φ Nil  =  x . Each element of the resulting array is a tuple of nine elements. The meaning of each element of the tuple is as follows: c is the center element; s is an array of the elements on the south of the element; similarly n , e and w are arrays of the elements on the north, east and west respectively; ne , nw , se and sw are arrays of the elements on the northeast, northwest, southeast and southwest. Note that this surrounds needs O( n 4 ) memory space for a matrix of n × n . We show some examples written with the general form. imagefilter ker = mconv (conv ker) shrink1 = iter n (mconv (conv ker) shrink 1 ) FDM n ker where shrink1 = id × B × T × L × R × BL × BR × TL × TR B =(||,),¢|), T =(||,(,¢|), L =(||,¢,(|), R =(||,¢,( BL=( | |, ≫, ≪ |, BR=( | |, ≫, ≫ |), TL=( | |, ≪, ≪ |), TR=( | |, ≪, The function imagefilter ker is an image filter with the coefficient matrix ker , which is used to compute weighted sum of the surrounding pixels. The shrink 1 reduces each part of the gathered surrounding elements to the element closest to the center, and the function conv ker calculates the weighted sum of them. The functions B and T take the bottom row and the top row of the input array respectively. Similarly, each of L, R, BL, BR , TL and TR takes corresponding part of the input array. Figure 3 shows an image of execution of the sharpen-filter by the above general program. The function FDM n ker performs the finite difference method, where iter is an iteration function and each iteration step is the same as image filters with specific coefficients. The following example calculates the array of which element at ( i, j ) is the maxi- mum in the i -th row and the j -th column, i.e. the maximum in the cross. The shrink max reduces each part of the gathered surrounding elements to the biggest element in the part, where the binary operator ↑ takes the bigger element. The function max 5 takes the maximum of the column and the row including the center element. crossmax = mconv max5 shrinkmax where shrink max = max × · · × max max = ( id , ↑ , ↑ ) max5 (c, n, s, e, w, -, -, -, -) = c ↑ n ↑ s ↑ e ↑ w As shown in this example, shrink is allowed not only to shrink the shape of the sur- roundings but to perform some calculation.
4.2  Surrounding Theorem
In this section, we give the theorem to optimize the general form by fusing shrink to surrounds. Image filters and difference methods usually have the shrink of the fixed size win- dow that takes the fixed-size rectangle region (window) of the surrounding elements. The function that takes a fixed number of columns (rows) can be written as a homo- morphism. For example, the function right = ( | |, ∞, ≫ |) takes the right-most column, which is used in the examples in the previous section. Thus, we here consider the gen- eral shrink that consists of homomorphisms. It is defined as follows. shrink = 9c × h n × h s × h e × h w × h n e × h n w × h se × h sw where h n = (19n, Bn, Bn) ,  h s = (19n, Bs, Bs) ,     h c = (19c, Bc, Bc)
h w = (19nw, Bnw, Bw) ,  h n c = (19ne, Bne, Bne) , h nw = (19nw, Bnw, Bnw) h se = (19sc, Bsc, Bsc) , h sw = (19sw, Bsw, Bsw) Here, ⊕ χ and ⊗ χ are extended to satisfy the following equations: Nil ⊕ χ x = x, x ⊕ x Nil = x , Nil ⊕ x x = x , and x ⊕ x Nil = x . The general form using this shrink uses O( n 4 ) operations for a two-dimensional array of n × n. Then, we give the result of the optimization by fusing shrink to surrounds. Theorem 1  (Surrounding). Let the function shrink be defined by homomorphisms as above. Then, there exist a projection function proj and operators ⊕ f, ஸ f, ⊕ f, and ஸ f, whose complexity is bounded by the largest of ⊕ χ and ⊕ χ , and the program mconv f shrink is optimized to the following program. map (f ◦ proj) ◦ scanr(⊕ ρ , ⑧ ρ ) ◦ map f , ′ ◦ scan(⊕ ρ , ⑧ ρ ) ◦ map f f Proof. The theorem is proved by the promotion of map shrink with extending the tuples. See the master's thesis [8] for details . The resulting program uses O( n 2 ) operations for a two-dimensional array of n × n , while the original general form uses O( n 4 ) operations. The parallel complexity of the resulting program is O((n 2 /P + √ π 2 /P log P)T( @x,@x)) for P processors, provided that the calculational complexity of ⊕ χ and ⊕ χ in the homomorphisms are T ( ⊕ χ . ⊕ χ ) All the examples shown in the previous section have the shrink functions described with homomorphisms. Thus, we can apply this theorem to all of them, and they are executed in O( n ² / P + √ n ² / P log P ) complexity using the skeletons. As mentioned above, the function that takes a fixed number of columns (rows) can be written as a homomorphism. Thus, this theorem holds for the shrink of the fixed size window that shrinks the surrounding elements to a fixed size, which is often seen in image filters and difference methods. Corollary 1  (Fixed Size Window). Let the function shrink be the fixed size window. Then, the program  mconv  f  shrink  is optimized to that of O(n 2 ) operations. Note that the homomorphism taking h × w subarray of a two-dimensional array has the operators of O( wh ) complexity. Thus, the total complexity of the program of fixed size window is O(n²wh). Finally, we note that we may perform more optimizations by using the shifting of the edges instead of butterfly computations for the global computations of scan and scanr, provided that the operators influence only a fixed number of elements [8]. This leads to the parallel complexity of O(( n 2 / P + √ n 2 / P ) T( ⊕ x, ⊕ x ) ) for P processors.
5
Experimental Results
We implemented the program¹ using our parallel skeleton library [14] and did our ex- periment on a cluster (distributed memory). Each of the nodes connected with Gigabit Ethernet has a CPU of Intel® Xeon® 2.80GHz and 2GB memory, with Linux 2.4.21 for the OS, gcc 2.96 for the compiler, and mpich 1.2.7 for the MPI. Figures 4 and 5 show the speedups and the calculation times of the sharpen-filter. The program is an optimized one from the general form (an equivalent of the program in Fig. 2). The inputs are images of 1000 × 1000 and 2000 × 2000. The computation times of the program on one processor are 0.70s and 3.85s respectively. The result shows programs described with skeletons can be executed efficiently in parallel, and proves the success of our framework. The program achieves almost linear speedups , and the total computational complexity of the optimized program is O( n 2 ) (thus, its parallel complexity is O( n ² / P ) for small P ). However, the serial performance is rather poor due to the overhead of using general skeletons (i.e. scan and scanr). We think this problem can be solved by replacing the general skeletons with those special- ized for this domain, and it can be automatically done by compilers (future work).
Related Work
6
SKiPPER [17] is a skeleton-based parallel programming environment for real-time im- age processing. It has skeletons specialized for image processing, while we use general skeletons on two-dimensional arrays. Thus, a program developed with SKiPPER may be faster than that written with our skeletons, but, the latter program can be easily composed with other programs and be optimized by fusion due to generality and solid foundation of our skeletons. There are several other skeletal parallel approaches (libraries), such as eSkel [1], Muesli [13] and P3L [6]. Their formalizations of skeletons on two-dimensional arrays are not enough (e.g. they have no scan skeletons, and the reduction takes only one op- erator) to deal with matrix-convolutions suitably. Our skeletons have a solid foundation, so that we can easily deal with matrix-convolutions and perform optimizations.
7
Conclusion
In this paper, we proposed a general theorem, called surrounding theorem , for optimiza- tion of a general skeleton program into an efficient application-specific program. It can deal with a wide class of matrix-convolution problems including image filters and dif- ference methods. The experimental results show that the optimized program can be ex- ecuted efficiently in parallel. We are now working on making an automatic mechanism for translating the sequential code to our general form with skeletons, and further an optimization mechanism for the application-specific program with respect to its global communication and sequential performance. 9
Acknowledgment
8
This work is partially supported by the Grant-in-Aid for Scientific Research (B), No. 17300005, Japan Society for the Promotion of Science. We are grateful to the ref- erees for their detailed and helpful comments.