NASA Contractor Report 178024

lCASE REPORT NO. 85-55

leASE

A PARTITIONING STRATEGY FOR NON-UNIFORM

PROBLEMS ON MULTIPROCESSORS

Marsha J. Berger

Shahid Bokhari

No0]{-Cr-rif?;CD C.)

NASA-CR-178024

19860010476

I .

Contract No.

November 1985

NAS1-17070 t/\[. --, I,_~: '\I,~,.-i c:-·q~;~

L::,"',: ;Y, : :.',~ \

r/;1r..':i',UlJ. 'v'1'{(~If!L'1.

INSTITUTE FOR CO~WUTER APPLICATIONS IN SCIENCE AND ENGINEERING

NASA Langley Research Center, Hampton, Virginia 23665

Operated by the Universities Space Research Association

NJ\SI\

National Aeronautics and

Space Administration

Langley Research cent.

Hampton. Virginia 23665

1

., ",.,

1 RN/NASA-CR-178024

_0._- _

~\lR:jA-t~K -1 /t:(J~q

ICASE-85-55 WAS i~2G:i78024· CNT*: NASl-17070 DE-AC02-76ER-030?7-V

85/11/00 30 PAGES UNCLASSIFIED DOCUMENT Ui~L: A partitioning strategy for nQnUn~f0rm problems on multiprocessors iL~~:

National Aer0naut~cs and Space

/*OATA PROCESSING EQVIPMENT/*MULTIPROCESSING (COMPUTERS)/*PARTITIONS

GKIDj/ MAPPING/ PARTIAL DIFFERENTIAL EQUATIONS/ ----- iRt.t.:j

A Partitioning Strategy for Non-Uniform Problems on

Multlproeessors

MarshaJ. Berger*

CourantInstitute of MathematicalSciences

New York University

251 Mercer St.

New York, NY 10012

Shahid H. Boldmrl**

Institute for ComputerApplicationsin Scienceand Engineering

and

University of Engineering and Technology

Lahore-31, Pakistan

ABSTRACT

We considerthe partitioningof a problemon a domain with unequal

workestimatesin differentsubdomainsin a waythatbalancesthe workload

acrossmultipleprocessors. Sucha problemarisesfor examplein solvingpar.

tial differentialequationsusing an adaptivemethod that places extragrid

pointsin certainsubregionsof the domain. We use a binarydecomposition

of the domainto partitionit into rectanglesrequiringequal computational

effort. We then studythe communicationcostsof mappingthispartitioning

ontodifferentmultiprocessors:a mesh-connectedarray,a treemachineand a

hypercube.Thecommunicationcostexpressionscanbe usedto determinethe

optimaldepthof the abovepartitioning.

*Supported in part by Department of Energy Contract No. DEAC0276ER03077.V, and by

th.e.Na.fionalAeronautics and Space Administration under NASA Contract No. NAS1-17070.

wml.e me. auraor was m reslaeenceat the Insti_te for Computer Applications in Science and

tmgmeenng, JNP_sALangley ttesearcn t:enter, rmmpton,VA 23665.

*Supported by NASA Contract No. NASI-17070 while the author was in residence at the

stitute for computer Applications in Science and Engineering (ICASE), NASA Langley Research Center.

L) q/?
-2-

I. Introduction

We consider the partitioning of a problem on a domain with unequal computational

work estimates in different subdomains, in a way that balances the work load across multiple

processors. Such a problem arises, for example, in solving hyperbolic partial differential

equations using an adaptive method that places extra grid points in certain subregions of the

domain (see e.g. [5] and [6]).Suchan approachhas also been proposed in the multigrid literature ([2],[3],[12],[15]). At a given instant of time a typical computational mesh for either of

these problems might look like Fig. 1.1.

F_. l.l Increased resolutions is obtained by superimposing fine grid patches over an

underlyingglobal coarse grid.

Hence, any simple partitioning scheme must accountfor the unequal amount of work to be

done in the left half versus the right half of the domain. H Fig. 1.1 representsthe grid for a

time dependent problem, the work load in the left and right halves is significantly different,

since the finest grids in space typically need a smaller step in time (if an explicit finite different scheme is used to integrate the solution), and so many time steps are taken on the fine

grid for every one on any coarser grid.

Other numerical examples which give rise to a problem with unequal work estimates

might come from solving a pde with different equation sets in different regions. For example, in calculatingtransonic flow around airfoils, the Navier Stokes equations may be used in

a boundary layer around the airfoil, and the Euler equations or even the potential equations

can be used in the farfield. These different sets of equations have very different costs associated with their corresponding difference schemes. Another possibility is that the work esti- "

mates come from different physics in different parts of the domain, for example in weather
-3-

calculations, depending on whether a region is over water or over land. Other examples of

unequal work include an iteration scheme such as SOR where a certain percentage of the

domain is relaxed a second time before iterating on the entire domain again. This ad hoc

procedure, applied to say the 10% of the grid with the largest residualcan greatly reduce the

cost of convergence.

In addition to the above static domains, we are interested in investigating the possibilities of solving adaptive mesh refinement problems on a multiprocessor system. A majordifficulty is that, no matter how portions of the mesh are initially assigned to processors, a

change in the mesh refinement will ultimately cause the computational load on the processors

to become unbalanced. Attempts at rebaiancing are complicated by the need to keep the

interprocessor communication overhead at a minimum. Since the adaptive mesh refinement

strategyin [6] is already based on a partitioning of the domain into rectangular gridpatches,

we can derive an approachpresented here which is simple and tractable. Other approaches

are given in [9], [17].

Most partitioning strategiesuse some type of domain decomposition to balance the work

load over many processors. Typically, these uniform mesh problems can be divided into

boxes (Fig. 1.2a) or strips (Fig. 1.2b).

(a) (b)

Fig. 1.2 Two common partitioning strategiesfor rectangular mesh problems.

The benefits of one over the other depend on the cost of transferring information around the

perimeter of a box to the neighboring partition/processor (which depends on the machine

architecture), and the order in which computations on such configurations can proceed.

Papadimitriou and Ullman [13] discuss communication/time tradeoffs for such partitionings,

_d obtainlower bounds on the rate such computations can proceed.
-4-

A different kind of partitioning is evident in the work of Adams and Jordan [1]. By

partitioning down to the grid point level, using a multi-color SOR iteration scheme, simultaneous updates can proceed for any given color grid point throughoutthe entire mesh. This

can be useful on processor arrays as well as vector computers.

This paper is organized as follows. We describe the binary decomposition used to partition the work load in section 2, and discuss some of its propertiesin section 3. In sections 4,

5 and 6, we study the communicationcosts of mapping this partitioning onto different types

of multiprocessors: a nearest neighbor array, a tree machine, and a hypercube. We derive

expressiom for the communicationversus computationcosts which can be used to determine

an optimal depth for the above partitioning. Section 7 summarizesthe results.

2. Binary Decomposition of the Domain

In this section we describe the strategy used to partition a domain into subunits requiring equal computational effort. In this presentationwe will assume that the number of available processorsis a power of 2, although many of our results generalize. Another underlying

assumption is that the numberof gridpoints N>>p, the number of processors. We will concentrate on the static case, and only say a few words about adaptivelyrebalancing the decomposition later in the section.

Suppose that work estimates on a given domain have already been obtained, through a

priori knowledge, or from an initial computationon a uniformmesh using a partitioning as in

Fig. 1.2. Given these work estimates, we can now make a vertical cut throughthe domain so

that the left and right segments each contain half the work (or as near as possible given the

constraint that the line is vertical, and the numberof gridpoints in each segment increases by

a finite amount on shifting the location of the cut by one column). If there are four processors available, the two segments are each partitioned using two different line segments of a

horizontal cut line next, into a total of four equally balancedwork loads. This procedure continues by re.cursivelypartitioning using first vertical then horizontal cut line segments, so that

the length of the longest side of any subregion is reduced every other step. A typical decomposition for the grids in Fig. 1.1 using 16 processorsis shown schematicallyin Fig. 2.1. The

idea for this decomposition was inspired by the similar looking rectangular regions used by

Bentley [4] in answeringtwo dimensional point domination questions.

We emphasize that the computationalwork of an iteration on any rectangular region in

Fig. 2.1 is identical. However, the communication requirements across the perimeters are

not. In particular, if the source of the problem is a grid such as Fig. 1.1, then the grid point

density along any given line segment varies, depending on whether a cut line intersects a

finer grid or not. In general, therefore, we will only obtain upper bounds for communication
-5-

2 8 10 0

9 11

3

1 14

12

4 6

-- 13 15

5 7

Fig. 2.1 A binary decomposition using 16 processors.

costs assuming the worst case gridpoint distribution. Instead, therefore, we will now assume

that the problem gives rise to different work requirementsin different regions, but is based

on an underlyinggrid of N2points which are uniformly distributed.

There are several points to note about the decomposition depicted in Fig. 2.1. First, by

restrictingthe subunits to be rectangular blocks, we avoid a messy problem with data structures. If more general L-shaped regions or diagonal lines were used, the specification of a

region would be more difficult. All that this approachrequiresto spedfy each block is the

four corners of the rectangle (2 will do). This is sufficiently low overheadthat every processor can keep a map of the entire domain with each processor/rectanglepair. A tree data

structure can easily be traversed for any neighbor information that is needed, for example,

for a non-sharedmemory machine.

Secondly, this approach does not suffer the drawbacksthat other decompositions suggested for this problem have. For example, Fig. 2.2 indicates a grid configuration with four

fine grids superimposed on a global coarse grid, and2 further refined grids nested in 2 of the

four. It is tempting to use these grids, which already form one type of decomposition of the

domain (and are each regular with a simple data structure) as the basis for assignment of

work to a processor. However, there is no attempt at load balancing in this approach. In

addition, if the mesh is later changed so that there are 8 subgrids, for example, either the

computation must request 8 processors, or some of them were idle beforehand.

We mention that the partitioning itself is easily accomplished. In principle, the partitioning can be done by summing the numbers of grid points (or work per gridpoint times the

number of grid points) first rowwise, then colurunwise. The partitioning cut is then made in
-6-

Flg. 2.2 A domain decomposition with 4 subgrids, and2 sub-subgrids.

the middle.

We point out one final advantage of the binary decomposition. As the computation

proceeds, if it turns out that one region gets more work (say a finer mesh is introduced), a

local rebalancing can be done without necessarily redoing the entire partitioning. For example, in Fig. 2.1 if region 14's work load increases by some amount B, the last cut line in the

partition, which separates regions 14 and 15, can be adjusted so that both regions have an

imbalance of only 6/2. If the k previous cut lines are adjusted, the imbalance is reduced to

B/'2k. The cost of rebalancing (which might include time to send data to the new governing

partition, for example), can be traded off against the lost time of having an imbalance in the

work load. This can determinehow high in the tree (the number k above) to rebalance.

3. Analysis of Partlflonlngs

In this section we present several definitions related to our partitionlngs and analyze

some of their important properties. These are essential to our discussion of mappings of

subregions onto various processor architectures, which we present in subsequentsections.

3.1. Deflnltlons

The depth of a partitioning is the number of times the domain has been partitioned.

This equals the depth of the correspondingbinarytree, with root node correspondingto the

entire domain, and leaf nodes corresponding to each rectangle in the final partitioning.
-7-

Each partition line is divided into a number of segments by the incidence of other partition lines.

The total number of segments of a partitioningis the sum of all such segments.

The depth of a segment is the depth of the partitionline to which it belongs.

For example,Fig. 3.1 showsa partitioningof depth2. The depth1 partitionline a-bis

dividedinto 3 segmentsby the inddenceof the two depth2 partitionlines c-d ande-f. The

totalnumberof segmentsin this partitioningare 5.

e f

c d

i

b

Fig. 3.1 A depth 2 partitioning with 5 segments.

We observe that if two adjacentregions of a partitioning are assigned to different processors then the segment between them represents a communication requirementbetween the

processors. The following definition makesthis easier to appreciate.

The graph of a partitioning is the dual graphobtained in the usual way [8] by representing each region by a node and connectingtwo nodes if and only if the correspondingregions

are adjacent(share a segment on their perimeter). Each edge in the graph of a partitioning

represents a communicationrequirement between the two regions representedby the nodes at

its end points. Fig. 3.2 shows a partitioning of depth 4 along with its graph. The assignment

of regions of a partitioning to the processors of a multiple computer system is now equivalent

to the mapping of the graph of a partitioning onto the graph of a multiple computer system

[71.
-8-

v

Fig. 3.2 A depth4 partitioningandits dualgraph.

3.2. Properties

Clearly, the total number of segments of a partitioning is an important number in our

analysis, since it represents the total number of different communication paths required. A

lower bound on this number occurs when the domain has uniform computational density in

which case the partitioning is made up of continuous horizontal and vertical lines extending

from one side of the region to the other, It is easy to show that in this case the total number

ofsegments TL is

L+i

Tl.(k=)-2k+i --22 , k even (3.1)

k-I

Tt.(k)=k+i-3 2 .2 2 , kodd,

fora depthk partitioning for2kprocessors.

To obtainan upperboundTt.,(onk)thenumberofsegmentswe, needtoinvestigafteurther

propertieofspartitionings.

3.2.1. Property1

Notethatthe graphof a partitioningremainsunchangedno matterhowmuchtheverti.

cal(horizontalpar ) titioninglinesare displacedaslongasthe sequenceof the vertical(horizontal)coordinatesof these lines are not disturbed.Becauseeachedge of this graph

correspondsto a segment,this propertyimpliesthat the total numberof segmentsis

unchangedundersuchdisplacements.
-9-

3.2.2. Property 2

A partitioning of depth k may be considered to be made up of the juxtaposition of two

partitioningsof depthk-1. Theline atwhichtheconstituentdepthk-1 partitioningsabut

eachotherbecomesthe newdepth! line,andthedepthof the remaininglinesincreasesby

& &-1

one. If k is even, the constituent depth k-1 partitionings will have 2 2 by 2 2 sides;in this

case the sides with the larger number of regions will face each other. If k is odd, the consfi.

k-1

tuent regions have 2 2 regions along either side.

Because of Property 1 above, we can displace the segments perpendicularto the interface of the constituent regions to any extent as long as we do not disturb their sequence.

This allows us to distort the constituent partitions so that no two depth k-1 segments on

either side of the interface are collinear, and thus obtain the maximum number of depth 1

segments, which turns out to be

Smax(k-) 22 -I, k even (3.2)

k+1

Smax(k-) 2"/--I, k odd.

Itispossiblfore a partitioninofgdepthk tohavethemaximum possiblenumberof

depthIsegmentsi,ndependenoftitsconstituedntepthk- !partitionings.

3.2.3. Property 3

It thus follows that the depth k partitioning with maximum total segments is made up of

two depth k-1 partitionings with maximum total segments. This leads to the following

recurrencefor maximum total number of segments,

Tu(k)--'2 z -I + 2 Tu(k-1), k even (3.3)

k+l

T_,(k)=22 -I+ 2Tu(k-1)k,odd

Thesolutiotnostheserecurrencearse

Tu(k)= 2k+2- 2k- 2_+2+ 1, k even (3.4)

Tu(k )=2 k+2-2 k-3.2 2 +1, kodd
- 10 -

3.2.4, Property 4

Another property of great interestis the maximum degree of any node in the graph of a

partitioning. A simple counting argument shows that for 2k processors, k>4, the maximum

degree is

22 +2 2 +3, keven (3.5)

k.l

2 2 +3, kodd.

The node with this maximum degree borderson the depth I partitioningline (see Fig. 3.3).

i

/ \

Flg. 3.3 Maximum degree possible for a 2kprocessor decomposition, k=4.

Infact,a simplecountinpgrocedursehowsthefollowingF.or2&processorks,even,

thepartition consisoftsk/2segmentsofverticalcutlineisnterlaced withk/2horizontal cut

line segments. If we number the cut lines from j = 1 to k in total, then the number of

incident edges a segmentfrom a cut line of depth j canhave is at most

e,U)==22 I- -- 1, k even (3.6)

ek(j)= 2 2 -- 1, kodd.

Usingthefactthat2/-Ilinesegmentsmakeup thejlhcutlinethe , expressioinns(3.6)canbe

summed toprovidaenalternadteerivatioofn(3.4for ) themaximum numberofedgesinthe
- 11-

dual graph.

4. MappingontoNearestNeighborArrays.

Inthissectionwe investigatehowthebinarypartitioningscanbemappedontonearest

neighborarrays. By nearestneighborarrayswemeanmulticomputersystemsin whichthe

processorsmaybe thoughtof aspointsonanintegerlattice,whereeachprocessorhascommunicationlinksto its 4, 8 ormorenearestneighbors.Forexamplethe , Illiacmachine(if

we ignoreits wraparoundconnections)is an exampleof a 4 nearestneighborarray(4nn

array). TheFEM[11]isanexampleof an8rmarray.

In this analysiswe willconsiderrectangulararraysof size2k. For sucharrays,there

existswhatwe callnaturalmappingsof depth/cpartitionedregionsontoprocessorsd,efined

asfollows.Whenpartitioningthe domainintotworegions,partitiontheprocessorarrayinto

twoequalhalves. Assignthe leftsubdomainto the lefthalfandthe rightsubdomainto the

righthalf. Repeatrecursivelyuntilthe processorpartitionhave s exactlyoneprocessorin

them.Atthispointeverypartitionedregionhasbeenassignedto aprocessor.

4.1. Cardinality of Natural Mappings

One way to measure the quality of a mapping is to compute its cardinaUty[7], defined

as the number of edges of the problem graph that fall on edges of the processor graph

divided by the total number of edges in the problem graph. Mappingswith a cardinality of

one have minimum interprocessorcommunicationoverhead since all processes that need to

communicate lie on processorsthat are adjacentto each other. One such extreme case occurs

when the domain being partitioned is uniform (as described in Section 3.2), and the graph of

the partitioning matches the graphof a 4nn arrayperfectly.

At the other extreme, decompositions of the type illustrated in Fig. 3.3 have graphsin

which some nodes have exponential degree. We could not possibly acommodate the edges

incident on such nodes using any fixed degree nearest neighbor array. The question then is

how low the cardinality can be over all possible decompositions.

A simple but important observation towards this goal is the following. In a natural

mapping of the graph of a partitioning onto a 4nn array, the edges on the perimeter of any

subdomaln are all mapped on the correspondingedges of the mesh. This can be seen in Fig.

2.1, where the 12 edges on the perimeter of the complete domain fall on the perimeter edges
- 12-

of the array. At the same time, the 8 perimeter edges in the left and right hand subdomains

also fall on perimeteredgesof the twohalves of the arrayand so on.

It follows that the edges of the partitioning graph that fail to fall on edges of the array

graphcannotexceed the number of edgesthat extend across partitionsand are notperimeter

edges. The number of such 'misses' that extend acrossthe depth1 partitionis preciselythe

number of depthI segments(eq. (3.2))less 2 (theperimeteredges),givingthe following

recurrences,

M(k)=-2_+_-3+2M(k-1), k_ (4.1)

k+!

M(k)= 2"T-3+2M(k-D,ko_, k>1.

Itisimportanttoappreciate thatfork= I,M(k)-O.

Restricting to the case of k even, these recurrencescan be solved to yield

M(k) = 3.2k-1 - 2_'+2 + 3. (4.2)

Combining (4.2) with the expressionfor the total number of edges yields the cardinality

C(k)=-T_,(k)-M(k)

Tu(k)

3.2k-1-2

= k_+2 (4.3)

3"2t - 2 2 +1

with a similar expression for k odd. Clearly, when k = 1, C(k) = 1, and it drops gradually

as k increases. It can be seen, however, that as k becomes large the cardinalityconverges to

0.5 and does not drop further. For example, for k=4, C(k)=.67, but for k=10, C(k)=.52.

A similar analysis for the 8nn array gives the recurrence relation

M(k) = 2_+_ - 5 + 2M(k- 1), which revealsthat C(k) convergesto .79.

One question to consider is whether these natural mappings are near optimal, or

whether there exists other mappings of regions to processors with higher cardinality. A

worst-case partitioning can have TUedges, where T_,is given by eq. (3.4). A 4nn array has

only TL edges, where TL is given by eq. (3.1). An optimal mapping is one which uses all 4nn

Tv-T,

edges, and thus has cardinality of at most --- ffi 75%. In the worst case, there are Tv
- 13 -

natural mappings that use all 4nn edges, and have the maximumnumber of additional edges

for a total of Tu, and so by this measure, natural mappings are within 2/3 of optimal. We

conjecture that in fact, no other mapping can do better for a 4nn array. A similar analysis

cannot be applied to 8nn arrays, since even in the worst case, the number of edges in these

arrays exceeds the number of edges in the partitioning graph. An 8nn array is not well utilized by our partitionings,since it has so many unusededges.

4.2. CommunicationRequirements

The cardinalityexpressionsof the previousectionareinteresting,butdo notgiveprecise expressionsforthe communicationoverhead.Inthissection,we obtainupperbounds

for the total cost of communicationwhenrunningour dissectionson 4nn meshes. We

assumethatthesolutionmethodis asfollows.Allprocessecos mputeinparalleland,byconstruction,takethesameamountof time. Thecommunicationstepproceedsasfollows.First

the informationto be communicatedacrossverticalboundarysegmentsis transmitted(first

leftthenright)byeachprocessor.Thenthisstep is repeatedforhorizontalsegments(top

thenbottom).Inspectionof Fig.3.2 revealsthatwhenthe dualgraphof abinarydissectionis

naturallymappedontoa 4rmmesh,adjacentnodesof the dualgrapharemappedontonodes

of themeshthatlieatleastin adjacenrows t oradjacencot lumns.Thismeansthateachcommunicationstepabovehastwophases.Inthe firstphase,datamakesatmostonehorizontal

(verticalmovemen ) to getto thecorreccot lumn(row). In the secondphase,it travelszero

ormorestepsin thevertical(horizontaldi) rectionto get to thecorrectrow(column).The

hardwarecommunicationmechanismateachnodeis assumedto besuchthata unidirectional

transmissionof datacanbe performedon onecommunicationlinkin onetimestep.Thusa

twowayexchangeof dataovera singlecommunicationlinetakestwotimesteps.

Bydefinitionall processorshaveequalamountsof computation.It remainsto evaluate

howmuchcommunicationoverheadisincurred.Noticethatinthe bestcaseof a simpleuniformpartition,thelengthof thesideof anysquaresubdomainis

L,,_¥o_(k,N) __ N

for a total communication time

= 8.2_-. (4.4)

L

for a problem with N points on a side mapped on to a square mesh with 2 2 processors on a
- 14-

side. In the case of non-uniform regions, the degree of distortion of a partitioning determines the time requiredduringthe second phase of a communicationstep. This time is slight

for mildly distorted partitionsbut can be a major factor in partitions with large distortion.

4.2.1. Skewness of a partitioning

To quantify the amount of distortion in a partitioning, we introduce the following concept.

The x-skewness(Sx) of a given partitioning with N x N points and depth k is the ratio of

the length of the longest horizontal side of any subdomain in that partitioning to the length of

the side of the square in the correspondinguniformpartitioning.

The y-skewness(Sy) is similarlydefined for verticalsides.

For example in Fig. 3.3, Sx is about2.5 and Sy about 4.

We work under the assumption that there is always at least one point per processor.

k t_

This constrainsthe skewness to lie between I and2]'(1 - (22 - 1) ). N

4.2.2. Dilation of edges in a dual graph.

Skewness itself does not completely deterraine the communication overhead. A partitioning can be highly skewed yet have a dual graph that precisely matches a 4nn array. On

the other hand, a partitioning can have this same skewness, but with a large mismatch

between the graphof the partitioning and a 4nn mesh.

To more accurately describe this mismatch we define the dilation of an edge in a dual

graph that has been naturallymapped onto a 4nn mesh to be the number of edges that data

passes through between two communicating processors during the second phase of communications.

The x-dilation, dx (y-dilation dy) of a partitioning is the maximum dilation over all edges

in the x (y) direction.

In uniform partitionings (with 5x-$y=l), the maximum dilation is zero. As skewness

increases, the maximum possible dilation also increases. The biggest change occurs as the

t.

22

skewness increases from ! to 2, since the maximum y-dilation increases from 0 to T "The
- 15 -

general expression for the worst case dilation is

[SyI

k

The maximumpossible y-dilation in a partitioningof depth k is 2 2_ I. Here we assume that

the first cut we make in our partitioning is always a vertical line. However it is still possibl€

for a partitioning with very high skewnessto have zero dilation.

4.2.3. Data Transmitted Per Step

The x and y dilations and skewnesses allow us to compute the time required for communication. Each data point can be transmittedto its destination row or column in one time

step in phase I of the communication.Assume the hardware first takes care of all data that

is to move in the x direction and then all data to be moved in the y direction. For a nonuniform partitioning, the maximum number of data points transmitted per communication

step is 281L.n_fo,,.(k,N)in the x direction and 2S>L_von.(k,N) in the y direction, since each

region hastwo vertical and two horizontalsides. The time required for phase one is then no

more than

Tl(k,N ) = 2($x+Sy) 2-_.

Phase 2 of communication is complicated by the fact that there may be several overlapping communication paths in a single row or column. For example if the first processor in a

row is transmitting to the 4th, then the 2nd might be transmitting to the 5th at the same time,

causing congestion. Furthermore, the total time for communication is influenced by the x

and y-dilations.

4.2.4. CommunicationStrategies

We propose two communication strategies for this troublesome phase two of commtmi.

cation. These are the permutation strategy and the pipelined strategy. Both strategies are useful over the range of values of problem size N, depth of partitioning k and skewness S.
- 16-

4.2.$. The permutation strategy.

We may view phase 2 of communicationas the permutation of a set of data on a chain

of processors. Each processor sends a data value to a processor at most dx (d_.)processors

away in the x (y) directions. All processors can send one data value out to its destination in

2*dx or 2*d_.time steps. The constant2 arises because it takes 2 time units for a processor to

receive and transmit one data point. The constant is not 4, which it would be for arbitrary

permutations, because planarityinsures no processorboth transmits and receives in the same

direction. The time for phase 2 of communicationsis thus

Tpenn= 2 * (Max points per side) * ( Max dilation).

for each side. This works out to be

N,

N,

The totaltimeforphase2 ofcommunications isthesum ofthetwoexpressions multipliebyd

2,sinceachregionhastwoverticand al twohorizontasidl es,

Flt_,,=r4"( m $x*xd+ Sy*dy)N

"_. (4.5)

4.2.6. The plpelining strategy

Instead of viewing phase 2 as a sequence of permutations, we can think of it as a

sequence of data transfersin which each processor transmits all of its data points to all processors to which it needs to transmit in a pipelined fashion. That is, if processor 3 needs to

send data to processors 5, 6 and 7, it pipelinesthis transfer so that as soon as it finishes sending off data intended for processor 7, it startssending data for processor 6 etc. In this case it

is impossible for, say, processor 5 to send data to processor 8 (should it need to do so) until

processor 3 has finished. This situation can arise from the partitioning in Fig. 4.1a, where

the configuration gives rise to two separate, overlappingchains of communication. A chain is

a contiguous sequence of processorsin a column(row) that all receive data from a single processor in an adjacentcolumn (row). Fig. 4.lb shows the chains from the column in Fig. 4. la.

Thus there exists the problem of congestion which we define to be the number of over.

lapping chains of communication in a given row or column. As might be expected,
- 17-

Fig, 4.1 Overlapping chains of communicationcause congestion.

congestion varies with dilation. The form of this relationship is somewhat unexpected. As

dilation increases from 1, the congestion increases but then reaches a maximum at half the

maxiumum dilation and starts decreasingagain. The exact expressions are

t t__2

Gx(k)= dx dx _ 2 2

£-1 £-= (4.6)

2 - d, d,,>22

G [d_ t'-I ,.(k)= ' dy< 22 • £-1

t2+-,, ,,>2,

The amount of time required to complete phase 2 using the pipelined strategy is the maximum time for one region to send out its datamultiplied by the maximumcongestion.

The time required is thus

Tptpe= 2*(Max points/side + Max dilation) * congestion, i.e.

Tp_p,,,,(k,N)=2*($=_+dx) * G,(k)

% _ , N * G_(k),.,_(k_v)-2(%2_---_+dy)
- 18-

The totaltimeforphase2 of communicationsis againthesumof thetwoexpressionsmulti.

plied by 2, giving

T,,_t=p.4(Sx2_-- + dx)*Gx()k+ 4($y2-_- +

dy)*G(_k). (4.7)

Comparisonof (4.5) and (4,7) showsthat the permutationstrategyis alwayspreferable

for partitionswith low skewness,whenS<2. For $>2 and low to moderatedepthof parti.

tioning, pipeliningis better, sincein this casethe congestionG is small. Fig. 4.2 shows

graphically the ratio of the respective costs of these communicationstrategies, for a problem

where N = 4096 is the numberof points on a side.

4 I_

I 2.00 F_l < 0.25

1.74 _ 0.25- 0.50

1.52 _I 0.51 - 1.00

1.32 I_ 1.01- 2.00

1.15 II 2.01 - 4.00

1.06

1.05 I >4.00

1.04

1.02

1.01

4 6 8 10 12 14 16 18 20 22

depth of partitioning,k --_

FII. 4.2 Ratio of cost of pipeline to permutation strategy for communicatingprocessors.

$. Mapping onto Trees

In this section we consider the costs ofmapping the partitions ontobinary trees. At

first sight, tree structured multiprocessors would appear unsuitable for grid problems,

because of potential traffic bottlenecks at the root. However, they are natural to consider in

this case since we use a binary decomposition of the domain. Our results show that in the

worst case, with 2k processors, the performance of binary trees is within a factor of a constant times k of the mesh performance.
- 19-

In our mode], the leaf nodes do all the computation, and the rest of the nodes are used

only for communication. Fig. 5.1 indicates how the partitions are matched with the leaf

nodes. Regions that are separated by the last depth k partitioning cut are mappedto adjacent

leaf nodes. Regions separated by the first partitioning cut are in different halves of the tree.

The solution algorithm starts with all leaf nodes computing on their respective subdomains. The communicationstep can be thought of as having k phases. In the first phase,

the leaf nodes send up all data that must pass throughthe root node. This includes all nodes

that border the first, depth 1 partition cut, and takes time proportionalto the length of that

line. In the second phase, it sends data that rises no higher than the two childrenof the root

node. This communication is between nodes sharing one of the two depth 2 boundary segments, and takes time proportional to the length of the depth 2 segment. In the kth and last

phase, leaf nodes sharing a depth k boundary segment swap data. If each phase proceeds to

completion before the next phase starts, the communication time for phase j is proportional

to the length of the maximum depth j segment. The total communication time is then proportional to the sum of these, and has latency k2 throughthe tree. Instead, a leaf node canstart

the next phase of communication as soon as it is ready. This pipelining gives both a smaller

communication time and a smaller latency through the tree of 2k-1. We define a hyperperk

imeter, (in analogy with the perimeter estimates for the meshes), as H(k) = _ l/, where

/-1

l! = any one line segment of depth j. Fig. 5.2 shows the maximum length hyperperimeter

in the given partitioning. If the communication is pipelined, the hyperperimeter is a connected series of line segments. If instead, each communication phase proceedsto completion

before the next phase starts, the maximum hyperperimeterneed not be connected. Instead of

the maximum of the sums, we would take the sum of the maximum length segments for each

depth j.

To have a fair comparison with the performance of meshes, we assume that each nonleaf node can only receive or transmit in one direction over one link in one unit of time.

Thus, a node can receive 3 values on its input lines and transmit 3 values on its outputlines

in 6 units of time. In the worst case therefore, buffers of size no greater than 2N are

required. The total communication cost is then 6"hyperperimeter+ latency. For example,

for a uniform domain decomposition the maximum hyperperimeter is

N N "" + k-'_+ k--_+-_- = N[3---_-]. Thus the total time

N+T+_+ • 2_'- 1 2_--_ 2_- 2_-
- 20 -

0 I 2 3 4 5 6 7 8 9 I0 II 12 13 14 5

Fig. S.I A binary decomposition mapped onto a tree of processors.

T.nb,or, n(k,N) = 18N (I-2 2) + 2k-1. (5.1)

Clearly, this is much worse than the corresponding expression (4.4) for a mesh.

Fig. 5.2 Darkened hyperperimeter shows the maximum communicationrequirement in

the tree.
- 21 -

For a non-uniform partition, in the worst case the hyperperimeter can grow like

N+(N-1)+(N-2) • • • , except there must be at least 1 point per region. For k = 4, for

example, the hyperperimeter is N+(N-2)+(N-2)+(N-3). In general, the expression is

k1

H(k) = kN - 2 2 (k-3) - 3 for a total time of

L

T(k,N) = 6/oi - 3.2 2(k-3)-9+(2k- 1) (5.2)

Notice that the latency is essentially irrelevant. In comparison with eq. (4.7) for nearest

neighbor meshes, where the leading term is 4N, the trees are a factor -_-

worse. However, a

naive analysis of trees gives a worst case bound of 6k2N, which is avoided here by using the

extremely ordered properties of the binary decomposition.

These communicationestimates can be used to determine the optimal number of partitions to use to solve a given sized problem, for a given efficiency. We do a sample calculation for the uniformly partitioned case. For anN by N square grid, the amount of computation to be done using 1 processor is W1= Cz.N2. For a rather simple method, C1 might

include 20 floating point operations per point. For 2k processors, the amount of work per

processor is

W2 = CI.N2/2t + C2.18N,

where we have droppedthe lower order termsin the communicationcost. The speedupis

Wz

W21

and the efficiency is

E= S__= 1

2k C2 18.2_

1+

C1 N

For example, if the communicationtime for one item C2 takes approximatelythe same time

as one floating point operation, then if N = 102, and k _- 6, the efficiency is only 63%. For

an efficiency of 87%, k = 4, or 16 processorsshould be used.
- 22 -

6. Communication Cost Analysis for Hypereubes.

An analysis similar to the one for nearest neighbor arrayscan be performed for hyper-

& k_

cubes. A nearest neighbor array of size 2 2 by 2 2 can be easily embedded in a hypercube of

dimension k using the well known Gray code mappingmethod. Fig. 6.1 shows a 16 x 16 processor array that has been mapped onto a hypercube of dimension 8. All the edges of the

mesh but only some edges of the hypercube have been shown. Specifically, only the edges of

the hypercube that connect nodes that lie along one row and column of the 4nn mesh are

shown. These edges demonstrate the richer interconnection of the hypercube in relation to

the mesh. We would expect this richer interconnect, ion to reduce the communication overhead when running our dissections,

\[\1

N/ Iltl

JJllllll ]/ /1111

,/111111

,/111111

Jlll

IJVI

L,/'II

Fill. 6.1 A 16 by 16 mesh mapped onto a hypercube of dimension 8. Only some hypercube edges are shown.
- 23 -

6.1. M_applngonto Hypercubes.

To map the dual graph of a partitioningonto a hypercube, first map the dual graph onto

a mesh, andthen use the Gray codes to embed the mesh in the hypercube. The cardinalityof

the mapping in this case is no better than that obtainedwhen a 4nn mesh is used. However,

the richer interconnection of the hypercube leads to a smaller communicationoverhead when

a detailed analysis is done.

6.2. Communkatlon strategies.

As before, we assume a processor can send or receive only one item on one ]ink at a

time. We believe this is an accurate model of a scalable multiprocessor. A hypercube that

can transmit on all links at one time is realizable only for a fixed dimension, but cannot be

extended to higher dimensions. We againuse two different communicationstrategies. In the

case of the permutation strategywe take advantage of the logarithmictime to permute data

in a hypercube. For the pipelining strategywe exploit the logarithmicdiameterof the hypercube.

6.2.1. Permutation Strategy.

A hypercube of dimension k can perform any permutationof 2k elements, one per processor, in 4k- 1 time, underour communication assumptions[16]. We can do better than this

for our problem by noting that the edges connecting all the nodes in a single row or column

of a mesh embedded in a hypercube (see Fig. 6.1) form a sub-hypercubeof dimension k/2.

Recall that data flows only along columns or rows duringphase 2 of the communication step.

Thus each permutation step takes at most 2k - 1 time. This is a worst case analysis, how.

ever, and may not be optimal, since it does not depend on the dilation of a given partitioning.

We obtainthe following expressions,

r_..._(kJV)=S.-_y_2 *(2k - 1)

N,

foratotalphaseITtimeof

= +sy) N (6.1)

Thesecorrespondto(4.5)formeshes.
- 24 -

6.2.2. Pipelinlng Strategy

In the case of hypercubes, the pipelining strategy utilizes the logarithmic diameter of

hypercubes to improve communication times. If we examine a contiguous subchain of dx

nodes in a row or column, it is easy to verify that there exists a tree of diameter no more

oo odovo.odo, ovooy o,odo h.'n than 1+

tree may be found by doing a breadth first search outwards from the desired node, staying

within the graph induced by the contiguous subchain. Thus the requirementof the pipelining

strategy that one node transmits to several other nodes in a subchain takes no more than

[21og(dx)] time. The problem of congestion remains. The amount of congestionis pre2+

cisely the same as with 4nn meshes. This is because the problem and not the multicomputer

architecture determines the congestion. In the hypercube, a number of trees are pumping

data out towards their leaves in parallel, instead of a number of chainslinearly pumping data

out towards their end points, as was the case for meshes. The time required is

Tplp,,,(k,N)f(Sx 2_ + 21o82(dx)+2) * Gx(k)

Tplp,_.(k,N)=($y 2_ + 21o82(dy)+2) *

Gy(k)

for a total time of

rIi,tp,

Fig. 6.2 shows the ratio of the pipelining to permutation communication costs in a

hypercube. For any fixed skewness, the pipelining strategy is initially cheaper, but as with

meshes, the permutationstrategy eventually wins with increasing depth of partitioning.

We compare the communication costs in a nearest neighbor array versus a hypercube in

Fig. 6.3. For each type of machine, the cost used is the cheaper of the permutationand pipelining costs, for the given parameters. For a large range of skewness and low depth of partitioning, the mesh is not much worse than a hypercube. For high depth of partitioning and

moderate skewness, the hypercube is much better, but its performance approachesthat of a

mesh for very low and very high skewness. For large skewness, there is a region with a

large perimeter which takes a long time to transmit. The cost of this dominates both the

square root andlogarithmiccommunication latencies of the mesh and hypercube respectively.
- 25 -

1024

512

256

1

64

32

8

l 2.00 ii@i I---1 < 0.25

1.74 _i_i_iii D 0.25 - 0.50

 #:iiiii

1.52 iiiiiiii m 0.51 - 1.00

i.32 i i.Ol-2.00

0

_:_- 1.15 i 2.01-4.00

__ 1.06

co > 4.00

I.05

I.04

1.02

1.01

4 6 8 10 12 14 16 18 20 22

depthofpartitioning,k

Fig. 6.2 Ratioof costof pipelineto permutationstrategyfor a hypercube.

7. Conclusions

We have shown how domains with non-uniform workloads can be partitioned to

equidistribute the computational load. For these types of grid problems, our binary decompositions can be mapped in a natural way onto trees, nearest neighbor meshes and hypercubes.

For 4nn arrays, it is interesting that the map from the problem graph onto the array

graph has at least a 50% hit rate of edges. For an 8nn array this is 79%.

We further evaluated the performance by analyzing the traffic through these networks.

For an N by N problem using 2k processors in the worst cases, the communicationoverhead

using the pipelined strategy, was approximately43/for hypercubes, 8N for meshes, and 6kN

for trees. The performance of trees was found to be better than a naive analysis would suggest. While these results are encouraging, certainly a better approachis to partition using a

weighted sum of computational effort and communicationcosts. In addition, the more difficult problem of adaptive load balancingwill have to confront the problems of modifying data
- 26·

8

4

i

2.00

1.74

1.52

(J)

(J) 1.32

OJ

c: 1.15 3=

0> 1.06

(J)

1.05

1.04

1.02

1.01

D 1.00

[] 1.01 - 2.00

• 2.01 - 4.00

m 4.01 - 8.00

III 8.01 -16.00

• >16.00

4 6 8 10 12 14 16 18 20 22

de ptho f P8 rtitioni ng, k •

Fig. 6.3 Ratio of the communication cost of a nearest neighbor array and a hypercube.

structures within each processor. It will clearly be more efficient to tolerate small amounts of

load imbalance than to change partitions with every perturbation.

We believe the decomposition technique presented here would be beneficial even on

shared memory machines such as the Ultracomputer [10] or the IDM RP3 [14]. Our partitionings allow efficient load balancing across processors, without the overhead of a finegrained queueing mechanism that would otherwise be necessary. They would also reduce

memory traffic and increase the cache hit rate.

Acknowledgements

It is a pleasure to acknowledge several helpful discussions with Bob Voigt, and to thank

Vijay Naik for pointing out reference [16].
- 27 -

References

[1] L. Adams and H. Jordan, "Is SOR Color-Blind?", ICASE Report No. 84-14, May,

1984.

[2] D. Bai and A. Brandt, "Local Mesh Refinement Multilevel Techniques", Weizmann

Institute of Science Report, 1984.

[3] R. Bankand A. Sherman, "AlgorithmicAspects of the Multi-level Solution o Finite Element Equations", CNA-144, Center for Numerical Analysis, University of Texas at

Austin, October 1978.

[4] .LBentley, "Multidimensional Divide-and-Conquer",Comm. ACM 23, (1980).

[5] M. Berger and A. Jameson, "Automatic Adaptive Mesh Refinement for the Euler

Equations", AIAA Journal23, (1985).

[6] M. Berger and J. Oliger, "Adaptive Mesh Refinement for Hyperbolic PartialDifferential Equations", J. Comp. Phys. 53, (1984).

[7] S. Bokhari, "On the Mapping Problem", IEEE Trans. Computers C-30,3 (1981).

[8] N. Deo, Graph Theory with Applications to Ensineerin8 and Computer Science, PrenticeHall, 1974.

[9] D. Gannon and J. Van Rosendale, "Parallel Architectures for Iterative Methods on

Adaptive Block Structured Grids", in Elliptic Problem Solvers I1, G. Birkhoff and A.

Schoeustadt, editors, AcademicPress. 1984.

[I0] A. Gottlieb, R. Grishman, C. Kruskal, K. McAuUffe, L. Rudolph and M. Snir, "The

NYU U]tracomputer. Designing an MIMD SharedMemory Parallel Machine", IEEE

Trans. ComputersC-32,2 (1983).
- 28 -

[11] H. Jordan,"A SpecialPurposeArchitectureforFiniteElementAnalysis",in Proc. 1978

Conf.onParallelProcessing,Aug. 1978.

[12] S. McCormickand J. Thomas, "The Fast Adaptive Composite Grid Method for Elliptic

Equations". To appear in Math. Comp., 1986.

[13] C. Papadimitriou and J. Ullman, "A Communication-TimeTradeoff", in IEEE Proc.

25_hAnnual Symp. on Foundations of Computer Science, (1984).

[14] G. Pfister, et al, "The IBM Research Parallel Processor Prototype (RP3): Introduction

and Architecture", Proc. 1985 Intl. Conf. Parallel Proc, (1985).

[15] J. Van Rosendale, "Rapid Solutionof FiniteElementEquationson LocallyRefined

GridsbyMulti-levelMethods",Ph.D. Thesis,Universityof IllinoisUIUC, May1980.

[16] A. Waksman,"APermutationNetwork",I. ACM15,1, (1968).

[17] P. Zave andW. Rheinboldt,"Designof an Adaptive,ParallelFinite-ElementSystem",

ACMTrans.Math.Software5, (1979).



1. Report No. NASA CR-178024

ICASE Report No. 85-55

I2. Government Acc8U,on No. 3. Recipient's Clotalog No.

4. TItle and Subtitle

A PARTITIONING STRATEGY FOR NON-UNIFORM

PROBLEMS ON MULTIPROCESSORS

5. Repon O.te

November 1985

6. Performing Orpanilltion Code

7. Author!s) 8. Performing Organization Report No.

Marsha J. Berger and Shahid H. Bokhari

10. Work Unit No.

9. Performing Organization Name and Address

Contractor Report

!:;f'lt; .':11 .Q':I.•f'l1

11. Contract or Grant No.

National Aeronautics and Space Administration 14. Sponsoring Agency Code

Washington, D.C. 20546

Institute for Computer Applications in Science

and Engineering

Mail Stop l32C, NASA Langley Research Center UAcl .17f'l7f'l ~--HtTT Htl1&E-eHc:-WTT~-i"l~''l~'l"&-'{,. l"'~ "":3"'~"l'f.z!"l'!-§- "------------------1 13. Type of Repon and Pefiod Covered

12. Spons~fing A';ncy Name and Addi'ess --

15. Supplementary Notes

Langley Technical Monitor:

J. C. South Jr.

Fin;:!l R~nf'lrt

Submitted to IEEE Trans. Comput.

16. Abstract

We consider the partitioning of a problem on a domain with unequal work

estimates in different subdomains in a way that balances the work load across

multiple processors. Such a problem arises for example in solving partial

differential equations using an adaptive method that places extra grid points

in certain subregions of the domain. We use a binary decomposition of the

domain to partition it into rectangles requiring equal computational effort.

We then study the communication costs of mapping this partitioning onto

different multiprocessors: a mesh-connected array, a tree machine and a

hypercube. The communication cost expressions can be used to determine the

optimal depth of the above partitioning.

17. Key Words (Suggested by Authorlsll

partitioning problem, multiprocessors, load balancing,

hypercubes, trees, meshes

18. Distribution Statement

59 - Mathematical & Computer Sciences

(General)

62 - Computer Systems

Unclassified - Unlimited

19. Security Oassif. (of this report)

.Unclassified

20. Security Oassif. (of this pI\IIl

Unclassified

21. No. of Pages

30

22. Price

A03

For sale by the National Technical Informal Ion Service. Springfield. Virginia 22161