Sec. 7 -5 Perfonnance in Practical Application 

127 

!..- III I 
0, 0 

log n 

FIgure 7-3 Experimental efficiencies for two-dimensional finite difference scheme with variable number of processors [Meier 84]. Although nm on a hypercube, these results depend only upon a two-dimensional mesh connection. 
When the ensemble consists of two processors, as shown in Fig. 7-4, communica- tion involves only those nodal points along one edge of each processor. A four-node configuration involves twice as much communication at fixed n. With eight processors, the most burdened (and therefore speed-Iimiting) processors have three edges involved in communication. Finally, for configurations of 16 processors and more, the communica- tion load becomes "saturated," with some processors communicating on all four boun- daries. The results of Fig. 7-3 confirm this expectation. 
This example demonstrates quantitatively the magnitude of the penalty imposed by communication overhead in a typical two-dimensional mesh problem. In this particular example all configurations achieve efficiencies in excess of 90% for n ?; 100. 
It is also interesting to consider the significance of the experimentally determined values for c. As shown in Fig. 7-3, the efficiency results for Mark II hypercubes with 16 or more processors are well predicted by Eq. (7-6) with c:::l. This implies that for the machine and algorithm described, the times required to update and to communicate a sin- gle nodal point are roughly equal (to within a factor of order unity). It follows that the results in Fig. 7-3 may be approximately applied to machines with different assumed communication overhead costs through an additive adjustment of magnitude -Iog(tcomm/tcalc ). 

128 

Elliptic Problems in Two Dimensions I: Chap. 7 

DEPENDENCE OF COMMUNICATION ,ON NUMBER OF PROCESSORS ~ = COMMUNICATION 

~ 

2 processors (one side communicoles) 

8 processors ( three sides communicote ) 

