Sec. 16-3 Potential Energies with a Short-Range Intefaction: A More Complicated Example 

301 

COOe 16.6 Routine for Finding the Contributioo of a Single Particle to the Interactioo Energy 16-3.4 Parallel Decomposition 
The parallel decomposition of local dynamics problems is very similar to that of 
local finite difference and finite element problems (Chaps. 7 and 8). Such problems are decomposed spatially, so that each processor controls either an equal volume of space or an equal number of degrees of freedom. The update algorithm then follows the general scheme: 
(I) exchange degrees of freedom near edges with neighbors. (2) compute exactly as in the sequential case, using the boundary conditions obtained from neighbors. 
The same approach works for local dynamics problems. Again, the processors are assigned to equally sized regions of space. The degrees of freedom are now the positions of the particles. The particles in cells on the edges of a processor's region will interact with particles in neighboring processors, and it will be necessary to communicate relevant information about them to all processors which might be affected. This com- munication is analogous to the exchange of boundary degrees of freedom in finite differ- ence problems. An important difference is that the number of particles near the edges is not known in advance, so that communication is only need-predictable. In fact, the number of such particles and even their location in memory will vary from processor to processor, as well as from one boundary to another within the same processor. Thus, the algorithm that exchanges particles across boundaries can no longer simply call vm-shift. 

302 

Particle Dynamics: Problems wid1 a Slightly Irregular Communication Structure Chap. 16 

The parallel decomposition divides the cells evenly among the processors and ini- tializes the guard layer in each processor. We will assume that there is available sufficient memory to store all of the additional particles contained in the guard layer. The algorithm for computing the short-range potential energy is given in Code 16- 7. 

proc-begin compute-energy ( evaluate the contribution, V. to the interaction energy arising from this processor ) 

