Sec. 20-4 Concurrent LU Algorithm 

377 

Table 20-2: Estimated Times for Stages in LU Decomposition 

Code 20.1 The split yipe Sub-cube Broadcast Algorithm 
This is larger than the value obtained for matrix multiplication in Eq. (10-19) for two main reasons. First 32-bit numbers are used in Chap. 10 and 64-bit numbers are used here: tcomm (32-bit) = 0.5 tcomm (64-bit), whereas the 8087 chip used on the Mark II hypercube performs single and double precision calculations in roughly comparable times. This accounts for a factor of two in the discrepancy with Chap. 10. The second 

378 

Matrix Algorithms ll: Banded Matrix LU Decomposition Chap. 20 

cause of additional overhead is the nonoptimal implementation of the crOS III routine vshift on the Mark II hypercube. The scattered decomposition leads naturally to the use of vshift to transfer rows and columns of the matrices. However, vshift on the Mark II hypercube is implemented in terms of cshift so that two additional memory-to-memory copies are necessary .In an efficient implementation this copying can be avoided on hypercubes like the Mark II where messages are packetized. The routine vshift will always be inefficient if the hardware is optimized for DMA (direct memory transfers) from contiguous memory to the communication system. An implementation of these algorithms without using the scattered decomposition and so with efficient communica- tion had t = tcomm/tcalc=2.7 .We can expect an optimal version of vshift to give a t value similar to this. 

- I ~ - 
'i ~ Q) 
> 0 

0.4 

0.3 

0.2 

0.1 

0.0 

LU Decomposition on the Caltech/JPL Mark" Hypercube ",".. (b) Nopivoting. 't=4.7 ,"' j( , (a) Load balanced ,"' X... 
(b) Maximal imbalance ,"' ",:f (a) , 'j' .t: 4 processors ".: ~... ..".: 
Models ," K' ,... 
64 processors} ..-1...'... ,.r ,tr ,".:... 
......~ '... 
..' " " " ...~,... ,' -,. ,'... ~' , :" '.:... ...:"' 
,,:;. ...:"' x 4 
processors 
;. ..- 
,.r~:"~-' + 16 processors 
?~ .64 processors ~ 

0.0 

0.01 

0.02 

0.03 

0.04 

0.05 

" l/m 

