Sec. 5-3 Example Software StructUJ'e 

101 

F1gure 5-2 The example of Fig. 5-1 drawn to emphasize the processor connection topol- ogy .As the blow-up of processors 2, 6 and 7 illustrates, the update of "boundary" nodal points in each processor requires communication with nearest neighbor processors in a one-dirnensional chain, which is a subset of the connections available in a hypercube (shown in inset). 
Software Supplement to this volume provides a full implementation of the illustrative text examples. 
We begin by presenting the control process program. The purpose of this program. which runs only in the control processor. is to act as coordinator of collective ensemble operations and to mediate such external input/output activity as is required. In this exam- pIe. the control process program obtains as input data the information needed to set up the decomposed domain and initial conditions. It then instructs the ensemble to com- plete the desired number of time steps. and finally receives the results. 
An examination of the pseudocode implementation in Code 5-1 (and the corresponding FORTRAN and C programs given in Appendix D) shows that the control process program initiates actions in the ensemble by issuing global communications to which the node programs respond. We may note in passing that if the number of nodal points in the problem domain is not evenly divisible by N. the number of processor nodes. the calculation cannot be perfectly load balanced. In such a case. the best that can be done is to distribute the load such that there exists no more than a one-point difference between the loads in any pair of nodes. This load imbalance condition will enter into our efficiency discussion in Sec. 5-4. During the main task of update computation. the con- trol process simply waits for the completion of the task and thus has no impact upon the concurrent efficiency. Finally. dumping and printing of the final solution occurs via controVensemble communication. Although this might not be the best way to solve the wave equation in practice. it provides a useful illustration of the important issues 

102 

Introduction to Decomposition and Concurrent Algoritluns: Wave Equation Chap.5 

involved in such a problem. 

proc-begin main (run the control process program) for-begin (interactive conunand loop) 
[read command from terminal or standard input] [ vm-bcastcp-broadcast command to ensemble] if-begin (conunand is INITIAUZE) then 
proc-call initialize yroblem (signal ensemble to set starting values) else-if (conunand is UPDATE) then 
[read number of steps, k ] 
[vm-bcastcp-broadcastk to ensemble] else-if (conunand is OUTPUT) then 
proc-call output-results (receive results from ensemble using vm-dwnpcp and direct it to output stream) 

else-if (command is STOP) then [terminate program] 
if end for end proc-end 

Code 5-1 Cootrol Process Pseudocode Program for Wave Equation Solver 
We next turn our attention to the ensemble node program. The node program has a structure which complements that of the control process program. During the setup and output phases of the programs, control process "write" calls are balanced by corresponding node "read" calls, and vice versa. The node program, however, contains as its most important part the section that performs the time step calculation. communi- cations calls in this part are for the mutual exchange of data between ensemble proces- sors, and do not involve the control process. 
The node program outlined in Codes 5-2 and 5-3 runs independently in each of the processors. The processors take requested action depending upon the values of the parameters passed to them from the control process and their neighbors. 
The node program begins by defining the communication links required for the linear decomposition topology of this problem. Following this preparatory step, the node program receives data from the control process necessary to generate the initial configuration of the vibrating string. Once all the starting values are obtained, the nodes begin their update steps. During this "working" phase, the control process is not involved and only communication between nodes occurs. After the specified number of updates, the node program may receive the request to forward its results through the control pro- cess. 

Essentially all of the interesting action associated with the problem solution occurs within this node program. Its flow centers on the central solution algorithm, and is punc- tuated by its points of contact with the control process, and with its "sibling" node processes. 

Sec. 5-3 Example Software Structure 

103 

Code 5.2 Main Node Program for Wave Equation Solver 
The above example program serves to illustrate the practical issues associated with implementing a particular problem decomposition within the virtual machine environ- ment defined in Chap. 4. Further practical details may be found in Appendix D which includes explicit C and FORTRAN control and node programs respectively. 
It may be seen in this example that the communication calls unique to the con- current computation environment are infrequent enough that the basic structure of the corresponding sequential algorithm is more or less preserved. On the other hand, the intricate "handshaking" required between the control process and node programs in order to accomplish global data manipulations presents something of an obstacle to the reader trying to trace the flow of information and control. 
The positive and negative aspects of this example serve as motivations for the developments to follow in subsequent chapters. In particular, we will expand upon the themes of "natural" concurrent decomposition while introducing programming extensions (in Chap. 6) which remove much of the clutter caused by explicit control process com- munication tasks. 

104 

Introduction to Decomposition and Concurrent Algoritluns: Wave Equation Chap. 5 

Code 5-3 Updating Subroutine for Wave Equation Solver (Node Program) 

5.4 Concurrent Efficiency in the Wave Equation Example 
We have already given a qualitative discussion of the factors influencing the con- current efficiency of the vibrating string calculation under a linear decomposition. In this section, we will quantify the treatment somewhat, in order to give a feeling for the kinds of performance to be expected in practice. 
We should recall the definition of efficiency within the context of the current prob- lem. Suppose that T 1 (M ) is the machine time required to perform an update cycle over M nodal points within a single (sequential) processor. This serves as our baseline for comparison with the concurrently computed case. We shall denote by TN(M) the time required to update the same M points within a concurrent ensemble of N processors. In this case the efficiency £ is given by: 

