Sec. 15-2 Specification of CUBIX 

277 

after the point where the previous one finished. Each call is completed in all processors before the next begins. Thus, if Code 15-3 were provided with the input: 
01234567 
the following variable bindings would result on a four-processor system: ;0= 0, jo=4 ;1=1, jl=5 ;2=2, j2=6 ;3=3, j3=7 
where subscripts refer to the different values of the variables in the different processors. Input functions must be called loosely synchronously, even if their stream argument is in multiple mode. The stream argument must be identical in all processors, although the other arguments need not be. 
The function ungetc presents a problem for multiple input. It is only possible to replace one character on the input stream, even if the stream is in multiple mode. Thus, ungetc replaces the last character read by the last processor, regardless of the multiplicity of its stream argument In practice, this is not a serious limitation. Both of ungetc's arguments must be identical in all nodes. 
Since output is buffered, it is possible to perform multiple output completely asyn- chronously (at least to the extent of putting the data into the buffer). Thus, putc , putchar , printf , fwrite , etc., may be executed completely asynchronously on multiple streams. Code 15-4 demonstrates a simple error-reporting mechanism which makes use of this feature. Even if the error condition occurs in only one processor, the call to fmulti places stderr in a state capable of handling such a situation. Of course, fmulti may be called at any time prior to the fprintf, and it need only be called once, allowing for numerous fprintfs , putcs, etc. 

Code 15-4 Example of Error Reporting Usingfmulti 
The example also includes a call to ffiush which is not necessary when using singular output Multiple output buffers are never flushed automatically. This has three important consequences: 
I) The user MUST callfsingl,fmulti ,ffiush ,fclose or exit to see anything written to a multiple output stream. These functions must be called loosely synchronously, even if the named output stream is in multiple mode. 

278 

CUBIX-lmplementation Issues Chap. 15 

2) In order to avoid errors, the user must be certain to call ifiush frequently enough to insure that the buffer does not overflow. Failure to obey this rule will result in the loss of data. 
3) At dIe time that ising/. .fmu/ti, fflush. idose, or exit is called, the buffers are flushed to the output device in order of increasing procnum .This fact may be used to control the interleaving of the output from different nodes, as was demonstrated in Code 6-2. 
Buffers are normally BUFSIZ bytes in length, and are allocated automatically by the system. BUFSIZ is a constant defined in the header file stdio.h, and is normally between 512 and 4096. The function setbuffer which is defined in Code 15-5 may be used to specify a new buffer of arbitrary size. It is similar to the function setvbui, which is defined in [ANSI 87], except that it does not change the buffering policy associated with the given stream. Normally, but is obtained by calling ma//oc , or by declaration as an external array. Setbuffer is useful for providing CUBIX with large user-allocated buffers which will not overflow at inopportune times. 

Code 15.5 Specification of setbuffer 
The Standard 1/0 Library provides for unbuffered output streams in which a read/write system call is generated for each byte of 1/0. The stream stderr is unbuffered by default, and stdout is unbuffered if it is connected to an interactive terminal. Nor- mally buffered streams may be forced into an unbuffered state by use of the routines set- but or setvbut. Unbuffered streams are supported by CUBIX with the following provi- sion: when an unbuffered stream is converted to multiple mode, it becomes a buffered stream with an automatically allocated buffer of the default size, BUFSIZ. Multiple streams revert to their unbuffered behavior upon calling fsingl , if they obtain their buffer by this mechanism. Thus, although stderr is normally unbuffered, fmulti (stderr) will work correctly, and automatically provide a buffer of BUFSIZ bytes. A subsequent exe- cution of fsingl (stderr ) will cause it to revert to its unbuffered state. 

15-3 Implementation 
The implementation of CUBIX is divided roughly into two parts. First, we describe a universal control processor program and a control processor/node protocol that provides the concurrent program access to all of the control processor's operating sys- tem. We must arrange that the concurrent program can make the system calls in Table 15-1. We also introduce the two new routines mread and mwrite at this level. 
Second, we must implement the C Standard 1/0 Library in terms of the system rou- tines just mentioned. Fortunately, the bulk of the work has already been done, because the C Standard 1/0 Library has, presumably, already been implemented in terms of the 