[initialize u = 0] 
for-begin ( each cell in this processor and the guard ring cells ) [initialize a linked list for this cell] for end 
for-begin ( each particle, p, in this processor ) [find out which cell particle p is in] [insert particle into list for this cell] for end 
proc-call bound-exchg ( exchange boundary cell data with neighboring processors, store incoming data in appropriate guard ring cells ) 
for -begin ( each cell resident in this processor ) 
for -begin ( each particle, p i' in the list for this cell) 
proc-call interact ( add contribution to interaction energy due to p i ) for end 
for end 
proc-call vm-combine ( compute sum of partial energies, U, 
accumulated in each processor) 

proc~end 

Code 16.7 Concurrent Code to Evaluate Interaction Energy 

The principal difference between the sequential program of Codes 16-5 and 16-6 and the parallel program of Code 16-7 is the use of the procedure bound-exchg. In the parallel case, the guard layer in each processor is initialized by bound-exchg which exchanges the information about particles near processor boundaries. 
Another addition in the parallel program is the use of the procedure vm-combine . As discussed in Sec. 4-4, vm-combine globally merges quantities across the ensemble, according to a specified combining function. In this case, the function, would perform floating-point addition. Thus, after the call to vm-combine, the energy, U, is the sum of all the partial energies computed concurrently in the processors. 

16-3.5 Communication Routines 
The subroutine bound-exchg is used to exchange all relevant information about particles in boundary cells. The communication arrangement is shown in Fig. 16-4. The guard strip on each boundary is filled with the contents of the true edge of its neighbors. This problem requires diagonal communication, after the manner of Sec. 16-2. As in that treatment, the diagonal moves will be composed of two orthogonal transfers, adding 

Sec. 16-3 Potential Energies with a Short-Range Interaction: A More Complicated Example 

303 

almost nothing to the complexity of the code. The main business of bound-exchg is to deal with each cell on each boundary in turn. The subroutine bound -exchg is outlined in Code 16-8. 

PATTERN FOR COMMUNICATION OF CELLS BETWEEN PROCESSORS 

F1gure 16.4 Rows of cells are mvved vertically and columns of cells are moved horizon- tally so that the extra (dotted) strip of cells on the boundary of each processor is filled with particles fr<Xn neighboring processors. In (a) the downward moving cells are highlighted with \\\ and the rightward moving cells with III. In (b) the upward moving cells are highlighted with \\\ and the leftward moving cells with 111. 
Figure 16-5 shows schematically which cells are communicated when edge-exchange is called for each of the four directions. To see how particles traverse the corners, consider the lower-left cell in the upper-nght processor. In Fig. 16-5(a) it moves leftward into the guard layer of its neighbor. In Fig. 16-5(b), both the original cell and its copy are transmitted down, so that copies of the cell now appear in the guard layers of each of the three neighboring processors. 
Forwarding around corners follows from the fact that the strips of cells communi- cated in Fig. 16-5 extend into the guard layer at both ends. Otherwise, in step (b) the 

304 

Particle Dynamics: Problems with a Slightly Irregular Communication Structure Chap. 16 

proc-begin bound-exchg ( exchange boundary cell data with neighboring processors ) for -begin ( each of the cardinal directions, d ) 
proc-call edge-exchange ( send data in appropriate strip of boundary cells to processor in direction d, store incoming data in opposing guard cell strip ) 
for end proc-end 

Code 16.8 The Routine bound-exchg 

F1gure 16.5 The cootents of lower left corner cell in processor p o is transmitted across: arc (a) when the cardinal direction is LEFT. in processor p 0. arc (b) when the cardinal direction is DOWN. in processor p 0' arc (c) when the cardinal direction is DOWN. in processor p 1. Thus. it reaches the correct cell 00 the boundary strip of processors p 1. p 2. and p 3. 
copy of the cell would not have been propagated downward. Thus, by carefully includ- ing newly received guard cells in subsequent communication, the simple algorithm of Code 16-8 correctly transfers cells around comers. The algorithm of Code 16-8 is a spe- cial case of a general purpose algorithm for long-distance communication in the hyper- cube, the crystal-router , which will be discussed in Chap. 22. In the general case, the loop over directions is not limited to the cardina! directions of a two-dimensional simula- tion, but rather includes all possible directions in the hypercube's hardware topology .As with Code 16-3, and implicitly in Code 16-8, incoming data is copied to an outgoing buffer, as well as local memory, according to whether it should be kept, kept and for- warded, or simply forwarded. 
The procedure edge-exchange is not a simple application of vm-shift because the cells which it must exchange do not reside in contiguous or regularly spaced memory locations. Edge-exchange, in Code 16-9 is implemented with an auxiliary buffer to hold particles scheduled for export, and/or insertion into cells. 

Sec. 16-4 Fully Unpredictable Problems 

305 

Code 16-9 The Routine edge-e~ha/1ge 

16-4 Fully Unpredictable Problems 
A significant class of scientific and engineering problems are need-predictable, and crystalline style "operating systems" are appropriate in such cases. Nevertheless, fully unpredictable problems are also important, and we will comment briefly on them here. Fully unpredictable problems could, in principle, be treated with crystalline communica- tions also. It would be necessary in such a case to poll sufficiently often to be sure that whenever something "interesting" has occurred in another processor, it would not take long to learn about it Interrupt-driven operating systems become useful when the time spent in fruitless polling becomes prohibitive. This is the case, for example, in certain types of Monte Carlo simulations, in which an interesting change in state in a neighbor- ing processor is a relatively improbable occurrence [Johnson 86]. The only way to guard against such an occurrence in a crystalline environment, however, requires that each pro- cessor frequently poll all of its neighbors. This is usually a waste of time, since by the nature of the problem interesting state changes happen only rarely. Thus, the communi- cation overhead due to polling is unacceptable. 
In contrast, an interrupt-based operating system can be used to eliminate the polling overhead. A processor can determine whether any of its neighbors have sustained any interesting state changes simply by asking whether any messages are waiting in a mail- box. The result is the elimination of a large number of null or empty messages. 

306 

Particle Dynamics: Problems with a Slightly Irregular Communication Structure Chap. 16 

16-5 Conclusion 
We have analyzed two types of need-predictable problems, position updates and short-range interactions. In both cases, the fact that the amount of communication could not be determined until run time presented no serious problems. We also found that com- munication between diagonally connected processors could be naturally achieved without direct connections. Forwarding around comers is accomplished by routing incoming data from one direction to the queue waiting to go out in another direction. It is impor- tant to stress that this is done completely within the crystalline environment There is no inherent need for complicated message routing to be performed by an operating system or by higher-Ievel system subroutines. 

17 

WaTor-An Irregular Distributed Simulation Problem 

17-1 Introduction 
This chapter describes a model problem which exhibits many features of more advanced applications on a parallel computer. The algorimm is WaTor [Dewdney 84]: a simple Monte Carlo dynamical simulation in which idealized fish and sharks live, move randomly, breed, and eat one anomer in a two-dimensional periodic ocean. The name WaTor derives from me toroidal topology of me imaginary watery planet. Almough me algorimm may not be relevant to realistic biological population dynamics studies, me program is intriguing and provides a rich learning environment for me parallel program- mer. It serves to illustrate many crucial ideas in dynamic time-and event-driven simula- tion. Our goal in me present discussion is to achieve an efficient optimal implementation on me concurrent machine. Our approach to concurrent decomposition will require me simulation objects to be organized in machine memory as doubly linked lists. This struc- ture, in turn, implies that me order in which me fish and sharks are updated is not predetermined, and mat contention between processors attempting to make mutually conflicting moves must be an issue. The contention problem can be handled by me roll- back strategy .Simply stated, rollback involves taking update steps normally under me assumption mat contention events occur relatively rarely. In me event of a conflict, me sim~ation is backed up in a such a manner as to undo me erroneous step. 
In addition to rollback, this algorimm displays anomer characteristic of many advanced applications: me computational load can be imbalanced because me computa- tional graph (and hence me load in each processor) is changing wim time. In mis exam- ple, me load imbalance comes about naturally because of me dynamics of the problem: the fish and sharks tend to aggregate in schools as mey breed, move, and eat each omer . If we perform a simple rectangular subdomain decomposition of me space (me ocean), some processors may contain very few if any fish and therefore will fall idle. Two memods for counteracting this dynamicalload imbalance will be described. 
In a more abstract vein, WaTor can be thought of as a two-dimensional cellular automaton [Wolfram 86]. By giving each fish and shark a "time to update" parameter during which time the update routine becomes inactive or sleeps, WaTor can be made to behave in a highly inhomogeneous fashion. This gives the simulation some of the characteristics of a wargame simulation, although involving only short-range interac- tions. 

We will begin the discussion of WaTor by specifying the dynamical rules govern- ing the moves of sharks and fish during the simulation. After briefly describing a 
307 

308 

WaTor-An Irregular Distributed Simulation Problem Chap. 17 

simple-minded approach to coding this problem, we will proceed to examine the superior linked list method. The resultant concurrent application (making use of the CUBIX 1/0 system) will turn out to be fully portable between concurrent and single-processor com- puters. We will conclude with an examination of the rollback method and a discussion of load balancing techniques. 

17-2 Dynamical Rules 
The dynamics of WaTor consist of the rules to which the fish and sharks adhere as they move about. The underlying space of the simulation is taken to be a two- dimensional ocean with the periodic boundary conditions of a torus. The domain is discrete in that fish and sharks are allowed to live only on points of a mesh. The simula- tion input data consist of five numbers: n./ish and nshark, the numbers of fish and sharks at the start of the simulation;jbreed and sbreed, the ages at which they breed; and finally, starve, the length of time a shark can survive without eating a fish. Here, time refers to the number of iterations the simulation has undergone, so that the simulation is tem- porally as well as spatially discrete. 
The fish move in a simple way. Each fish examines its four nearest-neighbor loca- tions. If empty locations are found, the fish is moved into one of these sites at random, making use of the prand routine described in Chap. 12. If no empty locations are discovered, the fish remains where it is. In either case, the age of the fish in incremented. In addition, if the age has reached jbreed and the fish has actually moved, a new fish of age O is placed at the old location. 
Sharks are similar to fish, except they are hungry (presumably, the fish are eating universally abundant plankton). This means that sharks look for fish at their four neigh- boring positions. If any are found, the shark picks one at random, eats it, and resets its starve parameter to zero. If no fish are nearby, the shark resorts to moving at random, just as the fish do. Again, the shark breeds if its age has reached sbreed. If a shark goes for starve time units without eating, it dies. 

17.3 The Software and Algorithm 

17-3.1 The Sequential Code 
Before delving into the parallel program, let us first discuss what constitutes an intelligent algorithm for the simulation of this dynamical system. A very simple pro- gramming strategy is immediately apparent whose parallel extension is quite straightfor- ward. However, we shall see that this simple method can be significantly slower than a clever alternative, which will be the ultimate goal of our present discussion. 
The simple method involves storage of all relevant information at the ocean sites (the mesh locations). One could imagine a few two-dimensional arrays: oceanx,y , agex,y , and starvex,y .A fish at location 17 ,26 could be represented by a "1 " in ocean 17,26, a shark by a "2," an empty spot by "0." The update step would then simply require a doubly nested set of loops which run over all locations. 

Sec. 17-3 The Software and Algorithm 

309 

The parallel extension of this method is sttaightforward. By analogy with the Laplace equation solver of Chap. 7, the ocean array is divided among the processors via the usual rectangular domain decomposition. The standard statement. "the parallel code is like the sequential code, but with special boundary conditions," again applies. The processors update their subdomains and communicate when necessary. Since the ocean consists of a regular homogeneous mesh, the problem of preventing processor contention is simply handled in the same manner that we accommodated nonsimultaneous relaxation in Chap. 7. An example of contention in the present context would be two processors trying to move fish to the same ocean location, resulting in an unforeseen collision. With our simple decomposition scheme, a four-color update ordering sttategy, such as that dis- cussed in Sec. 7 -6, or the segregation of edge points would suffice to forestall any such contention problems. 
The drawback of the simple method is that the time to perform an update scales with the size of the ocean. In this approach, the program can spend significant time loop- ing through empty regions of the ocean, doing nothing useful. The standard approach to overcoming this inefficiency is to store the fish and sharks as linked lists. This approach allows the program to step through the active objects, concenttating on useful work only. The result is that the time to update now scales with the number of fish. If there is a mass die-off of fish, the program will run proportionately more rapidly until the population again increases. 
Although the simple method has such a sb'aightforward parallel extension, it will be abandoned here in favor of the richer problem of devising a parallel version of the linked-Iist approach. 

17-3.2 Parallel Implementation 
We begin by defining the data structures that are centml to the linked-list approach. The attributes of a fish or shark are stored in a structure denoted fishy in Code 17-1. 

declare-struct fishy { declare-lnt type; (fish or shark) declare-lnt age; 
declare-lnt starve; (applies only to sharks) declare-lnt locx; (x-coordinate) declare-lnt locy; (y-coordinate) 
declare-lnt nupdate; (update number) 
declare-bur pnext; (pointer to next item in list) declare-bur pprev; (pointer to previous item in list) 

Code 17-1 Fish/Shalk Data Strod.ure 
The type member tells whether this is a fish or shark; members age and starve are the storage locations for the age and starvation factor of the fish; and lacx and lacy give the location of the fish in the ocean. Finally, pnext and pprev are pointers to the next and 

310 

WaTor-An Irregular Distributed Simulation Problem Chap. 17 

previous fish in the (doubly) linked list The other important data structure is that which defines the ocean mesh. Though the update loop of the program will run only through the actual fish, some sort of mesh is still needed. This facilitates the search for neighboring fish. In order to find out where a given fish can legally move, the ocean structure is queried for the existence of fish or sharks nearby. In the event of a positive result, ocean also contains a pointer to the fish or shark in question so that further queries may be addressed. For each location in the ocean, we have the structure defined in Code 17-2. 

Code 17-2 The Ocean Domain Data Stnlcture 
The parallel WaTor algorithm employs the familiar rectangular subdomain decom- position for the ocean. assigning to each node of the concurrent processor a subdomain of ocean and an independently linked list of fish (one which incorporates no links crossing processor boundaries). This decomposition is illustrated in Fig. 17-1. 
The heart of WaTor is the update routine. called move-.fish. It runs through the linked lists and attempts to move/breed each member of the list. The basic outline of the routine is given in Fig. 17-2. In order to update all the fish in its region. each node must know the contents of neighboring boundary strips of the ocean mesh. As a result, the first step. shown in Fig. 17-2(a). involves the interprocessor exchange of these strips of ocean. This array includes an extra set of sites all around the subdomain so as to be able to store the boundary strips. This closely parallels the guard ring approach to the Laplace equation decomposition in Chap. 7. This first stage is accomplished by the subroutine send-bound which. in turn. calls the crOS III routine vshift to do the communication. 
In Fig. 17-2(b). moveysh processes the linked lists and performs updates for every fish and shark in the region. A point to which we will return later is that fish and sharks are updated separately. By this we mean that first all the fish in the region are updated with the sharks held fixed. Then the sharks are moved. holding the fish fixed. Update order is not included as part of the canonical WaTor dynamical rules. so we shall neglect the possible influence of this special choice of update ordering. Different choices of update order will be discussed later. 
As the fish are moved. those which decide to leave the processor subdomain are not immediately communicated. but are instead removed from the linked lists and stored in a separate buffer. They are communicated at the end of the update step. The routine return-bound performs the communication (calling the ~rOS III routine cshift) and is sketched in Fig. 17-2(c). In addition to the actual communication. return-bound also handles the task of receiving the fish and inserting them into the appropriate linked lists of the destination node. It also detects contentions and rebuffers these fish. 

Sec. 17-3 The Software and Algorithm 

(a) 

311 

Typical WaTor Configuration Fish (do's) and Sharks (0) 

1'1gure 17-1 A typical WaTor configuration: (a) WaTor, 00 a 128 x 128 ocean. Fish are represented by dots, sharlts by diamonds, (b) The problem decomposed onto a 4 x 4 mesh of nodes with the fish linked lists drawn. No links cross any processor boundaries. 
The final phase of move jsh is a call to the routine named resolve, which is illus- trated in Fig. 17-2(d). This routine sends back those fish which encountered collisions and attempts to move the fish in a different direction. Contention resolution will be dis- cussed in more detail below. The pseudocode outline of move ysh is given in Code 17- 3. 

The routine send-bound simply consists of four ca11s to the crOS communication routines in order to transfer the boundary strips of the ocean mesh. The structure of update ysh is defined by the implemented dynamical rules outlined previously, with the 

312 

WaTor-An Irregular Distributed Simulation Problem Chap. 17 

Steps in WaTor Algorithm 

(a) Sendbaund Step t+y 

( b) Movefish Updote I "" "" 

(d) Resolution of Collisions 

F1gure 17.2 WaTor algorithm illustrated on a configuration with 5 sharlcs and 54 fish: (a) The sendbound step: a single node is shown. The vshift of +.£ is acbJally shown. Each node sends out its own +.£ strip of ocean and then turns around and receives the corresponding data into the extra padding of ocean. (b) Movefish updates everyone in its subdomain, (c) Sendfish sends and receives buffers, inserts incoming fish, (d) Reso/ve any collisiros by picking roe fish to go back and try again. 
few modifications necessary for the parallel case. While the details of the algorithm may be obtained from a detailed inspection of the Software Supplement, we will roughly out- line the workings here. 
When a shark eats a fish, the first thing that happens is that the fish is removed from its linked list by the routine re~sh .If the fish is in a boundary strip, and is therefore actually stored in a neighboring processor, this is unnecessary. For this reason, before actually calling remfish , we test to determine whether or not the fish is in the boundary strip. If not, only then is the fish removed. In either case, the shark moves to the new location. This movement, done by the routine movf , tests to determine if the fish or shark is leaving (that is, if it has entered a boundary strip). If the fish or shark is leaving, 

Sec. 17-3 The Software and Algorithm 

313 

it is removed from the linked list and buffered in an array of fishy structures, jbuff .The actual communication to neighboring processing nodes is done by return-bound. 

proc-begin move fth ( perfonn an update cycle ) for -begin ( each cell ) 
proc-call send-bound ( exchange boundary fish with other processors ) proc-call update fth ( update the fish in this cell ) 
proc-call return-bound ( return updated boundary fish to their home processor ) proc-call resolve ( resolve contention conflicts ) 
for end proc-end 

Code 17.3 The WaTor Update Cycle 
If a shark is eating an out-of-processor fish, the actual eating of the fish (that is, the rem./ish of it) will occur within the return-bound call in the neighboring processor. In this case, the shark is sent to its destination, the fish to be eaten is removed by rem./ish , and the shark is inserted by add./ish . 
So far, the program has been explained mainly in reference to the parallel machine. The code, however, continues to work on a sequential machine as long as the appropriate I/O subrou.tine library routines are defined In the limit of a single-processor machine, doc, procnum , and other members of the cubenv structure (Sec. 14-5.6) are still defined. These variables simply assume the appropriate values for a single processor: doc =0, procnum =0, etc. The gridmap routines still work, although trivially. The communica- tion routines, though not sending information to another processor, are not entirely trivial in the sequential case. This is because they still must perform memory-to-memory copies from one strip of ocean to another, as shown in Fig. 17-3. If we assume the existence of a dummy sequential crOS / CUBIX library , the very same code is applicable to both parallel and sequential machines. This code, when executed on the sequential computer , will not run at quite the optimal speed for the sequential machine. This is due to the fact that we are spending part of our time doing the memory-to-memory copies. These opera- tions need not be done on the sequential machine and therefore waste time. However, it is also clear that the wasted work is a surface effect and for large problems will be negli- gible. In fact, this asymptotic "immunity" holds for exactly the same reasons that the algorithm is efficient on the MIMD machine in the first place (see also Sec. 3-5 and the discussion of communication in other chapters). In general, we may well conjecture that for any algorithm which is reasonably efficient on the parallel computer, it will also be reasonable to run the same code on a sequential machine with a dummy library of con- current utility routines. 
A certain degree of code portability has been arrived at in this example. It is of more than passing interest to be able to run an application on various machines, some of which may have superior debuggers, output devices, discount prices, and so on. Clearly, we can make no claim to have solved the classic supercompiler problem, in which the ideal is to obtain a parallel application transparently, starting from a sequential program. 

314 

WaTor-An Irregular Distributed Simulation Problem Chap. 17 

Rather we have found that one correct way to write code is to do so with the parallel machine in mind, retaining the option to run on different types of machines (e.g., sequen- tial computers, multiprocessor, shared-memory machines) through the use of simple dummy routines. 

( 

F1gure 17-3 v shift 00 sequential canputer. The v shift of the +1. strip is shown. On the sequential machine, v shift does a memory- to-memory ~y. 
It is tempting to make the sttong assertion that some sort of parallel language / operating system / subroutine library standard needs to be established. Once this is done, codes are written according to this standard and run on all machines of MIMD architec- ture. We have attempted a modest start in this direction with our standard crOS / CUBIX communication system. 

17-4 The Resolution or Collisions and Conflict 

17-4.1 A Simple Technique 
Occasionally during the simulation, conflicts or collisions such as that shown in Fig. 17-4(a) will occur. The problem occurs along processor boundaries: since the boun- dary strips are communicated only at the beginning of an update sweep, the information contained in them becomes out of date. As shown in the diagram, a fish from each pro- cessor has spotted the same empty place in the ocean and both have attempted to occupy it Three-fish collisions, shown in Fig. 17-4(b), can also occur. 
A simple fix for the case of a collision may be incorporated in WaTor. The fish which has crossed the processor boundary and now has no place to go is simply returned back to the processor from which it came. Once returned, an attempt is made to reinsert the fish in the spot from which it came. This will usually be possible, although a double- conflict situation such as that shown in Fig. 17-5 can occur. In the current implementa- tion, if the reinsertion is permitted, it is perfonned. Otherwise, the fish in question is sim- ply thrown away. 
This simple version of rollback is not rigorously correct Due to the collisions and the consequent return of the fish (or shark), they will have a slightly more difficult time 

Sec. 17-4 The Resolution of Collisions and Conflict 

315 

Conflicts in Update 

(a) 2 fish 

... 

F1gure 17.4 (a) A two fish conflict. Both fish are moving to the same empty locatioo. (b ) A three fish conflict. 
crossing over processor boundaries than in making nonboundary moves. In effect, the ocean contains lines (corresponding to the processor boundaries) which exert on fish and sharks a slight repulsive force which tends to keep them in their current subdomain. A more correct rollback code must keep track of the old locations of the fish, so that the multiple rollback drawn in Fig. 17-6 can be performed. Correct rollback involves undo- ing the simulation in an unbiased manner sufficiently far back so that the update may be retried. It is important to note that rollback for this time-driven simulation can retreat onlya limited amount. It is possible (though very unlikely) to encounter the need to roll- back all the fish in the ocean. This is the largest possible rollback however, since the rollback doesn't extend back into previous time steps as is more generally the case for event-driven simulations [Jefferson 85]. 

17-4.2 Alternative Update Strategies 
We have seen that our insistence on writing an efficient, linked-list-based code has led to the somewhat difficult problem of rollback. In the present context, this provides us with an interesting new problem in parallel programming. Our primary motivation has not been to solve the WaTor problem per se, but rather to illustrate the issues common to 

316 

WaTor-An Irregular Distributed Simulation Problem Chap. 17 

F1gure 17.5 The frustrated fish of Fig. 17-4(a) is sent back to its original processor, ooly to discover that its spot has been taken by yet another fish. 

Conflict in Rollback 

, 

.fish ---move -rollback 

r Processor "",=---"Boundary 

Figure 17.6 The correct thing to do is to continue rolling back until everyone has a place togo. 
other important applications. Before pressing on, however, we should mention other pos- sible ways of writing the code. 
The first alternative has already been mentioned in the introduction. In this approach, we abandon the linked lists and store all information in the ocean structure. The program proceeds by looping through this array. In this strategy, the program can be so structured that processors avoid fish collision contention problems. Each node may first update all of its nonboundary domain. Then, all nodes concurrently update their left-hand boundary strips. The processors are kept synchronized by the crOS 

Sec. 17-5 Load Imbalance and Its Resolution 

317 

communication routines. They then proceed to the right-hand, the top, and the bottom boundaries. No contention can occur since the processors are always operating on boun- dary cells well removed from those currently being updated in any other processor. 
A second technique involves keeping the linked-list nature of the current program, but locking the neighboring ocean sites of each fish currently undergoing update. Once the fish is moved, the neighboring sites are unlocked. To be more explicit, consider a node attempting to update a particular fish. Before actually doing anything, flags are set at the neighboring ocean sites which have the meaning: "do not move anything into these sites- they are already under consideration." In general, messages will need to be sent to other nodes to request the current contents and to lock up these neighboring ocean sites. Since this message traffic will be random in nature, it must proceed within an asynchro- nous communications system. If a processor encounters previously locked sites, it notes the fact and moves on to the next fish in the linked list. These failures are later tried again until they are successfully updated. In this way, the processors are able to randomly update their fish. When they bump into one another, only one updates, while the other tries somewhere else. 
Finally, a third alternative approach to writing the WaTor simulation without hav- ing to resort to rollback avoids contention through checkerboard ordering. We again keep the linked-list nature of the program, but, as shown in Fig. 17- 7, we further subdi- vide the domain of each processor into four smaller domains. Each of these subdomains is associated with its own linked-list and the processors start by all updating their lower left-hand subdomains. Communications then occur, having the effect of synchronizing the processors. All processors next proceed to their lower right-hand subdomains and update these regions, and so on. In this way, the updated cells are segregated so that communicated ocean members never become out of date. Thus it is always safe to use a boundary site without concern for changes being made by another processor. 

17-5 Load Imbalance and Its Resolution 

17-5.1 Load Imbalance in WaTor 
Figure 17-8 shows some typical WaTor fish/shark configurations. A clear effect is the tendency of the fish to fonD large schools with the sharks eating around the edges of these areas. The asymptotic configurations typically bear little resemblance to the ran- dom starting configuration. The system develops a coherent structure. This structure results in a perfonDance impact, since with the simple rectangular decomposition so far described, the work load within different processor regions may vary widely. Processing nodes which happen to be responsible for domains containing few fish or sharks will have little to do. They will update their fish and then call the communication routines within return bound .At this point, they will wait until the neighboring nodes catch up to them. We-cannot relax the waiting condition. Otherwise, some processors could get far ahead of others (in tenDS of simulation time) and very costly rollbacks would become necessary . 

318 

subcells-< 

WaTor-An Irregular Distributed Simulation Problem Chap. 17 

Subcells within o Single Processor 

processor boundory 

\L independent linked lists 

FIgure 17-7 A single processor's subdOOlain which has been further subdivided into four smaller subcells. The fish and sharlts in each subcell are stored in their own linked list; i.e., the processor has four sets of linked lists. 
We here encounter an issue of fundamental importance for parallel computing: the case of a dynamically changing, fluctuating load and the problem of keeping the perfor- mance of the parallel computer high throughout such a computation. The problem is dynamical because the computational graph is continually changing with simulation time. Fish are eaten and sharks starve; the schools change shape and move about. If we resort to a more general decomposition than the simple rectangular case, (for instance, a decomposition like that shown in Fig. 17-9), a significant imbalance can quickly ensue after a few time steps. We touched upon the problem of load balancing earlier in this volume, but the present problem offers the best context for a more extensive discussion and treatment of the possible solutions. 
It has been found [Fox 86a] that the problem of finding decompositions can be frtJitfully recast in the form of an equivalent physical system. This is done by identifying a Hamiltonian or energy function which somehow parameterizes the properties of a desirable decomposition. Another name for this Hamiltonian is the objective function which was introduced in Chap. 13. A low-energy decomposition will be a good decom- position. In other words, a decomposition with a well-balanced load and relatively low communication overhead cost would correspond to a minimal Hamiltonian. The problem of decomposition now becomes that of finding a ground state (configuration of lowest energy) solution for the Hamiltonian. 
Before proposing a specific Hamiltonian, let us abstract a bit from the case of WaTor. Consider the fundamental computational objects or elements as particles which are free to move about in the space of the concurrent processor. These computational objects are, in the case of WaTor, individual fish or sharks. In a more general setting, the 

Sec. 17-5 Load ImbaIance and Its Resolution 

319 

F1gure 17-8 Two typical WaTor configurations. 
objects may be independent processes or tasks as viewed by the operating system of the concurrent machine. In this case we are mapping the complex system consisting of the WaTor ocean with fish as members into a complex system of interacting particles (the system studied in the previous chapter). The requirement of load balancing acts as a short-range repulsive force, causing the particles, and therefore the computation, to spread throughout the parallel computer in a homogeneous, balanced manner. The situa- tion is somewhat similar to a gas or fluid filling up a container. This particular analogy , however, is not precise. In a perfect gas, the repulsive pressure which fills the container is due to the random microscopic motion (velocity) of the particles, rather than to any actual repulsive force between them. In the case at hand, we do not want the particles (fish or processes) to retain a memory of velocity. Instead, we want them to move slowly 

320 

WaTor-An Irregular DistributOO Simulation Problem Chap. 17 

Approximately Load Balanced Domoin Decomposition 

FIgure 17.9 A more general decomposition of the WaTor configuration shown in Fig. 17-8. This deCOOlposition was drawn by hand and is meant to be approximately load bal- anced. 

so that they stay put in processors unless forced to move on. A better analogy, therefore, is that of particles interacting via a repulsive force with the system at a low temperature. 
The requirement in conflict with that of load balancing is interparticle communica- tions. The various elements (e.g., neighboring fish) of the computation need to commun- icate with one another. If the particles are far apart (distance being defined as the number of communication steps between them) large delays will occur, slowing down the com- putation. We therefore add to the physical model a long-range attractive force between those pairs of particles which need to communicate with one another. This force will be made proportional to the amount of communication traffic between the particles, so that heavily communicating parts of the computation will coalesce and remain close to one another in the concurrent computer. Figure 17-10 shows the communication require- ments for partial differential equation-based, finite difference algorithms as well as for irregular short-range simulations such as WaTor. 
There are many different ways to search for low-energy solutions to Hamiltonians. We will describe two approaches of general applicability below. 

17 -5.2 Simulated Annealing Approach to Load Balancing 
The method of simulated annealing [Kirkpatrick 83] has already been discussed within the context of the traveling salesman problem in Chap. 13. Here we will describe how an annealing Monte Carlo method can be used to find and adapt decompositions for irregular problems such as that posed by WaTor. This section describes a reasonable choice of the objective function (or Hamiltonian) for the annealing and mentions some 

Sec. 17-5 Load Imbalance and Its Resolution 

321 

(b) WoTor Interconnection Pattern ~ <>-,-e 

r 

0-1 

Interconnection becomes Communication when links cross processor boundaries 

FIgure 17.10 (a) The communication requirements for a two-dimensional partial dif- ferential equation-based problem. The dots represent degrees of freedom located at the grid sites. They interact in a nearest neighbor fashion and this is represented by the con- necting lines. These lines form the "communication pattern" of the algorithm. (b) The communication pattern for WaTor. 
additional issues pertinent to dynamical cases such as WaTor. The technique has an obvious motivation in at least one sense. We are seeking more general decompositions than the simple rectBngular case, allowing the decomposi- tion to dynamically adapt to the simulation by following the schools of fish as they develop. The crux then becomes: "How do we effectively and rapidly readapt the decomposition?" Simulated annealing provides a powerful general answer to this ques- tion. This point conceded, we are forced to ask, "Can we anneal rapidly enough to be able to closely track the computation?" Remember, the annealing steals computational cycles away from the principal computBtion tBsk of interest. 
Consider some general irregular problem executing on the MIMJ;> parallel machine with decomposition of the general type shown in Fig. 17-9. The Monte Carlo annealing procedure starts by quantifying the current effectiveness of the decomposition. If serious load imbalances or unnecessarily long-distance communications are in evidence, the sys- tem suspends the computBtion and proceeds to anneal the decomposition. This is accom- plished by the Metropolis method [Metropolis 53]. Trial moves are offered which are unconditionally accepted if they lower the energy. If the energy is raised (i.e., heading toward an even worse decomposition) the move is rejected except for a fraction of the attempted moves. The probability of any such move's acceptBnce decreases exponen- tially with the magnitude of the energy increase. It is important to note that it may be relatively expensive to actually move the fundamental particles or processes between processors. A significant improvement, therefore, involves moving tokens of the processes which consist of small datB packets sufficient to describe the essential features of the process. The annealing proceeds for some length of time, with the aim of 

322 

WaTor-An Irregular Distributed Simulation Problem Chap. 17 

substantially improving dIe decomposition. At dIe conclusion of annealing, dIe processes are actually moved to dIe final locations of dleir corresponding tokens and dIe computa- tion is restarted. 
To further take into account dIe costs associated widl moving processes, it is advis- able to add a third type of force to dIe Hamiltonian in addition to dIe repulsive load- balance term and dIe long-range attractive communication term. This dlird term is a force which attracts dIe process to dIe processor in which it is currently residing. We can picture it schematically as a spring connecting dIe process to its current location. The rationale for dlis kind of term is to prevent dIe system from wasting time making moves of only marginal gain. If it doesn't matter much whedler a process is in processor 21 or in processor 22 one step away, it is generally better to leave it where it is. 
In dIe simulated annealing approach, dIe annealing algoridlm acts analogously to a physical heat badl, widl dIe task of keeping dIe computation "cool." It is important to note dlat dIe annealing procedure can accomplish only so much. If annealing is allowed to continue for a very long time (attempting to reach dIe optimal solution) dIe cost is too many computational cycles stolen away from dIe actual calculation. This realization sug- gests dIe importance of dIe concept of temperature in dIe annealing computation. Given a specific algorithm, dIe annealing procedure can keep dIe computation only so cool. The minimum achievable temperature is characteristic of dIe underlying algoridlm and is independent of dIe parameters of dIe computer used. More detailed discussions of simu- lated annealing as applied to dIe problem of decomposition can be found in [Fox 85a, 86a,86g] and as applied to dIe finite element medlod in [Flower 86]. 

17-5.3 The Scattered Decomposition Approach to Load Balancing 
The scattered decomposition is an effective solution strategy for the decomposition problem posed by the large class of irregular problems of which WaTor is an example. The scattered decomposition has also been successfully applied to the parallel solution of irregular finite element problems [Morison 86] . 
The concept is illustrated in Fig. 17-11. As in the usual rectangular decomposition, we map the parallel machine onto a two-dimensional mesh. Instead of simply superpos- ing this r:.esh on the ocean of WaTor, however, we subdivide the ocean into many rec- tangular regions and subsequently distribute these subregions or templates among the processors of the concurrent machine. Each processing node ends up being responsible for many small mutually disjoint regions, which we will refer to as granules. Figure 17- 12 shows, for example, the granules for which processor O would be responsible in a typ- ical scattered decomposition. 
The scattered decomposition accomplishes load balancing for all but the most improbable of possible WaTor configurations. Neglecting any conspiratorial arrange- ment of fish, each processor has a piece of the action across the entire space of the prob- lem. In essentially any simulation run, each processor will end up with roughly the same work load. This spatial averaging approach to load balancing works best as the granule size decreases, although the communication overhead increases in this limit. 

Sec. 17-5 Load Imbalance and Its Resolution 

Scattered Decomposition for 16 Processors 

D Granule 

Template 

323 

F1gure 17-11 The scattered decanposition. Shown is a 4 x 4 array of processors over- laid over the problem 16 times. That is, there is a 4 x 4 array of templates, each of which consists of a 4 x 4 set of granules. 
Another appealing feature of the scattered decomposition is that the required com- munication pattern between the nodes retains a regular rectangular topology. Since this regular topology is often directly reflected in the machine hardware and the associated low-level message-passing system (e.g., crOS III), communications will generally be fast in this method. On the other hand, as mentioned above, small granule sizes yield a greater communication perimeter relative to the contained area of useful work (i.e., the surface-to-volume ratio is increased). As a result, it is not immediately clear whether scattered decomposition leads to higher or lower absolute communication overhead costs. The answer to this question is dependent upon the choices of both machine and algo- rithm. 

There remains one important point to be discussed with regard to coding of the scattered decomposition algorithm. A blunder in programming technique can actually prevent an algorithm which employs the scattered decomposition from achieving load balance. Care must be taken not to impose excessive synchronization conditions during update of a template. If communication-imposed synchronization occurs after each granule of the template, the program is self-defeating, in that at any given time different 

324 

WaTor-An Irregular DistributOO Simulation Problem Chap. 17 

Role of Processor O in Scattered Decomposition 

D Granule 
~ Granule associated with Processor O 
Te~late 
F1gure 17-12 The granules for whim processor O is respoosible. 
processors may be working on granules with very different loads. In fact, the result will be a worse load imbalance than without the scattered decomposition! In order to correctly achieve load balancing, we need to allow processors to move immediately to the next granule without waiting, so that the load balance is averaged over the whole problem. The processors run free without communication until reaching the last tem- plate. As a result, the communications between all of the granules throughout the ocean must be buffered and stored as the processors sweep through the ocean, and only at the very end of an update cycle are the communications actually performed. 
Throughout this book, we have emphasized the need for MIMD hardware; we must be careful not to shoot ourselves in the foot on the software level by mistakenly coding in a manner more appropriate to SIMD machines! 
As described so far, the scattered decomposition algorithm still must deal with the contention problem. One solution is rollback. Another, possibly simpler solution is to have processors skip over some of the templates. The update now proceeds as a sequence of "mini-cycles". As an example, in the first mini-cycle, processor ° would update its granule in template 0,8, etc., processor 1 would update its granule in template 1.9. and so on. The pattern is chosen so as to load balance and keep the processors far 

Sec. 17-5 Load ImbaIance and Its Resolution 

325 

apart during any mini-cycle. As we proceed through the mini-cycles, the patterns change and are chosen so as to eventually update all the granules. 
It seems that we have hit upon the solution to all of our load balance/conflict prob- lems. By virtue of the scattered decomposition, the computation is well load balanced and a straightforward method for resolving the conflict problem fits easily within this approach. In any practical application, however, the problem-dependent and machine- dependent factors influencing the solution of irregular problems may tilt the balance in favor of one or another of the alternative approaches.