E= NTN(M) ,- .~ 
As mentioned earlier, mere are distinct factors mat cause TN(M) to measurably exceed N-IT I(M) in an actual calculation. It is possible, however, to lump togemer the distinct sources of inefficiency into me net concurrent overhead,/which we encountered in Chap. 3. The overhead is defined such mat 
/ = 1--1 = NTN(M)-T I(M) (~-~) E 

T l(M) 

T l(M) 

Sec.5-4 Concurrent Efficiency in the Wave Equation Example 

105 

The concept of overhead is useful because components of the overhead arising from different effects simply combine additively to determine the overall efficiency. Thus, we can separately analyze the two important influences on efficiency in our partic- ular example. 
We first consider the inefficiency due to communication overhead. In order to quantify this effect, we must have values for tcalc and tcomm , the calculation and commun- ication times of a single member of our program. 
Although the actual magnitudes of these times may be of interest, for the discussion of efficiency we shall be principally interested in their ratio, tcomm/tcalc. As discussed in Secs. 2-9 and 4-7, the value of tcomm/tcalc for actual machines depends upon issues of hardware-resident software as well as message size. Extrapolation of current trends in concurrent computer development suggests that tcomm/tcalc values near unity will be typi- cal in the near future. In our illustrative calculations in the present chapter, we will con- sider the two values of 0.1 and 10, in order to bracket the probable range of actual machine characteristics. 
The contribution of communication overhead to concurrent inefficiency is simply expressed: 

