12 

Monte Carlo Methods and the Generation of Random Numbers 

12-1 Introduction 
Many important algorithms in engineering and the sciences are of the Monte Carlo type. This means that they employ pseudo-random number generators to simulate physi- cal systems which are inherently probabilistic or statistical in nature [Knuth 73b]. A typ- ical example of such an algorithm is the simulation of the neutrons emitted by a nuclear reactor. In the laboratory , one can measure the properties of the interactions of neutrons with various materials, such as those found in the walls of the reactor. These experimen- tal measurements are used to obtain a set of probabilities, which convey information such as: "If a neutron of velocity v hits material x at angle e, then there is a probability of y that the neutron will be absorbed (damaging the material) and there is a probability den- sity z «<I» of it being scattered by angle «I>." These probability functions can in turn be used to simulate processes in the nuclear reactor and to help in making engineering deci- sions regarding the design of the system. 
This type of Monte Carlo computation proceeds by directly simulating the histories of individual neutrons, from production in the core of the reactor to eventual absorption in some material. For each neutron, the effects on the reactor (damage to walls and so forth) are tabulated. Once a large number of neutrons has been simulated, a reasonably accurate picture of the simulated system emerges (assuming, of course, that the initial modeling of the system was correct). We conventionally call the history of an individual neutron an event or path. We will refer to this type of algorithm as an event-tracing Monte Carlo method. 
In this chapter we discuss a reliable method of generating concurrent random numbers. Chapter 13 will describe two application areas (simulated annealing, lattice gauge theories) which make use of the method. 

12-2 Overview of the Generation of Random Numbers 
At each decision step of the neub"on path problem introduced in Sec. 12-1, random numbers are employed to make the decision in accordance with the probability function representing whatever it is that the neutron has run into. Strictly speaking, the random numbers used in computer algorithms are not really random at all. They are generated by some definite method, and if the algorithm is run a second time, with identical starting conditions, exactly the same "random" numbers will be produced. These generators are therefore called pseudo-random number generators which produce a sequence of 

201 

202 

Monte Carlo Methods and the Generation of Random Numbers Chap. 12 

numbers that appears to be random. One could conceive of attaching some specialized device to a computer to give truly random numbers, such as a Geiger counter measuring radioactive decays of a piece of uranium. Even if we ignore the implications for the reli- ability of our computer, this random number generator would be undesirable for reasons of reproducibility .Suppose we make a modification to the code (speeding up a certain subroutine, for example), and we wish to demonstrate that the algorithm has survived this modification without any bugs introduced. With the pseudo-random number generator, this is easy: we start with some standard initial condition and then put the program through its paces, exercising all relevant parts of the program. If the results are identical, we can be reasonably confident that no bugs have been introduced into the program. 
When Monte Carlo algorithms with their associated random number generators are adapted for a parallel computer, questions immediately arise concerning correlations among the various parallel random number sequences (e.g., the random number sequence of each node of the parallel computer). If the sequences are correlated in some way, several severe problems present themselves. At the very least, if there are correlations, the effective speed of the algorithm may be much less than the naive speedup of the con- current processor. This is because if correlations are present, much of the information coming from the various processing nodes will be redundant and will not help to improve the statistical accuracy of the overall computation. Even more dangerous is the fact that in many types of Monte Carlo computations, correlations will lead to actual biases, and incorrect results will be obtained. In the nuclear reactor algorithm, for example, if more neutrons end up going to the north than the south, due to some correlation of the random numbers, one may falsely conclude that the south side of the plant is correctly or incorrectly engineered. 
It is important to realize that even in the sequential case, the question of correla- tions is a nontrivial one. A given pseudo-random number generator may be such that the nth and the n+18th number are correlated in some way, leading to the problems men- tioned above. In the normal sequential case, ways of generating random numbers which are largely free of correlations have been found. If we could somehow reduce the paral- lel case to something equivalent to the sequential method, we would be satisfied. For the most popular method of generating pseudo-random numbers, this is possible. There is a modified, parallel method which can be used so that the parallel computer exactly mimics what a sequential machine would do in the same computation. All questions regarding correlations among the parallel random number sequences are thereby reduced to those applying to the sequential case which has been extensively studied. In this approach it may be said, "If you were happy with the sequential algorithm and the correlations among the pseudo-random numbers in that case, you will be happy with the parallel algo- rithm." 
Another important feature of the method is that the parallel and sequential algo- rithms can be made to match exactly; therefore, debugging of complex parallel codes is greatly simplified. 
Before plunging into the details of the parallel method and examples, let us men- tion some references. The fundamental formula appears in [Knuth 73b] and ideas related 

Sec. 12-3 Algorithm for Parallel Random Nwnber Generation 

203 

to those discussed here are treated in [Frederickson 83]. The parallel method is explicitly put forth in the context of vector computers in [Brown 83, Barkai 84]. 

12-3 Algorithm for Parallel Random Number Generation 
We employ the most popular type of generator of pseudo-random numbers: the linear congruential algorithm. It is defined by: 
X"+l = ( aX" +c )mod m (12-1) 
where Xi is the i th member of the sequence and a, c, and m are chosen in order to make the members of the sequence appear as random as possible (see Knuth 73b, p. 102, for a table of good choices for a, c and m ). 
A fundamental relation is that, given Eq. Eq. (12-1), it is easy to write down the n +k th member of the sequence in terms of the n th: 
X"+k =(A X,,+C )mod m (12-2) 

with 

A =ak 

c(ak-1\ C = c(1+a+a2+ ...+ak-l) = ~ \.. ~I 
a-1 
Of course, (ak-1)/(a-1) is, typically, difficult to compute in practice (a is usually a very large number). It is usually simpler to use the more direct expression of Eq. (12-3), and to compute C by summing the powers of a directly using Homer's Rule. Equa- tion (12-2) shows that the n+kth member of the sequence is related to the n th through the same kind of relation as Eq. (12-1), with just a different value of a and c . 
Consider a generic concurrent processor, consisting of N processing nodes con- nected by some communication network. It does not matter here if the concurrent pro- cessor is MIMD or SIMD. The idea is to have each node compute random numbers using Eq. (12-2), with k=N (A and C in Eq. (12-2) are calculated only once and stored). The nodes are given a staggered start, so that their sequences don't overlap. Let the seed, or Oth random number, be denoted by y 0. Using this as a seed, a sequential computer would compute the sequence: 

Yo=Yo Y 1 = (aYo+c )mod m Y z=(aY l+C )mod m 

The seeds in the nodes of the concurrent ensemble are set as follows (subscripts denote position in the random number sequence, superscripts denote processing node): 

204 

Monte Carlo Methods and the Generation of Random Numbers Chap. 12 

x~O)=Yo x ~1) = (aYo+c )mod m = Y 1 

(12-5) 

x(N-I) - 
y0 -N-I 
This defines the staggered start. The nodes now use Eq. (12-2) to compute the next member of each of their sequences, with k=N .Therefore 
X 1(0) = y 0t-N (12-6) 
X(I) - 
yI -I+N 

XfN-l) =YN-l+N 

and, 

xiO) =y (}+-2N X ~1) = y 1+2N 

X'f-l) =YN-l+2N 

and so on. The method is illusb'ated in Fig. 12-1 for the case N=4. As is clear from the picture, the idea of the parallel method is to put the processing nodes into a staggered start and then to have them leapfrog one another using Eq. (12-2). The interleaved sequence corresponds exactly to the single sequence of random numbers that would be obtained on a sequential computer. 
The simple pseudocode for a sequential (31-bit) random number generator is given in Code 12-1. The routine rand-set initializes the constants a and c (as given in Knuth, for example) and starts the random number sequence off with seed ran-seed. 
The parallel random number generator, prand , and setup routine, prand-set , are outlined in Code 12-2. The routine prand is of the same form as rand, except that the multiplier and summand that it uses will be chosen to produce the large hops of N down the sequence of random numbers. The setup routine prand-set does most of the hard work of the parallel random number generator. It computes the correct multiplier and summand, A and C of Eq. (12-3), which produce the large hops of N steps. It also sets up the necessary staggered start determined by ran-seed. More specifically, in the 

Sec. 12-4 Example 1: Simple Event Tracing Monte Carlo 

YII YIO Y9 Ye Y7 Y6 Y5 Y4 Y3 Y2 YI YO 

serial computer sequence 

-seed 

205 

staggered start 

FIgure 12-1 
A comparisoo of the sequential and concurrent generatioo of random numbers. 