Figure 20.8 Measurements of the concurrent overllead 1/£ -1 as a function of the natural parameter -/N" Im for LU decomposition with no pivoting. We show results from the Caltech/JPL Mark II hypercube for 4, 16 and 64 node machines. The results (a) correspond to the balanced case where m is an exact multiple of ..IN .(b) corresponds to m values that are 1 more than such exact multiples. General values of m lie between these extreme cases. We also show the expectations for (a) and (b) C(Xlling from the model Eq. (20-12, 13) for the parameter value I It 
1=4.7 comm ca c 
It is apparent from Fig. 20-8 that the model agrees quite nicely with the measure- ments. The overhead is proportional to 11m and the model adequately describes the difference between the best and the worst case efficiencies (i.e., the effect of the load imbalance when -{Ff does not divide m exactly). A fortuitous cancellation occurs between the 1.5Im-{Ff and 't (-{Ff -2)lm2 terms in Eq. (20-12), so that the best case efficiency shows little dependence on N when plotted as a function of m This cancella- tion does not occur for Eq. (20-13), and the divergence of the N = 4 and N = 64 meas- urements is a direct consequence of the pipe startup term being proportional to (-{Ff -2). 

Sec. 20-5 Pivoting in the Concurrent LU Decomposition 

379 

In practical applications, we might expect that a typical value of m lies between the best- and worst-case results, and quite good performance is obtained in this intermediate case. For the case of a 20 x 20 matrix in each node, the concurrent overhead is about 35% on a 64-node hypercube. 
As long as ni= {ii" is large compared with (~ -2), the inefficiency follows the expected relation: 

E n corresponding in the language of Chap. 3 to a system of system dimension ds = 2. This is a general result; all full and banded matrix algorithms have dimension 2. This is an interesting result because matrices have no metric associated with them. It would appear that the system dimension ds can reflect the dimension of the world being decomposed 
lId rather than some local property of a metric space where l/n s is a natural formula for 
the generalized surface to volume ratio. Both matrix multiplication and LU decomposition have dimension 2, correspond- ing to an overhead proportional to l/n 112, with n equal to the grain size. In multiplica- tion, we find a basic calculational cost proportional to n 3/2 and communication costs linear in n. In LU decomposition, calculation is proportional to n and communication varies as -fn .The ratio of the communication to calculation costs has the same (1/{ii:) dependence in each case, corresponding to the same system dimension. It is interesting that this uniformity occurs even though the calculational complexity has a quite different dependence on grain size in the two cases. 

f= 

constant -1/2 

(20-15) 

20-5 Pivoting in the Concurrent LU Decomposition 
The algorithm described in Sec. 20-2 and 20-4 is appropriate when the diagonal elements Akk are dominant. This 'is often the case, but in general, attention must be given to numerical instability resulting from the division by Akk if its value is small or vanish- ing. This problem is handled by means of a pivoting rearrangement of rows and columns in order to obtain a larger element in the top comer of the window. We consider here the case of partial pivoting in which we search the first column in the window for the max- imum entry Ak+i. k (0 :$;i < m ). If this maximum corresponds to row number k +i max then the rows corresponding to i =0 and i =i max are swapped. Then the algorithm proceeds according to the description of Sec. 20-2. There arises one obvious complication in that the resultant matrix within the window is no longer square. As shown in Fig. 20-9, the window could be as big as m rows and 2m-2 columns before this step. The number of columns depends on the details of pivoting in the previous m -1 steps. Pivoting at step k mayor may not increase the window size. The A matrix window is at most m by 2m -1 in size after the pivoting at step k. This possible increase in window size affects the speed of the algorithm but does not alter in any essential way the steps outlined earlier in Section 20-2. 
The variable size of the window encountered here is handled naturally within the context of the earlier scattered decomposition. This decomposition ensures that the loads 

380 

Matrix Algorithms II: Banded Matrix LU Decomposition Chap. 20 

of different processors differ by at most one column. This illustrates the general utility of this simple decomposition approach. While originally motivated by the need to cope with the banded matrix structure, it represents a powerful method of dealing with other sources of imbalance. Despite the generality of this approach, we have simplified matters in the present example by assuming a constant window size of m by 2m at all steps. This restriction simplifies both the sequential and concurrent algorithms. A11 columns beyond the end of the significant window are assumed to be padded with zeros. This simplification was used in the example runs whose measured performance values are given in Fig. 20-11. 

Row k ~ 0 1 2 3 0 1 2 3 0 1 456745674567 8 9 1011 8 91011 8 91011 121314151213141512131415 012301230123 456745674567 8 910 11 8 910 11 8 910 11 121314151213141512131415 012301230123 456745674567 8 910 11 8 910 11 8 910 11 121314151213141512131415 

Window Edge if 
no pivoting """ 

Range of possible window edges after pivoting 

21 
Maximal 14 bandwidth 2 before 6 . 
d10 wID ow 
14 processed 2 6 10 14 

11 
" 
Lower edge 1 ~ / 
of band Diagonal 
Step k of LU Decomposition with Pivoting 

12301230123012301230123 , 
6 " 

, 

;;;1 

Original Bandwidth 23 16 Processors 

, 

4 

" 

9 

, 

14~pper edge of band 

FIgure 20.9 The window in the matrix A at step k of LU decomposition when partial pivoting is used. The example matrix has the same parameters as that in Fig. 20-5. 

Sec. 20-5 Pivoting in the Concurrent LU Decomposition 

381 

5 Initial Stages of 
LU Decomposition with Pivoting 
Shaded matrix elements L!] Active Unshaded D Inactive 

(a) Find Pivot in First Column 

Find maximum element P; in first column in each of processors i=O, 4,8, ]2 holding such elements 
FIgure 20-10(a} The first of five initial stages in the pivoting algorithm for the krA step 
of the dec<XDposition. These are followed by the stages in Fig. 20-6 just as in the case with no pivoting. 

The concurrent pivoting algorithm incorporates five additional steps. These sup- plementary implementation steps are illustrated in Fig. 20-10 and are described as fol- lows: 

(a) 

(b) 

Each processor containing elements of the first column determines the local maximum of I Ak+i .k I for those elements in the node. This step only involves ~ nodes. The other processors are idle. 
Next, determine the global maximum over these ~ local pivot candidates. At first glance it might seem appropriate to use the combine collective routine which was introduced in Chaps. 4 and 14. In this case, we would apply com- bine not over the full hypercube but over the nodes containing elements in the first column of the window. These nodes make up a subcube. The pseu- docode for performing a general subcube combine is given in Code 20-2. For historical reasons subcube combine is not part of the basic crOS III library , 

382 

(c) 

(d) 

(e) 

Matrix Algorithms II: Banded Matrix LU DecomIX>sition Chap. 20 

but it is given as a utility in the Software Supplement, Vol. 2 of this text. The subcube combine algorithm in Code 20-2 may be used to find the pivot ele- ment and the number of the row containing it, and return these values in a data structure. However, an alternative method is to use a "customized" ver- sion of the pipe algorithm given in Code 10-3. In this algorithm, shown in Code 20-3, each of the processors in the pipe, except the first one, receives a candidate pivot element and its row number from the processor below. The element received is compared with the candidate pivot element of the proces- sor, and the largest is passed on to the processor above, together with the corresponding row number. At the end of this process the processor at the end of the pipe contains the pivot element and its row number. 
Full ensemble broadcast utilities are used to distribute the current pivot loca- tion to all nodes. 
The new pivot row is swapped with the old top row. If they are contained in the same processors (the top row in the processor array), this step involves no communication. Otherwise, the exchange of rows is accomplished by coor- dinated communication between the appropriate processors. In view of the step (e) below which is needed when there was no pivoting, one only sends the old top row to its destination in this stage. 
Finally, the pivot row is sent to all processors. This completes the swap and also implements the second stage of the ordinary algorithm shown in Fig. 20-6. 

(b) Transmit Local Pivots to Processor O 

Processors 

Figure 20-10(b) The second of five initial stJlges in the pivoting algorithm for the kIJ. step of the decomposition. These are followed by the stJlges in Fig. 20-6 just as in the case with no pivoting. 

(c) Broadcast Pivot Information To All Processors 

Figure 20-10(c) The d1ird of five initial stages in the pivoting algorithm for the k~ step of the decomposition. These are followed by the stages in Fig. 20-6 just as in the case with no pivoting. 
(d) Send Old Top Row to Position of Pivot Row 

Processors 12, 13, 14, 15 now contain the same subsets 

of top row as processors 0,1,2,3 respectively. 
Figure 20-10{ d) The fourth of five initial stages in the pivoting algorithm for the kIlo step of the decomposition. These are followed by the stages in Fig. 20-6 just as in the case with no pivoting. 

383 

384 

Matrix Algorithms ll: Banded Matrix LU Decomposition Chap. 20 

(e) Send Pivot Row Vertically 

Processors 0 4 8 12 1 5913 2 6 1014 371115 

now contain same subset of row originally in 
12, 13, 14, 15 only. 

~ 

Figure 20-10(e) The fifth of five initial stages in the pivO(ing algorithm for theklk step of the dec<Xnposition. These are followed by the stages in Fig. 20-6 just as in the case with no pivoting. 

- I 
~ - 
'i ~ 0 

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 

LU Decomposi'ion on dte Cal'ech/JPL Mark II Hypercube With pivo'ing. 't= 4.7 
4 processors 
}Models 
64 processors ,"'. ,.. + ...+ ,.' + x.,..J<"" ,. + ~ , , , , '1.. t ..>s. .. -,. ... 
..: x 4 processors 
;... + 16 processors ~..-;. .64 processors ." ' 
-I I I I I 0.0 0.02 0.04 0.06 0.08 0.1 

,"' . ," + 
.,' . 
+ ~ >5, 

l/m 
F1gure 20-11 Measurements of the concurrent overhead 11£ -1 as a function of the natur- al parameter {N Im for LU decomposition with pivoting. We show results for the Caltech-JPL MarkII hypercube for 4, 16 and 64 node machines. We show the expected dependence of 11£ -1 00 .IN Im from the model of Table (20-3) for the parameter value III--, = 4.7. We note that as described in the text both the coocurrent and sequential 
c~ ~c 
algorithms have been simplified and are nonoptimal. 

Sec. 20-5 Pivoting in the Concurrent LU Decomposition 

385 

Code 20-2 Pseudocode for General Combine Over a Subcube of the Hypercube 

proc-begin find yivot ( find pivot row, p -row. and pivot element, pivot ) 
if begin ( this processor has elements in first colwnn of window) then 
-[find maximum element, max-el! in first column for this processor] [store row number of maximum element in max-row ] If end 
if-begin ( this processor contains the element in the last row of first colwnn of the window) then 
[send max-el and max Jow to processor above] 
else-lf ( this processor contains the element in the first row of the first colwnn ) then [receive caM-el and cand-row from processor below] If-begin ( max-el > cand-el] then [pivot (;- max-el] 
[p -row (;- max Jow] 

else 

[pivot ~ cand-el] 
[p -row ~ cand-row ] if end 
else-if ( this processor contains elements in the first column of the window) then [receive cand-el and cand-row from processor below] if-begin ( max-el > cand-e/) then [ cand-el ~ max-el] 
[ cand-row ~ ;- c' 1 

if end 
[send cand-el and cand-row to processor above] if end 
[broadcast p -row from processor containing element in first row and column of window to all other processors] 

proc-end 

Code 20.3 Detennination of the Pivot Element and Row 

386 

Matrix Algorithms ll: Banded Matrix LU Decomposition Chap. 20 

These pivoting steps are followed by those of the nonpivoting algorithm with the obvious changes required by the increased number of columns. It is apparent that while pivoting is an almost trivial operation on a sequential machine, it can require a fairly complicated algorithm in the concurrent implementation. The concurrent implementation is in this case significantly less transparent than the sequential equivalent. One approach described in Chap. 21 is that taken in Fox and Furmanski [Fox 86b] to build a library of general purpose concurrent matrix subroutines. A single such set of routines can be used to tackle essentially all matrix problems. By way of a general remark, we have found that matrix problems are in many ways the hardest for concurrent machines. The algo- rithms are so involved that although high speedup is possible, the software implementa- tion overhead cost is relatively greater than for the other problems discussed in this volume. 

Stage 

Sequenti~l 

~a) Find Pivot 
(b) Transmit to 
Corner Node 
fc) Inform All Nodes 
r .of Pivot Parameter 

I m~ tcalc 

0 

(d) Send Top Row in 

0 

Window to Position of Pivot Row 

(e) Send Pivot Row Everywhere 
~~rt au 

0 

tjnv 

fg) FonI1 Multipliers 

(m.JN -I) tcalc 

rh) Pipe Multipliers I Across 

0 

ri) Correct Matrix 

2(m-.[fl-l) (2m-.[fl- tcalc 

Time 

Concurrent 

I2(.;~C;)Ctcomm ilmtcalc 
2({Ff -l)tcomm 

2 log2 N tcomm 

({F[ -I) {F[ 

2In + .! ({FT -2)2 4 {FT -I 

tcomm 

tin" 

mtca/c 

A2 
4m tcalc 

Table 20.3 Estimated Times for Stages in LU Decomposition with Pivoting 

Sec. 20-6 Forward Reduction and Backsubstitution 

~R7 

In Table 20-3, we detail the simple model describing the performance of the con- current pivoting algorithm. We will only discuss the "best-case" scenario, when m is an exact multiple of ~ .The entry (d) is this table corresponds to a simple algorithm where the top row is sent by the shortest pipe between the source and destination nodes. In Chap. 21 we will formalize this utility as transfer and discuss optimal algorithms. 
This table contains more steps than Table 20-2 but the two share essentially the same structure. In Fig. 20-11 we compare measured concurrent overheads, 1/£-1, with the estimate obtained from Table 20-3 using the same parameter value tcomm/tcalc = 4.7 as in Sec. 20-4. There is now a more evident separation between the N = 4 and 64 measure- ments due to the increased importance of the subasymptotic (~ -2)/m2 and logzN /m2 terms in 1/£-1. 

20-6 Forward Reduction and Backsubstitution 
We now turn to the use of the factored form A=LU to complete the solution of ~ =Q or more generally AX =B .This task involves forward reduction Q f- L -1 Q , fol- lowed by backsubstitution ,!f-U-1Q , as described for the sequential case in Sec. 20-2. In the case of a single right-hand side, the discussion is largely academic, since the costs associated with this final phase are small (by a factor m) compared with those for LU decomposition. An optimized decomposition for the single right-hand side problem uses techniques discussed by Fox and Furmanski [Fox 86b] and discussed in Sec. 21-3.4. Instead, we will here suppose that Eq. (20-4) is to be solved for nb right-hand sides, and further that nb is large compared with {f[ . 
Forward reduction is illustrated in Figs. 20-12 and 20-13 and consists of three stages in the concurrent implementation of Eq. (20-9). The theoretical times are tabu- lated in Table 20-4 but we will not discuss them in great detail because they follow by analogy from the steps discussed in Sec. 20-4. 
In Table 20-4, we draw an analogy to Eq. (20-11) in making the definition: 
nb = nb/{f[ (20-16) We obtain the predicted concurrent overhead given by 
1'\ 1+~1 

1 +.!..1~+~ f = 1/£ -1 = ~ 2l m nb 

AA mnb 

(20-17) 

We have assumed here that {fl divides m and nb exactly. The modifications required to accommodate other problem sizes follow the treatment given in Sec. 20-4. 
Equation (20-17) is compared with measured performance on the Mark II hyper- cube in Fig. 20-14. The measured and theoretical overhead values are plotted as func- tions of l/nb at fixed m It is noteworthy that the overhead does not vanish as nb ~oo because of the presence of terms proportional to l/ni in Eq. (20-17). 
Backsubstitution is implemented in a fashion similar to forward reduction and is illustrated in Figs. 20-15 and 20-16. The expected performance of this algorithm is given in Table 20-5. 

388 

Matrix Algorithms II: Banded Matrix LU Decomposition Chap. 20 

Step k of Forward Reduction Form of Matrix L 

F1gure 20-12 The fonn of the lower triangular matrix L and its window at step k of for- ward reduction. 

3 Stages of Forward Reduction 

Shaded Matrix Elements Unshaded Matrix Elements 

~ Active D Inactive 

(a) Transmit L Window Horizontally to Right 

? Unit matrix element in processor O 

:;,. 5 :;,. 6 :;,. 7 9 10 11 13 14 15 1 2 3 5 6 7 
:;,. 9 :;,. 10 :;,. 11 13 14 15 1 2 3 5 6 7 9 10 11 
:;,. 13 :;,. 14 :;,. 15 

Figure 20-13(3) The first of three sll1ges involved in step k of forward reduction. We show this for the same size of matrix A as in Fig. 20-5 and a total nb = 8 right hand sides b. 

Sec. 20-6 Forward Reduction and Backsubstitution 
(b) Transmit Top Row in Q Window Downwards 

389 

Processors 0 4 8 12 now contain same 1 5 9 13 subset of top row 2 6 10 14 of current 3 7 11 15 !2. window 
FIgure 20-13(b) The second of three s~ges involved in step k of forward reduction. We show this for the same size of matrix A as in Fig. 20-5 and a total n" = 8 right hand sides b. 

( c ) Correct Elements of Q in Window 

-7j 

k 1 i 

B..~B.. l}~ I} i= k+l 
j= k 

-Bkj Lik to k+m-l to k+m-l 

Figure 20-13(c) The third of three stages involved in step k of forward reduction. We show this for the same size of matrix A as in Fig. 20-5 and a total nb = 8 right hand sides b. 

390 

Matrix Algorithms ll: Banded Mattix LU Decomposition Chap. 20 

Stage 

SeQuential 

(a) Transmit Active ~~umn of L to Right 

0 

(b) Transmit Top Row in B Window Downward 

0 

(c) Correct Elements of B in Window 

2(mVN -l)nb VN tcalc 

Time 

Concurrent 

2mnb tcalc 

- I 
~ - 'i 
! 

Table 20.4 Perfonl1ance of the Concurrent Forward Reduction Algorithm 

1.0 

0.8 

0.6 

0.4 

0.2 

0.0 

Forward Reduction on dte Caltech/JpL Mark II Hypercube 
7" .. 't=4. ,m=20 .. 

4 processors }16 processors Models 64 processors 

x 4 processors ..> + 16 processors 
.64 processors I I' I I 
0.0 0.05 0.1 0.15 0.2 0.25 

~ 

--. 

" l!nb 

FIgure 20-14 Measurements of the concurrent overhead 11E-1 as a function of the parameter v'N Inb for forward reductioo. The parameter ml.JN is fixed at 20. We show results from the Caltech/JPL Mark II hypercube for 4, 16 and 64 node machines. We show the expected dependence of 11E-1 of .IN Inb the model of Eq. (20-17) for the parameter value t It --, = 4.7. 
CDmm ~c 

Sec. 20-6 Forward Reduction and Backsubstitution 

391 

Step k of Backward Substitution Form of Matrix U 

FIgure 20-15 The fonD of the upper triangular matrix U and its window at step k of back substitution. 

4 Stages of Backward Substitution 
Shaded Matrix Elements El:J Active Unshaded Matrix Elements El:J Inactive 
(a) Transmit U Window Horizontally to Lefl 

~ 

t u stored -Ul stored in 
last row 

Figure 20-16(a) The first of four stages involved in step k of backsubstitution. The ex- ample problem has the same parameters as in Fig. 20-13. 

Matrix Algorithms ll: Banded Matrix LU Decomposition Chap. 20 

(b) Solve for Last Row in Window 

-7j 

0123012 4 5 6' 7 4 5 6 
8 9 1011 ~ 91011 1213141512131415 0123012 4567456 
8 910 11 8 9 1011 1213141512131415 0123012 
-E- last 
Flgure 20.16(b) The second of four stages involved in step k of backsubstitution. The example problem has the same parameters as in Fig. 20-13. 

X1ast,j = Blast,j ~ B1ast,j Ui~t 

(c) Transmit Last Row in Window Upwards 

Figure 20-16(c) The third of four stages involved in step k of backsubstitution. The ex- ample problem has the same parameters as in Fig. 20-13. 

Sec. 20-6 Forward Reduction and Backsubstitution 

393 

(d) Correct Elements of ~ in Window 

l 

last -7 

J 

Bij = Bij -UijXiast,j allj 
i <last 

Figure 20-16(d) The last of four stages involved in step k of backsubstitution. The ex- ample problem has the same parameters as in Fig. 20-13. 

Stage 

SeQuential 

(a) Transmit Active Column of U to Left 

0 

(b) Solve for X in Last Row of B, X Windows 

nb .[F[ tcalc 

(c) Transmit New Solution Upwards 

0 

(d) Correct Remaining Window Elements 

2(mVN -l)nb VN tcalc 

Time 

Concurrent 

nb tcalc 

AA 2mnb tcalc 

Table 20.5 Perfonnance of the Concun-ent Backsubstitutioo Algorithm 

394 

~ i ~ 0 

Matrix Algorithms II: Banded Matrix LU Decomposition Chap. 20 

1.0 

0.8 

0.6 

0.4 

0.2 

0.0 

Backsubstitution on the Caltech/JPL Mark n Hypercube " 
't=4.7,m=20 

+ 4ProceSSorS 
}"X " 16 processors Models ," " 64 
proccssors .. 
,p.." , 
...v,- 
...- " " ~.,' .~ 

,... 

x 4 processors + 16 processors .64 processors 

0.0 

0.05 

0.1 

0.15 

0.2 

0.25 

... l!nb 

Figure 20-17 Measurements of the concurrent overhead 1/£-1 as a function of the parameter ..IN Inb for back substitution. The parameter ml..JN is fixed at 20. We show results from the Caltech-JPL Malk n hypercube for 4, 16 and 64 node machines. We show the expected dependence of II£ -Ion ..IN Inb for the model of Eq. (20-18) with the parameter value t It--, = 4.7. 
~ ""c 
The inefficiency overhead associated with backsubstitution assumes the form: 0.5 °-3-+- 
2t r~+J-+~ m l m nb mnb 
which is compared with the measured performance in Fig. 20-17. The derivation ofEq. (20-18) neglects the mutual concurrency of communication and calculation which has the effect of nullifying the 0.5/ni term in Eq. (20-18). In this case, backsubstitution and for- ward reduction have identical predicted efficiencies. It should be noted that Eqs. (20-17) and (20-18) display the same l/ni behavior (proportional to 1/..[ii") obtained in Eqs. (20-12) and (20-13) for LU decomposition. 
Both forward reduction and backward substitution show excellent agreement in Figs. 20-14 and 20-17 between the measured efficiencies and the model Eqs. (20-17) and (20-18). We use the same value of t = tcomm/tcalc = 4.7 as given in Eq. (20-14). 

I-+ f = II£ -= ni-.[fl 

(20-18) 

20- 7 End Effects 
Referring to the earlier Figs. 20-1 through 4, it may be seen that all the foregoing discussion of this chapter applies to the situation in which the window is not near the beginning or the end of the matrix bands. At the beginning steps k <m or at the end k>M-m, the algorithm of Sec. 20-2, 4, 5, 6 is modified in correspondence with the reduced size of the matrix. This raises load-balancing issues as illustrated in Fig. 20-18. However, just as in the other situations involving potentialload-balance problems, the 

Sec. 20- 7 End Effects 

395 

scattered decomposition is effective. As shown in Fig. 20-18, the matrix fluctuates around exact balance but the number of rows and columns in individual nodes differs by at most one. As before, this imbalance leads to inefficiencies of order l/ni for the relevant steps, but for the total problem it becomes unimportant because it only arises in those steps at the beginning and end of the band. As a result, the contribution to over- head from end effects is given approximately by: 

The same load-imbalance issues arise in LU decomposition for full matrices but in the latter case it is even more important because the effect is not reduced by the factor m/M. 

k=M-m ~ 

Last row -7 k=M-l 

I II II I III I I I 

Successive Windows at End of Matrix 

Balanced Decomposition 

Imbalanced Decomposition Window sizem=12 
Figure 20-18 The window 00 the banded matrix for the last stages k =M -m . 
showing the decreasing size and varying load-balance. 

..M-I. 

396 

Matrix Algorithms II: Banded Matrix LU Decomposition Chap. 20 

At this point we can profitably return to the discussion of software overhead ini- tially given in Chap. 10. In Code 20-4 we illustrate three implementations of the back- substitution step which incorporate the necessary tests to cope with the end effects. 
Consider step (d) in Table 20-5, in which the elements of B in the current window are modified. In the naive implementation of Code 2O-4(a), the tests for end effects are contained in the inner loop over rows and columns of B .This leads to a contribution to inefficiency overhead that does not vanish as ni and nb become large. In this case: 
f software ~ a nonzero constant as ni, nb ~ 00 (20-20) 
In the better implementation of Code 20-4(b), the tests are moved outside the inner loop over rows of B .In this case we obtain: 

f software A \~V-~J.J m 
As discussed in Chap. 10, this overhead has the same dependence on problem size as those due to communication and load imbalance. It is therefore difficult to distinguish these effects in an analysis of performance measurements. 
In the best implementation of Code 20-4(c), the testing is moved outside the loops over both rows and columns. This yields the software overhead dependence: 
constant 

= 

constant 

J software = 

AA mnb 

(20-22) 

which is of negligible size compared with other sources of inefficiency. The foregoing discussion typifies the care required in designing efficient concurrent matrix algorithms. 

Code 20-4(a) Naive Backsubstitutioo Algorithm 

Sec. 20-8 Conclusions 

397 

COOe 20-4(b) A Better Backsubstitutioo Algorithm 

Code 20-4(c) The Best Backsubstitution Algorithm 

20-8 Conclusions 
We have demonstrated the manner in which the techniques introduced in Chap. 10 may be generalized to a more complicated matrix algorithm involving matrices of non- trivial structure. Among the key issues identified as important in making the transition from matrix multiplication to the richer problems are the following: 
(i) The importance of the scattered decomposition in reducing load imbalance. (ii) The use of the square subblock decomposition. 

Matrix Algorithms ll: Banded Matrix LU Decomposition Chap. 20 

(iii) The increased importance of short messages with the corresponding greater demand on both the hardware for message transmission and the algorithm used for the pipe. 
(iv) Communication inefficiency retains the dependence on 1/..[ii as in matrix multiplication. 
(v) The presence of components of the calculation for which many nodes are idle is irrelevant in a correctly designed algorithm. This conclusion represents an illustration of the irrelevance of Amdahl's law for large grain size machines applied to this class of decompositions. 
(vi) Concurrent implementation of matrix algorithms requires great care. We have found that this class of problems leads to greater ratio of implementa- tion time on concurrent versus sequential machines than for most other algo- rithms. 

21 

21-1 Introduction 
This chapter generalizes and builds upon the methods introduced in Chaps. 19 and 20. In this chapter, we present a set of communication utilities that use algorithms simi- lar to addvec (Chap. 19). Where appropriate, we illustrate these algorithms by applying them to matrix problems. In particular, Sec. 21-3.4 contains an extensive treatment of matrix-vector multiplication. A summary of the issues involved in the performance of matrix algorithms on the hypercube is given in Sec. 21-5. 
It should be emphasized that these advanced algorithms are nearly all hypercube- specific and are not immediately applicable to machines of lower connectivity .This res- triction stands in contrast with the discussion of Chaps. 10 and 20, in which a two- dimensional mesh topology was adequate. Also noteworthy is the fact that the utilities developed in this chapter are equally applicable to a full hypercube and to any subcube of lower dimensionality contained within the ensemble. For clarity we will present only the full hypercube versions of each algorithm. 

21-2 Cube Geodesics 
The problem of multiple scalar products discussed in Chap. 19 was solved by the introduction of the routine addvec. We note that Chap. 19 made use of a rather general communication strategy which we will call the method of ordered cube geodesics. In fact this strategy is also used in Chap. 18 for the parallel shellsort. This general approach will apply to several other interesting utilities. Consider the task of sending a message between any two nodes n 1 and n2 of the hypercube. In the case of addvec, we accumu- lated the product, altering the message during transmission. However, in the utilities dis- cussed in the present chapter, we will primarily concern ourselves with algorithms which leave the message untouched after transmission. 
A shortest path between nl and n2 is a geodesic, and ifnl and n2 are h hops apart then there are h ! geodesics in the hypercube connecting them. This point is illustrated in Fig.21-1. In Chap. 19 all the messages traveled along geodesics. Now consider a gen- era! problem in which a collection of J messages, labeled by the index j, are to be sent concurrently from nl(j) to n2(j). It is natural to use a geodesic for each message but this leaves two problems to be addressed. Firstly, there may be many geodesics to choose from, and secondly, in order to minimize the execution time on the hypercube, it is not strictly sufficient to make messages travel along a geodesic. This approach could lead to 

399 

400 

Communication Strategies and General Matrix Algorithms Chap. 21 

communication imbalance or "hot spots" and resultant poor concurrent performance. However, we note that if every message travels along a geodesic and the communication is balanced (equal) among processors, then we are guaranteed an optimal implementa- tion. The method of ordered cube geodesics yields a particular choice of geodesics that ensures load balance for many important applications. The algorithms in Chap. 19 used this met.'od as do most of those described subsequently in this chapter. This particular solution is specific to the crystalline communication environment. An alternative approach to this message-routing problem would allow each message to "find" an optimal path dynamically, based upon the changing state of the communication channels in the hypercube. The dynamical approach is inappropriate in the present context since we do not have unpredictable or chaotic message traffic. Rather we have a large number of messages whose movements are highly correlated, so it makes sense to consider the message traffic as a whole and optimize it. This task is straightforward within the syn- chronous environment provided by the crystalline communication system. Generally speaking, the problems discussed in this book exhibit message traffic that is highly corre- lated in space and time, so optimizing the overall message traffic, rather than the routing of individual messages, represents a viable approach. 
In the next section we present specific examples of the ordered cube geodesic method and relate them to Chap. 19. In Sec. 21-4 and Chap. 22 the discussion is extended to the important class of problems in which a message is to be sent from one node to several other nodes. 
The following auxiliary crystalline utilities will be introduced in this chapter: 
index fold 
expand 
transpose transfer scatter 
broadcast 
In addition, the discussion will include the following routines which have been intro- duced in Chaps. 11 and 19. 
global-reverse addvec 
All these routines are described in the Software Supplement with full implementations in C that employ the crOS III primitives. In practical applications it may be advisable to optimize these basic implementations for a specific machine environment. In the case of the Caltech Mark II hypercubes, the optimized code runs about twice as fast as the port- able routines that are written in C and call crOS III primitives [Fox 86b]. 
Figure 21-2 schematically summarizes the action of all the utility routines dis- cussed here and in ChaD. 22 on a simDle 4-node machine. 

Sec.21-2 Cube Geodesics 

a) 

b) 

401 

Geodesics from nlto n2 

One hop 

Two hops 

Three hops 

l'1gure 21.1 Sample geodesics for ooe, two, and three h~ distances between nodes nl and n2" 

21.3 Homogeneous Communication Algorithms 

21-3.1 Data Structures 
Many of the crystalline utilities can be usefully applied to convolutions of the gen- era! mathematical fonn: 

g<!) =F<!,f(!» 

(21-1) 

Communication Strategies and General Matrix Algorithms Chap. 21 

t-o k=l I--j ~ j~ 

i~O o 1 2 3 4 , 6 7 
i-1 o 1 2 3 4 , 6 7 fold . 
2--7 1- O I 2 3 4 , 6 7 
i~3 o I 2 3 4 , 6 7 

k=O k-t f--j ~ j~ 

transfer -- 

i=O o 1 2 3 4 , 6 7 
i-l o 1 2 3 4 , 6 7 addvec i=2 o I 2 3 4 , 6 7 ~ i=3 o 1 2 3 4 , 6 7 

i=O 
Bll11234567 i=1 a b c d e f 9 h crystal- i-2 a p y 8 £ .X 1( ~ i=3 A B C D E F G H 

Any Configuration 

i=O i=l i=2 i-3 

Any Distribution of~~~~ ~~t6l~ 

~ = identity element 

F1gure 21-2 A summary for the case of four nodes of the transforn1ations described in this and the following chapter. The indices i,j , and k are described in Sec. 21-3.1. 

where F is a functional and f and 9 functions. The variables ,! and y may both be appropriate for decomposition over the nodes of a concurrent computer. As an example, consider the problem of Chap. 19. Dropping the vector notation, and letting y label the distinct sums and x designate the components of each sum, the problem of Chap. 19 can be written: 

S(y) = r.A(Y,x)B(y,x) (21-2) x 
We start with the range of the variable x fully decomposed and with each node containing all of the y values. However, after the summation stage of the algorithm (called/old in Sec. 21-3.3), the variable y is fully decomposed. Problems of this type are intrinsically "Iong-range" in the sense that data must travel allover the cube. Contrast 

Sec.21-2 Cube Geodesics 

a) 

b) 

401 

Geodesics from n1to n2 

... 

~ 

j 

One hop 

n2 

~ 

Two hops 

Three hops 

Flgure 21.1 Sample geodesics for ale, two, and three hop distances between nodes nl and n2" 

21.3 Homogeneous Communication Algorithms 

21-3.1 Data Structures 
Many of the crystalline utilities can be usefully applied to convolutions of the gen- eral mathematical fonD : 

g~)=F~.f(!» 

(21-1) 

402 

k=O k=l f--j ~ j--+ 

i=O 
-I 2 3 4 , 6 7 
i=1 o I 2 3 4 , 6 7 fold 
. 
2~ 1= O 1 2 3 4 , 6 7 
i=3 o I 2 3 4 , 6 7 

Communication Strategies and General Matrix Algorithms Chap. 21 

k=O k=1 f---j ~ j~ 

i=O O I 2 3 4 5 6 7 
i-l O I 2 3 4 5 6 7 addvec i=2 O1234567~ i=3 O I 2 3 4 5 6 7 

i=O i=l i=2 i=3 

-aaA4£eE ~2ycC6XgG reverse 1 bBS.fF 
36dD71ChH 

transfer ~ 

i~O [OIII213141516l7J i=1 I 
Iscatter 
. 
2~ 
1~ i=3 

Any Configuration 

i=O i~l i=2 i=3 

Any Distribution of~~~~ ~~~~ 

~ = identity element 

Figure 21.2 A summary for the case of four nodes of the transfonnations described in this and the following chapter. The indices i. j. and k are described in Sec. 21-3.1. 

where F is a functional and f and 9 functions. The variables ,! and y may both be appropriate for decomposition over the nodes of a concurrent computer. As an example, consider the problem of Chap. 19. Dropping the vector notation, and letting y label the distinct sums and x designate the components of each sum, the problem of Chap. 19 can be written: 

S(y)=~A(y,x)B(y,x) (21-2) % 
We start with the range of the variable x fully decomposed and with each node containing all of the y values. However, after the summation stage of the algorithm (called/old in Sec. 21-3.3), the variable y is fully decomposed. Problems of this type are intrinsically "Iong-range" in the sense that data must travel allover the cube. Contrast 

Sec.21-3 Homogeneous Comm1mication Algorithms 

403 

this example with the standard wave equation problem of Chap. 5. In that case. y and :! are members of the same domain. where f (,!) is the displacement at time t. and-g Cy ) is the same field at time t+dt .In the wave equation problem and. in fact. in the majority of problems in this book. one is not faced with two variables with inconsistent decomposi- tions. Usually f (,!) and 9 Cy ) naturally share the same decomposition. As a result. we usually only need nearest-neighbor rather than long-range communication. As we will see in the next subsection. the FFf of Chap. II is a good example of an exception. in which:! and y actually do have different decompositions. 
For the-most general problem class described by Eq. (21-1). it is natural to label the data (e.g.. the A .B Cy .x) in Eq. (21-2» by three indices i .j. and k. The integer i runs from O to N-I and labels the subsets into which:! is decomposed; j runs from O to N-I and labels the subsets into which y is to be decomposed. As usual. N is the number of nodes. The index k supplies any necessary extra labelings. We will typically ignore the k label in discussion since it will usually be used to index only an inner loop of no immediate interest. The situation is illustrated in Fig. 21-3 for the caseN=8. Each node i holds the data structure labeled by j and k. We will allow degenerate cases in which j does not assume all N values; in general we let j run from O to J-I. 

~ 

6 
1Internal 5 variable j running over sal 4 values as i 

We are now prepared to describe various transformations in the i-j plane which comprise the crystalline communication utilities which are the subject of this chapter. 

21-3.2 index 

This utility swaps the indices i and j and is illustrated in Fig. 21-4. To take a par- ticular example, after a can to index, the original contents of node 3 are spread over the whole hypercube, and stored in the location j=3 of each node. Figure 21-5 shows the paths taken by this data with those going to node 3 highlighted. The appropriateness of the name index for this utility can be appreciated by considering a database example. 

404 

Communication Strategies and General Matrix Algorithms Chap. 21 

Using the fonnalism of Sec. 21-3.1, let,! designate the records in a database and y represent the entries in an index to this database. Consider the case in which,! is decoffi- posed and the contributions to the database index are calculated concurrently. Each node contains contributions to the complete index that need to be assembled together. This assembly action is precisely that perfonned by index, which decomposes the actual index over the hypercube and yields a data structure in which all contributions to a particular index entry are held in a single node. 

INDEX 

7 

~~ I :::! I I I I~Remaining 
indices k 
01234567 Node Number i 
Figure 21.4 All illustration of index for the data structure introduced in Fig. 21-3. The routine index reflects across the diagooal in the i, j two-dimensional plane. 

\ 

:2f---l- 

~( ...; 

x 

x 

Reflect across diagonal 

The routine index uses essentially the same algorithm as addvec in Chap. 19 and is illustrated in Code 21-1. Note that as in Chap. 19, we systematically run through each channel starting with 0. This procedure is characteristic of the method of ordered cube geodesics. 
Note that index is well-defined even if j does not assume the full number of values N. Load imbalanced cases will not be treated here; a full discussion will be found in [Fox 86b] and follows reasoning similar to that already given for addvec in Chap. 19. 
An interesting application of index is to a revised concurrent FFf algorithm, as illustrated in Fig. 21-6. This problem represents a convolution of the type given in Eq. (21-1). Consider a 32-node decomposition of a 212 = 4096 point one-dimensional FFf. As described in Chap. 11, the basic sequential algorithm involves twelve steps, one for each of the binary digits in the representation of a typical ,! or y value. In the concurrent FFf problem of Chap. 11, one needs nearest-neighbor commUnication (in a hypercube) for the first five steps, while the remainder involve no communication. In the revised algorithm, as described in [Fox 86c], index is used to swap first the top and bottom five digits in the decomposition of,! and y; this operation is shown in Fig. 21-6. Now the first five steps of the conventional FFf are local to each node of the hypercube and may 

406 

Communication Strategies and General Matrix Algorithms Chap. 21 

be performed without communication. We now re-apply index and complete the last seven steps as before. The relative performance of the two algorithms (Chaps. 11 and 21) is analyzed in [Fox 86b], showing that the communication burdens are identical, while the revised algorithm is actually better because it avoids the calculational load imbalance present in the Chap. 11 algorithm. Defining the ratio, 't = tco".,../tca/c , we have: 
f = .!. -1 = ~ logN (FFT-Chap.ll) (21-3) C e 5 log(Nn) 
2't lo"M . 
= -.vb' .(revzsed FFT -Chap. 21) (21-4) 5 log(Nn) 

a) INDEX for the FFf 

b) GLOBAL REVERSE 

c) DEGENERATE INDEX 

-?> 

d) DEGENERATE GLOBAL REVERSE 

-?-i2IiJli41 
1'1gure 21-6(a) The use of index in the concurrent FFr algorithm given in Sec. 21-3.2. This is compared with global-reverse in (b) while (c) and (d) illustrate degenerate cases. 

It seems remarkable that trading nearest-neighbor communication in the Chap. 11 algorithm for long-distance communication in index yields a net gain in efficiency! 

Sec.21'-3 Homogeneous Commmlication Algorithms 

An7 

The index algorithm is very similar to that of global-reverse introduced in Chap. 11 to reverse the bits in the decomposition of y. In fact, when j assumes all N values, (1 =N), index and global-reverse differ by onlya simple permutation in the local data structure. If I is smaller than N then the two algorithms are completely distinct. The use of index in matrix algorithms will be described in Sec. 21-4.2 and 21-5. 

21-3.3 fold and expand 
Expanding upon our general set of utilities, we will find it useful to divide addvec as described in Chap. 19, into two separate functions. The first,fold, is illustrated in Fig. 21- 7 and combines all data for a given location onto the diagonal of Fig. 21-3 (in the node j =i ). The routine fold can be used with any operation -in Chap. 19 we used simple addition of the partial scalar products of two vectors. The routine expand is in some sense the inverse of fold, since it takes data stored on the diagonal and broadcasts them horizontally in Fig. 21-3. As illustrated in Fig. 21-8, for a given j all nodes hold the same data values after calling expand. In Chap. 19 fold was used to accumulate scalar products in a unique node and then the expand algorithm was used to ensure that all nodes contained the full set of scalar products. We need not repeat the algorithms for fold and expand given in Codes 19-1 and 19-2. However, we note that their communica- tion paths are complementary to one another and that they involve the same cycling through channels (characteristic of the method of ordered cube geodesics) shown for index in Code 21-1. The performance offold and expand has been presented previously in Fig. 19-5, in which the parameter p is the number of values taken by the indices j and k in the notation of Sec. 21-3.1. The routine expand can be considered as a special case of concat, introduced in Chap. 14, where each node contains the same amount of infor- mation. This simplification allows an efficient implementation on the hypercube. 

FOLD 

3 t ..7 ~.. .. 2 

~ol ~~... Rem~~n~ng""J/ I III; II I 
indiceskO 1 2 3 4 S 6 7 Node Number i 

F1gure 21.7 An illustratioo off old for the data structure introduced in Fig. 21-3, 

408 

Communication Strategies and General Matrix Algorithms Chap. 21 

EXPAND 

7 
fj 6 Copy from diagonal. 
5 All nodes have 
identical entries for 4 
h. 
eac J. 

~ 

!~ I I I I' j I 
01234567 Node Number j 
F1gure 21.8 An illustration of expand for the data strocture introduced in Fig. 21-3. 

Re~ 
indices k 

21-3.4 transpose and Matrix- Vector Multiplication 
We next illustrate the use offold and expand in matrix algorithms by discussing the problem of matrix-vector multiplication: 

y =A!. (21-5) 
where y, and ,! are vectors of length M , and A is an M by M matrix [Fox 86b, 86c, 86e] . We consider the case in which A is a full matrix (few zero entries). We will employ the same decomposition of A as in Chap. 10, with a square (111= M IN) x 111 submatrix stored in each node of the hypercube. We can write Eq. (21-5) in block form: 
yr=r.A~vi J (21-6) vJ 
where 1.1 and v are block indices running from 0 to N -1. The vectors yl1 and i v are 111- dimensional vectors obtained by dividing l: and ,! respectively into N separate sub- vectors: 

xmv+j =x; (21-7) 
y.-yA 11 mJL+i -, 
Figure 2l-9(a) illusb"ates the storage of both A and,! .In the two-dimensional decomposition of processors, the nodes are labeled by the index pair ~, v). Let us sup- pose we start with the initial state shown in Fig. 2l-9(a), where the node ~,v) holds the block matrix A JLv and all elements of x v. 
The algorithm for matrix-vector multiplication consists of four stages. In the first, shown in Fig. 21-9(a), we form: 

Sec.21-3 Homogeneous CommWlication Algorithms 

409 

Matrix-vector Decomposition (16-node) 1 ( a) Stage I: Local Calculation 
11 =0 11=1 11 =2 11 =3 

f-L =0 

1".1", 1".1".1 ~ 

1"0 ix, ix,lx, i 

~~I.~l.~ 

f'.lx.lx.lx, 

1..1..1..01."1 

~I 

my(:~y!' 
yr 

jx,lxsiX.lx,j 

l~rL 

fL = 

j.L =2 1 
1(0 
I(0 '0 i~ 

Ix.i x, I ~ 

~loJX11!1 

}.I. =31~1 

~i 
~I f1. =0 

~I 

(x.fi;-rx;j .71 I 

r"r"IX,4IX,,11 

~I 

}.I. 

E;EJ 

J-L =2 

11 fL =3 

=a 

!I= 

1/=2 

II =3 

Figure 21-9(a) The first of four stages of matrix-vector multiplication discussed in the text and illustrated for a 16 x 16 matrix A decomposed on a 16 node hypercube. 

yI1<V)=r.AijvxjV (21-8a) j 
with no summation over the index v. The y I1<v) are partial contributions to the vector y IL . This stage is local to each processor and is load-balanced up to differences of at most one row or column in the size of the blocks in each processor. 
We now usefold as illustrated in Fig. 21-9(b) in each subcube formed by a row of processors. This step completes the sum over v: 
yIL= r.yI1<V) (21-8b) v 

Next, we call expand over the same row of subcubes in order to broadcast y IL , so that each processor in a given row holds the vector y IL .This step is shown in Fig.21-9(c). The last step depends on the particular way one wants to use Eq. 21-5. Some applications may involve iterating y = ~ a number of times, replacing :! at each iteration with the new ~ on the right-haild side. In such a case, the implementation is 

410 

Communication Strategies and General Matrix Algorithms Chap. 21 

Matrix-vector Decomposition (16-node) 1(b) Stage ll: Fold 
II =0 11=1 II =2 II =3 

Figure 21.9(b) The secood of four stages of matrix-vector multiplication discussed in the text and illustrated for a 16 x 16 matrix A decomposed on a 16 node hypercube. 
completed by ensuring that y is stored in exactly the same manner as :! .This is achieved by adding a stage before the-application of expand, which b'ansposes the storage of y 11 to lie along columns rather than rows [Fox 86b, 86c, Ho 86b]. This introduces the next util- ity to be discussed. The routine transpose transposes matrix data stored in a two- dimensional mesh. After fold the vector yl1 is distributed so that it forms a matrix yl1v, as illustrated in Fig.21-9(b). The matrix is initially decomposed so that subblock yl1v is stored in processor number 1.1 + vN. After application of transpose yl1V is stored in pro- cessor v + 1.1.jN. The algorid1m for d1is is given in Code 21-2. We use the method of ordered cube geodesics in transpose with d1e usual cycle over communication channels, although in this case it is necessary to run over only half of d1em. 
Not surprising is the form of the communication overhead f c' which is given via Eq. (5-5) as: 

fc = .!. -I = .!. £ A m 

~ 
]2~ 

tcomm tcalc 

(21-9) 

Sec. 21-3 Homogeneous Communication Algorithms 

411 

Matrix-vector Decomposition (16-node) 1 ( c) Stage ill: Expand in Rows 
II =0 II =1 II =2 II =3 

FIgure 21-9(c) The third of four stages of matrix-vector multiplicatioo discussed in the text and illustrated for a 16 x 16 matrix A decomposed on a 16 node hypercube. 

with m= ..fn representing the square root of the number of matrix elements stored in each node. This efficiency formula has the usual form for an algorithm of system dimension d = 2, already seen in Chaps. 10 and 20 for other matrix problems. 
Note particularly the timing of transpose in Table 21-1. The ratio logzN/..[N can be thought of as the ratio of distances in a hypercube to those in a two-dimensional mesh. This ratio is characteristic of the underlying hypercube. It would be replaced by a number of order unity on a two-dimensional mesh architecture in which the transpose routine would be required to transmit data through a distance of order ..[N hops. Figure 21-10 shows the actual speedup obtained for the multiplication of a 16 x 16 matrix by a vector, as a function of the number of nodes, N. 

412 

Communication Strategies and General Matrix Algorithms Chap. 21 

Code 21-2 The transpose Algorithm for a d-Dimensional Hypercube 

Action 

(a) Local Block Calculation 

(b)fold 

(c) transpose 

(d) expand 

Time 

(2Jn2-m)tcalC 

[l-l/N]m(tcalc + tcanm ) 

.~t m~ 

comm 

(l-l/-1N)mtc~ 

Table 21-1 Timing for Steps of Concurrent Matrix-vector Multiplication 21-4 Inhomogeneous Communication Algorithms 
We will complete our discussion of the auxiliary crysmlline communication rou- 
tines by considering problems with intrinsic imbalance. Such problems might include the use of broadcast to send information from one to all other nodes. The routine scatter will be used to distribute information in one node among all others, while transfer will send information from one node to a single destination node. 

21-4.1 broadcast and forest 
We have already encountered the need for broadcast in Chaps. 10 and 20, and have commented on its implementation as a tree and pipe algorithm there and in Sec. 14-5.3. The action of broadcast on the data structure of Sec. 21-3.1 is shown in Fig. 21-11. In [Fox 85b, Fox 86b and Ho 86a], a detailed discussion of the optimal algorithms for broadcasting on machines without special hardware assistance is given. The optimal 

Sec.21-4 Inhomogeneous Comm\mication Algorithms 

413 

50 

20 
Speedup 10 

16 by 16 Matrix-Vector Multiplication 

t- O 

~--~ 
Direct Decomposition r-- 
't~ 1.5 

5-1 

~ 
10 100 1000 Number of Processors N 
Flgure 21-10 The speedup S as function of dte number of nodes N for a 16 by 16 matrix A. We show results for both a perfect hypercube ('t = t CDmm = 0) and a machine similar to dte Caltech Malic L n widt 't =1.5. 
algorithm depends greatly upon the ensemble hardware so it is not appropriate to give a detailed discussion here. On the Caltech Mark I and II hypercubes, the best broadcast algorithm depends on the value of m, the number of words to be sent, and N, the number of nodes. In Fig. 21-12 three "phases" are illustrated with different algorithms dominat- ing in different parameter ranges. The pipe is typically the best algorithm except for short messages. If the hardware imposes significant startup time for messages, as shown in Fig. 4-2 for the initial INTEL iPSC system, then the pipe is no longer attractive, since it requires the message to be broken into several segments, each having its own startup time. Since broadcasts over cubes and subcubes are important. we expect that future hypercubes will include special purpose hardware to support broadcasts. As a result of such modification, the analysis of Fig. 21-12 would be significantly changed. 
In Fig. 21-12, one of the phases is dominated by an algorithm referred to asforest (see [Fox 86b] and personal communication in [Ho 86a]). The routineforest distributes broadcast trees over the hypercube in a way that incurs less startup time than pipe. How- ever, the cost per word is lower for the pipe routine, which therefore becomes preferable as the message length increases. The routine forest is a generalization of the simple broadcast algorithm described in Chap. 14 that uses a simple tree decomposed over the hypercube. The routine forest takes the several trees formed by different parts (words/packets) of the message to be broadcast and obtains better load balancing by overlapping the different trees. The third phase, scatter followed by expand, uses an algorithm which is described in the next subsection. 

414 

Communication Strategies and General Matrix Algoritluns Chap. 21 

BROADCAST 

Remaini~ indicesk 111111 234567 
Node Number i 
FIgure 21.11 An illustration of broadcast for the data stmcture introduced in Fig. 21-3. 

~ 

-+--+ O 1 

1if 3 
Number of Points 
2 
-=::. pipe phase ====- 
- 
~ 
~ 
I 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Dimension of Cube 
Figure 21.U A sketch of the domain spanned by 111, the number of words to be broad- cast, and d = 1ogzN ' the dimensioo of the hypercube. We show three regioos where dif- ferent algorithms dominate. These are the pipe discussed in Sec. 20-4; scalIer followed by expand discussed in Sec. 21-4.2 and 21-3.3;foresl discussed in Sec 21-4.1. 

Broadcast Phase Diagram 

100 

Sec.21-4 Inhomogeneous Communication Algorithms 

415 

21-4.2 scatter An important utility, scatter, is illustrated in Fig. 21-13. The routine scatter takes data stored in a single node and distributes it homogeneously throughout a subcube, which may be the full hypercube, forming the set of destination nodes. The correspond- ing algorithms given in [Fox 86b ] use a combination of ordered cube geodesics and pipes. The scatter algorithm is outlined in pseudocode in Code 21-3. 

SCATrER 

v 
+-+ 
O 1 2 3 4 
Node Number i 
l'lgure 21.13 An illustration of scatter for the data structure introduced in Fig. 21-3. 
An obvious application of scatter is to distribute a data set stored in a given node, 
for example node 0, or perhaps in a disk drive attached directly or indirectly to the node. The effect of scatter is to divide this data set into N parts and deliver equal portions to each node. Figure 21-14 illustrates another possible use of scatter. It will be recalled from Chap. 20 that many sources of load imbalance are associated with banded matrix LU decomposition. In particular, finding the multipliers from the first column of the current window occupied only the ~ processors in the first column of the node array. Figure 21-14 shows how we can apply scatter to subcubes consisting of rows of nodes and equalize the load on the nodes. Unfortunately, this method is not immediately practi- cal as stated because scatter takes time of order mtcomm , where we distribute m numbers as m/~ packets in ~ nodes. The load imbalance that scatter is intended to remove involves a calculational cost of order mtca/c .Thus, in the usual situation where tcomm ?; tca/c , no net gain is realized. This approach may be improved, however, by simultaneously applying the idea of Fig. 21-14 to the ~ columns involved as the "first column" in the next ~ steps. This actually corresponds to applying the index algorithm described in Sec. 21-3.2. The routine index can be thought of as the homogeneous ver- sion of scatter, just as expand is the homogeneous version of broadcast. The index cal- culation requires a time of order mlogzN tcomm , which is to be compared with the calcula- tional imbalance given by m~ tca/c for the total of ~ steps. The use of index formally reduces the calculational load imbalance to zero in the limit of large ~. Note that the reduction factor of logzN /~ is characteristic of the hypercube, and that these special 

j 

"" 

)/ 

Rerna~ of indic~~1ng""l 

5 

+ 6 

416 

Communication Strategies and General Matrix Algorithms Chap. 21 

techniques do not apply to simpler architectures such as the two-dimensional mesh. Another interesting application of index to a problem of this general type is the conver- sion of matrices from the local square subblock to the scattered square subblock decom- position (i.e.. from the decomposition of Fig. 20-3(a) to that of 20-3 (b)). 
Returning to Fig. 21-12 for broadcast. we note that for a single transmitted word or a packet of minimum size. scatter has no action. and expand reduces to thc standard tree algorithm used in Chap. 14 for broadcast. Thus. the low m phase of broadcast is dom- inated by a simple generalization of the tree algorithm; scatter followed by expand. 

proc-begin scatter ( scatter vector in processor O among other processors ) declare-buf V; ( a vector; in processor O contains data to be scattered ) declare-int M; ( number of items in vector V) [ let N = the number of processors ] 
[ let procnum be the processor number] [ nloop = minimum of N and M ] 
if-begin( this processor is processor O ) then 
[ divide V into blocks of N elements (last block may contain less than N) ] for -begin ( loop index, i, from nloop -1 to 1 ) 
[ m = number of blocks having an il11 data element ] [ CL = bit number of lowest 1-bit of i ] 
proc-caU vwrite ( send il11 item from each of the m blocks on channel CL ) for end 
else-if (procnum < M ) then 
[ CH = bit number of highest 1-bit of procnum ] c +1 
[ npipes = N 12 H -1 = number of forwarding pipes ] 
for-begin (i from 1 to npipes) 
[ CL = bit number of lowest 1-bit of i ] 
proc-caU vread ( receive data on channel CH and forward on channel CH + CL + 1 ) 

for end 
proc-call vread ( receive data items on channel CH and store 
in V(procnum), V(procnum + N). etc) 

if end proc-end 

Code 21-3 The scatter Algorithm 

21-4.3 transfer The final utility, transfer, is perhaps less generally useful, at least in application to matrix algorithms. As illustrated in Fig.21-15, it transfers information between two given nodes. This is needed, as shown in Sec. 20-5, when employing pivoting in LU decomposition; transfer is used to move the original first row into the position of the chosen pivot row. The latter needs to be broadcast as described in Sec. 20-5, so transfer is not required to put the pivot row in the place of the first. 

Sec.21-4 Inhomogeneous Comm\mication Algorithms 

Before 

417 

After SCA TfER 

--7! 

, 
-~"' I.' »;",! i:t!; 

i!;r1i 

i Imbalanced Balanced 
First column 
FIgure 21.14 A possible use of scatter to distribute the first column over a matrix unifonnly over an nodes. 

TRANSFER 

j 

, I I I I 
01234567 Node Number i 
Figure 21-15 An illustration of transfer forthe data structure introduced in Fig. 21-3. 

~ 

Rem~ 
indices k 

.E; "' 
t.Move information between two distinct nodes 

Communication Strategies and General Matrix Algorithms Chap. 21 

proc-begin transfer ( transfer a vector, V, of m elements between processors S and D ) define-function CL (i) = bit number of lowest I-bit of i define-function cH(i) = bit number of highest I-bit of i [ let route = bit-wise XOR of S and D ] 
[ let r = bit-wise XOR of S and procnum ] 
if-begin ( r = 0, i.e., this is source processor) then 
for-begin ( all elements, Vi, of V, i = 0, I,..., M -1 ) if-begin (i is even ) then 
proc-caU cwrite (send Vi on channel CL (route» else-if (i is odd) then 
proc-caU cwrite (send Vi on channel CH (route» if end for end 
else-if (r = route, i.e., this is destination processor) then 
for-begin ( all elements, Vi, ofV , i = 0, I,..., M -I ) if-begin (i is even ) then 
proc-caU cread ( receive Vi on channel CH (route» else-if (i is odd ) then 
proc-caU cread ( receive Vi on channel CL (route) ) if end for end 
else-if (r = lowest (CH (r ) + I) bits of route ) then 
for -begin (all even-numbered elements of V) 
proc-caU cread ( receive element on channel CH (r ) and forward on channel CL (route -r) ) 

for end 
else-if ( r = highest (d-CL (r » bits of route) then for -begin (all odd-nUmbered elements of V) 
proc-call cread (receive element on channel CL (r ) and forward on channel CH (route -r) ) 

for end if end proc-end 

Code 21.4 The transfer Algorithm 

The routine transfer can be implemented optimally using two geodesic pipes. In Sec. 20-4, we discussed the best pipe algorithm for broadcast, which we can call a linear (split) pipe. The pipe connects all nodes in a subcube by decomposing them into a one- dimensional line or ring. If we just wish to connect two nodes, then rather than connect- ing the nodes by a one-dimensional route, it is better to use a geodesic, and we obtain the dual geodesic pipe outlined in Code 21-4. In an optimized version of this algorithm, pieces of data alternate between the two geodesic pipes. 
The optimization used in Code 21-4 is hardware-dependent, and the algorithm assumes a machine for which reading and writing by the intermediate nodes of the pipe 

Sec.21-4 Inhomogeneous Communication Algorithms 

419 

take roughly equal time and cannot be overlapped. If message latency is high a better implementation of transfer would send all the even and odd elements of the V through their respective pipes together using vread and v write, rather than sending one element at a time with cread and cwrite. 

21-5 Comments on General Matrix Algorithms 
A great deal of work on matrix algorithms for the hypercube has already been per- formed; in particular by groups at Caltech [Fox 85b], Yale [Johnsson 85], Oak Ridge [Heath 87] and INTEL [Moler 86]. In this section we extend the results of Chaps. 10 and 20 to other algorithms, and exploit the communication utilities described earlier in this chapter. 
We can classify matrix algorithms into three broad classes: (1) those used for full matrices with few zero elements, (2) special techniques that exploit features such as banded or tri-diagonal structure, and (3) iterative algorithms. The latter tend to be used for very sparse matrices since they require less storage and are often more efficient than the direct solution methods. The algorithms used in Chaps. 7 and 8 of this book are of the iterative type. 
Let us first consider full matrices and banded matrices where the bandwidth is large compared with N, where N is the number of nodes. These classes of matrices require techniques which are quite similar to those discussed in this book. Three important algo- rithms are: 
Matrix-Vector Multiplication, which is often used as a part of an iterative algo- rithm. It has been discussed in Sec. 21-3.4. 
Matrix Inversion and its generalization, LU decomposition, for solving sets of equations. Some of the relevant algorithms have been described in Chap. 20. 
Matrix Eigenvalue and Eigenvector determination, which typically use a technique known as Householder's method to reduce a general symmetric matrix to tri- diagonal form. The fast QR algorithm is then used to diagonalize the resultant tri- diagonal matrix. 
We now discuss some of the issues that come up for these algorithms. 1. The decomposition of the matrix into square subblocks is again optimal and to be preferred over row, column, or rectangular decompositions on hypercube hardware without a large message startup time. 
2. The scattered decomposition is crucial in both LU decomposition and Householder's method since as these algorithms progress, rows and columns are eliminated or deactivated. We can see this effect in Fig. 20-18 because the full matrix algorithms are structurally similar to the end steps of the banded algorithms. As shown in the examples of Sec. 20-5, approximate load balance is maintained even after pivoting, as long as one uses the scattered decomposition. 
3. In [Hipes 86] a particular modification of Gaussian elimination, called the Gauss- Jordan method, is shown to be optimal for full matrix inversion. It naturally preserves load-balance. In Fig. 21-16 we measure the performance of this algo- rithm on the Caltech/JPL Mark II hypercube. For instance, we find an efficiency of 

420 

4. 

5. 

Communication Strategies and General Matrix Algorithms Chap. 21 

Gauss Jordan Inversion With Pivoting On The Mark II Hypercube 0.8 f~ :::=:=:=:=:=: 

.0 

Efficiency 

0.6 

1/; 

0.4 

r 

/ 

0 50 100 150 200 250 300 Global Matrix Dimension 
FIgure 21.16 Measured efficiencies for the Gauss Jordan matrix inversioo algorithm with pivoting. Results are shown on 2. 4. 8. 16. and 32 node hypercubes of the CaltecbflPL Mark II design. The algorithm uses 64-bit arithmetic and the results are plotted against the matrix order M. Each node has a grain size 11 =M2IN. Even hyper- cubes use a square and odd hypercubes a rectangular dec<Xnposition. 
over 80% on a 16-node hypercube with a grain size corresponding to a 25 x 25 submatrix in each node, where the full matrix is 100 x 100. Figure 21-16 shows that the measured efficiency is generally high and indicates the suitability of the hypercube for these problems. 
The straightforward subblock d~mposition needs another generalization when dealing with symmetric matrices, which are usually the type for which eigenvalues are needed. A square matrix naturally decomposes onto a concurrent processor with the topology of a two-dimensional mesh. One does not need a new triangular topology to deal with symmetric matrices. Figure 21-17 shows one of many simi- lar ways that one may decompose a symmetric matrix (or effectively any triangular structure) onto a two-dimensional mesh [Fox 84d]. There is a load imbalance, which is of order 1/-rn-, where -rn- is approximately the number of extra matrix ele- ments in diagonal compared to off-diagonal nodes. For clarity, we show in Fig. 21-17 the local decomposition that is analogous to that of Fig. 2O-3(a) for unsym- metric matrices. There is a natural scattered version of this decomposition. 
The last stage of the eigenvalue determination, the QR step, only involves a small number of matrix elements since the original M (M + 1)/2 independent elements have been reduced to 2M -1 in the tri-diagonalization step. It is interesting that this reduction in the number of variables, which aids the sequential code, makes the concurrent implementation much harder, since there are not enough matrix ele- ments to spread over the nodes. A recent observation by Furmanski has led to a good concurrent algorithm for the QR step [Furmanski 86]. Recent work 

p 

0.2 

0 

D : 1 row 2 cols + : 2 rows 2 cols O : 2 rows 4cols x :4rows 4cols " .4 rows 8 cols 

Sec. 21-5 Comments on General Matrix Algorithms 

421 

6. 

7. 

[Dongarra 85] has also led to new recursive algorithm called Cuppens Method that may be better than the QR method on sequential computers, and is much easier in its concurrent implementation [Ipsen 87] . 
We always find the typical two-dimensional efficiency form for these matrix algo- rithms to be: 

.! -1- ~ (21-10) £ ;In 
where n is the number of mabix elements stored in each processor and: 
tcQmm const = a+b ( -= t) tcalc 
In Eq. (21-11), the nonzero a is due to load-imbalance and was not present in the mabix multiplication algorithm of Chap. 10. However, we see the effect of load- imbalance in Eqs. (20-12, 13) for LU decomposition. One might expect b to depend logarithmically on N for eigenvalue determination because of the need to calculate many scalar products. However, the new algorithm for global scalar pro- ducts described in Chap. 19 leads to values for a and b that are asymptotically independent of N . 

const 

(21~11) 

Figure 21.17 A local decomposition appropriate for a symmetric matrix. The nodes are labelled as appropriate for a two dimensiooal mesh and would be straightforwardly changed 00 a hypercube. 

422 

8. 

9. 

10. 

Communication Strategies and General Matrix Algorithms Chap. 21 

Sparse matrices are surprisingly interesting even for the simple matrix-vector multiplication algorithm described for full matrices in Sec. 21-3.4. In the case of sparse matrices, one is faced with a load-balancing problem that can be tackled with the sophisticated optimization techniques outlined in Chaps. 3 and 17 [Fox 86d]. 
The utilities described in Sec. 21-2, 3, and 4 can be used directly to improve both the methods in Chap. 20 and those of similar algorithms. In particular, using the tricks outlined in Sec. 21-4.2, one can often use index to remove the term in a in Eq. (21-11), which is converted into a term proportional to (logzN I~) tcomm .The results of [Furmanski 86] are presented in Table 21-2. The communication overheads are never large and are compar- able, for similar size data sets, to those obtained in simple nearest neighbor problems, such as those in Chaps. 7 and 8. 
Finally, we note that a key problem in the field of linear programming is an optimization technique involving the solution of many linear inequalities. These can be formulated as a matrix of inequalities. Important strategies for solving these inequalities are the simplex method and an improvement due to Karmarkar [Dantzig 68, Karmarkar 84]. Both of these methods can be implemented efficiently on a hypercube using similar methods to those we have already discussed for other matrix algorithms [Felten 86b]. 

Sec.21-5 Comments on General Matrix Algorithms 

423 

Algorithm 

MultipJication -Optimal (Chap. 10) 

LU Decomposition + forward reduction 

Naive Square Subblock 
with Scattered Decomposition and No Pivoting 

Optimal -No Pivoting 

Optimal -Pivoting 

Gauss-Jordan Inversion 

Optimal 

Householder Tridiagonalization 

Optimal 

Banded LU Decomposition + forward reduction 

Naive -No Pivoting but 
with Scattered Decomposition 

Optimal -No Pivoting 

Optimal- Pivoting 

Const (See Eq, 21-11) 

't 

9 3't -+- 4 2 

t[ 

t[ 

1 + 1.- IOgi¥ 
]2-W--'t 2 + 1.- IOgi¥ 
]4-W--'t 

1 + .!. IOg-JV 
]4~'t 

3 
[SlOglv' 
]4 3+sw- 't 

1 't+- 
2 

't[1 + .J::.-~ ]2{N 
't[2 + .J::.-~ 
]43{N 

Table 21-2 Communication Overheads in Optimal Matrix Algorithms 

22 

22- 

Introduction 
In the previous chapters of this volume we have seen that the crOS III communica- tion system can be used to solve a wide range of scientific and engineering problems on concurrent processors. crOS III is best-suited to problems in which the nature of the communication is known a priori, since provision can then be made for any necessary message forwarding. In this chapter we describe the crystal-router algorithm, which can be used to pass messages between arbitrary nodes of a hypercube when the need for com- munication, but not its nature, is known beforehand. The crystal-router is therefore use- ful in irregular problems in which the message traffic changes dynamically and/or involves message passing between non-nearest neighbors. Section 22-3 discusses exam- pIes of the use of the crystal-router in circuit simulation and the evaluation of non-binary FFfs. The crystal-accumulator, introduced in Sec. 22-4, is a modification of the crystal-router which allows messages destined for the same processor to be merged en route according to some user-specified function. 
In Sec. 22-5 we describe the implementation of a virtual concurrent processor, similar to that introduced in Chap. 4, known as the virtual machine loosely synchronous commzmication system (Vw..$CS). The crOS III communication routines are blocking. However, in the Vw..$CS this constraint is removed, pennitting non-blocking, loosely synchronous communication between arbitrary processors. Details are given of how to implement the VMLSCS by means of the crystal-router on a hypercube with a synchro- nous message passing protocol. 

22-2 crystal-router and General Message Traffic in the Crystalline Environment 
In this section we shall see that the crystalline environment is quite capable of deal- ing with rather arbitrary message traffic. The only essential requirement is that the com- munication calls occur in loose synchronization. In other words, the receiving nodes must know to expect some sort of data. The aim is to be able to transmit messages between arbitrary nodes of the hypercube when the size of the messages is not known a priori by the receiving node. All that should be required is a large enough buffer for incoming messages. It seems likely that such a general crystalline router will, in many circumstances, offer equivalent functionality and greater performance over more general message-passing (e.g., interrupt-based) operating systems. In particular, crystal-router is useful in irregular problems of the type discussed in Chaps. 8, 16, and 17, where most 

425 

426 

General Message Passing in a Loosely Synchronous Environment Chap. 22 

communication is local and the accommodation of occasional irregularities (both in size and destination) of messages is needed. 
Consider first the ordered cube geodesic method in the case that a given message must be sent from one node to several others. We have already encountered such cases in Chaps. 9, 10, and 19, where information was broadcast to all nodes of either the cube or a subcube. As shown in these chapters, it is usually not optimal to send one message for each destination node. Rather, it is better to construct a tl grand tour" such that the message takes a single path visiting each destination. The routing of multidestination messages is not an easy problem. In general, it involves the solution of a kind of traveling salesman problem, as discussed in Chap. 13, to minimize the path on the hypercube. This point illustrates the fact that it is not simple to extend the methods of Chap. 21 to more general optimal routing methods. Not only is the best route for a given message unknown, but in irregular problems, it is not always obvious how to balance the message traffic. In spite of these reservations, the following discussion presents a natural generali- zation of the methods of Chap. 21 in the form of the crystal-router algorithm. 
Consider a concurrent application whose stages of calculation and inteIprocessor communication are constrained to remain in loose synchronization by the periodic occurrence of expected communication. We will assume that our generic application's communication steps may be characterized as follows: 
At a given stage of the calculation. each node of the hypercube has a set of messages to be sent to other nodes. The destination nodes expect messages but know neither how many messages will arrive nor the identities of the ori- ginating nodes. 
This requirement is termed need-predictable in Chap. 16, and is characteristic of loosely synchronous problems. This situation is typical of the message traffic in many problems discussed in this book, but we consider particularly the cases in which the mes- sages are not just sent to nearest-neighbors but also to more distant nodes. Such irregu- lar, but time-synchronized, message traffic occurs in any simulation of an irregular sys- tem. Examples follow directly from the problems of Chaps. 7, 8, 16, and 17, when applied to irregular geometries. In addition, the applicability of crystal-router to simula- tions of circuits (Sec. 22-3), and neural networks is of considerable interest. 
It is useful to distinguish between two broad categories of the crystal-router algo- rithm. The dynamic crystal-router, as described above, routes messages of arbitrary length between arbitrary nodes of the hypercube. The generality of such an algorithm is useful in many applications in which the communication is need-predictable, but of an unknown nature. However, there are applications, such as the circuit simulator described in Sec. 22-3, in which the routing of messages is the same for every call to crystal-router. In such cases the communication overhead can be reduced by using the less general static crystal-router. Unlike the dynamic crystal-router, the static version does not determine ab initio whether each message must be sent over each channel every time crystal-router is called. Instead the static crystal-router has an initialization, or "dry-run", phase during which each node constructs a routing table, which determines how messages are to be routed on subsequent calls to crystal-router. The columns of the