fc =[~J~ (5-6) M tcalc 
This expression assumes that the decomposition places equal numbers of points in each processor (M IN is an integer) and that every processor contains two "communi- cated" points. Thus for every M IN update calculations, there occur two communication operations within a given node. 
The treatment of cases for which the processors have unequal loads requires modification of Eq. (5-6). The example implementation of the last section yields a load imbalance if M IN is noninteger. In the imbalanced case, there exists a one point differ- ence between the loads of the processors. In this case, we add an overhead term equal to N IM , which is the fraction of ideal run time spent waiting for the overloaded processor(s) to catch up. The net overhead now consists of communication and load bal- ance components: 

[2N~tcomm f=fc+fL= ---+ M tcalc 
Here, me quantity within me ceiling function rl in me latter tenn is zero when M IN is an integer and unity omerwise. 
Now me concurrent efficiency is obtained directly from Eq. (5- 7) via me definition: 
E = (I + f )-1 (5-8) 
The result of evaluating mese expressions for me two example values of tcommltcalc is shown in me plot of Fig. 5-3. It is apparent that wim increasing grain size (n = M IN), me load-balanced efficiency asymptotes to unity. Furthennore, me "sawtoom" effect caused by load imbalance when M is not evenly divisible by N (8 in mis case) 

N M 

MmodN N 

106 

Introduction to Decomposition and Concurrent Algorithms: Wave Equation Chap. 5 

diminishes as grain size grows. The upper envelope of each sawtooth curve corresponds to the case of perfect load balance, while the lower envelope delineates the imperfectly balanced case. 

1.0 

0.8 

c 
6~ 0. a c: 
8 0.4 

C4.) 

0.2 

0.0 

~ 

= 0.1 

- 

8 PROCESSORS (N=8) 

0 

200 

n 

400 600 number of updated points 

800 

FIgure 5.3 Plots of concurrent efficiency as a function of total problem sim, for an 8- processor ensemble. Shown are cases for two values of the ratio t It 
/c. 
comm ca Figure 5-4 presents the same calculated perfonnance trends as Fig. 5-3, but in a manner which makes clearer the functional relationship between efficiency and grain size. As stated above, the quantity E-l -1 represents the net communication/imbalance overhead fraction. It may be seen from Eq.(5-6) that for a fixed value of tcommltca/c , this overhead is inversely proportional to the grain size, (n = M IN). As a result, a plot of E-l -1 versus M-1 should fall between envelopes that are straight lines passing through the origin. The slopes of the envelopes marking best- and worst-case load balance are detennined by the magnitudes of tcommltca/c and N . 
Superimposed on the calculated plots of Fig. 5-4 are measured values for the wave equation problem solved on eight processors of the Caltech Mark II hypercube. It is apparent from these results that our efficiency calculations are valid, with a tcommltca/c ratio of approximately 3.7. 

5-5 Summary 

The example put forward in this chapter is quite simple in terms of its root equa- tions and implem~entation. It does, however, serve to illustrate many of the characteris- tics common to applications written for the ensemble/control process architecture. The homogeneous problem exemplified by our discrete update algorithm serves to highlight the importance of the factors of communication overhead, topology, and load balance in determining performance. From this point we may move onto the consideration of prob- lems that are more complex from one or more of these standpoints. 

Sec. 5-5 Summary 

107 

n 

reciprocal problem size 

F1gure 5-4 Plots of inefficiency overhead as a function of reciprocal problem size. Machine parameters are as in Fig. 5-3. Points represent actual timings of calculations on eight nodes of the Mark II hypercube. The measured values are for load balanced configurations which correspond to the lower envel~ of the sawtooth curve. 
We emphasize that the domain decomposition used in this chapter and discussed in Chap. 1 and 3 is very simple but general. A preview of the range of applications for this method to be discussed in subsequent chapters is illustrated in Fig. 5-5. 
Chapters 7 through 13 will deal with physically motivated computational problems that introduce the issues of higher dimensional decomposition topologies as well as more diverse solution methods. In later chapters, further problems associated with load balancing and concurrent algorithms are pursued. 
The next chapter, however, will deal with another consideration in practical imple- mentation. We have seen that while the structure of the ensemble program in our exam- ple follows in a straightforward manner from its sequential equivalent, the overall clarity is obscured somewhat by the necessity for two enmeshed application programs to handle the control and node processes respectively. The following chapter introduces an aug- mentation of the programming environment which precludes the necessity for rewriting a control program for each new application. This environment, built up from the already described operating system primitives, provides a clearer programming style for most problems. At the cost of imposing a slightly thicker layer of "management" software, we will realize a relatively great gain in algorithmic clarity and transportability. 

108 

Introduction to Decomposition and Concurrent Algorithms: Wave Equation Chap. 5 

(b) 

F1gure 5.5 (a) a simplified game tree for computer chess 00 a four-processor machine. Realistic trees have larger branching factors (- 35), (b) a typical finite difference mesh with d<Xnain decompositioo for four processors (Otap. 7), (c) spatial dec<Xnposition (16 processors) and load imbalance illustrated for an irregular particle dynamics problem (Otap. 16). 

6 

6-1 Introduction 
In Chap. 5, we used the example of a vibrating string to illustrate some of the gen- era! concepts involved in constructing a concurrent application. When designing a paral- lel program to treat a specific problem, we encountered the necessity to write two pro- grams. The first of these, running in the control processor, is the control program, while the other, running in the nodes of the parallel machine, is the node program. The con- trol program manages the program's interaction with the outside world and coordinates the global actions of the parallel program. This style of programming is a natural one to apply to distributed memory, large node, message-passing parallel computers because they frequently provide no direct facility for parallel I/O, except indirectly through the control processor. This I/O isolation of the node processes represents a serious impedi- ment to robust programming. Application programs in this model are considerably longer and more error-prone than would be the case if direct I/O involving node processes were possible. A related difficulty arises from the fact that the application con- sists of two tightly coupled programs, which must communicate in a rigidly coordinated manner. In the indirect I/O model, unanticipated I/O activity on the part of any or all processes most often results in the complete breakdown of the program. Maintenance, modification, and debugging of programs is clearly very difficult in such an environment. 
In this chapter, we present an alternative to the control program/node program methodology, in which the control program become$ an unchanging part of the system, and the only thing that changes from one applicatIon to the next is the node program. This approach becomes possible if we provide a library of routines for node programs which permits each application node process to engage in I/O activities analogous to those associated with an ordinary sequential process. 

6.2 Introductory Input/Output Formatting 
The purpose of this chapter is to introduce the general principles and framework of a concurrent 110 subroutine library .This generic approach to the concurrent 1/0 problem was first implemented in the C language on the Caltech Mark II hypercube, under the name of CUBIX. A detailed discussion of CUBIX within the CrOS III Crystalline Operating System, will be given in Chap. 15. For the present, however, we will discuss 

I11 

112 

An Approach to the Concurrent Input/Output Problem' 

Chap. 6 

CUBIX in only the degree of detail necessary to illustrate its basic utility. A problem arises here in attempting to present this environment in a language- independent way. Whereas the C language handles 1/0 through subroutine calls, FOR- mAN incorporates 1/0 tasks and their associated formatting functions as parts of the language definition. This means that if we propose a subroutine-based CUBIX system, it will give the appearance of being biased in favor of C over FORTRAN. While the C programmer can simply make use of CUBIX versions of the familiar print/ and scanf functions, the WRIlE and READ statements familiar to the FORTRAN programmer are useless without language/compiler extensions. Of course, the FORTRAN programmer may simply adopt the subroutine-based 1/0 approach with CUBIX, but it will be at the cost of sacrificing some of the language independence of our concurrent programming approach. 
In this chapter, we will adopt the perhaps inelegant expedient of casting all CUBIX 1/0 operations in a C-like syntax. While this choice will be comforting to the ingrained C programmer, the FORTRAN devotee may feel quite abandoned. For this reason, we devote this section to the definition of CUBIX 1/0 formatting conventions. In most par- ticulars, the conventions follow the standard C library definitions [ANSI 87]. As described in Appendix C, we will apply these ideas to FORmAN as extensions to the conventional WRIlE and READ statements. 
Formatted data is read from the standard input device, usually a terminal, through calls of the scan/(format, argl, ...) function. In the simplest usage, the first argument is a character string specifying the format of the data to be read. Among the valid possible contents of the format string are: 
the substring "%d", which indicates that an integer should be read, and the substring "%r', which indicates that a floating-point number should be read. A number of such tokens may be concatenated within a single format string. The format string is scanned from left to right. As the data are read and interpreted, they are placed at the addresses specified by subsequent arguments argl, arg2, Since these arguments must be addresses, it is common practice in the C language to use the address-of operator, ' &,' prepended to a variable name. 
Formatted output is produced with the printflformat, argl, ...) function. The first argument is also a character string in which the same two character tokens may be used to identify integer and floating-point numbers in the argument list The format string can also contain other characters, which are printed without modification in the output stream. The string "\I}" may be used to insert ~ new-line character in the output at the corresponding point 
Finally, we will need to refer to specific output streams. The name stdin can be used to refer to the standard input stream, i.e., the place from which scant reads, and the name stdout can be used to refer to the standard output stream, i.e., the place to which 
print/writes. 
Code 6-1 illustrates the basic use of scanf and print/ in a C program. 

Sec. 6-3 The CUBIX Library Functions and Modes 