16 or more processors (011 four sides communicoteJ 

F1gure 7-4 mustmtion of the dependence of the net amount of c«nmunication on the nDnlber of processors. The limit is reached with 16 processors, after which at least some processors communicate on all four edges. 
7-6 Higher-Order Differencing Schemes 
One obvious generalization of the simple example discussed in the earlier sections of this chapter is to consider a higher-orderfinite difference scheme. Operationally, this means that more than nearest-neighbor nodal points are used to update the potential at a point. An extended update stencil clearly increases the number of arithmetic operations required, whether the update is performed sequentially or concurrently. 
In the case of a concurrent algorithm, the larger numbers of boundary values com- municated between processors also impacts the communication cost. However, the overall concurrent efficiency of the higher-order update scheme may actually be increased at fixed problem size, since the computational operations increase by a larger factor than the number of communicated data items. The actual effect in any specific case must depend upon the ratio of calculation to communication times and the scaling of communication overhead costs with message size. We can quantify these remarks using the communication overhead f c introduced in Chap. 3. The efficiency discussed in the last section, takes the form: 

e= (1 + fC)-l = 1-fc (7-7) 
The basic issues are summarized in Fig. 7-5. Consider first Fig. 7-5(a), which illus- trates the first-order problem described in the opening sections of this chapter. Using n as before to denote the maximum number of nodal points stored in each processor, we 

Sec.7-6 Higher-Order Differencing Schemes 

see that there are 

129 

4n 

and 

floating-point calculations (cf. Eq. (7-3» 

4{ii" words to be communicated 
Here we have considered the general case of 16 or more processors in which com- munication is required on all four edges (cf. Fig. 7-4). This leads to a communication overhead given by 

fc = -r t'-lS) vn tca/c 
Combining Eqs., (7-7) and (7-8), we observe dtat dte Mark II hypercube results of Fig. 7-3 are consistent widt dte result c=1 in Eq. (7-6) if 
tcomm 
-I 

1 tcomm -- 

tcalc 
This result falls in the expected range for machines and applications of this kind, for which it has been noted earlier that tcomm and tcalc can vary by factors of up to two. 
Now consider a higher-order differencing scheme shown in Fig. 7-5(b) for which Eq. (7-3) is replaced by: 

+ 11:"'. (k-l) + 16"'. (k-l) "'1'.+% 'I'.+y 
"' (k-l) "' (k-l) -'l'i-~ -'l'i-2y 
-"'. (k-l) -"'. (k-l) 
]'I'.+2x 'I'.+2y 
Here, we encounter 9n calculations and approximately 4(2.{ii") = 8.{ii" words to be com- municated. Even though communication has increased, f c decreases by a factor of 1.1 from Eq. (7-8). In Fig. 7-5(c), we show what happens when off-diagonal terms are added to the stencil. Now f c is decreased by a factor of 2 cCi>mpared with Eq. (7-8). One way of understanding this decrease is to note that in Fig. 7-5(a) each communicated word is only used in a single update; for Fig. 7-5(b) some are ~sed once and some twice; while in Fig. 7-5(c), each communicated word is typically used in three updates. 
In Fig. 7-5(d) it is shown that as the update stencil increases in size, we obtain a value of f c = c /.{ii" where c steadily decreases as the stencil increases in size. There is an interesting limit when the stencil grows so big that it covers the full domain. Now f c takes the form: 

const tcomm (7-10) 
n tcalc 
Expressed in words, c /-.rn (c decreasing) tends to the limit 1/n as the stencil increases in size. 

fc=- 

130 

Elliptic Problems in Two Dimensions I: Chap. 7 

STENCIL DEPENDENCE OF CALCULATION AND COMMUNICATION illustrated far n= 16 

dl 

INCREASING SIZE OF UPDATE STENCIL 

. 

, 

f -communicof"n I( -1- c -colculofion ,/0 

e) 

--'"" --v 

-1- 2,/ii 

-1- 3.;;; 

-1-- ~ 4.;n x t'O. 

. 

==> 

STENCIL COVERS FULL DOMAN 
fc " 1 !;""" n loo" 

Figure 7-5 mustratioo of the stencil dependence of the communication overllead f c .(a) Shows the simple fonnula corresponding to equatioo (7-2) while (b) through (e) show stencils of increasing complexity. 
The overhead described by Eq. (7-10) is identical to that of the long range force problem discussed in detail in Chapter 9. It is worth noting that one may think of Figs. 7-5(d) and (e) as applying to particle dynamics problems with a force of successively increasing range. 

Sec.7-6 Higher-Order Differencing Schemes 

The discussion of higher-order differencing in this section has dealt only with simultaneous update schemes, rather than the successive relaxation methods illustrated earlier in this chapter. As mentioned above, the extended stencil of higher-order update schemes enlarges the border overlap between processors. In Sec. 7-3, we saw that a suc- cessive update scheme like SOR could be performed concurrently on a decomposed domain if each processor region was subdivided into two colors, so that no two adjacent regions in different processors were updated simultaneously. 
For the extended stencil of Fig. 7-5(b) this analysis still holds, since no communi- cation occurs across processor diagonals. Stencils like Fig. 7-5(c), however, with diago- nal connections, require more than two colors to preserve a consistent successive update scheme. In this case, the four subrectangles of each processor region (Fig. 7-2) are assigned four distinct colors, and four phases are required to update all points in a proces- sor. Figure 7-6 illustrates this modification. 

131 

132 

Elliptic Problems in Two Dimensions I: Chap. 7 

neighboring processors, as well as those displaced only horizontally or vertically. In order to update a comer nodal point. a given processor will need to communicate with diagonally neighboring nodal points, as well as with its Cartesian nearest neighbors. This presents no particular difficulty for our virtual concurrent processor as set forth in Chap. 4, since any processor in the ensemble may be specified as a partner in a communication call. In an actual machine and operating system (e.g., a hypercube), however, the diago- nal next-nearest neighbor connection between processors may not be directly available. In this case, comer nodal points must be treated specially, being communicated once in the vertical direction and once horizontally in order to arrive at the desired destination. While this small number of "special" nodal points should have no serious impact on per- formance, it can increase the complexity of writing the program, since the operations needed to pass a value around the corner must be coordinated among all involved proces- sors. 

The additional programming necessary to handle diagonal passing is actually minimal in the rectangular domain illustrated in this chapter. Since any comer grid point in a subregion which has a diagonal neighbor always has rectilinear neighbors as well, the diagonal communication is achieved "for free" by performing the up/down and left/right communication steps in succession. While this approach is completely robust for domains whose decomposition is topologically equivalent to a rectangular tiling, it is not applicable, in general, to less regular decompositions such as might be dictated by irregularly shaped or multiply connected domains. The issues associated with diagonal passing are discussed further in Chap. 8, in which the finite element method, which virtu- allyalways involves non-nearest neighbor communication, is discussed. 

7.7 Scientific Applications 
Finite difference algorithms of the kind described in this chapter have proven use- ful for solving a variety of differential equations that arise in scientific problems. While many combinations of update procedures, stencils, and equations remain to be explored, some experience has been gained in the solution of systems with nontrivial scientific interest. We present here a couple of examples from the work with which we are fami- liar on the Caltech hypercube. A great deal of other exciting science is being or could be addressed on concurrent computers using the techniques described in this chapter. 
An early application which goes beyond the simple elliptical Laplacian problem is the solution by Meier [Meier 84] of the equations governing two-dimensional hydro- dynamic flow. In the absence of magnetic, gravitational, and radiation fields, the prob- lem becomes one of solving dIe coupled system: 
~+v.(py)=o (continuity) (7-11a) ~+v.(P~)=-vP (momentum) (7-11b) 

oE ot 

+ V-(Ey) = -V.(py) 

(energy) 

(7-11c) 

Sec. 7- 7 Scientific Applications 

133 

Here, p, p, v, and E represent mass density, pressure, velocity, and energy density, respectively. This system has been solved on the Mark II hypercube as a means of simu- lating flow phenomena in astrophysics. 
In a manner similar to the wave equation example of Chap. 5, the spatial domain serves as the basis for a two-dimensional decomposition, with centered spatial differenc- ing used for mesh updating. Temporal differencing again occurs synchronously in all processors, so that the concurrent efficiency is governed by the same considerations which governed the Laplace equation example. 
Another field in the applied sciences to which concurrent finite difference algo- rithms have been profitably applied is that of acoustic and elastic wave propagation. Clayton [Clayton 84] and Frankel and Clayton [Frankel 86] have also employed the Mark II hypercube to solve the acoustic wave equation in inhomogeneous media for seismological investigations. The second-order finite difference update scheme in two dimensions again requires a spatial update stencil which encompasses only Cartesian nearest neighbors. As in the above example, the efficiency analysis for the time- dependent problem carries through from the Laplace case, and the results of Clayton and Frankel have confirmed this expectation. The earlier use of numerical simulations to study seismic propagation and scattering employing sequential computers has been gen- erally limited to relatively small two-dimensional domains. The encouraging results with concurrent processors may portend an extension of such work to larger and finer meshes, as well as to three dimensions. 

8 

Elliptic Problems in Two Dimensions II: 

The Finite Element Method 

8-1 Introduction 
The previous chapter was concerned with the numerical solution of an elliptic boundary value problem through domain discretization. It was this discretization in the sequential problem which formed the basis for a concurrent formulation of the method. In the present chapter we will consider the same basic problem, the solution of Laplace's equation, but make use of an alternative technique, the finite element method. 
To the reader unfamiliar with the distinction between finite elements and finite differences, it may be unclear why both methods should be discussed. In fact, some who are familiar with the subject may raise the same question! Our reasons are twofold, hav- ing to do with the different contexts of applicability of the two methods, and with differ- ences in the issues of concurrency involved. 
In most cases, the finite element formulation of a boundary value problem entails a greater amount of programming and computational effort than a comparable finite differ- ence application. The justification for this greater complexity lies principally in the facil- ity with which the finite element method can be used to handle discrete grids with an irregular or complicated distribution of points. Such domain discretizations, which arise commonly from the need to model "real world" domains or to resolve phenomena over a range of spatial scales, are addressed in a natural and general way by the finite element method. The fact that finite elements are central to many modern scientific and engineer- ing simulations therefore constitutes the first half of our justification for their considera- tion here. 
The second rationale is based upon some of the new algorithmic issues raised by extending our treatment to finite elements. In the previous chapter we obtained a very simple finite difference matrix structure from the regular spacing of nodal points in the discretized domain. This simple structure allowed the implementation of stationary itera- tive solvers without the need for any explicit reference to global matrix indices. 
In the case of finite elements, the analogous matrix of linear equation coefficients lacks the regular,/predictable structure characteristic of the finite difference problem. This circumstance will force us to consider solution schemes that accommodate more general matrices. While the stationary methods discussed in the last chapter may be for- mally generalized to handle finite element matrices, such an implementation would prove to be awkward and inefficient in practice. It is therefore opportune at this point to discuss 

135 

136 

Elliptic Problems in Two Dimensions II: The Finite Element Method Chap. 8 

a different class of iterative solvers, which handle such mabices by design. The so-called gradient methods to be discussed below are proposed as an effective approach to con- current finite element analysis, although it should be noted that nothing prevents their effective use in the finite difference setting of Chap. 7. [Ortega 85] discusses the applica- tion of various iterative methods including the gradient methods to a range of vector and parallel computer architectures. Our specific problem of iterative solution of finite ele- ment systems in a MIMD environment was pioneered in the work of [Adams 82,84] using the NASA ICASE Finite Element Machine. 

8-2 The Finite Element Formulation 
The method of finite elements is the subject of a large number of textbooks and research papers, ranging in their approaches from mathematically theoretical to prag- matic and heavily applied. It will be impractical in this brief discussion to do justice to a description of the method, but we will attempt to assemble the most basic concepts needed to discuss our Laplace equation example. It is hoped that this will provide the reader unfamiliar with finite elements an appreciation of the basic issues as they relate to concurrency, rather than a complete introduction to the subject. 
We begin by considering a spatial domain on which we seek approximate solutions of the original equation 

V2<!1 = O (8-1) 
subject to boundary conditions as in Chap. 7. The method of approaching this problem again calls for the determination of an approximate solution for <!1 at a finite number of discrete points in the domain. Figure 8-1 shows a typical finite element discretization of a two-dimensional domain. 
It is apparent from Fig. 8-1 that the typical finite element grid does not have an evenly spaced or rectilinear structure. As in the finite difference case, we seek to con- struct a system of linear equations whose unknowns <!Ij comprise the desired approximate solution. The following discussion outlines the construction of this system. 
Besides the geometric complexity, the important new concept illustrated in Fig. 8-1 is that of the element, a domain area (or volume) delineated by line segments (or curves) between grid points. Every grid point is a part of at least one element, and in the finite element formulation there exists a nonzero matrix element connecting any two unknowns <!Ij and <!1 j whose points share a common element. From this it may be seen that the task of computing matrix elements may be broken up into the computation of element matrix contributions which may be subsequently assembled into the glabal matrix system, given a knowledge of the geometric arrangement of elements. 
It is the proper computation of the element matrices that is one of the principal topics of finite element treatises, and determines in part the accuracy of the resulting solution. At the heart of this problem is the recasting of the governing differential equa- tion in the form of a variational equation or so-called "weak form." Variational formula- tions of this type are often referred to as virtual work principles, particularly in mechan- ics. 

Sec. 8-2 The Finite Element FOmlulation 

137 

FIgure 8.1 The finite element method involves solving for nodal degrees of freedom whose interactioo is based upon element intelpOlation functions. 
Standard finite element references such as [Hughes 87] or [Bathe 82] describe the application of this variational principle to the discrete problem. Broadly speaking, the space of acceptable solutions is assumed to be adequately represented by an M - dimensional space of basis functions (sometimes called "shape functions") where M is equal to the number of grid points. In the discrete problem (under the Bubnov-Galerkin approximation), we use these basis functions 11; to construct a trial solution. 
The function 11; has the useful property of vanishing at all grid points except the i tIt point, at which it equals unity .This property facilitates the element-wise breaking up of computations related to the matrix elements. In practice, basis functions are usually chosen which have a simple linear or polynomial variation within an element. Our approximate solution '(j) is assumed to be given by the expansion: 
M 
'(j)=L<1>;11; (8-2) ;=1 

where as before, the coefficients cl>i are the approximate nodal values. When the varia- tional principle is applied to this representation of the solution space, a linear system of equations or equivalently, a matrix problem is the result. 
If we use <1> as shorthand notation for the vector of unknowns cl>i , the matrix prob- lem may be expressed in the form 

K <1> = F (8-3) 
In this equation, the matrix of coefficients K is commonly referred to as the "stiffness" matrix, and the right-hand side vector F contains "driving" terms originating in the boun- dary conditions. 

138 

Elliptic Problems in Two Dimensions II: The Finite Element Method Chap. 8 

The form of a general entry in dle stiffness matrix is determined by dle governing variational principle of dle problem. In general, a matrix element Kjj consists of an integral over dle problem domain involving dle basis functions Tlj and Tlj : 
J' {}Jlj OT)j {}Jlj OT)j 
Kjj = !(Tlj,Tlj,a;-' ~' -ay' -ay'...)dx dy (8-4) 
The fact dlat dle shape function Tlj is nonzero only widlin dlose elements of which nodal point i is a part allows us to restrict dle above integral to dle area of a single ele- ment to obtain dle element stiffness contribution kij .In dle specific case of dle Laplacian equation (8-1), dle element stiffness matrix can be shown to have dle form: 
k.e.=J(~~+~~)dxdy (8-5) I} ax ax ay ay 
The total matrix element connecting dle ith and jth nodal points is obtained by summing over all elements of which dle i th and th points are bodl members: 
K..=~k.e. 
(8-6 
)I} ~ I} e 

We have now sketched out the construction of the left-hand side of the matrix equation (8-3), and it remains to define the generalized force vector, F. In general, the treatment of driving terms in the finite element method is also connected with the varia- tional form of the problem, but in the case of simple specified value (Neumann) boundary conditions, the form of F is relatively easy to demonstrate. 
Consider the system of Eq. (8-3) to include equations for the specified boundary values, as well as for the bona fide unknowns. These "fixed" degrees of freedom are removed from the linear system by moving all terms involving such values to the right- hand side. Thus, if the fixed boundary potential at point i is denoted by gi, we may write the contribution from a given element to the right-hand side vector F as: 
fie=-L kij gi (8-7) i-i 
Here, kij is the element stiffness of Eq. (8-5), and the summation runs over all nodal points j in element e at which Neumann boundary conditions apply. As with the stiffness, the global right-hand side F is obtained through assembly of element contribu- tions after the manner of Eq. (8-6). 
All of the foregoing discussion has been concerned with the matter of constructing the linear system (8-3). While it is this task that is most specialized and specific to the finite element method, it is really somewhat tangential to our principal topic of con- currency. We will be concerned with the setup and structure of the stiffness matrix and right-hand side of (8-3) mainlyas these influence the task of solving the system. 
Our relative emphasis on solving Eq. (8-3) rather than the details of constructing it is based upon two observations. First, although calculating the required matrix elements from Eq. (8-5) by numerical integration is not an entirely trivial computational cost, it is in many practical cases a small cost compared with that of inverting K. Although this is 

Sec. 8-3 The Equation Solver-Conjugate Gradient Method 

not universally true, the solution step represents the main obstacle to extending finite ele- ment methods to large domains in many applications. 
The second reason for skirting the issue of matrix setup is that the concurrent decomposition of this task is comparatively trivial to implement Since both the matrix K and vector F naturally separate into independent element contributions, the processors of a concurrent ensemble can construct the element quantities for different elements in parallel, without the need for any interprocessor communication. Therefore, this step of equation setup has a concurrent implementation identical to that for a sequential machine. One potential complication would appear to arise from the need to combine element con- tributions from different processors into a global matrix equation. As we shall see in the next section, however, our choice of solution technique will allow us to do without this global assembly step. 

8-3 The Equation Solver-Conjugate Gradient Method 
As mentioned in the introduction to this chapter, iterative equation solvers called gradient methods, exist which are applicable in the case of finite element equations. It will be useful to start by briefly describing the basic approach of gradient methods before detailing a particular one. 
Consider the task of iteratively solving M simultaneous equations. We may regard this problem as equivalent to finding the M-dimensional vector which minimizes some residual error quantity defined on the M -dimensional space. A gradient method makes use of trial values for the variable at step i to generate new values at step i + I correspond- ing to a reduced value of the error function. By successively moving "downhill" toward the error minimum, the method converges toward the desired solution. Different gradient methods differ in their technique for choosing descent directions. For the discussion in this chapter, we will examine the conjugate gradient method [Jennings 77], which employs descent direction vectors that are mutually orthogonal relative to an inner pro- duct weighted by the stiffness matrix. 
We begin by defining the quantities used in the conjugate gradient algorithm. The stiffness matrix K and "force" vector F are assumed to have been already constructed, and will remain unchanged by the algorithm. The remaining quantities are three M - vectors and two scalars whose contents are changed with each iteration. The iteration number to which a given M-vector applies will be indicated by a superscript. Thus, <l>(O) denotes the starting trial solution, which in many practical cases is simply zero. 
The other two M-vectors required are r (k), the residual vector, and p (k), the conju- gate search direction. The conjugate gradient algorithm begins by initializing these M- vectors: 

r (0) = p (0) = F -K <l>(0) Then, dle following iterative loop labelled by k is performed: 
(r (k) .r (k» 
I. a = 
k (p (k) .Kp (k» 

(8-8) 

140 

2. 

Elliptic Problems in Two Dimensions ll: The Finite Element Method Chap. 8 

cI»(t+l) = cI»(t) + at p (t) 

3. r(k+l)=r(k)-<xkKp(k) 

4. 

(r (1+1) .r (1+1') ~1 = 
(r (1) .r (1') 

5. p (1+1) = r (1+1) +J31 p (1) 

(8-9) 

6. k = k+ I; go to I. (continue until converged) 

The iterations are terminated when cI>(k) has converged to within some desired accuracy, as determined by the magnitude of the vector of residuals, r . 
This algorithm does not depend upon any particular column or band structure, but it does require the matrix K to be symmetric and positive definite. Both of these conditions are satisfied by the finite element matrix of the present problem. It is significant that the only use that is made of K in this algorithm is to form its product with an M-vector. In order to perform this product, it is not necessary to assemble the element contribution matrices into a global K, since the element contributions to any such product may be computed separately and summed afterward This approach to computing a stiffness pro- duct is not of great importance for a conventional sequential computer in which the entire matrix K is readily available. In the concurrent implementation, however, a given pro- cessor may need access to parts of K residing in another processor. In this context, the element-wise modular approach to calculating products is vital. 

8-4 Concurrent Decomposition 
In order to see how best to implement this algorithm concurrently, let us look again at the operations carried out in Eq. (8-9). The iterative loop includes two basic kinds of operations: the vector inner product (such as r (k) .r (k,>, and the matrix-vector product Kp (k). Each of these operations may be cast in a parallel form, as outlined below. 
In the case of the vector inner product, the decomposition is easily seen. Each pro- cessor of the concurrent ensemble is given responsibility for a subset of the degrees of freedom which comprise a vector such as r or p. Each processor concurrently calculates its contribution to the dot product, and the final result is o"btained by summing the proces- sor contributions via a collective combination step. If we neglect the single collective communication needed to assemble the inner product, this calculation involves no exchange of data between processors, and is limited in its concurrent efficiency essen- tially only by load balancing. As long as the processors are given equal numbers of degrees of freedom, the concurrent efficiency of this part of the calculation is close to unity. 

We might assume at this point that the manner in which degrees of freedom are grouped for partitioning among processors is irrelevant While this would be true for the 

Sec. 8-4 Concurrent Decomposition 

141 

vector dot product alone, the need to work with the stiffness matrix K dictates an element-based spatial decomposition, similar to that employed in the finite difference case. As stated above, a given row (or column) of this matrix contains nonzero elements only between degrees of freedom corresponding to points that share a common element Consequently, a given entry of the Kp product vector only makes use of information that is spatially local to the corresponding grid point Thus, it makes sense to group the degrees of freedom in processors according to contiguous spatial groupings of elements. 
At this point. we can see the importance of the fact that it is not necessary to assemble the matrix K in its global form. The elements of the Kp vector may be con- structed from the individual element "pieces" of K given in Eq. (8-5). Special considera- tion is required in dealing with degrees of freedom which fall on the boundary of a processor's region of responsibility. Since such an "internal boundary" point requires stiffness contributions from elements residing in more than one processor, it becomes necessary to exchange data between neighboring processors. 
The calculation of the Kp product is handled by first performing computations for all elements internal to a given processor. This is followed by a communication cycle which brings in the boundary Kp contributions. The communication cycle must be designed in such a way as to additively accumulate all element contributions to a given degree of freedom and then to distribute the result to all processors which share that nodal point We satisfy this requirement by invoking a two-pass sequence of communica- tion between processors. 
According to our definition of a processor region, each finite element belongs unambiguously to one processor, but the nodal points along outside edges are shared. We shall designate the edge points along the -x and -y boundaries of a given processor to be accumulated within that processor (Fig. 8-2). Conversely, the remaining +x and +y edges are designated as sent degrees of freedom, which are to have their contributions to the Kp product accumulated within the neighboring processor. 
In this way each nodal point is accumulated in one and only one processor. The +xy corner nodal point is special in that its accumulation destination is in the processor located diagonally to the upper right In the present discussion, this contribution reaches its destination through two successive communication steps. First, a left-to-right shift and accumulation is performed. Then, a down-to-up step is performed which leaves complete values of the Kp result in the destination slots. Since each processor needs a copy of the Kp result for subsequent calculation, the communication cycle is completed by reversing the direction of passing. An up-to-down pass followed by a right-to-Ieft pass distributes the result to all processors. It is important that in the reverse pass the communication steps deliver copies, as opposed to performing additive accumulation as in the forward pass. The passing of corner values could be accomplished directly within the capabilities assumed for our virtual concurrent processor environment In the present example, however, we use the two-step process to illustrate a means of implementing th.e method on machines such as hypercubes with restricted interprocessor links. 
Figure 8-2 might suggest that this scheme for performing boundary communica- tions requires a regular rectilinear arrangement of elements and nodal points in order to 

142 

Elliptic Problems in Two Dimensions ll: The Finite Element Method Chap. 8 

l'lgure 8.2 Processor regions overlap on the outer layer of nodal points. Accumulated degrees of freedom (shaded) are a processor's local responsibility, while the other edges and comers are under neighboring processors ' jurisdictioo. 
work. In fact, it is irrelevant in this scheme whether the processors contain regular or congruent arrangements of elements. The only requirement in this regard is that neigh- boring processors agree on the number and kind of shared boundary nodes along their common edges. Aside from this consideration, it is only necessary to ensure that roughly equal numbers of degrees of freedom are contained within each processor to achieve a good load balance. 
The generality of problems which can be handled by the above described scheme depends on machine-specific communication considerations. A grid such as that in Fig. 8-3(a), which is easily mapped onto a simple two-dimensional processor mesh, will work under even the restrictive assumption of a machine that supports only nearest-neighbor mesh communications. On the other hand, a grid which is used to model multiply con- nected or "dendritic" regions as in Fig. 8-3(b), may be difficult to map efficiently onto a machine with only hardwired mesh connections. Such problems may be more effectively handled by concurrent processors with the capability of more general inte1]>rocessor communication. In the final section of this chapter, we briefly discuss the use of crystal- line communication routing utilities in this regard. In such cases, the details of the pass- ing scheme described above for a simple mesh would no longer hold, but the concepts of accumulated and sent degrees of freedom within the conjugate gradient algorithm would remain valid. An alternative approach to handling irregular meshes within a restricted communication topology is the scattered decomposition [Morison 86] which is used in a related context in Chapter 17 for decomposing an inherently irregular problem. 

Sec. 8-5 Example Implementation 

143 

simply connected multiply connected (a) ( b ) 
FIgure 8.3 Schematic illustration of contrasting dOOlain decOOlposition problems. (a) A 
simply connected region is easily accommodated by mesh processor communication. (b) Multiply coonected dOOlain requires more COOlffiunication flexibility or the scattered decomposition. 
8-5 Example Implementation 
The following examples of pseudocode serve to illustrate the concurrent conjugate gradient algorithm outlined in the previous section. Much of the programming needed to accomplish such tasks as intemctive control, grid setup, and input/output may be similar in form to that already given in the case of finite differences. For this reason, we focus the following examples on the conjugate gmdient algorithm and its associated communi- cation opemtions. 
We will examine here the structure of the node subroutine conj-grad, which per- forms the central task of solving the linear system of Eq. (8-3). This subroutine would normally be called from the main-Ievel program following the steps of data input and matrix genemtion. The overall structure of the routine follows that of the itemtive loop in Eq. (8-9). The loop is punctuated at intervals by communication calls needed to exchange neighboring data or to combine results across the ensemble. The modified con- jugate gmdient algorithm has the following form: 

1. 

a) b) c) 

perform r .r product within nodes (calc.) 
sum r .r contributions across ensemble (comm.) perform Kp product within nodes (calc.) 

144 

2. 3. 4. 

5. 6. 

Elliptic Problems in Two Dimensions II: The Finite Element Method Chap. 8 

d) accumulate neighbor contributions to Kp (comm.) e) perform p .Kp product within nodes (calc.) 
t) sum p .Kp contributions across ensemble (comm.) g) calculate a (calc.) update <1» (calc.) update r (calc.) 

a) perform new r .r product (calc.) 
b) sum new r .r across ensemble (comm.) c) calculate 13 (calc.) update p (calc.) 
increment iteration counter and continue 

In the above concurrent algorithm, we have labeled the individual steps according to their function as either a calculation (calc.) or communication (comm.) task. The sequential counterpart of this algorithm is regained by simply deleting the steps labeled as pure communication steps. Thus, a relatively clean separation between functions specific to the concurrent algorithm and those of the conjugate gradient method proper is apparent. Code 8-1 contains the pseudocode implementation of this loop, and is indica- tive of the programming necessary to implement the conjugate gradient solver within a typical hypercube or other message-passing machine environment 
In this example, the iterative loop is "primed" by initially calculating the norm of the starting residual vector r (0). Concurrent calling of the subroutine dot-global results in the calculation of the global scalar product of vectors whose elements are distributed among the ensemble processors. The details of this subroutine will be discussed below in connection with Code 8-2 and Fig. 8-4. 
The iterative loop itself makes use of similar constructs. The intermediate scalar products of p .Kp and r .r are also obtained through calls to dot-global. The remainder of the loop simply fqllows the sequential algorithm of Eq. (8-9). 
The routine Ky yroduct is not described in detail here, because its structure does not involve any explicit issues of concurrency. As indicated in our earlier discussion, all that is important for our example is that K y yroduct is designed to obtain the product of the stiffness matrix K with an M -vector, using an element-wise additive decomposition of K. In calling this routine, each processor obtains its contributions to the product vec- tor as if the elements contained in other processors do not exist The overlapping contri- butions to Kp from neighboring processors are subsequently accumulated in the global- ize call. 
Before dropping the subject of the formation of Kp , it should be noted that in an application on a given concurrent processor, the choice of an optimal strategy for han- dling the stiffness matrix is quite important, being impacted by such machine-dependent factors as the amounts of fast memory available to each processor. Alternate strategies may be formulated to minimize either the storage requirements or arithmetic operation count in obtaining Kp, depending upon the machine-imposed constraints. 

Sec. 8-5 Example Implementation 

145 

proc-begin conj-.Krad( perfonn k conjugate gradient iterationS ) 
proc-call dot-global ( obtain ensemble-wide scalar product of r with r place result in old-norm) 

for-begin (a total of k cycles) 
proc-caU K y yroduct (multiply current p by K in each processor; place result in vector Kp ) 
proc-caU globalize (bring in and combine neighboring processor parts of Kp; leave result in vector Kp ) 
proc-caU dot -Blobal ( obtain ensemble-wide scalar product of p with Kp ; place result in p -dot-Kp ) 
[ 11 is replaced by old-normlp -dot-Kp ] 

proc-call add-scale-vector ( multiply vector p by a and add to solution vector x; leave result in vector x ) 
proc-call add-scale-vector ( multiply vector Kp by -(a) and add to vector r; leave result in vector r ) 
proc-call dot -Klobal ( obtain ensemble-wide scalar product of r with r ; place result in new-norm ) [ 13 is replaced by new-normlold-norm] 

proc-caU add-scale-vector ( multiply vector p by (3 and add to vector Kp; leave result in vector p) 
[ old-norm is replaced by new-norm] for end proc-end 

Code 8-1 Conjugate Gradient Iterative Solver 
We next examine the two routines dot-global and globalize. which incorporate the interprocessor communication step vital to the concurrent implementation. The imple- mentations of both of these routines make use of an assumed ordering of the degrees of freedom within a given M -vector. This structure is illustrated schematically in Fig. 8-4. 
In the construction of intermediate vectors such as r or Kp. we have chosen to order the entries so that the sent and accumulated degrees of freedom are grouped separately in memory .Such an ordering simplifies the task of selectively manipulating quantities of different communication type. As an example. the vector product operation of dot global must take care not to "double-count" the contributions to the product which come from shared nodal points. In the illustrated scheme. the first ndof -sends entries of the M-vector correspond to those degrees of freedom for which some other processor has responsibility .Therefore. the dot-global routine is simply required to perform a dot pro- duct on the last ndof -total-ndof -sends entries of the multiplicand vectors. This step is followed by a global sum of the node contributions. Code 8-2 contains the pseudocode implementation of this routine. 

146 

Elliptic Problems in Two Dimensions ll: The Finite Element Method Chap. 8 

M = ndof -total (all degrees of freedom in processor) 

~ -~ ndof sends (degrees offreedom to be communicated) 
~ W ;:/ ;:/ ;:/ ;:/ d , 

." ' 
1sent" 

~ccumulate~ I 

"local" degrees of freedom ~ 

F1gure 8.4 The nwnbering order of degrees of freedom within each processor is con- trived to simplify the handling of communication. Accumulated and sent nodal points are stored as contiguous blocks so that they may be moved en masse during interprocessor communication. Shaded blocks represent local processor responsibility. 
Also illustrated in Code 8-2 is the implementation of globalize. Whereas dot-global is used to obtain and distribute a glDbal scalar quantity, the globalize subrou- tine is required to obtain neighboring processor contributions to a glDbal M-vector built up from element pieces. 
The globalize routine utilizes the two-pass cyclic communication scheme discussed 
above in connection with Fig. 8-2. Neighboring processors exchange and accumulate shared edge values via the communication buffers right-buf, left-buf, up -buf, and down-b~ Implicit in this treatment is the assumption that adjoining processors agree on the mapping between items communicated through these buffers and their ultimate desti- nations in the M-vector. In practice, such a mapping can be accomplished by indexing of the buffers, or in complicated geometries by explicitly passing a glDbal index with each value. The latter alternative entails a greater communication cost, but may be necessary in cases where it is impractical to enforce a matching of local indices across processor boundaries, as might easily occur with complicated grids. 

Sec. 8-6 Perfonnance of the Finite Element Exantple 

147 

proc-begin dot -Elobal ( obtain ensemble-wide scalar product of vector a with vector b ) for-begin (loop over degrees of freedom in processor's jurisdiction ; 
k running from ndof -sends to ndof -total ) 
[ multiply kth element of a by kth element of b and sum in result ] for end 
[ vm-combine -add values of result across ensemble and replace result ] proc-end ( return result ) 

proc-begin globalize (bring in and combine neighboring processor contributions to vector a ) [ place a entries to be sent to right neighbor in right-butbuffer ] 
[ vm-shift -send right-bufto right neighbor; receive from left neighbor into left -but] [ add entries of left -but to corresponding entries of a, replacing them ] 

[ place a entries to be sent to upper neighbor in up -bufbuffer ] 
[ vm-shift -send up -bufto upper neighbor; r~eive from lower neighbor into down-buf] [ add entries of down-bufto corresponding entries of a, replacing them ] 

[ replace summed a entries from lower neighbor in down-butbuffer ] 
[ vm-shift -send down-butto lower neighbor; receive from upper neighbor into up -but] [ replace entries of a with corresponding entries of up -but just received ] 

[ replace summed a entries from left neighbor in left -bufbuffer ] 
[ vm-shift -send left -bufto left neighbor; receive from right neighbor into right-buf] [ replace entries of a with corresponding entries of right-buf just received ] 
proc-end 

Code 8.2 C<Xnmunicatioo Routines for Conjugate Gradient Solver 
The routines outlined in Code 8-1 and Code 8-2 comprise the core of our conjugate gradient finite element solver. This particular example presents a somewhat special case- in that we have treated only two-dimensional problems and limited our discussion of complicated domains to general remarks. Nevertheless, the above implementation should serve to illustrate the central concepts common to more sophisticated programs, and to indicate some of the directions in which this work might be readily extended. 

8-6 Performance of the Finite Element Example 
As in the discussion of a finite difference application in the previous chapter, we now present the results of implementing our finite element/conjugate gradient algorithm on the Caltech Mark II hypercube. Much of the discussion in Chap. 7 relating to the theoretical efficiency of the concurrent algorithm is directly applicable here. The com- munication load varies as the number of perimeter nodal points per processor, while the computational load (which is dominated by matrix-vector products) is approximately pro- portional to the total number, nel' of elements per processor. As a result, we would expect the two-dimensionallattice efficiency dependence given in Chap. 7: 

148 

Elliptic Problems in Two Dinlensions II: The Finite Element Method Chap. 8 

£ ::: 1- cne/l12 (8-10) We should expect detailed differences in the communication and computation tasks 
to be reflected in different values of c for the finite element and finite difference exam- pies. Coincidentally, these differences approximately cancel to yield a value of c ::: 1 in Eq. (8-10), as obtained earlier for finite differences. In Fig. 8-5, the concurrent efficiencies are plotted on a logarithmic scale to illustrate the dependence of efficiency on problem size and number of processors. 

Measured Efficiency E of 'wo Dimensional Finite Element Code 

lognel 

JO 100 1000 10,000 grain size nel 
F1gure 8.5 Experimental efficiency results fran a two dimensional finite element pro- gram nmning on Caltech Mark II hypercube. Solid line indicates theoretical efficiency trend Eq. (8-10). 
As was noted in Chap. 7, configurations of fewer than 16 processors fall on higher parallel efficiency trends. The 32-node configuration, which involves processors com- municating on all four sides, yields efficiencies lying close to the c = 11ine. We may conclude that the finite element implementation presented in this chapter shows good agreement with theoretical efficiency predictions and yields useful concurrent speedup on actual machines for more than about 10 finite elements per processor. A performance analysis of finite element problems on a shared memory ELXSI 6400 and comparison with the hypercube willbe found in [Walker 86] and [Montry 87]. 

Sec. 8- 7 Extensions and Scientific Applications 

149 

8.7 Extensions and Scientific Applications 
The basic approach to solving finite element systems on concurrent processors out- lined in this chapter has been employed to solve a number of scientifically interesting problems in continuum mechanics. This approach also forms the basis for ongoing work with the goal of solving more challenging problems [Nour-Omid 87a, b]. 
The concurrent conjugate gradient algorithm described above has been used [Lyzenga 85] to solve finite element systems arising from linear elastostatic boundary value problems. The form of the operator and resulting stiffness matrix is somewhat more complicated than those in the Laplacian example. The governing equation is given in tensor index notation by: 

aijj = -ii (8-11) 
Here, i i denotes a vector body force (such as gravitation) and aij is the Cauchy stress tensor which is related to material displacements, Ui , by: 
a.. =A.uLL~.. + 
Jl(u. .+u. . 
)(8-12 
).1 1 I.} I.. where A. and Jl are the Lame elastic constants for an isotropic medium. We employ the usual conventions here of commas for differentiation, repeated indices for summation, and ~ij for the Kronecker delta unit tensor. 
This equation gives rise to a problem similar to the Laplacian case, in that it involves only second derivatives of the dependent variables Ui, and yields a symmetric, positive-definite matrix. In work performed on the Mark II hypercubes, the conjugate gradient solver has been adopted unchanged, except for the addition of a simple precon- ditioner based upon a rescaling of the matrix by its diagonal elements. This problem has been solved for both the two-dimensional domain discussed above, and for three- dimensional domains. 
The matter of preconditioning for iterative equation solvers is an important one if such techniques are to find broad practical application. While iterative methods like con- jugate gradient are characterized by comparatively modest computational costs per itera- tion, the number of iterations required to obtain a solution to within a desired accuracy can be highly problem-dependent 
In general, matrix systems with widely varying magnitudes of eigenvalues present a problem for gradient methods. In such cases, when working with finite numerical pre- cision, the solver may become "lost" in trying to find a residual minimum in M-space, since the shapes of level surfaces in this space are extremely nonspherical. Precondition- ing amounts to the application of some transformation to the original matrix system which is intended to render it more tractable for the iterative solver. 
In the finite element matrix, the diagonal elements are proportional to the areas (or volumes) of the elements associated with each degree of freedom. Therefore, one soun:e of ill-conditioning in finite elements is any strong variation in grid refinement Such ill- conditioning may be alleviated to some extent by rescaling the system so that the magni- tudes of the diagonal elements are constant. In practice, this kind of diagonal precondi- tioner can substantially accelerate the convergence of the conjugate gradient method. 
~ 

150 

Elliptic Problems in Two Dimensions ll: The Finite Element Method Chap. 8 

However, it may be of limited usefulness with other sources of ill-conditioning. For example, large contrasts in matrix elements within the same row or column of the stiff- ness will not be remedied by diagonal scaling. 
In general, the more closely a preconditioning transformation approximates the actual inverse of the matrix, the more successful it will be in accelerating convergence. A number of promising approaches to this problem are applicable to the concurrent case. One such approach represents the approximate matrix inverse by a product of local ele- ment matrix inverses [Muller 85]. This element-by-element preconditioner can offer good performance for some classes of problems, but it also involves a sequential matrix factoring scheme which must be performed in a consistent manner on a concurrent pro- cessor . 
By analogy with our discussion of the inherently sequential SOR update scheme in Chap. 7, it becomes necessary to apply the element-by-element preconditioning step in spatially disjoint subregions of adjoining processor regions. Since the finite element method yields diagonal connections between grid points as well as nearest neighbors, the four-color isolation scheme illustrated in Fig. 7-6 is needed. 
Although not formally regarded as a preconditioner, the multigrid method offers yet another approach to increasing the robustness of iterative finite elements. In the mul- tigrid method, a succession of meshes of progressively greater refinement is used. The solution obtained on a coarse grid at modest computational cost is projected onto the next finer grid as input to the next approximation to the ultimate solution. 
Under some circumstances, this process of projecting the "low-frequency" solu- tions onto progressively finer meshes has the effect of giving the ultimate iterative solu- tion step an easier system to solve. The application of this approach to large granularity concurrent processing is straightforward, as it involves the same domain decomposition schemes already discussed. The projection and interpolation of solutions between grids of different refinement involves concurrent vector/matrix operations completely analo- gous to those treated above in the discussion of simple conjugate gradients. A basic mul- tigrid scheme has been implemented on the Caltech hypercube machines [Cisneros 87; Van de Velde 87a, b]. 
At first glance, it might seem that the multigrid method is not easily implemented on concurrent computers because of the involvement of steps using a coarse grid which allow little concurrency and hence constitutes a sequential bottleneck. This is not true, however, and the point provides an interesting comment on the relevance of Amdahl's law to concurrent computing [Amdahl 67, Chan 87, Fox 84c]. Let us consider a mul- tigrid method in which the finest grid incorporates n nodal points per processor node. In its simplest form, the method is comprised of successive stages labelled by k, each involving nk = n/2k members. We see that the total calculation time is: 
T calc -L nk tcalc (8-13) k 

while the communication time is given by: 
T comm -I. nkl12 tcomm "