Sec. 15-3 Implementation 

279 

operating system primitives. Of course we must modify the sequential implementation so it distinguishes between multiple and singular I/O, and uses mread and mwrite when appropriate. We will restrict our attention to the extensions. 

15-3.1 Implementation of System Calls 
The first step in the implementation of CUBIX is to get the system started. Code 15-6 is an outline of the general approach. A CUBIX process is begun on the con- trol processor by running the main program of Code 15-6 with command line arguments indicating the desired cube dimensionality, the file to be loaded into the cube, and argu- ments to be passed to the cube program. After parsing the command line, the program of Code 15-6 calls cubeldve which starts the application's main routine in the cube, and passes the arguments requested on the command line along with the environment of the control process. 
The control processor now effectively gives up control and enters a loop which serves requests from the nodes. The loop terminates when the nodes request the exit sys- tem call, which terminates the control program. 

Code 15.6 Main Rootine for Universal Control Processor Program 
The subroutine serve-request reads one request from the concurrent processor, identifies what system service is being requested, performs that service and provides appropriate data, and return codes to the concurrent process. Such a routine is a standard programming exercise and can be written in many ways. A very elegant technique uses an array of pointers to functions, each of which is responsible for collecting and dissem- inating the data for one of the system calls. An outline is given in Code 15- 7. 

280 

CUBIX-Implementation Issues Chap. 15 

CcKIe 15-7 Control Processor Subroutine serve Jequest 
The various functions in the table are responsible for obtaining their arguments from the nodes and returning their results. Node programs will identify the control pro- cess program to execute by providing serve-request with the appropriate index into the array functbl. A header file, e.g., funcnums.h , can be used to consolidate this informa- tion. 

Code 15.8 The Header File funcnwns.h 
An implementation of the system call creat is shown in Codes 15-9 and 15-10. All of the routines listed in Table 15-1 can be implemented according to the tem- plate of Codes 15-9 and 15-10. The two new routines mread and mwrite are more com- plicated because each node must communicate with the control processor to either read or write distinct data. The crOS III routines fdumpcp, mdumpcp. dumpelt. ftoadcp. mloadcp. and loadelt are useful in implementing mread and mwrite. An implementation of mwrite is shown in Codes 15-11 and 15-12. 

Sec.15-3 Implementation 

281 

Code 15.9 Implementation of creat in Control Processor 

15-3.2 Implementation of the C Standard 1/0 Library 
The C standard I/O library is a collection of routines designed to provide the C pro- grammer with a fast, portable, and reliable input/output facility. The library.itself is, of course, operating-system dependent, but the interface seen by the C program is not. Thus, programs written using the standard 1/0 library are potentially operating-system independent. In practice, the vast majority of C programs use the standard I/O library , so it is of the utmost importance to implement it efficiently. 
In this section, we are interested in the implementation of the library , and espe- cially how it can be modified to handle the extensions of Sec. 15-2.2. Rather than embarking on a detailed exploration of the issues involved in a sequential implementation of the library , we will assume that the reader already has such an implementation. If the sequential library has a clear and simple design, only a very small amount of code will have to be changed to support the parallel CUBIX extensions. 

282 

CUBIX-Implementation Issues Chap. 15 

COOe 15.10 Implementation of creat for Nodes 
It is reasonable to assume that the sequential library makes use of system calls in only a few places, i.e., in a handful of common subroutines used by all the functions accessible to the user. Typically, the write system call is made only within fftush and a common routine Jsbuf, which is called by the putc macro. Similarly, the read system call is made only within ~lbuf , which is called by the getc macro. Hopefully, the only routines needing non-trivial modifications will be Jsbuf, fftush, and ~lbuf. Of course, we will also need to write .fmulti,fsingl, and ismulti from scratch. 
The first thing we must add is a way to distinguish between multiple and singular streams. All that is required is that a single bit in the stream structure be reserved for this purpose. The structure (in UNIX implementations) already has a field called Jag, which is treated as a bit-field and stores infomlation such as whether the file has been opened for reading or writing, whether 1/0 is to be buffered, whether any errors have occurred accessing this stream, etc. Adding a new bit-field to Jag is trivial. For example, ismulti can be implemented by placing Code 15-13 in stdio.h . 
The functions fmulti and fsingl are not as simple. Complications arise because of the need to allocate a new buffer when an unbuffered singular stream is switched into 