113 

Code 6-1 Example of C-Language Formatted I/O 
The following tmnscript, in which italics indicate user input, shows the result of an invocation of Code 6-1 : 
Enter an integer and a floating-point number. 423.14 
You entered 42 and 3.140000 Bye. 

The analogous FORTRAN program would follow directly from this example. For- mat strings would be passed via character variables, and because FORTRAN arguments are passed by address, variable names would replace the addresses of the C example. 

6-3 The CUBIX Library Functions and Modes 
The aim of the CUBIX system is to provide a reasonably general concurrent 1/0 environment which in its "sequential limit" reverts to an appearance and function as close as possible to conventional sequential 1/0 functionality. Code 6-1 illustrates a set of gen- eric 1/0 tasks whose meaning and aim is clear for a single-processor machine. What then is the meaning of the program in Code 6-1 when executed under CUBIX in the proces- sors of a concurrent ensemble? The answer to this question depends upon the program- ming context, and is tied up with the new concept of concurrent 1/0 modes . 
We can imagine a context in which all processors of the ensemble have identical values for the variables i and x to send to an output stream. In this case, the sensible action for the system to take would be to print a single output stream corresponding to the identical outputs of the node processes. This mode of concurrent 1/0 handling will be referred to as the singular mode, and represents the default mode in CUBIX. 
The use of singular mode 1/0 is subject to a few simple but important restrictions. Firstly, it should be clear that each node of the ensemble must concurrently execute identical 1/0 operations. A second restriction has to do with the meaning of concurrent 1/0 operations. This restriction follows by close analogy with the restrictions imposed upon crystalline communications in Chap. 4. In fact, the similarity stems directly from the fact that CUBIX routines are built up from the crystalline communication primitives. As with those primitives, the CUBIX routines must be called in loose synchronization. 

114 

An Approach to the Concurrent Input/Output Problem: ...Chap. 6 

While the calls need not occur perfectly simultaneously in real time, the node processes cannot proceed beyond the CUBIX subroutine call until all processors have arrived at the rendezvous point. Since communication takes place within the CUBIX subroutine, no node may have read or write operations pending at the time of the call. 
Concurrent 1/0 would clearly be rather trivial and uninteresting if the singular mode were the only option. We instead recognize the need to exchange distinct data with different processors of the ensemble, and introduce the multiple 1/0 mode to this end. 
When an input or output stream is operating in multiple mode, the individual processes may invoke CUBIX library routines with different arguments or data. The dis- tinct data are in this case routed sequentially to or from the 1/0 device in order of increas- ing processor number. The CUBIX system provides for the switching of a stream between singular and multiple modes through a pair of loosely synchronous subroutine calls. 

An I/O sb"eam can be switched between multiple and singular modes with the func- tions fmulti and fsingl. These functions require a single argument which identifies the sb"eam to be switched, e.g., stdin or stdout. Code 6-2 illustrates the use of multiple and singular I/O as a means of manipulating "common" and "private" processor data. The structure ce, the file cros.h and the function cparam are part of the crOS III communica- tion system. They will be discussed further in Chap. 14. They are used here so that the processor can identify itself. 

Code 6-2 Program illustrating CUBIX Singular and Multiple Modes 

Sec. 6-3 The CUBIX Library Functions and Modes 

115 

An interactive session involving the program of Code 6-2 on a concurrent proces- sor with eight nodes would have the following pattern: 

CUBIX is now in singular mode. Please enter a number. 42 
CUBIX now switching to multiple mode. Please enter 8 numbers. 235711 131719 
You gave 42 and 2 to processor 0. 
This unique message came from processor zero. You gave 42 and 3 to processor 1. You gave 42 and 5 to processor 2. You gave 42 and 7 to processor 3. 
You gave 42 and 11 to processor 4. You gave 42 and 13 to processor 5. You gave 42 and 17 to processor 6. You gave 42 and 19 to processor 7. 
We can see from this simple example how in multiple mode some of the restric- tions associated with singular mode are relaxed. As already noted, there is no longer the need for all processors to issue identical calls. Furthermore, in multiple mode there is no need for processors to produce output in loose synchronization. This is because output is buffered separately in the nodes. Writing into these distinct processor buffers with printf may proceed completely asynchronously. Flushing of the buffers, by calling ffiush , Idose or exit, however, must be accomplished with loosely synchronous calls. In multi- pie output mode, buffers are never flushed automatically, so the application program must explicitly call Iflush in order for the user to see the output. When flushed, the con- tents of the buffers belonging to each node appear in the output stream in order of increasing processor number. 
To demonstrate the advantage of using the CUBIX library, we will consider a modification of the wave equation solver of Chap. 5, using the CUBIX system. This time there is no need for a control program. Code 6-3 outlines in pseudocode all of the user- written node programming needed in order to accomplish the same tasks as the Chap. 5 example. Appendix D contains complete implementations of this same example problem in both C and FORTRAN. 

116 

An Approach to the Concurrent Input/Output Problem' 

Chap. 6 