Code 12.1 The sequential random number generator 
processing node whose procnum is equal to i, R will be set to the i th random number after the seed (Yj of Eq. (12-4». Note that the order of the computation of r and R are reversed in rand and prand. This is intentional and is done to get exact matching of the parallel and sequential sequences. 

12.4 Example I: Simple Event Tracing Monte Carlo 
Our first example will be computing a one-dirnensional integral by the Monte Carlo method. We will evaluate the definite integral: 
1 
I = J f (x) dx (12-8) o 

206 

Monte Carlo Methods and the Generation of Random Numbers Chap. 12 

Code 12-2 The parallel randan number generator 
where in this case, f (x) = X2. Evaluating this integral by Monte Carlo amounts to gen- emting random numbers, ri' in the range O to I, and for each one computing the score, f (ri ). If this is done many times, the average score tends to the value of the integral. More precisely: 

1 L I = Jim -r. f (rj) (12-9) L-- L j=l 
This example seems rather trivial, and the reader may be wondering why we would bother doing a Monte Carlo for a one-dimensional integral. It turns out that in higher dimensions (that is, a multiple integral over many variables) the Monte Carlo method actually becomes an efficient means of evaluating integrals. In statistical mechanics and quantum field theory, where multiple integrals over thousands of variables (at least) are encountered, Monte Carlo is the only efficient method available. The one-dimensional integral also forms a simple example of the event tracing Monte Carlo. The events are just the generation of the random numbers themselves. The score (the summation of Eq. (12-9» is analogous to keeping track of the interaction of each neutron with the reactor walls. The pseudocode which correlates with Eq. (12-9) (for the case f (x) = X2 ) is listed in Code 12-3. 
The parallel code for this problem is simple. We have already resolved the correla- tion problem among the random numbers. In addition, the events are completely independent of one another. The parallel algorithm, therefore, involves having each node generate events and keep its own estimate of the average score. At the end, the N scores are averaged and printed. The individual scores of each event could be accumulated as the algorithm runs, but this would generate much needless communication in the con- current processor, leading to a less efficient implementation. By combining scores only at the end, the nodes are completely decoupled during the actual computation. If each 

Sec. 12-5 Example 2: Calculation of 1t via Monte Carlo 

207 

1 Code 12-3 Code to compute fx2 dx 
o 

node generates many events (the typical case), the algorithm will achieve a speedup very close to the perfect value of N. Code 12-4 contains the pseudocode for the concurrent version of the Monte Carlo integrator. 

Code 12.4 Parallel version of 12-3 

12-5 Example 2: Calculation of 1t via Monte Carlo 
In this example we will compute an estimate for 1t via the "dartboard" algorithm. Consider a circle of unit radius centered at the origin with a square circumscribed about it (the comers of the square are at the points (-1,-1), (1,-1), (1,1), and (-1,1) ). This circle forms the dartboard. We now throw darts at the board, but in a rather special way. They are thrown so that they always hit the square, but have an equal probability of hitting any 

208 

Monte Carlo Methods and the Generation of Random Numbers Chap. 12 

part of the square. During this process, we keep a score which is incremented by 1 for each dart inside the unit circle, and O for darts outside the circle. If we then divide the total score by the total number of darts thrown, we obtain an estimate of the ratio of the area of the circle to the area of the square, since the darts have equal probability of hitting any part of the square. The ratio of the areas is 7t/4, so that the total score divided by the total number of trials (x4) yields an estimate of 1&. 
As in Example 1, the events (the throws of the darts) are independent and so we will have each node of the concurrent processor run its own dartboard simulation and make its own estimate of 1&. At the end, the N estimates will be globally combined. There is, however, a somewhat nontrivial aspect of this example as compared with the previous one. In fact, this is our main motivation for presenting it-there are, after all, better ways of computing 1&! 
The complication arises from the fact that each event depends on two random numbers. As each event is processed in one of the nodes, two consecutive calls to prand will give numbers N steps apart in the fundamental sequence of random numbers, while the sequential program calls rand and so gives two consecutive random numbers, Yi and Yi+l. We have therefore lost the one-for-one matching of events between the parallel and sequential codes. Fixing this is easy: in a simple generalization of the previous prescrip- tion for parallel random numbers, we merely need to generate the numbers (on the con- current processor) in pairs, (Yi.Yi+J and then do a large hop down the sequence of2N-l steps. This will restore the exact ~atch with the sequential computation. 
We present in Code 12-5 the modified random number generator and setup routine. 

Code 12.5 The initia1izer for dte "big jump" random number generator 
We see that prand-setB produces two sets of multipliers and summands, one set- A ,C for the big jumps of 2N -I steps, and the other set, a ,c the usual pair for steps of size lor the small jumps. As a result, we need two different random number routines, prandB and prandS which follow the pseudocode of Code 12-6. With this issue resolved, the rest of the parallel code (Code 12- 7) is simple. The structure is very similar to that of Example I. The nodes execute the same algorithm to generate the events as in the sequential case, with the only difference being to first call prandS , producing the small jump, then prandB , giving the big jump. 

Sec. 12-6 Remarks on the Use of prand 

209 

Code 12.6 The big and small jump random number generators 

proc-begin main ( run parallel "dartboard" estimation of 1t in processor ) [ receive number of darts to be thrown by this processor ] [ receive random number seed, ran-seed ] 
proc-call prand-setB ( find initial random number, R and multipliers and summands for this processor ) 
[ initialize score = 0 ] 
for -begin ( each dart thrown by this processor ) 
proc-ca1l prandS ( evaluate next random number, r , in range O to 1 using small skip ) [ set x-coord = 2r -1 ] 
proc-ca1l prandB ( evaluate next random number, r. in range O to 1 using big skip ) [ set y-coord= 2r -1 ] 
if-begin ( dart has landed inside circle) then [ increment score by 1 ] if end for end 
[ sum score over all processors ] 
[ 1t = 4 * score/ ( total number of darts thrown by all processors)] [ output 1t ] proc-end 

Code 12.7 A parallel program to C<Xnpute 1t 

12.6 Remarks on the Use of prand 
In the examples, we have seen how to employ the parallel set of random numbers to robustly perform Monte Carlo calculations on the concurrent computer. Before con- cluding this chapter, it will be useful to say a few words about more realistic, complex algorithms and the issue of exact matching between the sequential and concurrent pro- gram s . 

210 

Monte Carlo Methods and the Generation of Random Numbers Chap. 12 

In Example 2, we saw that even if the events of the Monte Carlo problem depended upon multiple random numbers, we could modify the concurrent code to provide exact matching. In realistic cases, this strategy can break down. To illustrate, let us go back to the example of neutrons flowing out of a nuclear reactor and interacting with the sur- rounding environment. An event consists of the entire history of a single neutron. One possible way of decomposing this problem onto the concurrent processor is to have each node run its own simulation of the entire reactor and at the very end combine the results. The fundamental problem here is that the number of calls to the random number genera- tor made by each event is inherently unpredictable. The number of calls to prand, in fact, depends upon the values of previous random numbers. A certain value for the first call may send a neutron into some lead, absorbing it and so ending that event. However, a different value for the first call to prand may cause the neutron to miss the lead, and hence have a much longer history, involving many more calls to prand .This results in a problem for the concurrent machine: each node, when it wants to process its next event, and therefore make a big jump in the random number sequence to avoid the random numbers already used by the other nodes, needs to know how large of a jump to make. This would mean that at the end of each event, all the nodes would have to communicate to determine how many random numbers each of them used, to find the proper size of the big jump, and to find the correct values of A and C to produce this jump. 
Besides involving much overhead (recomputing new values for A and C each time) such a solution could give rise to a severe load imbalance on the concurrent machine. Some events will take much longer to process than others. If we resynchronize at every event, we will virtually guarantee this imbalance, but if we let the nodes stay decoupled through the processing of many events, computational load fluctuations will tend to average out, producing a much more load-balanced computation. A much more attractive alternative therefore suggests itself. Instead of modifying the concurrent code in some drastic way, we should slightly modify the sequential code. The philosophy is to have the sequential code simulate the concurrent machine insofar as the random numbers are concerned. 
Suppose we wish to have the sequential algorithm exactly match the same compu- tation on a concurrent machine of N processing nodes. On the concurrent machine, we use the simple parallel set of random numbers as shown in Example 1. This insures that no random numbers will overlap on the concurrent processor. On the sequential machine, we run the usual sequential program, but with a modified set of random numbers. N independent streams of random numbers are set up ( e.g., turn the R of Code 12-1 into an array of length N) and the multiplier and summand for all N streams are the values A and C given in Eq. (12-3). These streams are given the same staggered start as in the concurrent case. Then, in the simulation of the first event, the first stream is used, in the simulation of the second event the second stream is used, and so on. This represents a fairly simple change in many sequential codes, and causes it to exactly match the concurrent result 
Finally, we should say that this is not the end of the story .Linear congruential pseudo-random numbers are often not random enough, and one may wish to use more 

Sec. 12-6 Remarks on the Use of prand 

211 

powerful algorithms such as linear shift register methods. A parallel version of these seems problematical and to our knowledge has not been done. Recent research in paral- leI random numbers appears in [percus 87] . 

13 

13-1 Introduction 
A reliable technique for generating pseudo-random numbers on a parallel machine was described in the previous chapter. Here we will describe some applications of this method. Two examples will be discussed: the Traveling Salesman optimization problem and the use of Monte Carlo methods for studying quantum field theories. 
The Traveling Salesman application may seem at first to be an esoteric toy prob- lem, but it serves to introduce the important subject of simulated annealing whose use in concurrent decomposition and load balancing of complex problems will be discussed in Chap. 17. Quantum field theory simulations are important in physics, and we will also find that they lead to an interesting approach to general finite difference algorithms on a MIMD concurrent computer. 
The applications in this chapter will be described within the restricted context of hypercube architectures, because it is with this class of machines that the authors have experience. This specialization should not, however, significantly limit the applicability of the algorithms of this chapter to other parallel machine architectures. 

13-2 The Traveling Salesman Problem 
The classic Traveling Salesman Problem (TSP) is that of finding the shortest tour for a traveling salesman so that he visits each of a set of cities exactly once, and then returns to his initial city. The difficulty of this problem is illustrated by the immense number of possible tours: (N-1)!/2 for N cities. As this number grows faster than any finite power of N , one is limited to small problems if applying a brute force search of all candidates. 
The TSP is a simple example of a difficult optimization problem and is NP- complete. Examples with more practical importance include circuit layout and wire placement, resource allocation, and logistical problems. In this chapter we will discuss the simulated annealing approach to the solution of these problems. Graph theoretic and linear programming techniques are covered (for sequential computation) in [Luenberger 84] and [papadimitriou 82]. Important early work on the TSP appears in [Lin 65, 73], the branch and bound approach is covered in [Held 70, 71] and a parallel version of branch and bound is given in [Felten 87]. A promising approach to optimization, that of neural nets, is first given in [Hopfield 86]. 

213 

214 

Applications of Monte Carlo Methods: Traveling Salesman and QCD Chap. 13 

The simulated annealing method [Kirkpatrick 83] is a Monte Carlo approach to solving the TSP. The algorithm samples a small number of the (N-l)!/2 possible tours to find a good solution. The method does not in general find the optimal solution, but rather an approximation to the best solution. In many practical contexts this near-optimal solu- tion is quite adequate, so that simulated annealing represents a viable approach. 
As mentioned above, we will discuss the TSP in terms of its implementation on a hypercube concurrent processor. This is to provide some specificity to our discussion and also to reflect work that has actually been done [Felten 85b]. As will be seen, the technique generalizes to MIMD machines of many other possible topologies. The hyper- cube is, as we've seen, well suited for mesh problems-those of a regular geometric nature. The TSP does not fall into this class of problems and this provides further motivation for studying it. It is of interest to learn how effectively concurrent processors of regular architecture can be applied to such a seemingly irregular class of problems. 
The cities of the TSP are assumed to be arbitrarily positioned on a two-dimensional plane. An attempted physical or equal area decomposition of the problem would likely lead to load imbalance, with some processors of the ensemble having significantly more work to do than others. This load imbalance would be a major source of inefficiency. The correct decomposition is found by realizing that the topology of the TSP matches that of a ring. Since the traveling salesman visits each city exactly once before returning to the initial city, there exists a ringlike connectivity between cities in a given tour. Hypercubes include rings (see Fig. 13-1) as topological subsets, so that the problem maps nicely onto such a machine. The same is true of architectures based upon periodic meshes. By assigning each processing node an equal number of cities, we can minimize load imbalance. Each processor is assumed to store a list of cities. We begin by assum- ing there are 2d processing nodes and N cities, and that the tour starts with the first city on the list in node 0. The candidate tour follows in sequence the list of N /2d cities in that node, and then continues with the list in node 1. After going through each of the nodes the tour returns to the first city in node 0. This is a decomposition based upon city label- ing. As the optimization algorithm runs, the cities will be shuffled between nodes of the concurrent processor, eventually reaching their proper location in the optimal tour. 

FIgure 13.1 A three-dirnensional cube with the ring emphasized. 

Sec. 13-3 Simulated Annealing as an Optimization Technique 

215 

13-3 Simulated Annealing as an Optimization Technique 
A possible approach to solving dIe TSP is dIat of iterative improvement. This medIod starts widI some initial configuration (initial tour) of dIe system and dIen changes are made to dIe system. The simplest kind of basic change dIat can be applied to perturb a tour is to swap pairs of cities. These swaps are not necessarily limited to adjacent pairs. Changes dIat decrease dIe tour lengdI are accepted, while dIose dIat increase it are rejected as shown in Fig. 13-2. The process is repeated in some sequence until no further changes are accepted. 

(a) 

B 

N 

/ 

~ 

~ 

(b) 

FIgure 13.2 Adjacent swap algorithm. (a) initial tour (b) tour obtained by swapping B and C in the tour order. 
Let us call the length of the tour the objective function. This function is defined over the space of all possible tours. It goes by this name since the objective of the TSP is to find the location of the minimum of this function. At this point, it is easy to see why iterative improvement is not a viable method for solving the TSP. Changes that decrease the objective function are always accepted, while those which increase the objective function are always rejected. For almost all initial tours, therefore, iterative improvement will rapidly get trapped in a local minimum, far above the true global minimum. To avoid getting caught in local minima, the algorithm should be allowed to conditionally accept some increases in the objective function. One achieves this by using the simulated annealing algorithm, which employs an analogy with statistical mechanics. 
As with iterative improvement, we start with some initial tour and make changes to it As before, changes which lower the objective function are always accepted, but instead of rejecting all moves which raise the objective function, we accept them with conditional probability e-AgIT , where L1g is the amount of increase in the objective func- tion (tour length), and T is a parameter that is subject to control. In the statistical mechanics analogy the objective function is the energy of the system and the parameter T is the temperature. By occasionally accepting increases in the objective function the system can climb out of local minima. In the physical analogy, iterative improvement corresponds to the algorithm with T=O and freezes the system very rapidly. This pro- duces a metastable configuration usually far above the ground state (optimal tour) of the 

216 

Applications of Monte Carlo Methods: Traveling Salesman and QCD Chap. 13 

system. By lowering the temperature slowly, one anneals the system-the system is given time to escape from metastable configurations. 
One problem with this approach is in developing adequate annealing schedules. An annealing schedule is the prescription for the number of iterations required at each temperature as the temperature is gradually lowered. Typically, the system is susceptible to freezing at a certain temperature (or temperatures). As a result, the temperature must be lowered very slowly through these points so as to avoid trapping in metastable states. 

13-4 Concurrent Implementation of the Traveling Salesman Problem 
We will now discuss the parallel implementation of simulated annealing for the TSP. A crucial aspect of the algorithm is the method of inb"oducing random changes into the tour. Initially, we might consider the use of simple swaps of adjacent cities as our move sb"ategy .Using adjacent swaps as the move method is reasonable at low tempera- tures, but not at high temperatures. This is because at high temperatures the goal is to find the correct globallarge-scale structure of the optimal tour. It takes many adjacent swaps to generate large-scale structure changes in the tour, so this does not work efficiently. We will ignore this point for now, but return to it later. 
Using adjacent swaps has the advantage that Ag can be computed locally. That is, for most of the cities in a given node of the ensemble, the information necessary to com- pute Ag is already resident on that node. If a node has cities A , B , C and D and we are considering swapping B and C (see Fig. 13-2) then: 
Ag =(AC +BD)-(AB +CD) (13-1) 
where AC is the distance between A and C .Each node stores a part of the tour. A cycle of the algorithm begins by having each node perform adjacent swaps of the cities away from the boundaries of the tour pieces. Once these are completed, the nodes communi- cate in a ring in order to complete adjacent swaps of the cities on ends of the tour seg- ments. This sequence comprises one iteration. The process is repeated, allowing cities to move throughout the concurrent processor, eventually finding their proper locations in the optimal tour . 
This algorithm is very efficient, but only when compared with the same (adjacent swap) algorithm on a sequential machine (see Table 13-1). The main difficulty, as men- tioned above, is that this algorithm is inefficient in making large-scale improvements. When the system reaches a local minimum, it sometimes requires a large conspiracy of simple adjacent swaps to climb back out of the minimum. In order to speed up this pro- cess, we will at this point exploit the richer set of communication channels available in the hypercube architecture. 
The improved algorithm is the same as before, except that it uses a more complex method of moving cities around the tour in addition to the adjacent swap. By using the exb"a communication channels we can move a city a longer distance through the tour dur- ing a single step. In Fig. 13-3, a 25 cube (32 processors) is shown in such a manner as to make the embedded ring obvious. In the adjacent swap algorithm, only those communi- cation channels along the circumference are used. Clearly, jf a city needs to move from 

Sec. 13-4 Concurrent Implementation of the Traveling Salesman Problem 

217 

node 16 to node 24, it would be much faster to move it directly. Table 13-1 

Figure 13-3 Even though it is not obvious, this is a 5 cube drawn in such a way so as to show the embedded ring. The small circles represent processors; the lines are communi- cation channels. The channels going clear across the ring are those used during "hy- perswaps" and allow cities to make very large moves in the tour. 
To describe this nonadjacent swap between two nodes. let us designate the involved nodes by x and y. Assume that node x has cities labeled Ax. Bx. and Cx and node y has cities labeled Ay. By. and Cy. If we propose to swap Bx with By (see Fig. 13-4) then 

L\g = 

AJi; + B;c; + A;ii: + B:c;J A;ii; + ii;c: + A;B; + B;c; ] 

(13-2) 

218 

Applications of Monte Carlo Methods: Traveling Salesman and QCD Chap. 13 

A reasonable implementation of the required calculation involves node x sending Ax, Bx and Cx to node y. Node y then computes L\g, decides whether or not to accept the swap, and so sends back to node x either Bx or By. This "hyperswap" algorithm gives the k'h city in each node the option of swapping with the k'h city in any of that node's "hypemeighbors" within the hypercube. 

" 

A~ 

" 

"' 

), 

By.( 

, 

/ 

¥ 

>B% 

/ 

" 

~y 

/ 

~ 

\ 

\{ 

(a) 

(b) 

FIgure 13.4 Hyperswap algorithm. (a) initial tour (b) tour obtained by swapping B. and By in the tour order. 

The hyperswap algorithm is less "efficient" than the pure adjacent swap method because every city is on the edge-that is, requires communications. On the other hand, the hyperswap is much better from a standpoint of obtaining the solution quickly, because it is less sensitive to the annealing schedule chosen and is several times more effective than using only adjacent swaps. The best algorithm may be a mixture of adja- cent and hyperswilps, with fewer hyperswaps at low temperatures, since most of the large-scale structure of the system has locked in. 
There is nothing special about the hypercube topology of the hyperswap algorithm. The aim is simply to allow cities to make large moves within the tour in single steps. On a parallel machine of lower connectivity than a hypercube, the large steps could be implemented via message forwarding. Conversely, on a shared-memory machine, large- step swaps could be implemented with even more freedom than in the hypercube exam- ple. Other types of moves which can be more effective than long-range swaps are dis- cussed in [Kirkpatrick 83] . 
Because the optimal sequential annealing algorithm is not known, we cannot mean- ingfully discuss the efficiency of the concurrent program as we have done for other prob- lems. In using (5-4}, one should in fairness always use the best sequential program as a baseline for efficiency calculation, rather than a degenerate sequential version of the parallel algorithm. 

Sec. 13-6 Quantum Field Theories 

219 

13-5 An Example Traveling Salesman Problem 
A 64-city problem, organized as 4 groups of 16 cities each positioned on a unit square, is shown in Fig. 13-5. The illustrated solution was found using a 64-node Caltech Malic I hypercube. The cities were grouped in order to observe the hierarchical nature of the solution process. After running at high temperature T, the four links connecting the groups were discovered by the optimizer and frozen in. The small-scale structure within the groups was then discovered by annealing at lower temperature. Since there still remain a few crossings within the groups, we surmise that we have not found the abso- lutely optimal path-further running at low temperature would be necessary to remove these defects. 

r 

r' 

Yl 

\;1 

~ 

Figure 13.5 An approximate solution to a 64 city problem, with the 64 cities grouped as 4 clusters of 16 each. 

13-6 Quantum Field Theories 
Quantum chromodynamics (QCD) is the theory of quarks and gluons and how they interact. One would (well, at least a theoretical physicist would!) like to solve this theory and derive its predictions for such things as the force law between quarks. Furthermore, it is of interest to compute the theoretical properties of bound states of quarks and glue for comparison with the real world-protons, neutrons, and the like. If one can actually match all the observed properties in a simple, elegant way, then this will provide some indication that QCD is, in some sense, the correct theory . 
QCD seems to be impossible to solve by analytic means. Not only is the relevant physics inherently quantum mechanical (so that semiclassical approximations are use- less), but the fundamental interactions of the theory are nonlinear. The approach to the solution of QCD, and other quantum field theories by numerical means [Creutz 83] has 

220 

Applications of Monte Carlo Methods: Traveling Salesman and QCD Chap. 13 

only been realized in recent years. This approach (lattice gauge theories) is actually a numerical realization of Feynman ' s path integral formulation of quantum mechanics [Feynman 65]. 

13- 7 The Monte Carlo Approach to Lattice Gauge Theory 
The reader interested in a detailed description of the Monte Carlo technique in lat- tice gauge theories should consult the references; here we will simply sketch what goes on and point out those algorithmic features of interest in this discussion. We begin by setting up a four-dimensional grid or lattice. This lattice represents a volume of space- time. Associated with the sites of the lattice are variables which represent the quarks. On the connecting bonds or links between the sites reside other variables which represent the glue fields. In the case of QCD, these terms are actually SU(3) matrices. 
The Monte Carlo algorithm proceeds by cycling through all the links and randomly changing their values as well as those of the nearby link variables (the interaction is local). Once the link variables have settled into a physically correct pattern of random walks, various statistical correlations between the link variables are evaluated, which are then used to obtain the desired physical properties of QCD. As an example, the force between quarks is related to the average value of large loops of link matrices. B y aver - age value is meant an average over many cycles of the Monte Carlo algorithm; a loop is the product of link matrices along some closed path through the lattice. 
In practice, one is not really interested in the properties of QCD in a finite box, so the spacetime volume that the lattice represents should be made as large as possible. This requirement, coupled with the four-dimensional nature of the lattice and the large number of degrees of freedom at each site and link, make for large memory and compu- tational demands on any computer. Lattice gauge Monte Carlo simulations often contain 106-107 degrees of freedom. As a result of the Monte Carlo statistical approach, the answers have statistical uncertainties. Interesting observables, such as the closed loops mentioned above, tend to have very unfavorable signal-to-noise ratios. The very long runs needed to resolve vanishingly small correlations further increase the demands upon the computer. It would appear that even current supercomputers are not up to this task in general. As a result, lattice gauge theories have provided strong motivation for the development of extremely high-performance, special-purpose parallel supercomputers [Christ 86, Beetem 86, Marinari 86, Gaines 86, 87, Schilling 87] . 
Figure 13-6 shows the result of an early lattice gauge computation on a hypercube computer [Otto 84a, Brooks 84]. The plotted result is the mutual potential energy of a pair of quarks, as a function of their separation distance. This is the QCD analog of the Coulomb force law between electrons. The error bars represent statistical errors due to the Monte Carlo method. The computed force law turns out to yield agreement with the expected qualitative behavior of QCD. At short distances, we observe Coulomb-like behavior (potential -l/r ), which accords well with the limit of asymptotic freedom. At long distances, a different term dominates the potential, which grows linearly with dis- tance, r. This again agrees with previous expectations of the behavior of QCD at large distances. The linearly increasing potential corresponds to the formation of a flux tube 

Sec. 13- 7 The Monte Carlo Approach to Lattice Gauge Theory 

between the quarks, leading to quark confinement. 

221 

,0 

2.0 

0.0 

v~ 

.2.0 

-4.0 

I 

Heavy Quark Potential 

--n 

/ 

~ 

- 

- 

~ 

'r 

-6.0 I. '; , 0.0 0.4 

0.8 x 

1.2 

1.6 

Figure 13.6 The heavy quark potential result as described in [Otto 84a]. The ordinate la- bel V~ refers to the potential energy of the quark, anti-quark system: the force law between quarks is given by the negative derivative of the potential. The abscissa label x refers to the distance between the quarks. In absolute units, x = 0.4 roughly corresponds 
10-14 to cm. 
The above described computation was perfonned on a 123 x 16 lattice mapped onto a six-dimensional hypercube (64 nodes). The hypercube was mapped onto a 4 x 4 x 4 three-dimensional mesh. As a result, each node was responsible for the update or evolu- tion of a 33 x 16 sublattice of the entire problem. The concurrent efficiency of the com- putation was measured as 95%. The 5% inefficiency is accounted for by a 2% loss due to load imbalance (not all processors doing the same amount of work) and a 3% loss due to communication overhead. It is significant that this low communication overhead is not due to the short range of interaction between the degrees of freedom, as shown in the traditional lattice problems of Chaps. 7 and 8. In fact. most of the variables in the 33 x 16 sublattice are communicated at least once during an update cycle. The overhead is low because of the proportionally greater amount of computation required for each variable. 
We can quantify these remarks by noting that the communication component of the concurrent overhead for this problem on the Caltech/JPL Mark I, n hypercubes follows the relation: 

f = l/e- 1 = 

0.25 n 1/3 

(13-3) 

where n=33 X 16 is me grain size. The coefficient of n-1/3 is somewhat lower man the 

222 

Applications of Monte Carlo Methods: Traveling Salesman and QCD Chap. 13 

value of unity found in Chaps. 7 and 8 for iteratively solved elliptic problems. The lower coefficient reflects the increased calculational complexity of the QCD algorithm com- pared with the finite-difference and finite-element approximations to the Laplacian. In order to further elaborate this point, let us consider a QCD calculation for SU(m) gauge theory .The calculation involves multiplication of m x m matrices and has an associated computational cost proportional to m 3; the communication cost is proportional to the number of matrix elements or m2. Thus, we find that: 
1 1 tcomm 
(3-4)l/e- 1 °c ---1 mnl/3 t 
calc 
and decreases as m increases. The measured value in Eq. (13-3) corresponds to the case m=3. 

An important lesson from the issue of calculational complexity is that one need not be limited to short-ranged interactions; more nonlocal versions (which do occur in prac- tice) are also efficient. This latter point has been experimentally verified-see [Otto 84b]. 
An important algorithmic constraint on these Monte Carlo calculations is the so- called detailed balance condition. This condition places an important restriction on the update cycle of the concurrent algorithm. In order to satisfy detailed balance, any two points that are "connected" through the fundamental interaction (Lagrangian) must not be updated simultaneously [Otto 84b]. The basic decomposition involves division of the lat- tice into equal subregions which are assigned to individual nodes. Concurrency is real- ized by performing parallel Monte Carlo changes of the field within each processor so that N points are being updated at anyone time. During this concurrent update, care must be taken to ensure that none of these points are linked in the interaction. This is easily accomplished for QCD by requiring that any two updated points are separated by at least one lattice point that is unchanged at the current step. There arises an analogous constraint in condensed matter Monte Carlo simulations and these require a somewhat more subtle algorithm [Johnson 86] to ensure detailed balance. The simulation of a liquid with irregularly distributed particles affords no obvious method of assuring that adjacent processors do not simultaneously try to move two mutually interacting molecules. In Monte Carlo calculations involving long-range force fields, detailed bal- ance precludes the update of any two particles simultaneously. In such a case, con- currency must be sought not in simultaneous updates but through decomposition of the calculations associated with a single update. This approach has been applied to Coulomb force fields, and has proven quite successful [Fucito 85]. 

13.8 The "Walking Program" 
Before concluding this chapter, it will be useful to discuss an important method of structuring the QCD algorithm in order to simplify the application code. As is the case in much practical programming, there is a right way and a wrong way to write the code. 
The issue of concern relates to the previously mentioned loops. The fundamental variables of the calculation are matrices and a loop is a product of these matrices along some closed path through the lattice. On a distributed memory machine such as the 

Sec. 13-8 The "Walking Program' 

223 

hypercube, the encoding of these can become quite difficult since the loop can overlap more than one processor subdomain. Explicitly keeping ttack of where a given loop crosses processor boundaries is a daunting task, especially in four dimensions! 
An effective way to deal with this problem is to devise a walking program for the calculation of loops. Figure 13- 7 shows how to describe the shapes of arbitrary loops. A loop is specified by a string of symbols, each of which gives the direction for a single step along the loop. The loop shown is specified by the string: (+.i, +j, +i, +j, -.i , -.i , -j, -j). To compute the loop, we feed the walker routine a starting location and this string as input. The walker reads the first entry of the string, takes a step through the lattice in that direction, and sets a matrix, answer, equal to the matrix it finds there. It then reads the next string entry, steps in that direction, and multi- plies the matrix it finds there by answer. This process continues until the end of the string is reached. At this point, the walker will have walked completely around the loop and answer will finally contain the desired result The walker travels around the loop, carrying the partial answer along with it 

y 

x 

.x 

-x 

-5> 

-5> 

" +x 

+51 

f Start 

+x 

+5> 

F1gure 13.7 The loop corresponding to the walker string (+x. +.9. +i, +.9. -x. -x. -.9. -.9). 
The beauty of this technique is twofold. Firstly, given a general walker routine, loops of all sizes and shapes are handled with equal ease. All that is required is to con- struct the appropriate string for the loop. Secondly, this algorithm translates quite easily into a parallel environment. All that needs to be added is a check, at each step, to deter- mine whether a processor boundary is being crossed. If so, the partial result and the remaining string are passed to the appropriate processor. Once a processor passes along its work in this manner, it checks to see if any work requests (consisting of partial answers and strings) have entered its domain. If so, it processes them. Since the loops close, we are insured that once the dust settles, the loop answers will reside in their correct processors. That is, a final answer will land in the processor which contains the 

224 

Applications of Monte Carlo Methods: Traveling Salesman and QCD Chap. 13 

start of the loop. 
The parallel walker routine allows the computation of loops of arbitrary shape and size to proceed, simply by specifying loop strings. Since the entire QCD algorithm is based on loops, the routine allows the high-level coding of the algorithm to be tran- sparent to issues involving processor boundaries. 

13-9 Finite Differences and the Walker Method 

It turns out that the walker methodology helps tremendously in writing the QCD code. Can similar ideas be applied elsewhere? The answer seems to be yes. An exam- pIe is found in finite difference codes, typified by the Laplace relaxation algorithm, described in Chap. 7. In that case, the fundamental step replaces a field variable with the average of its neighbors: 

[ +.i , fetch, +i, -.i , add, -.i , -i, add, -i, +.i , add, +yJ 
The entries "fetch," "add" indicate the action to be taken on the field variable(s) at the current site. Note that, in this case, the walker simply passes over some of the sites without using the variable residing there. 
The use of this approach would allow coding of the Laplace solver without explicit reference to boundary points. Admittedly, this may represent overkill for the simple finite difference Laplacian problem, and entails added overhead associated with process- ing of the walker string. The potential importance of the method, however, becomes clear when one realizes that it allows modification of the differential operator and/or the dif- ferencing scheme (to arbitrarily complex templates) during execution of the program. For more complex operators and differencing schemes, the overhead of interpreting the string becomes a progressively less important consideration. 
The programming methodology described in this section has quite general applica- bility. By properly designing the parallel implementation we have hidden the processor subdomain boundaries from view (allowing the solver to be easily coded). An added benefit is the enhanced adaptability and potential functionality of the overall program. 

Sec. 13-9 Finite Differences and the Walker Method 

. 

225 

. 

. 

. 

. 

(a) 

(b) (c) 
Figure 13.8 (a) The template for the simplest Laplace solver. (b) and (c) are equivalent walker strings for this template. 

14 

14-1 IntrOOuction 
In Chap. 4 we presented a virtual concurrent processor so that we could discuss concurrent algorithms in concrete terms. The virtual concurrent processor was made very general to avoid restricting the discussion of the algorithms. However, the virtual concurrent processor was not designed to be implemented efficiently as a real machine, so the implementation issues were intentionally ignored. In this chapter we resume the discussion of communications in a concurrent processor and examine in detail a specific implementation of a crystalline communication system on a practical machine. 
We shall give implementations in the C programming language of many of the rou- tines discussed in Chap. 4. The C routines are meant to illustrate the algorithms used to implement the collective routing of data through the concurrent processor. However, the routines as presented are not necessarily the optimal or fullest possible implementations on a real machine. Some features, such as the handling of errors, are treated superficially since a more detailed treatment would detract from the discussion of the algorithms without providing any important benefits. In particular, no effort is made to correctly set the external integers e"no and errcnt that are introduced in Sec. 14-3. 
The crystalline operating system that we present in this chapter is called crOS III. It was designed at Caltech to run efficiently on a specific concurrent processor -the Mark III hypercube described in Chap. 2 -and so will have some features and routines that were not discussed in Chap. 4. In addition to the discussion in this chapter, Appendix E provides a concise definition of crOS III in the form of manual pages for each of the rou- tines. We begin this chapter with a detailed description of a hardware implementation of the communication channels. Following the discussion of the channels, we discuss in detail the implementation of the basic crystalline communication routines. The imple- mentation of the basic routines will not refer to any particular interconnection topology. Although crOS III was motivated by a particular design of hypercube, it has been designed to be easily implemented on other machines. Currently it runs on essentially all hypercubes and the shared-memory BBN Butterfly. 
Before we can proceed to the implementations of the collective communication routines, we must select an interconnection topology for the ensemble. The fully inter- connected topology is impractical for a real machine and is not needed for the applica- tions we have examined in previous chapters. Mter comparing the features of several topologies, we conclude that the hypercube topology is both adequate for all the applica- tions we have examined and practical to build. In addition, we assume that the control 
227 

228 

ClOS m: A Specific Implementation of the Progranuning Environment Chap. 14 

process, whose location has been left unspecified until now, will reside in a separate pro- cessing element that will be known as the control processor. Having specified the topol- ogy encompassing the nodes of the ensemble and the control processor, we proceed to the discussion of routines that implement collective communications on a hypercube. Emphasis is placed on the algorithms used by the routines rather than other important, but less interesting, details of the implementations; 
We then return to the issue raised in Sec. 4-6 of assigning nodes in the ensemble to regions of the decomposition topology. The issue is much more important now, since the hypercube does not have channels connecting all pairs of nodes and we would like adja- cent regions of the decomposition topology to be connected by communication channels. We describe an algorithm for assigning nodes to a restricted class of Cartesian grids. 
The chapter concludes with some comments on the limitations of the crystalline communication system and introduces a new method of implementing communications, the interrupt-driven approach, which overcomes some of the limitations of the crystalline system. Since none of the applications we examine in this book require an interrupt- driven communication system, a more detailed treatment of the subject is not given here. Finally, a few comments are made on the application of the topics discussed in this chapter to concurrent processors operating in a shared-memory rather than a message- passing environment. 

14-2 Implementation of Communication Channels 
Before we can discuss the implementation of the crystalline communication sys- tem, we must examine some of the details of the communication channels that will be used by the system. The channels provide the only means of communication between the nodes of the ensemble. Only two nodes can be connected by a single channel, so all communications are point-to-point. In contrast, if many nodes resided on a common bus, broadcast communications could be a primitive function. The channels that we consider are implemented in the node hardware and consist of a transmitting component and a receiving component. The transmitting component of one node is connected to the receiving component of the other node and vice versa. 
While the transmitting and receiving components of a channel can be implemented in many ways, we will illustrate the important features of a channel by considering one particular implementation based on that used in the Caltech!JpL Mark I, II, and III hyper- cubes. In this implementation, each receiving component of a channel is associated with a buffer. The transmitting node can write into the channel buffer and the receiving node can read from it However, the transmitting node must not write into the buffer if it con- tains data that has not yet been read. Likewise, the receiving node should not attempt to read from a buffer that does not contain valid data. Thus, the channel hardware must provide information on the status of the channel buffer to both the transmitting and receiving nodes. 
In the hardware implementation that we are considering, two interrupt lines provide the channel's status information, where one of the lines is connected to the transmitting node and the other to the receiving node. A channel's interrupt lines can either be polled 

Sec. 14-2 Implementation of Communication Channels 

229 

by the processor or cause the processor to be interrupted; in the current implementation we will always use the lines in polled mode. The channel hardware sets the transmit interrupt line when its associated channel buffer is empty, and clears the line otherwise. Similarly, the channel hardware sets the receive interrupt line when its channel buffer is full, and clears the line otherwise. In this implementation of the commlli'lication chan- nels, all data must be sent in fixed-Iength packets whose length is determined by the size of the channel buffer. The data must be sent in packets because the receive interrupt can- not indicate that the buffer is partially full and the transmit interrupt cannot indicate that the buffer is partially empty. 
As mentioned previously, a node must have a transmitting and a receiving com- ponent for each channel, so that it can both send and receive messages. If the transmit- ting component and the receiving component of the channel are completely separate, each node connected by a channel can write into the other's channel buffer without interference. Such a channel is known as fu/1-dup/ex since both sides can transmit independently. However, two sets of wires are needed to connect a full-duplex channel since the two sides of the channel are independent. Fewer wires are needed to connect a channel if a half-dup/ex circuit is used. In a half-duplex circuit the same wires are used by both sides of the channel, so only one of the nodes can be transmitting at a time. In the context of the crystalline communication system, using a half-duplex channel does not restrict communication, since in the crystalline system it is required that a write on one side of the channel be complemented by a read on the other. 
In general, each node has several such communication channels, each consisting of a channel buffer and receive and transmit interrupt lines. The processor must have a means of addressing which of the channels it wants to access, and specifying whether it wants to transmit or receive. When a packet is read from a channel, the channel is first selected for reading. When the processor attempts to read the packet, only the selected channel sends data to the processor. Likewise, when a packet is to be written to a chan- nel, the channel is first selected for transmitting, so the packet is transmitted by only the selected channel. 
Clearly, having only one channel selected for reading is essential because other- wise several channels could try to send different packets to the processor. However, allowing several channels to be selected for transmitting at the same time is a useful feature. No interference can occur as it could if several channels were selected for read- ing, since the identical packet is being transmitted simultaneously on all selected chan- nels. When we discussed matrix algorithms in Chap. 10, we saw an instance in which being able to write the same packet to several channels at once is desirable. In Sec. 14-5, we implement a broadcast communication routine similar to that introduced in Chap. 4, which writes the same packet to several channels at once. One additional generalization that will prove useful is to allow a packet to be transmitted on selected channels as it is being read from a different channel. We simply note here that allowing simultaneous transmits and receives does not substantially increase the complexity of the nodes' hardware. 

230 

crOS m: A Specific Implementation of the Progranuning Environment Chap. 14 

While having the ability to perform simultaneous b"ansmits and receives on several channels is a desirable feature, not all concurrent processors have this ability. Even if the processors only allow single channels to b"ansmit or receive at once, the discussion in the following sections will still be useful. However, the implementations of the basic com- munication routines will be slightly more complicated, since they will have to loop over the indicated channels to simulate the simultaneous channel operations. The basic com- mu.lication routines will then have the same functionality as those implemented with channels that allow simultaneous transmits and receives. The only difference to keep in mind is that the broadcast algorithm presented in Sec. 14-5 will not always be optimal. 
An important issue in the design of the channel that we have overlooked is the size of the channel buffer. Since only full packets can be sent by the channel, if the packet size is much larger than a message, time will be wasted padding the buffer and sending more data on the channel than necessary. However, if the buffer is too small, the over- head of polling the channel is proportionally greater, and the maximum data b"ansfer rate of the channel is reduced. 
It has been found in practice that the software overhead associated with polling the channel and transferring data from memory to the channel is substantially greater than the time the channel hardware takes to transfer the data across the channel. Thus, send- ing a moderate anlount of unnecessary data on the channel is not very important com- pared with performing extra software operations. If data can be b"ansferred into and out of the channel buffer through a direct memory access (DMA) circuit, the major cost becomes that of testing the status of the channel and initializing the DMA conb"oller. In such a system, the packet size can be quite large before it becomes a dominating factor. The Mark III hypercube at Caltech, which has a DMA conb"oller connecting the node memory and channel hardware, uses a channel buffer of 64 bytes. Our present estimates indicate that the buffers could be considerably larger without substantially slowing the transmission of short messages. 

14-3 Basic Crystalline Routines 
Having defined a specific hardware implementation of the communication chan- nels' we are ready to proceed with the implementation of the CrOS III crystalline com- munication system. The implementation that we consider will reflect many of the features of the virtual concurrent processor presented in Chap. 4. However, several new features will be added, since the system we describe is meant to be a real system that gives good performance on a real concurrent processor. In addition, we will pay atten- tion to details that were ignored in Chap. 4, such as the detection of communication errors. 

Another difference from the discussion of the virtual concurrent processor involves restrictions on the interconnection topology. In the case of the virtual concurrent proces- sor the communication channels are numbered by the processor number of the node on the other side of the channel. Thus, communication is said to be node-addressed. In our crOS III crystalline communication system the basic routines cread and cwrite imple- ment communication only between nodes directly connected by a communication 

Sec. 14-3 Basic Crystalline Routines 

231 

channel. Since we are not assuming a fully interconnected topology , node-addressed communication is no longer desirable. Instead, we simply number the channels with con- secutive numbers starting at 0. This type of communication system is said to be channel-addressed. 
One of the features we add here to cwrite is the ability to send a message to several channels simultaneously, instead of only one channel. The channels on which the mes- sage is to be sent can be efficiently specified by a bit field. A bit field is a data object whose individual binary digits represent quantities that can assume only two states. Thus, a bit field consisting of one bit for each channel can specify the channels on which a mes- sage is to be sent The bits corresponding to the desired channels are set while the rest are cleared. The bit field specifying the desired channels is known as a channel mask, which replaces the channel number as an argument to cwrite in the current implementa- tion. Thus, cwrite is defined with the calling syntax given in Code 14-1. 

Code 14.1 The Routine cwrile 
The channel mask, chan, is an unsigned integer since a bit field has no sign associ- ated with it We define the following special cases for the arguments of cwrite. If the channel mask is 0, indicating the message is to be sent on no channels, the call to cwrite is ignored and it should not be complemented by a cread on the other side of the channel. However, if the channel mask is nonzero and the byte count is 0, the call to cwrite must still be complemented by a call to cread since a message whose length is ° is sent to the destination node. Sending messages of length ° is sometimes useful if some nodes have nothing to write, since all nodes typically read and write together when crystalline com- munications are used. 
The data sent by cwrite consist of a header that specifies the length of the message, the message itself, and a trailer that consists of the checksum of the header and the data. A header is necessary so that the corresponding cread on the other side of the channel can determine the length of the message. The checksum is useful because it allows the complimentary cread routine to detect communication errors when the data are read. While cwrite can easily calculate the checksum in software, generating the checksum in the channel hardware is relatively simple and allows a faster implementation of both cwrite and cread. If the hardware calculates the checksum, cwrite first writes the mes- sage except for the checksum and then instructs the channel hardware to append the checksum. Since at the user level, the length of the message is specified in terms of bytes rather than packets, and a header and trailer are added to the message, the number of bytes that cwrite transmits may not be a whole number of packets. However, at the sys- tem level the channel requires that data be transmitted only as whole packets, so if the data does not fill the last packet in the message, cwrite must pad it with additional bytes 

232 

crOS m: A Specific Implementation of the Programming Environment Chap. 14 

to fill the packet. Having discussed the fornlat of data to be sent by cwrite, we now describe how cwrite interacts with the channel hardware in order to send the packets containing the data. The process of sending packets on the selected channels consists of two basic steps. First, cwrite writes one packet, which the channel hardware simultaneously transmits on each of the channels indicated by the channel mask chan. After writing the packet, cwrite polls the transmit interrupt line for each of the indicated channels until an the channel buffers are empty .The length of time that cwrite spends polling the interrupt lines is indefinite, leading to the possibility of deadlock if the destination node never calls cread to empty its channel buffer by reading its contents. 
After deternlining that the buffers of the selected channels are empty , cwrite repeats the cycle of sending a packet and polling the transmit interrupt lines until all the data are transmitted; cwrite will not return until the whole message has been transmitted and read by the destination node. Unless an error occurs, cwrite returns the number of bytes in the message. Communication errors cannot be detected by cwrite since it receives no acknowledgement from the receiving node, but other types of errors may occur. If the buffer refers to data residing outside of the user's data space and the node has memory-management facilities, cwrite should indicate the error. The channel mask is invalid if it specifies nonexistent channels, and the number of bytes in the message could be too large for the system to handle. 
The routine cwrite indicates an error condition by returning -1 and identifies the specific type of error by setting the external integer errno. All the routines in crOS III use the convention of returning -1 to indicate an error and identifying the type of error by setting errno. Since errno is not reset after each error. it should only be examined after an error occurs, as indicated by a return value of -1. 
In order to incorporate the enhancements in the channel hardware that we dis- cussed in the previous section, we add an additional argument to the version of cread dis- cussed In Chap. 4. The new argument is a channel mask indicating the channels on which the message is to be forwarded as it is being read from the incoming channel buffer. To avoid problems with deadlock, the m:sJage must not be forwarded to the channel from which it is being read. In addition, ~ incoming channel is also specified as a channel mask, so that the channels are always labeled the same way. Thus, the capa- bility of the channel hardware to read a message and forward it in one operation is reflected in crOS III. We note that allowing .cread to forward messages as it is reading them requires relaxing the requirement that cread always be complemented by cwrite. A call to cread may now also be complemented by another call to cread that forwards the message. Of course, cwrite must still be caned at the source node since cread cannot ori- ginate a message. The usage of cread is defined in Code 14-2. 

Sec. 14-3 Basic Crystalline Routines 

233 

Code 14-2 The Routine cread 
The routine cread performs a sequence of operations complementary to those per- formed by cwrite. In order to read a packet, cread first polls the receive interrupt line of the channel specified by inchan until its channel buffer is full. As with cwrite, the length of time that cread spends polling the interrupt line is indefinite, leading to the possibility that cread will never return if the source node does not fill the buffer by sending a packet. When the incoming channel is ready, a packet is read from the channel buffer of the channel specified by inchan. As the packet is being read, the channel hardware simul- taneously transmits the packet on any channels indicated by outchan. The transmit inter- rupt lines of any forwarding channels are then polled until all of their channel buffers are empty as in cwrite. 
The number of bytes and, hence, the number of packets to be read by cread is determined by the message length, which is prepended to the message sent by cwrite. After reading the first packet of the message, cread knows how many additional packets it must read. The header containing the message length is kept by cread and the remain- ing bytes of the first packet are stored in buffer. The argument nbytes specifies the max- imum number of bytes that cread will store in memory .If the message length exceeds the maximum, cread truncates the message and returns -I to indicate the error. In addi- tion, the number of bytes in the untruncated message is saved in the external variable errcnt. In any case, the whole message is always read from the incoming channel and transmitted on any forwarding channels so that the complementary call to cwrite can finish sending the whole message and each of the receiving nodes can always receive the whole message. 
When cread receives the last packet in the message, it obtains the checksum appended to the data by cwrite. Any additional bytes following the checksum that were used to fill the last packet are discarded, though they must be read from the channel buffer in order to empty the buffer, allowing the next transmission to occur. By keeping its own checksum of the data in the message, cread can compare the two checksums to determine whether any communication errors occurred. If the checksums match, no errors occurred, but if they are different, cread returns -I to indicate the error. The appropriate value of errno is set to identify the communication error. Thus, even though the transmitting node cannot detect the communication error, the error can be detected before the data are used in subsequent computations. Since the transmitting node does not know an error occurred, it cannot simply retransmit the data without defining a higher-level protocol. However, in the present context, the error rate on the channels can usually be made low enough that error correction is not an important need. In addition to the errors just described, cread can report errors in any of its arguments, as was the case 

234 

CJOS ill: A Specific Implementation of the Programming Environment Chap. 14 

with cwrite. In particular, cread repom an error if the channel mask inchan indicates more than one channel. 
In some cases, an application program may want to forward the incoming message without keeping a copy in local memory .If buffer has the special value defined by the macro NULLPTR, the data is only written to the channels indicated by outchan. In addi- tion, we define several other special cases for the arguments of cread. If outchan is 0, the message is read into memory without forwarding on any channels. If inchan is 0, the call to cread is ignored. and it should not be be complemented by a cwrite on the other side of the channel. However, if inchan is not 0, and the byte count is 0, a call to cwrite must still complement the call to cread, even though cread will not write any bytes into memory. 
We are now ready to examine the synchronizing effects of the crystalline routines cwrite and cread in more detail than in Sec. 4-3. In particular, we want to justify our decision to prevent cwrite from returning until the destination node calls cread. While the implementation of cwrite would be more analogous to that of cread if it first polled the channel status and then transmitted a packet, such an implementation leads to an undesirable feature in the synchronizing condition that cwrite imposes. We illustrate the problem by considering the effect of writing a message that is short enough to be sent in a single packet If cwrite polled the channel status before sending the message, it could return before the destination node called cread, since the message is sent in a single packet. However, if the message required more than one packet to be sent, cwrite could not return until the destination node called cread, since the second packet could not be sent until the first was read. Thus, the synchronizing condition imposed by cwrite would depend on the message length and on the size of the channel buffer. 
A more severe problem than the dependence of the synchronizing conditions on the size of the channel buffer, is that the response of cwrite to a programming error also depends on the size of the channel buffer. To illustrate the problem, we consider a pro- gramming error in which two nodes write to each other at the same time. Such a situa- tion is an obvious error since we have stated that a call to cwrite must be complemented by a call to cread in the destination node. If both nodes attempt to write messages that require more than one packet to send, deadlock occurs after they both write their first packet, since neither node can call cread to read the incoming message until it has finished writing its message. However, if both nodes wrote a single packet message, deadlock would not occur because each node could write the single packet into the chan- nel buffer without waiting for the destination node to call cread. Then after calling cwrite, both nodes could call cread to read the messages. Thus, the effect of making a programming error would depend on the length of the message and the size of the chan- nel buffer. Such a situation is a serious problem because simply changing the length of a message or moving the prognun to a concurrent processor with a different packet size could cause a working program to break, or conversely, cause a formally incorrect pro- gram to work. 
A preferable choice is to ensure that cwrite always causes deadlock when it is used incorrectly, regardless of the message length or the packet size. By polling the status of 

Sec. 14-3 Basic Crystalline Routines 

235 

the channel buffer after sending a packet, rather than before, we guarantee that cwrite can never return before the destination node calls its complementary cread. The channel status does not need to be tested before writing a packet since cwrite always leaves the channel ready for the next packet to be written. With such an implementation of cwrite, the synchronizing condition that it imposes is exactly analogous to that of cread: neither routine can ever return until its complement is called. 
Since we are inttoducing a crystalline communication system that is meant to pro- vide communication facilities for real application programs, we expand our list of basic routines to include three additional routines. Two of the routines, v write and vread, are closely related to cwrite and cread but are more flexible in the types of data that they can handle. While cwrite and cread can only write and read contiguous blocks of memory, vwrite and vread write and read vectors whose elements are regularly spaced in memory. The ability to write such vectors is very useful in certain problems in which the data to be communicated are not contiguous. A simple example is sending the first row of a matrix when the matrix is stored by columns. Without the ability to send vectors, the user would have to first gather the data into a contiguous block and then use cwrite to transmit the data, a process that is cumbersome and inefficient 
Instead of accepting only the total number of bytes to transmit as cwrite does, vwrite requires as input the size of the data items, the separation between data items, and the total number of such items. The usage of vwrite is defined in Code 14-3. 

Code 14.3 The Routine vwrite 
Any vector whose elements are regularly spaced can be specified as an argument. As with cwrite, a call to vwrite must be complemented by a call to vread in the destina- tion node. The implementation of vwrite is identical to that of cwrite, except that when it writes packets to the channel, it writes size bytes, increments the memory pointer by offset bytes, and writes another size bytes, until the packet is full. If the packet becomes full before all size bytes are written, vwrite finishes writing the data item in the next packet. We note that offset can be negative so that data objects can be sent in reverse order from buffer. Another distinction of v write is that the value of nitems is used as the message header instead of the byte count as in the case of cwrite. The polling of the channels and the resulting synchronizing condition that vwrite imposes is identical to that of cwrite. A checksum is also sent with the message in the same way as in cwrite and the error conditions that vwrite can report are similar. 
The routine vread incorporates the same conventions in its arguments as vwrite. The argument nbytes in cread is replaced by the three arguments, size, offset, and nitems. 

236 

crOS ill: A Specific Implementation of the Programming Environment Chap. 14 

The new argumen~ have the same meaning as in vwrite, except that nitems is the max- imum number of data items that vread will read. The actual number of data items read is determined by the number sent by vwrite as indicated in the message header. The usage of vread is defined in Code 14-4. 

Code !4.4 The Routine vread 
The implementation of vread is again identical to that of cread, except that the data items are not written into memory in a contiguous block. Instead, the first item is written at the start of buffer, the second item is written at byte number offset of buffer, and so on. The procedure continues until either all the data items in the message have been read from the channel and stored in memory , or the maximum number of data items specified by nitems has been reached. The routine vread and its complementary v write must agree on the size of the data items, but they do not have to use the same value of offset. A com- munication error can be detected by vread by comparing the checksum it calculates with the one that vwrite appended to the message. The errors that vread can report are similar to those associated with cread. 
The final routine that we introduce in this section is rdstat, which provides informa- tion about the status of the channels. Since cwrite and vwrite always leave a channel ready to write another packet, checking the transmit interrupt lines is pointless. How- ever, checking whether a message is ready to be received may be useful in some cases. For instance, an application may increment a counter while waiting for a message to arrive in order to measure the amount of time it had to wait By calling rdstat, the appli- cation can determine when the message arrives so it can stop incrementing the counter and call cread or vread to read the message. 

int rdstat( chan) unsigned int chan; 

1* mask of channels to check *1 

Code 14.5 The Routine rd\'tal 
The routine rdstat, defined in Code 14-5, tests the status of the channels indicated by the channel mask chan. If messages are ready to be received on all of the indicated channels, rdstat returns I. If the channel mask is invalid because it specifies a nonex- istent channel, rdstat returns -I to flag the error. Otherwise, it returns 0. 

Sec. 14-4 The Hypercube Interconnection TolX>logy 

237 

14-4 The Hypercube Interconnection Topology 
In the previous section we discussed the implementation of the basic crystalline communication routines without explicit reference to the interconnection topology of the concurrent processor. However, before we continue with the discussion of implementing the collective communication routines, we must select an interconnection topology. Without a specific topology we cannot describe algorithms that implement the collective routines. The topology we select should be sufficiently general to be useful on a rich variety of problems, yet practical enough that a real concurrent processor can be built with that topology . 
First, we shall briefly examine the fully interconnected topology that we assumed for our virtual concurrent processor in Chap. 4. Having an pairs of nodes directly con- nected by a communication channel is clearly the most general interconnection scheme, which is why it was adopted for the virtual concurrent processor. However, building a real machine with so many connections is completely impractical unless the size of the machine is limited to relatively few nodes. If the ensemble consists of N nodes, each node must have N -1 communication channels so it can communicate with all the other nodes in the system. The total number of channels in the ensemble is then N (N -1)/2. The number of channels required per node and the total number of channels, which are usually connected by some type of backplane, increases so rapidly with the number of nodes that building a large system using the fully interconnected topology is out of the question. 
In Chap. 9, we saw that only a ring topology is required for problems with long- range interactions. For an ensemble to be connected as a ring, each node needs two com- munication channels, regardless of the total number of nodes in the ensemble. Since the number of channels per node is fixed, a very large system can be connected as a ring. The ability to scale a concurrent processor architecture easily to larger systems is very desirable since this allows application prognuns to be written in a manner which is independent of ensemble size. However, the ring topology has a serious drawback because the number of channels is so restricted. Efficiently programming a problem that requires a richer interconnection topology than the ring may not be possible on such machines. Many problems require that data be sent between nodes that are not directly connected by a communication channel, and the maximum distance between nodes grows rapidly in a ring topology. In contrast to the fully interconnected topology in which all nodes are directly connected, the maximum distance between nodes in the ring is N /2. Such a rapid growth in the maximum distance between nodes motivates us to consider other connection topologies. 
By examining the extremes of the fully interconnected topology and the ring topol- ogy, we can see the trade-offs involved in selecting an interconnection topology for a real machine. An ensemble with few channels per node, such as in the ring case, is easily built and scales easily as the size of the ensemble grows, but the maximum distance between nodes grows rapidly, restricting its usefulness. On the other hand, a machine with many channels per node, such as one using the fully interconnected topology, is very general since all nodes are directly connected by a communication channel, but it is 

238 

crOS ill: A Specific Implementation of the Programming Environment Chap. 14 

not practical for a large system since the number of channels required becomes prohibi- tive with increasing ensemble size. 
Keeping these trade-offs in mind, we consider the communication requirements of some of the algorithms presented in previous chapters. Several of the applications, such as the matrix problems described in Chap. 10, require only a two-dimensional grid con- nection. An ensemble with such an interconnection topology can easily be built since the number of channels per node is fixed, as it is in the ring topology. A two-dimensional grid connection with periodic boundaries requires four channels per node, regardless of the total number of nodes in the ensemble. Since the number of channels per node is not affected by the system size, the two-dimensional grid can be used to build large systems. The maximum distance between nodes in a square grid with periodic boundaries, i.e., a torus, is {f[ , where N is the total number of nodes in the ensemble. Thus, we see that the two-dimensional grid topology requires twice as many channels per node as the ring topology but is more flexible because more nodes are directly connected by a channel, and the maximum distance between nodes grows more slowly as the size of the ensemble is increased. 
Another important feature of the two-dimensional grid connection is that the ring topology is contained within the grid. By ignoring certain communication channels, a ring topology can be obtained for applications that require only a ring connection. The results for the two-dimensional grid can be genemlized to any number of dimensions. The number of channels required for a d -dimensional grid is 2d and is independent of N , the number of nodes in the ensemble. The maximum distance between nodes in a d- dimensional grid with periodic boundaries is N lid. All grid connections scale nicely with the ensemble size, though the nodes become more complicated as the dimensionality of the grid increases since the number of communication channels per node also increases. Another desirable feature of grid connections is that grids of each lower dimensionality are contained within the topology . 
Though grid topologies are convenient for solving many problems, a grid whose dimensionality is fixed is not always adequate. Consider the case of the Fast Fourier Tmnsform. The algorithm given in Chap. 11 requires that the dimensionality of the grid increase as the number of nodes in the ensemble increases. An interconnection structure whose dimensionality depends on the number of nodes is the hypercube interconnection topology .The hypercube topology can easily be described inductively. A zero- dimensional hypercube, a O-cube for short, consists of a single node with no communica- tion channels, i.e., a standard sequential computer. A l-cube is constructed from two 0- cubes by connecting them with a single communication channel. A 2-cube is constructed from two l-cubes by connecting their corresponding nodes with an additional channel, forming a square with the nodes at the vertices and the communication channels on the edges. In general, a d -cube is constructed by connecting the corresponding nodes of two (d-l)-cubes with an additional channel. 
The number of nodes in a hypercube is N = 2d, and the number of communication channels per node is d. Thus, we see that the number of channels required per node increases as the size of the ensemble increases, unlike the fixed-dimensional grid 

Sec. 14-4 The Hypercube Interconnection TolX>logy 

239 

topologies. In tenus of N, the number of channels is logzN .Although the number of channels per node increases with N, the increase is slow enough that comparatively large ensembles are practical. The maximum distance between nodes of a hypercube grows more slowly with increasing N than any fixed-dimensional grid; the maximum distance is logzN .As in the case of the fixed-dimensional grid, the hypercube contains grids of each lower dimensionality. 
By now the direction of our discussion should be apparent The hypercube is ade- quate for all the applications we have discussed so far, yet it is practical for fairly large ensembles. A balance between the number of channels per node and the maximum dis- tance between nodes is achieved, since both increase logarithmically with the total number of nodes in the ensemble. In order to discuss the implementation of the crystal- line communication system in concrete tenus, we have chosen to use the hypercube inter- connection topology in the remainder of the chapter. While much of the discussion will be specific to the hypercube, many of the general ideas can be applied to other topolo- gies. 

Before resuming discussion of the implementation of the crystalline communica- tion system, we need a specific numbering of the nodes and channels in the hypercube. A d -dimensional hypercube consists of 2d nodes, each of which is connected to one other node in each dimension. Thus each node requires d channels. The processor number of a node is conveniently chosen to be a d -bit binary number where the i th bit of the processor number represents the coordinate of the node in the ith dimension of the hypercube. This immediately suggests a convenient convention for the local numbering of channels connected to a given node, since any two neighboring nodes in the hypercube have processor numbers that differ in exactly one bit position. Thus, in a 3-cube, chan- nels 0, 1, and 2 connect node ° to nodes 1, 2, and 4, respectively. Fig. 14-1 illustrates this numbering scheme for a hypercube of dimensionality 3. 
The important point here is that all nodes connected by a communication channel differ by a single bit in their processor numbers. The number of bits that differ between the processor numbers of two nodes gives the "distance" between the nodes. In program- ming examples using the hypercube topology, we will assume that in addition to the vari- ables procnum and nproc, which contain the processor number and the total number of nodes, respectively, the variable doc (dimension of cube} will contain the hypercube dimension. 
Because of the symmetry of the hypercube, any node could be assigned the proces- sor number ° arbitrarily, using the hypercube topology and the numbering scheme to determine the remaining processor numbers. However, we do have one asymmetry in the concurrent processor, namely the control process. Until now, we have said nothing about where the control process resides. In order to facilitate the discussion of the collec- tive routines involving the control process, we now assign the control process to its own processing element, which is distinct from the nodes composing the hypercube. We allow only one of the nodes in the hypercube to communicate directly with the control processor through a communication channel that is identical to those connecting the nodes in the ensemble. The node connected to the control processor will be assigned a 

240 

crOS ill: A Specific Implementation of the Programming Environment Chap. 14 

001 

011 

Control Jrocessor.. --v ~ 

L 

y101 

/ 

1 

J 

11 

100 110 
l'lgure 14.1 Labeling of processor numbers and channels for a three-dimensional hyper- cube. 

processor number of 0. A variable cpmask is defined which contains a channel mask to specify the channel connecting node ° and the control processor. Of course, the number of dIe channel connecting node ° and the control processor must be at least as large as dIe largest available hypercube dimension. In dIe nodes odIer dIan 0, cpmask contains ° so dIat any attempted communication widI dIe control processor is ignored. The variable cpmask will be used in Sec. 14-6, which describes dIe collective routines implementing communications between dIe nodes and the control process. 

14.5 Collective Communication Routines 
Having decided to resbict our discussion to a concurrent processor with a hyper- cube interconnection topology, we see more clearly the second reason given in Chap. 4 for using collective communication routines. In addition to providing a conceptual grouping to guide our thinking, the collecti,,~ routines allow machine dependencies in their implementation to be hidden from the application program. In this section we will examine implementations on a hypercube of the collective routines introduced in Sec. 4-4, using the basic communication routines described in Sec. 14-3. In addition, we introduce two new routines in this section that are part of CrOS III. Chapters 21 and 22 have other useful communication utilities which we consider as part of a user library rather than the basic crOS III system. 

Sec.14-5 Collective Communication Routines 

14-5.1 cshift 

241 

The first collective routine we will implement is the simplest one, cshift. The pur- pose of cshift is to allow a group of nodes to both write and read messages among them- selves without the application program specifying which nodes write first and which nodes read first For instance, a pair of nodes can exchange messages, and a group of nodes can pass messages one step around a ring using cshift. The essential feature is that cshift is complemented on the other side of a channel by itself, rather than by a different routine. Since all the nodes in applications for which the crystalline communication rou- tines are useful typically both read and write together, using cshift instead of cwrite and cread makes programming such applications much easier. In fact, cshift and the other collective communication routines we discuss in this section should eliminate the need to explicitly use cwrite and cread in many application programs. 
Stating the purpose of cshift makes most of its implementation obvious. It must call both cwrite and cread, but the order in which it calls them must be different for dif- ferent nodes. Since the only available means of distinguishing the nodes is their proces- sor numbers, the order of reading and writing must be based on only a node's processor number and the hypercube topology. 
As discussed above, an important feature of the scheme that we use for assigning processor numbers to the nodes of the hypercube is that nodes connected by a communi- cation channel differ by a single bit in their processor numbers. Thus, nodes with an even number of one bits in their processor numbers are only connected to nodes with an odd number of one bits in their processor numbers, and vice versa. Thus, we can use the number of one bits in a node's processor number to determine whether the node reads or writes first in cshift. 
Since a node's processor number does not change, we can calculate the number of one bits in an initialization routine to avoid recalculating it every time cshift is called. Instead of storing the number of one bits, we should just store a flag indicating whether the number of one bits is even. Let us assume that the initialization has been performed and the result stored in the variable iseven. If the number of one bits in a node's proces- sor number is even, then iseven contains I; otherwise, it contains 0. 
Given the initialized variable iseven, the implementation of cshift is almost trivial, and is given in Code 14-6. The implementation shows that we have changed two of the arguments of cshift from the calling sequence of vm-shift introduced in Sec. 4-4. Instead of using the channel number, we specify the channels on which the outgoing data are to be written with the channel mask outchan. Likewise, the channel from which the mes- sage is to be read is indicated by a channel mask for consistency with the way the output channels are specified. The value returned by cshift is the number of bytes read unless an error occurs, in which case -I is returned. 

242 

crOS m: A Specific Implementation of the Programming Environment Chap. 14 

Code 14.6 The Routine cshift 
In the current implementation of cshift, the buffers provided by the application for the incoming and outgoing messages must not overlap. Otherwise, cshift corrupts the outgoing data in the nodes in which it calls cread before calling cwrite. However, in many applications having the incoming data overwrite the outgoing data is desirable, since the outgoing data are no longer needed after being transmitted. Thus, a useful enhancement to cshift as presented here would allow outbuf and inbuf to be identical. Such an enhancement could be implemented in several ways. The most obvious imple- mentation involves copying the contents of outbuf into a temporary buffer in those nodes where cread is called' first. After calling cread, cshift calls cwrite to transmit the data in the temporary buffer. In order to avoid problems with the availability of a temporary buffer of adequate size, cshift could declare a buffer of some reasonable size and transfer the message in pieces no larger than the temporary buffer. 
If the communication channels of a particular machine make use of full-duplex communication, the implementation could use the channel buffers as the temporary buffer by alternating between writing a single packet into the channel buffer and then reading an incoming packet. Of course, cshift would directly access the communication channels in such an implementation, rather than using cwrite and cread. Using the chan- nel buffer as the temporary buffer is faster than copying the outgoing message into a buffer in memory, which involves copying the message and additional time. However, such an implementation may not be optimal in a concurrent processor with half -duplex communication channels, because both sides of the channel cannot write simultaneously 

Sec. 14-5 Collective Communication Routines 

243 

in all cases. To prevent problems arising from both sides of a channel attempting to write simultaneously, a node must test to determine whether the node on the other side of the channel is writing a packet before writing its own packet. Such a test slows the data transfer rate, and is not needed in an implementation that requires writes on a channel to be complemented by reads, such as cwrite and cread. Thus, we see that when inbuf and outbuf are the same in cshift, the communication speed depends on the details of the hardware implementation of the channels. In order to obtain maximum communication speed on any system, inbuf and outbuf should refer to distinct regions of application memory. 
Several additional features are built into the standard crOS III definition of cshift that are not included in the simplified definition of Code 14-6. In addition to allowing inbuf and outbuf to be identical, the special cases already discussed for cwrite and cread are supported, such as the suppression of the read or write when the channel mask is 0. One additional special case is provided to allow for easier programming of applications that use periodic boundaries. An important feature of a good implementation of a hyper- cube application is that a specific dimensionality of the hypercube not be programmed into the implementation. If the algorithm uses a two-dimensional grid topology and requires that each node writes a message to its left neighbor and reads a message from its right neighbor, periodic boundaries do not cause any problems for cshift. However, if the application is run with only one node in the grid direction in which the communication occurs, cshift no longer works because the node has no neighbors with which to trade messages, so both inchan and outchan are 0. The application will still work in such a situation if the contents of outbuf are copied to inbuf, as if the node had sent a message to itself. Since calling cshift with both inchan and outchan set to 0 has no other meaning, no ambiguity arises from adding such a copying feature to cshift. In particular, if non- periodic boundaries were used, cshift might also be called with inchan and outchan set to 0 to avoid unnecessary conditionals in the application program. However, in such a case the data in inbuf is ignored anyway, since the application would not try to use data from outside the boundary of the decomposition. 
~ 

14-5.2 vshift 

One of the new routines we introduce in this section is the collective communica- tion routine, vshift, which is similar to cshift except that it uses vwrite and vread instead of cwrite and cread in its implementation. All the comments on the use of cshift and its special features apply to vshift, except that vshift will work correctly for identical values of inbuf and outbuf only if insize and outsize and inoffset and outoffset are also identical. If all three pairs of arguments are not identical, the special feature does not apply and errors may result. Since the implementations of cshift and vshift are so similar, we will not repeat the previous discussion. As in vwrite and vread, the byte counts are replaced with the size of the data items, their offset from each other, and their number. The cal- ling syntax is given in Code 14-7. 

244 

crOS ill: A Specific Implementation of the Programming Environment Chap. 14 

Code 14.7 The Routine vshift 

14-5.3 broadcast 
The next collective communication routine we examine is broadcast. Instead of implementing a version of broadcast that is as general as the one described in Sec. 4-4, we impose a restriction on which nodes can be included in a broadcast. Since specifying a list of all nodes involved in the broadcast is rather cumbersome, we restrict broadcasts to subcubes within the hypercube topology. A subcube is a hypercube whose dimen- sionality is less than or equal to the dimensionality of the full hypercube being used by the application. An important feature of subcubes is that when the hypercube is mapped onto decomposition topologies using the specific gridmap routines that we discuss in Sec. 14-7, all of the nodes along a dimension of the topology comprise a subcube. Thus, broadcast is useful in matrix algorithms for sending a submatrix from one node to all other nodes along a row of the matrix as described in Chap. 10. Since a subset of the available channels span any subcube, a node can uniquely specify a subcube by indicat- ing the desired subset of channels. A channel mask is an efficient means of specifying the subset of channels involved in the broadcast. 
All nodes involved in the subcube broadcast must specify the same subset of chan- nels and the same source node. The source node then writes the message in its buffer on the indicated channels. Since all the nodes in the subcube are not connected to the source node by a communication channel, some of the nodes must retransmit the message received from the source node. The broadcast algorithm is a tree algorithm that results in each node in the subcube receiving a single copy of the message sent by the source node. 
Before treating the more general case, we will examine the algorithm for broad- casting to a hypercube when the source node is node 0. The procedure for broadcasting to a O-cube is trivial since node 0 is the only node in the system. The situation is barely more complicated for a l-cube since the only possibility is for node 0 to send the mes- sage to node 1. The procedure for broadcasting to a 2-cube is similar to that for a l-cube since a 2-cube consists of two l-cubes whose corresponding nodes are connected by an additional communication channel. Node 2 corresponds to node 0 and node 3 corresponds to node 1. Thus, node 2 sends the message to node 3 and node 0 sends the 

Sec. 14-5 Collective Communication Routines 

245 

message to node 1. Since node 2 is connected to node 0, node ° can send the message to node 2. 
The broadcast algorithm continues by repeating the same procedure for each higher dimension. The broadcast tree of a d-cube contains broadcast trees for two (d -1)- cubes with node ° sending the message on the one additional channel connecting it to the second (d -1)-cube. The broadcast tree is illustrated in Fig. 14-2 for a 4-cube. 

1A 

tIIO tIlt 

Broodcast Tree for 4-cube 

0000 
~~ooo 00010100 0010 
1001 
010 0110 0101 0011 
1011 0111 

'01101 

F1gure 14-2 broadcast tree for a four-dimensiooal hypercube. 
Having described the broadcast tree, we must now design an algorithm that allows each node in the tree to determine on which channel to read and on which channels, if any, to retransmit the message. As with cshift, the determination can be based on a knowledge only of the hypercube topology and the node's processor number. The chan- nel on which to receive the broadcast message will be indicated by the channel mask rchan, and the channels on which to retransmit the message will be indicated by the channel mask wchan. By examining Fig. 14-2, we see that rchan consists of only the rightmost one bit in the node's processor number, and wchan is the mask of all channel numbers less than rchan. In the special case of the source node (node 0), rchan is not used and wchan contains all the channels in the hypercube. 
Now we are ready to genemlize the algorithm for broadcasting from node ° to a full hypercube to allow it to broadcast to a subcube using any node as the source node. The method of genemlizing the previous algorithm consists of one transformation of each node's processor number, and an additional condition for determining rchan and wchan. The transformation consists of combining the node's processor number and the source node's processor number with the exclusive-OR operation. The exclusive-OR operation leaves a one bit in the bit locations where the two processor numbers differ and a zero bit where they are identical. The effect of the exclusive-OR operation is to temporarily renumber the nodes of the subcube so that the source node becomes node 0. Thus. we 

246 

crOS III: A Specific Implementation of the Programming Environment Chap. 14 

can use the previous algorithm, which specifies how to broadcast from node 0, to broad- cast to the whole hypercube using any node as the source node. 
In order to restrict the broadcast to a subcube, we impose the additional condition that only those channels spanning the subcube are to be used in determining rchan and wchan. With the additional condition, rchan specifies the channel corresponding to the rightmost one bit in the transformed processor number that is also one of the channels spanning the subcube. Likewise, wchan is the mask of all channel numbers less than rchan that are also contained in the mask of channels spanning the subcube. For the spe- cial case of the source node, wchan is the same as the mask of channels spanning the sub- cube. 

We have now taken the algorithm for broadcasting from node 0 to the whole hyper- cube and generalized it so that any node can broadcast a message to any subcube. We present the routine that implements the general broadcast algorithm in Code 14-8. 
From this implementation we see that the value returned by broadcast is the value returned by cwrite in the source node and the value returned by cread in all other nodes. We also note that the broadcast algorithm is implemented with only one call to a basic communication routine in all of the nodes, since we have allowed cwrite to write to mul- tiple channels at once and cread to both read and forward the message in one operation. 
If the hardware implementation of the communications is as described in Sec. 14-2, so that writing on multiple channels and reading and forwarding on multiple channels are actually single operations, the broadcast algorithm is optimal for all subcube dimensions and all message lengths. The first packet of a broadcast to a d -dimensional subcube arrives at the nodes farthest from the source node in d steps, which is the minimum number of steps possible with the hypercube topology .Each additional packet in the message takes one unit of time since each node in the broadcast tree forwards it in one operation. Thus, the algorithm must be optimal since each packet cannot require less than a single operation. 
However, if the channels were implemented so that each reading and writing operation was done separately, the broadcast algorithm would not be optimal for all sub- cube dimensions and all message lengths. In broadcasting to a d -dimensional subcube, the source node would always need d operations to transmit a packet, since it has d channels on which it must write the packet. If the message were broadcast in a chain with the source node at the middle of the chain, only two operations would be needed to broadcast each packet While the first packet would take longer to arrive at the furthest nodes when d > 2, if the message were many packets long, the smaller number of steps per packet would make the chain broadcast a better algorithm. This type of broadcast was called a pipe in Sec. 10-3.2, and will be elaborated on in Sec. 20-4 as the better split yipe. Thus, we see that properly designed channel hardware makes optimal sub- cube broadcasts easier to obtain. A more detailed discussion of broadcast algorithms for hypercubes whose nodes cannot write to multiple channels and read and forward in a sin- gle operation is given in [Johnson 85], [Ho 86a], and [Fox 86b]. This topic will be dis- cussed further in Sec. 21-4.1. 

Sec. 14-5 Collective Communication Routines 

247 

Code 14.8 The Routine broodcast 

14.5.4 combine 
The fourth collective communication routine that we will implement in this section is combine. Unlike vshift and broadcast, we will use exactly the same arguments for combine as we used in Sec. 4-4. In contrast to broadcast, combine will be implemented with a different type of tree algorithm. In fact, the tree that is utilized by the combine algorithm is actually a collection of binary trees with each node in the hypercube at the top of one of the trees. Fig. 14.3 illustrates the collection of binary trees. 
Each node starts with a block of data items that are to be combined with the similar blocks from all other nodes. First, all nodes exchange their blocks across channel 0. Each node then uses the combining function provided by the application program to com- bine the items in their original block with those in the block received on channel 0. Next. all nodes exchange their partial results from the first step across channel 1. The blocks received on channell are then combined with the current partial result to give a new par- rial result. The process of exchanging current partial results and combining the two 

248 

ctOs ill: A Specific Implementation of the Programming Environment Chap. 14 

partial results to obtain a new partial result is repeated for each channel in the hypercube. The resulting block is the combination of the original blocks from each of the nodes. This communication strategy is generalized to the method of ordered cube geodesics in Chap.21. The example implementation of combine is given in Code 14-9. 

Combine Tree for 3 dimensional hypercube 

000 001 010 011 100 101 110 FIgure 14.3 combine tree for a three-dirnensional hypercube. 

III 

int combine(bur, func, size, nitems) 
char *bur; 1* pointer to data buffer *1 
int (*func)(); 1* pointer to the combining function *1 int size; 1* size of data items being combined *1 int nitems; 1* number of data items to combine *1 

int i, j, temp[BIG], status, error; if(size > BIG) return(-l); error = 0; 
for(i = 0; i < doc; ++i) { 
status = cshift(temp, 1«i, BIG, buf,l«i, size*nitems); if(status = -1) error = -1; 
for(j = 0; j < nitems*size; j += size) 
if«*fimc)(buf + j, temp + j, size) = -1) error = -1; 

} 
if(error = -1) remm(error) else remm(status); 

} 

Code 14.9 The Routine combine 

Sec. 14-5 Collective Communication Routines 

249 

In this implementation, a temporary buffer is used to store the block of data items obtained from cshift. We have avoided the issue of how large the temporary buffer should be by defining its size with the macro BIG. In a more sophisticated application, all nitems would not need to fit into the temporary buffer since each item could be exchanged separately and combined, rather than exchanging all the items together. Still, the problem remains that at least one of the data items must be buffered temporarily, and the data items might be very large. One way of dealing with the problem is to have com- bine return -I to indicate an error if size is larger than its temporary buffer, as we have done. 

The combining function,tunc, has three arguments because in some cases allowing a single combining function to operate on variable-sized data objects may be useful. A generalized combining function could then obviate the need for having distinct functions for every instance of a general group of functions that operate on different sizes of data objects. For instance, the combining function could implement a bitwise OR of the two data items. If the combining function handles only single byte objects, combine must call tunc once for every byte to be combined. Although such a method works, it is slow because of the large number of calls totunc. A more efficient implementation would use a combining function that would perform the bitwise OR of two data objects that are size bytes in length, so that tunc would only need to be called once. In any case, the combin- ing function can always ignore the third argument, so the increased flexibility costs very little, even if it is not used. 
The algorithm for combining the blocks of data items from each node is redundant since each node calculates the same result However, the result for a single data item is still calculated in a minimum number of steps, since all of the redundant calculations are performed in parallel. The number of steps needed to combine the blocks of data items from all the nodes in a d -dimensional hypercube is d, which is the same as the maximum distance between nodes in the hypercube measured in communication steps. However, if the number of data objects being combined is large, the algorithm used in combine may not be optimal. Allowing different nodes to combine different objects and then distribute the results at the end of the calculation may be faster. A more sophisticated combining algorithm that distributes the calculation and collects the resulting combined items will be discussed in Chap. 19. 

l4-5.5 concat 

The last collective communication routine that we present in this section uses an algorithm that is similar to combine. The routine concat broadcasts blocks of data from each of the nodes to each of the other nodes in N -I steps, where N is the number of nodes. The algorithm uses the same set of binary trees illustrated in Fig. 14-3, except that a concatenation replaces each combining operation. 
Each node starts with a block of data that is to be broadcast to each of the other nodes. First, all nodes copy the block of data into the buffer that will contain a block from each of the nodes in order of processor numbers. The buffer contains nproc blocks of size maxsz bytes, since maxsz is the maximum size of any data block. After copying 

250 

crOS m: A Specific Implementation of the Programming Environment Chap. 14 

its own block into the buffer, each node writes its block on channel ° and reads a similar block from channel 0, storing the block in the appropriate location in the buffer. Next, each node exchanges the two blocks in its buffer with those of the node connected by channel 1. The process repeats for each of the doc channels in the hypercube, with twice as many blocks being exchanged in each successive step. When the process is complete, the buffer in each of the nodes contains blocks of data from each of the nodes. Since twice as many blocks of data are exchanged on each successive channel, a total of N -1 blocks of data must be exchanged by each node. 

int concat (fromptr, nbytes, toptr, maxsz, szarray) char .fromptr; 1. pointer to data block .1 
int nbytes; I. number of bytes in data block .1 char .toptr; I. pointer to buffer for all blocks .1 int maxsz; I. maximum size of any block .1 int .szarray; I. pointer to array of block sizes .1 

int i, j, send, recv, chan, status; char *in, *out; 

for (i=O;i<nbytes;++i) top1r(procnum*maxsz+i]=fromp1r(i]; szarray (procnum]=nbytes; status = 0; 

for (i=O; i < doc; ++i){ chan = 1«i; 
send = procnum/chan; 
for (j = send; j < (send + chan); ++j){ recv = j.chan; 
in = toplr + (recv*maxsz); out = toplr + (j*maxsz); 
szarray [recv] = cshift (in,chan,maxsz,out,chan,szarray[j]); if (szarray[recv]==-l) status= -1; 

if (status =:0) status=nproc; return (status); 

} 

Code 14-10 The Routine concat 
The array, szarray, contains the length of each block of data since a node may send a block that is smaller than maxsz. However, the buffer containing the concatenated blocks is indexed as if all of the blocks contained maxsz bytes. The routine concat calls cshift once for each block it receives, giving a total of N -I calls. Since N -I blocks must be received by each node in order to receive a block from every other node, and the com- munication occurs concurrently in all of the nodes, concat is the optimal method of