Sec. 15-3 Implementation 

283 

#include <cros.h> 
#define MAXPROCS 256 1* maximum number of processors in cube *1 #define HUGEINT Ox7fff 1* a large integer, for 16-bit machines *1 

Xmwrite() 

int fd, i; 
int nload, ndump; 
int szarray[MAXPROCS+l]; 
int array-of-intsz[MAXPROCS+l]~ struct{ 
int ret; 
int errno; } reply; 
struct cubenv ce; 

cparam( &ce ); 
cread(&fd, ce.cpmask, 0, sizeof(fd»; 

ndurnp = fdurnpcp(fd, HUGEINT, szarray); 

array-of-intsz[O] = ce.nproc; for(i=l; i<=ce.nproc; i++) 
array-of-intsz[i] = sizeof(int); 

nload = mloadcp(&szarray[I], array-of-intsz); if(nload < 0 II ndump < o){ 
reply.ret = -I; 
reply.ermo = errno; }else{ 
reply.ret = 0; 

} 
bcastcp(&reply, sizeof(reply»; 

COOe 15-11 Control Processor Implementation of mwrite multiple mode. Correctly programming tl1is feature requires careful study of how and when buffers are allocated and freed in tl1e sequential implementation. It is necessary to add additional bit-fields to -.flag to record tl1e ownership and status of such buffers. 
The only difference between multiple and singular output streams is manifested when tl1e buffers are flushed. Witl1 tl1e mwrite and write system calls, tl1is difference is trivial to include. If tl1e stream is multiple, tl1en call mwrite ; otl1erwise, call write. A small additional complication is tl1at fjlush does not necessarily call write at all if tl1ere is notl1ing to be written. This is a perfectly acceptable timesaving measure for singular files, because all processors are guaranteed to have notl1ing to write if anyone of tl1em has nothing to write. Unfortunately, it is not acceptable for multiple files, because even if 

284 

CUBIX-Implementation Issues Chap. 15 

Code 15-12 Implementation of mwrite for Nodes 

Code 15.13 Definitions in stdio.h to Implement ismulti 
one processor has nothing to write, other processors might still wish to write a nonzero number of bytes. In order to satisfy loose synchronization, the processor with nothing to write must also call mwrite , with cnt set to 0. 

Sec. 15-3 Implementation 

285 

15-3.3 Standard Input Functions 
Input streams are significantly more difficult than output streams. The difficulty arises because multiple input involves global behavior. What is read in one processor depends in detail on what has been read already by lower-numbered processors. Nevertheless, we would like to use as much code as possible from a sequential imple- mentation of the Standard I/O Library. We can use most of the Standard I/O Library input routines, with little modification, if we add a subroutine to the beginning of each one which waits for all the processors with lower processor numbers to finish reading their input When they are done, we adjust the stream's internal pointers in the next pro- cessor so that the next byte read from the stream is correct Code 15-14 shows how we propose to change the function fscanf so it performs properly with multiple streams. Within the function fscanf the macro getc should behave exactly as in the sequential library .This is in contrast to its behavior in an application program, which depends on the multiplicity of its argument To accomplish this, we define the symbol-LIBRARY in the code for fscanf and all other library routines, and make the change to stdio.h indi- cated by Code 15-15. The function mgetc, referred to by the getc macro in Code 15-15 is shown in Code 15-16. The constants -DONE and -REFll..L, and the structure cmd-struct will be used shortly. 

Code 15.14 Implementation of CUBIX scant 
The function wait, shown in Code 15-17, consis~ of a loop which tenninates when each of the processors with processor number lower than procnum has executed a call to release. During execution of this loop, one of the lower-numbered processors might exhaust the buffer associated with the stream and call ylbuf to refill it. Since ylbuf calls read, a system call in Table 15-1 to fill the buffer, we must arrange that all of the other processors also call read, or we would violate the rule of loose synchronization. 

286 

CUBIX-lmplementation Issues Chap. 15 

Code 15-15 Code added to stdio.h to Facilitate Multiple Input 

#define -LIBRARY #include <stdio.h> 

int -mgetc(stream) FILE *stream; { 

mtc; 

wait(stream); 
c = getc(stream); release(stream); return(c); 

} 

Code 15-16 Irnplernentation of Multimode getc 
This is the purpose of the broadcast and the conditional call to refill in the do-wbile loop in wait. After releasing control to higher-numbered processors, the function release , shown in Code 15-18, waits for the rest of the processors to call release, calling refill, in Code 15-19, whenever necessary. When the last processor has announced completion, release repositions the stream so all processors agree on the next character to be read. The function -.filbuf is outlined in Code 15-20. In several places, the code is strongly dependent on the details of the implementation of the Standard I/O Library. In these places, we simply indicate in italics what should be done. With these routines and an implementation of the Standard I/O routines, multiple standard input is accomplished. All of the input routines can be implemented by analogy with the template for fscanf in Code 15-14. 
The use of wait and release minimizes the modifications necessary to convert a sequential Standard I/O Library into a concurrent one, but it does so at the expense of speed. All the data is copied into all of the processors by refill, even though only one processor intends to use the data. In the case of scant this is probably of little 

Sec. 15-3 Implementation 

287 

consequence because the computation required to unformat each byte of input data is substantial. The performance of other functions, e.g., fread, can be severely degraded, however. An alternative to the use of wait and release is to use mread whenever the number of bytes to be read is known in advance, i.e., infread, getw, and getc. The use of wait and release are still necessary for library routines which read an amount of data that depends on the data itself, e.g., scant and gets. 

288 

CUBIX-Implementation Issues Chap. 15 

#include <cros.h> #include <stdio.h> 

wait(stream) FILE *stream; { 

slruct cmd-str cmd; int mask; 
slruct cubenv ce; cparam(&ce); 
mask = ce.nproc -1; 1* broadcast on all channels *1 for(source=O; source<ce.procnum; source++){ 
do{ 
broadcast(&cmd, source, mask, sizeof(cmd»; if(cmd.type = -REFILL) refill(stream); 
}while(cmd.type 1= -DONE); 

reposition the read pointer in the stream to cmd.newposition. 

} 

Code 15.17 Implementation of wait 

Code 15-18 Implementation of release 

Sec. 15-4 Conclusion 

289 

refill(stream) Fll..E *stream; { 

issue a call to read which fills the buffer associated with stream 

} 

Code 15.19 Implementatioo of refill 

Code 15-20 Irnplemen!Jition of JIlbuf 

15.4 Conclusion 
It is not difficult to implement a universal server process to run in the control pro- cessor of a concurrent computer which transparently provides system services to the nodes. Once this is done, libraries that provide a standard interface to the operating sys- tem can be executed in the nodes. The resulting environment is easy to use, familiar to the experienced programmer, and a significant improvement over the control processor/node environment described in Chap. 14. 

16 

Particle Dynamics: Problems with a 
Slightly Irregular Communication Structure 

16-1 IntrOOuction 
The problems discussed so far have involved predictable overhead, in the sense that the programmer could determine a priori when and how much interprocessor com- munication would be necessary .We have further seen that the class of predictable prob- lems is surprisingly rich. Nevertheless, there are many problems in which the interpro- cessor communication is not so easily predicted. We shall refer to the simplest class of nonpredictable problems as need-predictable. In need-predictable problems, the need for communication is predictable at compile time, but the details, such as exactly which data will be transferred, are not known until run time. This case stands in contrast to fully unpredictable problems in which even the necessity of communication is not known beforehand 

16-2 Position Updates: A Simple Example 
A simple example of a need-predictable problem is the position update phase of an N-body particle dynamics simulation. Suppose that the positions and velocities of the particles are known at some time and the positions are to be obtained at some later time. In general, we can relate the updated position, XIleW' to the previous position, Xold by Eq. (16-1): 

~ 
Xnew=Xo/d+~(X,V,...) (16-1) 
UsuallJ' ' I ~ I is small in comparison with other relevant length scales (in the simplest case, ~ = v dl, with dl a small time step) , so the particles don' t move very much in each time interval. Now, let us suppose that a spatial decomposition has been adopted in which each processor controls a fixed region of space, and all the data connected with particles in that region are stored in the corresponding processor's memory .Note that this is not necessarily a good decomposition for all types of dynamical simulations (see, for example, Chap. 9 on long-range interactions). 
In order to maintain the decomposition after the position update, some particles will have to move between processors. The precise identities of the movers is not known until run time. Nevertheless, the programmer knows that some communication must be done after each update, even if it is only to say, "I have no particle to be moved." Thus, the need for communication is predictable. 

291 

292 

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

Before discussing the position update problem in the context of crOS Ill, we may see how it is solved with a more general operating system. For this particular problem, the general operating system utilities, such as message queuing and routing, make the problem somewhat easier to program. As mentioned in Chap. 14, these utilities are costly in terms of performance, however, so it is important to ask whether crOS III (crystalline communications) can be used instead. The pseudocode fragment in Code 16-1 shows how one might update particle positions using a general operating system which automat- ically routes and queues messages for subsequent processing. 

Code 16.1 Updating Particles 00 an Advanced Operating System 
Code 16-1 contains perhaps the most straightforward solution of the problem as seen from the programmer's point of view, and we present it here as a point of departure for the rest of the discussion in this chapter. The algorithm could be improved in several ways. However, it is not our intention here to develop techniques for using general operating systems. Message forwarding and queuing are handled entirely by the operat- ing system. The program is divided very clearly into two phases: 
1) an update/scan phase in which particle positions are updated, and those which need to be moved are extracted from the local data structures and forwarded via the operating system 
2) a receive/cleanup phase in which incoming messages are interpreted and incoming particles are placed in local data structures. 
Unfortunately, simplicity has its price: general operating systems are slow because they must handle general message routing and queuing operations. For a need- predictable problem, the power of a general operating system is wasted. In the following section, we solve the same position update problem using the substantially faster 

Sec. 16-2 Position Updates: A Simple Example 

293 

crOS III polled communication scheme. 

16-2.1 A Polled Algorithm 
For concreteness, we will restrict the problem under consideration to two dimen- sions; the generalization to three or more dimensions is not difficult We will also assume that I~I is always smaller than the smallest linear dimension of any processor's region. Thus, a particle can never jump completely over a processor in a single step. This assumption is not necessary for the correctness of the algorithm that uses automatic routing and queuing, but relaxing the condition would substantially complicate the cry- stalline program. With this restriction, communication is only necessary between nearest neighbors in a rectangular two-dimensional grid. Of course, diagonal moves are also possible, and we must either provide diagonal communication paths in hardware or route particles around comers in software. We elect to route particles in software because the hypercube topology, commonly available in hardware, contains the Cartesian links, but not the diagonal links. 
The overall approach is similar to that taken earlier with the advanced operating system. The problem is divided into an update/scan phase and a receive/cleanup phase. An important difference is that we don't simply give the particles to the operating system during the update/scan phase, and rely on the system to manage sufficient queue space and the like. Instead, we transfer the items from the local list to one of several export queues. We will have one export queue for each communication channel, for a total of four in the two-dimensional problem. The application process manages the queues, and communication is only done at predetermined times, when both sides of the channel expect it The algorithm is given in Code 16-2. 
This algorithm is incomplete, however. The routing strategy necessary for moving particles around corners is hinted at, but not described. Forwarding is accomplished by conditionally placing incoming particles (those from the import queue) into a new export queue. This operation is in error, however, if the new queue corresponds to a direction in which communication is already complete. We must design the algorithm to make sure that this never happens. Codes 16-3 and 16-4 illustrate a correct routing strategy. 
We may consider a particular case to see how this algorithm guarantees that parti- cles will be forwarded properly. In Fig. 16-1, the position of particle p has changed so that it must move from processor p 1 to processor p 3, which is up and to the left from its initial location. We could move it there in two ways: up first and then left, or left first and then up. Since the constant LEFT, whose value is 1, is less than the constant UP, whose value is 2, loop 1 will export particles to the left before they are ex~orted up. Thus, it is imperative that p be placed on the left export queue first This is guaranteed by loop 2 in next-dir, which tests for export to the left before it tests for export up. Correctness is guaranteed for all directions because loop 1 and loop 2 cycle through the directions in exactly the same way. If a particle must go around a comer, it is always placed in the queue which will be exported first 
To see how routing works, let's follow particle p after it is exported to the left Based on the assumption that I~I is small, we can be sure that p is now in the correct 

294 

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

proc-begin update yarticles ( update particle data ) 
declare-buf export ( four buffers for storing outgoing particles ) declare-buf import ( four buffers for storing incoming particles ) for -begin ( each particle, p ) 
[update particle data for p] 
if-begin ( p is no longer in this processor) then [remove p from the local list of particles] 
[identify first direction, dir, in which p must move] [place data for p in buffer export[ dir ]] if end for end 
for-begin ( each direction, dir ) 
[write buffer export[ dir ] to neighboring processor in direction dir , while receiving incoming particles into buffer import[dir]] 
for -begin ( each particle, p, in import[ dir ] ) 
if-begin ( p belongs to this processor) then [copy p to the local list] 

else 

[identify next direction in whichp must move, nextdir] [copy p to export[nextdir]] 
if end for end for end proc-end 

Code 16-2 Outline of Particle Update Procedure Using a Crystalline Operating System 
column of processors, and that it does not need to move farther left Thus, when next-dir is called and p is in the import queue of processor p 2 in Fig. 16-1, UP is the first direc- tion for which test-export returns 1RUE. So p is placed in the up export queue, and in the next iteration of loop 1, in Code 16-4, it makes its way to its final destination. 

16-3 Potential Energies with a Short-Range Interaction: A More Complicated Example 
In the previous section, we saw how irregular short-range communication could be handled with crOS III. We chose a very simple problem to illustrate the salient points of message forwarding. In this section, a more demanding example is considered. We compute the potential energy of interaction of a collection of particles which interact via a short-range potential, such as the Lennard-Jones potential: 
U I! 12 I Io ~ 7J (16-2) 
This particular interaction is a semiempirical representation of the force between neutral noble gas atoms. The interatomic forces are attractive at large distances due to dipole fluctuations (see for example [Kittel 76]), and are strongly repulsive at short 

U(r)= ~ 

4 

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

295 

proc-begin update yarticles ( update particle data ) declare-buf export ( storage for outgoing particles ) declare-buf import ( storage for incoming particles ) [let RIGHT = 0] [let LEFf = I ] [let DOWN = 2] [let UP = 3] 
[let STAY-HERE = -I] 
for -begIn ( each particle, p ) 
[update particle data for p ] 
proc-caU next-dir ( check if p must move by returning next ) if-begln ( next * STAY -HERE) then 
[remove p from the local list of particles] [place data for p in buffer export[next]] if end for end 
for-begln ( each direction, dir = 0 to 3) 1* LOOP I *1 
[write buffer export[ dir ] to neighboring processor in direction dir , while receiving incoming particles into buffer import[dir]] for -begIn ( each particle, p, in import[dir] ) 
proc-caU next-dir ( find next direction, next, in which p must move ) if-begln ( next = STAY -HERE) then 
[copy p to the local list] 

else 

[copy p to export[next]] if end 
for end for end proc-end 

Code 16-3 Detailed Code for Updating Particle Positioos Using a Crystalline Operating System 
distances due to the "hard core" electrostatic and exchange interactions. The Lennard- Jones potential is short ranged because it falls off sufficiently rapidly with increasing r that to a gOOd approximation 

u (r ) = O for r > 30" (16-3) 
In this example 30" is the cutoff length or the separation beyond which pairwise interactions are considered negligible. It is the existence of such a cutoff that character- izes short-range problems. We can justify the introduction of a cutoff by imagining parti- cles distributed uniformly in space with density p. The average number of particles in a shell of thickness dr, a distance of r from the origin is: 
dN = 41tr2dr p (16-4) and their cumulative affect on the potential energy of a particle at the origin is: 

296 

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

Code 16.4 Procedures to Find the Next Direction to Send a Particle 

Five Steps (a), (b), (c), (d), (e) in TRANSIT of p from 1} to PJ 

(b) ~ 

0 Import Queue 

[i?] Down , , 
[g Left 

Gup 

I Export rQueue 

0 Right J 
Rest of processor is LOCAL LIST 

Figure 16-1 The route taken by a particle, P. which is required to move from processor p I to p 3' a diagonal connection. The steps in the route are: (a) from local list in p I to LEFf export queue in p l' (b) from LEFf export queue in p 1 to the import queue of p 2' (c) fr<Xn the import queue in p 2 to the UP export queue of p 2' (d) from the UP export queue in p 2 to the import queue in p 3' (e) from the import queue of p 3 to the local list in P.. 

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

12 

[7] l dr 

For large r, i.e., r > > 0', we find that 

297 

(16-5) 

;;] -4 d[ ~ 
SO the contribution from shells of increasing radii is a rapidly decreasing function of r /0". If we are interested in the potential energy of a particle at the origin, we introduce a negligible error by restricting our attention to shells smaller than a fixed cutoff, for exam- pIe r < 30". 
Short-range forces appear on a more macroscopic scale in granular flow problems [Werner 85]. In such problems, irregularly shaped objects which are touching interact via friction and elastic deformation. There is no interaction at all between particles which are not touching. The interaction is short-range, with a cutoff given by the diame- ter of the largest particle, but it can be extremely complicated over that range. 
Finally, short-range forces playa role in "Particle-Particle Particle-Mesh" N-body calculations, where they are used as a correction to a long-range force computed by other methods [Hackney 81a]. The Lennard-Jones potential treated in this chapter is intended as an illustrative example. More complicated interactions can be treated in much the same way. 

(16-6) 

16-3.1 The Sequential Problem: Cells and Linked Lists 
Before further considering local particle dynamics problems in parnllel, it is impor- tant to understand the problem in a sequential context. An important characteristic of the algorithm is how well it takes advantage of the locality of the interaction. Clearly, one could solve these problems by the methods discussed in Chap. 9 (Long-Range Interac- tions), but at the expense of a huge amount of wasted computation. The techniques of Chap. 9 explicitly compute the interaction between each pair of particles. Since the number of pairs grows as the square of the number of particles, these methods quickly become impractical as the number of particles increases. Thus, the methods of Chap. 9 are inappropriate for short-range problems. 
One approach might be to abandon the computation of the potential when the separation is found to be greater than cutoff in length. However, computing the separa- tions for each pair is still too much work because the computation of the separation must be performed for each pair of particles. What is needed is a means to quickly identify the neighbors of a particle i.e., those within the cutoff. To be useful, the time spent identify- ing neighbors of a given particle should be proportional to the number of neighbors rather than to the total number of particles. 
A simple way to identify neighbors is to presort the particles into cells of linear size equal to the cutoff length. Then the search for particles near a given particle can be res- tricted to those cells adjacent to the cell in which the given particle is found. This tech- nique, known as geometric hashing, is illustrated in Fig. 16-2. 

298 

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

PARTICLES. IN A TWO DIMENSIONAL DECOMPOSITION INTO CELLS OF SIDE 30- 

FIgure 16.2 To find all particles within 3<1 for particle p, it is necessary to search ooly in the shaded cells. Note that the cells are chosen to have linear size equal to the 3<1. 
We are interested in finding all the neighbors of particle p. The particles of interest are those within the cutoff indicated by the circle. Obviously, it is only necessary to search within the nearby shaded cells. The interaction does not extend into any other cells, and we would be wasting our time by computing the separation between p and par- ticles anywhere else. In Fig. 16-2, the savings realized by limiting our search is not espe- cially great We have managed to avoid searching about half of the cells. Fig. 16-2, however, might represent only a tiny part of a much larger problem such as that shown in Fig. 16-3. Searching only the particles within the shaded area there represents a huge sav- ings ove;r the full O (N2) search. 

16-3.2 Linked Lists 
Now consider the operations we will need to perform on the cells in our simulation. Firstly, and most importantly, given a cell, we must be able to quickly find all the parti- cles contained within it. Secondly, we must be able to fill the cells in a reasonable time. It must be possible to insert a new particle into a cell so that the time required to generate the entire cell structure is proportional to the number of particles. Thirdly, cells must be 

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

299 

..I I I " ., ' , 
'1fI1'r"j'1Ifjj 
c., - " 
, " 't . 
. ~- - 
F1gure 16.3 A larger scale versioo of Fig. 16-2 which shows that by searching only the nine shaded cells, we can immediately eliminate about 99% of the particles from those that could interact with our chosen particle. 
capable of holding different numbers of particles. This suggests that simple arrays can- not be used for the cells because an array cannot grow dynamically as a cell becomes full. If we made all the arrays big enough to hold the maximum number of particles, we would waste a tremendous amount of memory, since most of the arrays would be almost empty .A simple data structure which meets the requirements is a linked list. In a simple implementation, a linked list consists of a collection of items. Each item consists of a pointer to the contents of that item (the particle) as well as a pointer to the next item. Linked lists require a fixed additional amount of memory per particle (the space for the pointers) and, therefore, are not too costly. A detailed discussion of linked lists at this point would take us rather far afield, so the interested reader is referred to any standard treatise on algorithms and data structures, e.g., [Knuth 73c]. 

16-3.3 A Sequential Implementation 
We are now ready to present a sequential algorithm to compute the potential energy of a collection of particles interacting through a short-range potential. This algo- rithm is outlined in Codes 16-5 and 16-6. 
The program in Code 16-5 is quite sirnp1e. First, the array of cells is initialized. Next, the list of particles is scanned and the individual particles are inserted into the appropriate cells. Once the particles are sorted into cells, it is possible to find the neigh- bors of a particle by searching only its own cell and neighboring cells. This is the purpose of the procedure, interact, illustrated in Code 16-6, which computes the interaction between a particle and all particles in neighboring cells. In order to count each pair only 

300 

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

proc-begin compute-energy ( evaluate interaction energy. U ) declare int ncells x ( number of cells in the x direction ) declare = int ncells J ( number of cells in the y direction ) [initialize U = 0] 
for-begin ( i = oto (ncells-x + 1» 
for-begin (j = Oto (ncellsy + 1) ) 
[initialize a linked list for cell number (i ,j )] for end for end 
for-begin ( each particle, p ) 
[find out which cell particle p is in] [insert particle in list for this cell] for end 
for-begin ( i = 1 to ncells-x ) 
for -begin ( j = 1 to ncells y ) 
for -begin ( each particle, P, in the list for this cell ) proc-call interact ( add contribution of interaction 
energy due to p to U ) 

for end for end 
for end proc-end 

Code 16.5 Sequential Evaluation of the Potential Energy Using Linked Lists 
once, it computes the interaction between particles P 1 and P 2 only if the x -coordinate of P 1 is greater than or equal to that of P2. If the x-coordinates are identical, it computes the potential only if the y-coordinate of P 1 is greater than that of P2. This method is com- pletely robust, and any other unambiguous criterion would be equally acceptable. The particular test we have adopted here has the advantageous property that a particle will never interact with itself. Such a self-interaction would have disastrous consequences in the computation of energies which contain terms varying as l/r12. 
We have simplified the code in interact, at the cost of some memory expenditure, by introducing a guard layer of empty cells at the edge of our simulation. In a sequential program, this guard layer is always empty and serves no pUlJ>ose other than to provide a uniform neighborhood for all the internal cells. In the parallel version of the program, however, it will be used to store data received from neighboring processors, so that the neighborhood of a cell on the edge of one processor is consistent with the state of cells in neighboring processors.