Code 6-3 Main CUBIX Program for Wave Equation Solver 
The user-written code required to implement the CUBIX example of Code 6-3 is somewhat more lengthy than the corresponding element program of the previous chapter. It is significant, however, that the program is very much shorter than the combined ele- ment and control programs of Chap. 5. Furthermore, the clarity and ease of modifying the program is greatly enhanced using CUBIX, since there is no longer the need for the user to explicitly coordinate the control and element processes. 
A further advantage of the CUBIX approach stems from its previously mentioned compatibility with standard sequential 1/0 operations. The programs presented in Appendix D could be compiled and run on a conventional computer without modification. The specialized concurrent library routines such as shift, gridinit. jmulti, and so on need only be implemented in such a manner as to give the correct (or null) results in the limit of a single processor. This backward compatibility of CUBIX-style applications with sequential computers offers many advantages in the task of program development Initial development and debugging of the nonconcurrent body of a pro- gram is invariably easier within a mature sequential environment, and the applicability of a single programming style to a variety of machine types makes it all the more attractive as a general purpose approach. 
The applications to be developed in the remainder of this book will make use of the CUBIX style of concurrent programming for the reasons cited above. The detailed implementation of CUBIX within the crOS III crystalline communication system is deferred until Chap. 15. While it is good to retain an awareness of the importance (and limitations) of the control process as an intermediary for communication with the nodes of our concurrent computer, we will find that concealing this dependence as much as pos- sible proves to be a fruitful approach to a wide variety of computational problems. 

7 

Elliptic Problems in Two Dimensions I: Solution of Laplace's Equation by Finite Difference , 

7-1 Introduction 
In this chapter we extend our treatment of concurrent computational problems to the case of simple elliptic partial differential equations. The examples discussed in previ- ous chapters generally involve the use of concurrency in sufficiently natural ways that the choice of decomposition strategy is fairly obvious. In the current discussion, we will explore the concurrent implementation of algorithms with a somewhat less clean-cut parallel structure, while also considering higher-qimensional topologies for decomposi- tion. 

An attempt is made here to use relatively simple illustrative methods, and to avoid the complications of more sophisticated numerical approaches. As elsewhere in this volume, the general principles of concurrent decomposition and algorithm structure sketched out below should serve to guide the further development of more generally use- ful and powerful programs. Insofar as the main issues of this chapter have little to do with the basic control process tasks discussed in the earlier examples, the examples of this and succeeding chapters are presented in the form of CUBIX applications. This choice of example structure improves the clarity of the concurrent algorithms. 
For the purpose of this discussion we shall consider techniques for solving the two-dimensional Laplace equation: 

V2(j)=~+~=O (7-1) OX2 Oy2 
While we might as easily tteat other elliptic operators or three-dimensional cases, the present problem is sufficiently challenging to demonstrate a number of salient points regarding concurrent decomposition. 
We seek solutions for the unknown potential, (j), within some finite spatial domain, given fixed boundary solution values. In pursuing an approximate numerical solution to (7-1), one common approach is to consider a spatial discretization of the problem domain, so that the problem becomes one of finding values of (j) at a finite number of nodal locations. Note that here, as in Chap. 5, the term "nodal point" refers to the domain discretization, and should not be confused with the processor node of the concurrent ensemble. This approach is an extension of the one-dimensional spatial discretization employed in Chap. 5 to solve the wave equation. In the one-dimensional hyperbolic 

117 

118 

Elliptic Problems in Two Dimensions I: Chap. 7 

case, it was a straightforward matter to algebraically solve the finite difference expres- sion to obtain a simple update prescription for the solution. In the case of an elliptic problem such as (7-1), the simultaneous solution of the difference equations connecting the discrete cl>i leads formally to a matrix problem. The difficulty of solving this matrix system is closely tied to the structure and sparsity of the matrix. 
While a one-dimensional elliptic problem would give rise to a simple tridiagonal matrix, higher-dimensional operators yield more complicated matrix structures. We might attack any of these matrices with direct methods derived from Gaussian elimina- tion. In this and in the next chapter, however, we will emphasize iterative solution methods. This choice is made for both pedagogic and practical reasons. Iterative methods are preferable for the illustration in this chapter because they involve concurrent decompositions with a comparatively obvious correspondence to the physical dimen- sionality of the problem domain. In contrast, direct matrix methods generally require a decomposition based upon the two-dimensional row/column world of the matrix. These row/column-oriented concurrent matrix operations will be covered in Chaps. 10 and 20. On the practical side, iterative methods assume importance in large higher-dimensional problems, where the cost of direct techniques increases with problem size much faster than the cost (per iteration) of iterative methods. It is likely that iterative techniques will predominate for the large three-dimensional problems to which concurrent processing is increasingly being applied 
Both this and the following chapter deal with methods of solving the matrix sys- tems arising from discretized elliptic problems. At the level of the matrix problem, the difference between this chapter (finite difference) and the next chapter (finite elements) lies in the regularity of the matrix structure. 
Because the method of finite differences yields matrices with a more predictable structure than the finite element method, this chapter will deal with iterative solvers that do not require the full system matrix to ever be assembled. This will contrast with the finite element equation solver discussed in Chap. 8, which assumes less about the matrix structure, but might be applied to finite difference systems as well. 

7.2 The Finite Difference Formulation 
In the simplest application of the finite difference method, the problem domain is covered with a uniformly-spaced square mesh of nodal points, so that all points, except those on boundaries, have four equidistant nearest neighbors. We may designate a node's neighbors by the subscripts :tx and :ty. The first-order finite difference equation corresponding to Laplace's equation (7-1) is: 
4 <l>i -<l>i-x -<l>i-y -<l>i+x -<l>i+y = O (7-2) 
where i is an index running over all M nodal points. This leads to M simultaneous equa- tions for the <l>i' and a matrix problem characterized by a symmetric, diagonally dominant sparse system. 
Perhaps the simplest iterative technique applicable to this system is Jacobi iteration, also known as the method of simultaneous corrections. In this method, the trial value cl>i(k) 

Sec. 7-2 The Finite Difference Fornlulation 

119 

at the k th iteration is obtained by solving the i th equation for cl>i, using the values from the previous iteration at the other nodes. Inspection of (7-2) shows that this corresponds intuitively to replacing the unknown value cl>i at each node with the average of its neigh- bars' values from the previous iteration. The Jacobi update procedure for step k may be written: 

(;.-1) + "'. (;.-1) + "'. (;.-1) + "'. (;.-1) 
]'fl-y 'f1+X 'fl+y 

The Jacobi iteration scheme represents one of the class of stationary iterative methods [Jennings 77], the members of which are distinguished by the correction tech- nique applied at each update. In the case of Jacobi iteration, trial solution values are updated simultaneously. In another stationary method known as Gauss-Seidel iteration, the trial values <l>i are updated in sequence. For example, the value <l>i (k) is obtained from the i-l newly obtained values for step k and the remaining M-i "old" values from step k-l. 

It should be clear that when considering concurrent implementations of these methods, the Jacobi scheme is the most straightforward because it imposes no sequential order on corrections to nodal values. Therefore, neglecting communication steps, the Jacobi algorithm can be carried out completely in parallel without modification. This contrasts with the Gauss-Seidel case, which assumes a particular update sequence, which if observed in the concurrent case would require processors to work on regions of the domain sequentially instead of in parallel. 
Unfortunately, the Jacobi method is of much less practical importance than the 80- called nonsimultaneous methods because of its relatively slow convergence. At this point we might go forward with development of the concurrent Jacobi scheme despite its comparatively poor performance. However, this would not serve to illustrate much that has not already been discussed, since the update method simply amounts to a two- dimensional extension of the simultaneous updating used to evolve the wave equation forward in time in Chap. 5. Instead, we shall discuss the problem of concurrent succes- sive over-relaxation, of which the Gauss-Seidel method is a special case. 
In the case of successive over-relaxation (SOR) we replace the update prescription of (7-3) with: 

In the above formula, 0> is the over-relaxation parameter, which generally lies in the range 1<00<2 for optimum convergence. In the case that 0}=1, the SOR method is identi- cal with Gauss~Seidel [Jennings 77]. 
It should be noted that Eq. (7-4) has been written assuming that the cl>i (k) are updated in order of increasing x and y coordinates. In fact, for a sequential algorithm, any particular update order might be used, and newly updated values could be simply written over their predecessor values in memory. While alternative update orders might affect the convergence rate, there would be no fundamental difference in the necessary programming aside from indexing. When we turn to a concurrent implementation of 

120 

Elliptic Problems in Two Dimensions I: Chap. 7 

SOR however, update order becomes a more important consideration. In order for two (or more) subsets of points to be updated concurrently, neither subset may have an influence on the result of the update step in the other. More formally, this means that any two points may be updated concurrently only if their mutual matrix element is zero. In the context of our first-order finite difference problem, it means that no two regions with adjacent edges may be updated in parallel. This restriction forces us to adopt the modified update order described in the following section. It is interesting to note that if our coefficient matrix were full (i.e., all nodal points interacting) there would be no possi- bility of a concurrent decomposition based upon spatial regions, and an alternative "Iong-range" decomposition like that of Chap. 9 would be required. 

7-3 Concurrent Decomposition 
The task of obtaining a concurrent decomposition for our finite difference/SOR problem centers about two related considerations: domain partitioning and update order- ing. The first of these considerations determines to some extent the second, so we shall consider it first 
The finite difference formulation selected for this chapter restricts us to considera- tion of uniformly-spaced meshes of nodal points. While we could easily imagine the use of such a locally uniform mesh to approximately cover a nonrectangular domain, it will be convenient for the present example to consider the simple square domain illustrated in Fig.7-l. 
As was the case in Chap. 5, intuition provides an adequate means of selecting a good domain decomposition for multiple processors. In the two-dimensional case, we are required to select processor regions which balance the workload and minimize the communication overhead. Since communication occurs along interprocessor boundaries, the latter requirement seeks to minimize the perimeters of the subdomain boundaries. A convenient approximation to the best decomposition is obtained by breaking the square domain into congruent square subdomains. While intuition serves well here, the optimal solution for more convoluted meshes could be much less obvious. In such cases, we may choose to employ optimization methods of the kind discussed in Chaps. 13 and 17 in order to obtain a decomposition. 
In Fig. 7-1, one particular processor region of a 16-node decomposition has been highlighted for the purpose of illustrating the flow of data within and between regions. A grid point which does not lie within the shaded boundary strip of the region requires only data "local" to its own processor node in order to be updated. In contrast, the shaded edge points require one or two values from neighboring processors in order to be updated. Thus these shaded points constitute the two-dimensional analog of the com- municated end points given in the wave equation example of Chap. 5. In three dimen- sions, these communicated points would occupy the six faces and twelve edges of a parallelepiped. 
If we were doing a simple simultaneous update scheme (Jacobi iteration), the algo- rithm would simply consist of steps that concurrently update the processor subregions, separated by communication steps that involve the swapping of edge values between 
~ 

7-3 Concurrent Decomposition 

)::.:~~;~~: .~~ 11;;:~:: 

'-- processor region boundory 

... 

I "" 
, " I. ~ 

... ... ... ... ... .. .. .. .. 

... r:lill ;I;!... 
Ilill1lii ~ ~ : 

... 

... 

~ 

... 

Figure 7.1 D<Xnain decomposition for two-dimensional finite difference scheme. Nodal points within shaded region require informatioo fr<Xn neighboring processors. 
neighboring processors. An additional complication is introduced, however, by the con- sideration of consistent update ordering for nonsimultaneous corrections. 
Given the rectangular decomposition shown in Fig. 7-1, the requirement of nonad- jacent concurrent updating may be satisfied by invoking the so-ca1led redlblack or check- erboard ordering scheme of Fig. 7 -2 (see for example [Ortega 85]). Each square (or rec- tangular) subregion is in turn subdivided into four regions of two kinds. For the sake of definiteness, we refer to these partitions as "red" and "black," by analogy with the squares of a checkerboard. 
Inspection of Fig. 7-2 shows that no red or black partition of any subregion is adja- cent to a partition of the same color in a neighboring processor. Thus, a fully parallel SOR update may be realized by first updating all red points, then all black points, and so on. A sufficient condition to guarantee the consistency of this update approach is to make certain that all processors have completed one color update cycle before any start on the next color. . 
Strictly speaking, this ordering is different from that of standard sequential SOR, in which the update sweep proceeds unidirectionally across the global grid. The definition of SOR, however, only requires that there is some update order, but cannot say what that 

122 

Elliptic Problems in Two Dimensions I: Chap. 7 

processor boundaries 
.-""/~ 

"red" 

"black" 

succesive update in non-adjoining red regions 

FIgure 7.2 "Red/black" checkeJboard domain partitioning required for concun-ent imple- mentation of the SOR algorithm. Points in red regions are updated first, and following processor synchronization, black points are updated. This partitioo allows sequential SOR cycles to run concurrently in the processor nodes. 
order is, since it depends in detail on the specific problem. The reader may object at this point that since the concurrent algorithm presented is no longer strictly equivalent to the original sequential scheme, we have little justification for calling it true SOR. Indeed, from a practical standpoint, the programmer may wish to dispense altogether with the bother of maintaining a consistent update sequence and simply perform parallel SOR in the subregions without regard for the updates occurring in neighboring processors. This approach would amount to a kind of hybrid between the Jacobi and SOR methods, and would scarcely introduce a greater corruption of the "pure" SOR scheme than we have already inflicted. We will stay, however, with the red/black SOR approach in order to illustrate the issues involved in adapting an inherently sequential update for concurrent application. 
In the next section, we will detail the implementation of this example in our virtual concurrent processor environment. The main structure of the parallel algorithm we will 

Sec.7-4 Example Implementation 

123 

implement is briefly summarized by the following steps: 
1. exchange edge values with neighboring processors. 
2. perform SOR update on all red points in sequence to obtain cl>i (k). 3. exchange edge values with neighboring processors. 
4. perform SOR update on all black points in sequence to obtain cl>i (k). 5. advance stepk--+k+l. 
6. if finished, stop. Otherwise, go to 1. 
As was the case in the example of Chap. 5, the algorithm involves a highly predict- able and homogeneous allocation of computation and communication tasks, so that load balancing is straightforward. Also, the individual processors, while not running strictly in lockstep, are kept in the loose synchronization required for update consistency by the communication calls interspersed between update cycles. 

7-4 Example Implementation 
We have mus far presented me concurrent finite difference algorimm in a some- what abstract manner. In mis section, we will attempt to clarify me implementation of this scheme wim a pseudocode example. As elsewhere in me book, mis example of pseudocode correlates wim real examples in me Software Supplement. 
We will present selected portions of me node programs designed to implement me algorimm in me CUBIX environment (cf. Chap. 6). For me sake of brevity and clarity, not all of me subroutines called are fully described, particularly when meir function is tangential to me principal algorimmic issues. 
The main module of mis example program is given in Code 7-1. The main pro- cedure begins by using grid-mapping routines (Chap. 4) to establish communication pams wim up to four nearest neighbors in a two-dimensional processor mesh. Next, a loop is entered which accepts user commands via me CUBIX singular input functions. 
Within me command loop, me appropriate option is executed through a subroutine call. In mis very rudimentary example, me only options provided are mose required to read in me starting data, to perform a specified number of updates, and to dump me result via CUBIX multiple I/O routines. Clearly, any really useful program of this kind would incorporate more features, such as mose necessary to check for convergence. 
The initialization and dumping portions of me program are not central to me present discussion, and have merefore not been described in detail. Their principal func- tions are to manage me allocation of domain points and fixed boundary points to subre- gions, and to set up me necessary bookkeeping to identify neighboring values for each relaxed point. 
Of more immediate interest is me implementation of me update cycle, which con- stitutes me real heart of the program. Code 7-2 contains me pseudocode outline of mis update cycle. As described in me previous section, each full cycle consists of a commun- ication and computation phase for each of two colors of points in a processor. Strictly speaking, during me "red" phase of a cycle, edge data need only be received from neigh- boring black points, and need only be sent to neighboring red points. The reverse is true during me "black" phase. In me illustrated implementation, we have skirted me issue mat 

124 

Elliptic Problems in Two Dimensions I: Chap. 7 

Code 7.! Main Driver Module for Finite Difference Example 
only half of the edge values need be communicated during any call to get-edge. An intelligent implementation of get-edge should be designed to recognize whether the current phase is red or black, and communicate only the needed values. 
Glossing over this detail, the function get-edge in Code 7-2 illustrates a method for systematically exchanging edge values among the processor nodes. Each node program requires storage space for the values at each of its own domain points, plus a single layer or "guard ring" of points needed to update the interior values. These slots are to be filled with values supplied by interprocessor communication. 
In our rectangular decomposition geometry , the exchange is accomplished by four successive calls to vm-shift. First, a processor's values falling along the right edge are copied into a communication buffer and sent to the right neighbor while complementary left edge data arrives from the left neighbor. These data are stored in the appropriate guard ring locations, and the "right-shift" pass is ended. This is followed by analogous "left-," "up-," and "down-shift" passes, after which the guard rings of all processors have been filled with up-to-date values. We have made use here of the ability of vm-shift to transparently handle nonexistent neighbors in the case of processors with fewer than four two-dimensional neighbors. 
Following the completion of the get-edge communication step, each processor independently performs a sequential SOR update on its red points, according to the prescription of Eq. (7-4). In this sequential process, each nodal value is immediately overwritten by its updated value. The communication step is then repeated, followed by an update of all black points. 
This example serves to illustrate the relative simplicity with which a concurrent implementation of a finite difference solver may be designed. A more generally useful concurrent SOR program could be expected to follow the general pattern outlined above, 

Sec.7-4 Example Implementation 

125 

although the adaptation of this decomposition approach to irregularly shaped domains or higher-order differencing schemes may present significant additional programming chal- lenges. 

proc-begin update-SOR ( perfonn k update cycles ) for -begin (a total of k cycles) 
proc-caU get-edge (exchange edge values with neighbors ) for-begin (each red point, i) 
[locate four nearest-neighbor values for Jl 
[replace potential i with updated value (cf. Eq. 7-4) ] for end 
proc-caU get-edge (exchange edge values with neighbors) for-begin (each black point, i) 
[locate four nearest-neighbor values for Jl [replace potential i with updated value ] for end for-end proc-end 

proc-begin get-edge (exchange edge values with neighbors) [load edge values into right-edge array] 
[ vm-shift -send right-edge to right neighbor while receiving from left neighbor into left -edge ] [ store left -edge values for use in update] 

[load edge values into left -edge array] 
[ vm-shift -send left -edge to left neighbor while receiving from right neighbor into right-edge] [store right-edge values for use in update] 

[load OOge values into up -edge array] 
[ vm-shift -send up -edge to upper neighbor while receiving from lower neighbor into down-edge] [store down-edge values for use in update] 

[load edge values into down-edge array] 
[ vm-shift -send down-edge to lower neighbor while receiving from upper neighbor into up -edge] [store up -edge values for use in update] proc-end 

Code 7-2 SOR Update Module for Finite Difference Example 

126 

Elliptic Problems in Two Dimensions I: Chap. 7 

7-5 Performance in Practical Applications 
In designing and evaluating any practical application similar to the example just presented, a number of competing computational cost factors will determine the overall concurrent efficiency of the implementation on a given machine. Among the factors influencing performance are the order of the finite difference approximation (as discussed in the next section) and the dimensionality of the mesh. . 
While the relative performance relations may be derived theoretically, as was done in Chap. 5, it is useful at this point to consider actual performance figures that emerge from implementing a two-dimensional finite difference solver on a real machine. We present below the results obtained from a first-order finite difference program run on the Caltech Mark II hypercube in various configurations. The specific example cited here is described further in Sec. 7-7. Unlike our illustrative example, it employs a simultaneous spatial update scheme. However, its similar communication structure makes it a useful case study in program performance. Timing results were obtained as problem size and number of processors were varied, allowing calculation of experimental efficiencies. 
We can easily predict the load-balanced efficiency dependence on problem size from considerations discussed in the previous chapters. The ratio of computation to interprocessor communication for a given update step is very nearly proportional to the ratio of total number of nodal points per processor, n, to the number of points per proces- sor involved in communication. In other words, this ratio varies, in the two-dimensional case, ~ the ratio of the area to the perimeter of the processor region, or as n 1/2 for large n [Fox 85c, Reed 87]. 
The n 112 ratio immediately yields the predicted communication overhead 
Ic =cn-1/2 (7-5) which in turn leads to the expected efficiency 
£= (I + IC)-l ::: 1- cn-1/2 (7-6) In the above expressions, c is a constant factor dependent in part upon the machine- 
specific ratio tcomm/tca/c .This constant is of some interest in the present experimental efficiency analysis, since it also depends upon the number of processors in the ensemble, which is variable in the benchmark results presented below. 
In presenting the experimental efficiency results, it is useful to plot the quantity -log(l-£) as a function of log n. According to Eq. (7-6) such a plot should yield a family of parallel lines each with a constant slope of 1/2. 
In Fig. 7-3, the experimental points [Meier 84] are plotted along with the line corresponding to c = I. The experimental points group into families on lines of slope 1/2. The families of timings employing fewer processors fall on systematically higher efficiency trends than those for larger numbers of processors. The systematics of the vari - ations in efficiency may be understood by considering the varying amounts of communi- cation overhead for various numbers of processors at constant n. All ensemble configurations show an asymptotic trend to unit efficiency with increasing n. However, at constant n, different numbers of processors entail different communication overheads.