L 

NASA Contractor Reptt 181626 

ICASE REPORT NO. 88-14 

.a 

ICASE 

DATA TRAFFIC REDUCTION SCHEMES FOR 

SPARSE CHOLESKY FACTORIZATIONS 

Vijay K. Naik 

Mertell L. Patrick ~ ~ ~ - 

(liASA-CR- 18 1656) CAT A TRAFFIC EECUCTIO N N88-2083 1 

SCfiElES FOR SPARS€ CHCLESRY ?LCIC612ATICNS 

Eical BegCKt (&$SA) 35 CSCL 12A 

Unclas 

G3/59 01352S8 

Contract No. NAS1-18107 

February 1988 

INSTITUTE FOR COMPUTER APPLICATIONS IN SCIENCE MID ElZINBEBfNG 

NASA Langley Research Canter, hpton, Vlrglnla 23665 

Operated by the Universities Space Research Aeeociation 

.. 
Data Traffic Reduction Schemes for Sparse Cholesky Factorizations 

Vijay K. Naik 

and 

Merrell L. Patrick 

Institute for Computer Applications in Science & Engineering 

NASA Langley Research Center, Hampton, VA 23665 

and 

Computer Science Dept, Duke University 

Durham, NC 27706. 

ABSTRACT 

Load distribution schemes are presented which minimize the total data traffic in 

the Cholesky factorization of dense and sparse, symmemc, positive definite matrices 

on multiprocessor systems with local and shared memory. The total data traffic in factoring an n x n sparse, symmetric, positive definite mamx representing an n-vertex 

regular 2-D grid graph using nu, a< 1, processors is shown to be O(n 2). It is O(nT), 

when no, a 2 I, processors are used. Under the conditions of uniform load distribution 

these results are shown to be asymptotically optimal. The schemes allow efficient use 

of up to O(n) processors before the total data traffic reaches the maximum value of 

O(nY). The partitioning employed within the scheme, allows a better utilization of the 

data accessed from shared memory than those of previously published methods. 

l+S 3 

3 

Kcrearch supported by the National Aeronautics and !3paa Administration under NASA mtnd No. NAS1-18107 while the 

authors were in residence at ICASE. NASA Langley Research Center. Hamp. VA 23665. 

i 
1 

1. Introduction 

Consider the problem of solving a system of linear equations 

Ax=b 

where A is an n x n sparse, symmetric, positive definite mamx, x is an n x 1 vector of 

variables, and b is an n x 1 vector of constants on a multiprocessor system with both 

shared and local memory. In this paper the total data traffic involved in factoring the 

matrix A using the column version of Cholesky decomposition is analyzed. The 

analysis is restricted to systems for 2-D regular grid graphs. The results developed 

here, however, can be applied to a wider class of problems.' 

In [2], George et al. have given a parallel sparse factorization scheme for local 

memory multiprocessor systems that has a total data traffic of O(n' +a logzn) using na 

processors. This result was improved to O(n'+O) in 

1 +E 

has a total data traffic of O(n 2, using na, as 1, 

[4]. In this paper a scheme that 

processors is presented; with na, 

3 

a 2 I, processors the resulting data traffic is O(n'). Although a multiprocessor system 

with both shared and local memory is assumed as the model of computation, the 

results developed here are equally applicable to systems with only local memories. 

In Section 2, some background work is discussed and the model of computation 

is specified. In Section 3, expressions for the total data traffic required in factoring 

dense matrices are developed. Using these results the total data traffic for sparse systems is analyzed in Section 4. Conclusions are given in the final section. 

' V. K. Naik and M. L. Patrick. Communicotion requirements ofparallel sparse Cholesky foctwizotions. m prcpantion. 
2 

2. Some Preliminaries 

The model of computation used in our analysis is that of a multiprocessor system 

with a two level memory hierarchy such that each processor has local memory and all 

processors have access to a shared memory. The access cost per nonzero element in 

the shared memory is assumed to be unity for any processor. The total number of 

shared memory ilccesses from the beginning to the end of an algorithm is defined as 

the communication requirement or datu fraffrc of that algorithm implemented on the 

multiprocessor system. 

, 

The n x n matrix A considered here is assumed to represent a system 

corresponding to a 2-D regular grid graph, the vertices of which are ordered using the 

nested dissection method [3]. Figure 1 shows an ordering scheme for a 7 x 7 grid. 

1. 7- 4- 43 22 20 25 

3. 8. 6- 44 24 29 2? 

29 9. 5. 45 23 30 26 

19 2@ 21- 46 40 42 42 

1D 16 13 43 33 33 39 

12 l? 13 40 33 30 36 

12 18 19 49 32 39 35 

Figure 1: A 7 x 7 grid with nested dissection ordering. 

.ONGINAL PAGE IS 

EOOR QUALITY 
3 

Such an ordering has an optimal sequential operation count and fill-in. 

The basic algebraic scheme used in this paper to factor A is the column version of 

the Cholesky decomposition method [6]. An outline of this algorithm for factoring the 

n x n matrix A is as follows. Let A =UT; uij E A and fij E L. 

for j = 1 until n do 

begin 

Initialize lij = uij, i = j, - . ,n. 

for k = 1 until j-1 do 

for i = j until n do 

1.. = I..-/. IJ IJ 1.k * 4,k ; 

l,j = .lisJ ; 

for k =j+1 until n do 

ly = 1kj 1 1jj; 

end 

In the above algorithm for clarity, the values of fij are shown separately from those of 

aij. In practice lij may overwrite uij. 

In the next two sections, parallel schemes for the above algorithm are developed 

for the distribution of work among processors of the multiprocessor system and their 

communication requirements are analyzed. The analysis assumes some familiarity 

with graph theory concepts related to matrix representations of systems of equations, in 

particular, the notion of vertices, edges, separators, subgraphs of a graph, and the 

correspondences between the vertices and the rows and columns of the matrix, the 

edges and the nonzero elements, and addition of edges and fill-in during the factoriza- 

. tion of A. For details on these concepts see [3] and [7] and the references therein. 

Finally, the standard notations "O", "o", and 'Q" are used to characterize the 

asymptotic growth rates of computation and communication requirements. Precise 
4 

definitions of these notations can be found in elementary textbooks on the analysis of 

algorithms, e.g. 191. 

3. Communication Requirements of a SPOCC Factorization Scheme for Dense 

Matrices 

First, a scheme based on the column version of the Cholesky decomposition [6] 

for factoring the lower triangular part of an m x m symmetric, positive definite, dense 

matrix using p processors is described. Without loss of generality, further assume that 

p = - (2 + r) where r is an integer. The lower triangular matrix is partitioned into p 

1 

2 

m 

r partitions such that all, except r partitions, are s x s square blocks, where s = -. The 

remaining r partitions which lie on the diagonal of the matrix are s x s triangular 

blocks. Each of these partitions is assigned to a unique processor. Initially, each processor reads the data for its partition from shared memory into its local memory. The 

r processors in charge of the partitions containing the left most s x s blocks of the 

matrix commence the computations of their part of the factorization. As soon as an 

element of the factor is computed, it is written into shared memory for access by other 

processors. As the necessary data becomes available, the remaining processors initiate 

computations on the blocks assigned to them. Each processor accesses from shared 

memory only those elements that are needed for local computation and each element is 

read no more than once by any processor. This parallel factorization scheme is termed 

as the submatrix-partition oriented column Cholesky (SPOCC) factorization scheme. 

Next the communication requirements of the above scheme are analyzed. In the 

following two theorems it shown that the total data traffic involved is O(m2 $) and that 
5 

this result is asymptotically optima First ConsiGcr the data traffic associated with a 

generic square block f shown in Figure 2. In that figure, the darkened area represents 

thc (lata clcnicnts that ilrc rcqtrircd for the con1puti1tions in block 1. 1.ct block / hc 

bounded by columns (i-1)s + 1 and is, and by rows (i-1)s + 1 and Js where 1 s i 5 r and 

1 sj 5; r. The following lemma provides bounds for the communication cost of block I. 

Let ak.1 be any element in block f. During the factorization this element undergoes 

a sequence of modifications before it arrives at a final value. The computations leading to the final value of this element require the final resultant values of certain other 

elements and its own initial value. For determining the communication requirements 

the actual value, either final or initial, of an element accessed is not important. Only 

an account of this access to the shared memory is necessary. To keep the discussion 

brief, in the analysis of the data traffic, the term appropriate or required value of an 

(i - I).s + 1 i.s 

Figure 2: The data traffic associated with block I. 
6 

element is used. An appropriute or required value of an element at a particular step 

of the computations is defined to be the value of that element that is required by the 

data dependencies of the algorithm at that step of computation. In the proof given 

below, al; represents all the elements in row 1 of the lower triangular part of the matrix 

under consideration. 

1 1 

2 2 LEMMA 1 : A total data trafic of (2 - -).2 + -as is necessary and sufJicient for factoring the square block I; it is the same for all square blocks bounded by the columns 

(i-l)s + 1 and i.s. The data traffic associated with a s x s triangular diagonal block 

1 1 bounded by the columns (i-1)s + 1 and i-s, is (i - -)a2 + -.s. 2 2 

PROOF: Consider the computations at any element in block f. For these computations the appropriate values of the elements a1,. and ak,i for all i I 1 are required; i.e., 

the data required are all the elements in row f and all the elements in row k on and to 

the left of column 1. Now ak,/ is any element in block f and hence, 0-1)s + 1 I k 5i.s 

and (i-1)s + 1 I 1 I i-s. Thus, the computations corresponding to any element in block f 

depend on all the elements in some row 1 of the factor, (i-1)s + 1 I 1 I i-s, and on the 

elements to the left of column i.s+l in some row k, (j-1)s + 1 5 k sis. Furthermore, no 

other information is needed in performing the factorization at that element. Therefore 

all the information in the above mentioned rows is sufficient for performing the computations at all the elements of block f. Conversely, the block f contains all the elements at,l such that 0-1)s + 1 5 k IJS and (i-1)s + 1 I f I i.s. The completion of the computations at these elements requires all the elements in rows 1 of the factor such that 

(i-I)s + 1 I f 5 i.s and all the elements to the left of column i.s+l in all the rows k such 
7 

that ~-1)s + 1 5 k <is. Thus, it is necessary as well as sufficient to read the values of 

these specific elements for completing the factorization of the square block I. A count 

of the number of such elements gives a bound on the data traffic involved in factoring 

block I. It should be noted that this count includes the elements in block I which are 

locally computed, but the initial value of each of these elements has to be read prior to 

performing computations at that point and thus the communication cost remains the 

same. 

Thus, the data traffic that is both necessary and sufficient to factor block I is 

3 1 

2 2 = 2(i-1)*2 + --s2 + --'s 

1 1 

2 2 = (2i - -)2 + -*s. 

The above bound on the data traffic is independent of j and so it is the same for all the 

square blocks within the columns (i-l)s+ 1 and i-s. 

Similarly, we can show that the data traffic associated with a s x s triangular diag1 1 

2 2 onal block bounded by the columns (i-1)s + 1 and i.s is, (i - -)a? + -.s. 

Using these results we get a bound on the total data traffic as stated in the following 

Theorem. 

THEOREM 1 The total data traffic involved in the SPOCC factorization of an m x m 

dense matrix using p processors is O(m2 6). 
8 

PROOF: The total data traffic associated with all the blocks bounded by columns 

(i-l)s -t 1 and i.s, 1 4 i 5 t, is given by, 

(r-i) x dura trufic associated with a square block 

+ data trufic associated with a triangular block bounded by the given columns. 

From Lemma 1 and the fact that there are t such column partitions, we get the total 

data traffic involved in factoring the m x m dense matrix using the SPOCC factorization as 

r . qr-w 2 '1 - -).? 2 + -3) 2 + ((i - -).P 2 + -s) . 1 1 1 

i=l 

Ignoring the lower order terms, the total data traffic is 

i=l 

1 m 

2 r 

Since p = - (3 + r) and s = -, the total data traffic 

0 

The above theorem gives an upper bound on 

is 0(m2@. 

the data traffic involved in factoring 

the m x m matrix using the submatrix assignment scheme. In the following theorem 

we prove that this result is optimal in the order of magnitude sense. Before giving the 

proof, we again consider the data dependencies involved in computing an element of 

the factor. Consider the computations at any element ui2 To compute the value of the 

corresponding element in the factor, values at all the elements in row j and the values 

of elements in column 1 through j of row i are needed. Thus, any off-diagonal ele- 
ment uij requires values of 29 elements for completing the computations and any diagonal element (when i = j) requires values of j elements. We now have the following two 

obvious observations: 

(1) The values at all the elements in any row i can be used to complete the computations at exactly one element, namely, the diagonal element ai,i. 

(2) If i and j are any two rows, such that i>j, then the values of the elements in 

these two rows can be used to complete computations at exactly one off-diagonal 

element uij. No other information is needed to complete the computations at that 

element. 

These observations lead to the following lemma on the number of elements at which 

the computations can be completed by acquiring information from exactly k rows. 

LEMMA 2 : If the required values of all the elements in any k rows are available then 

k2 k 

22 the computations of at most - + - elements can be completed; conversely, the compk 

22 pututions corresponding to any - + - elements depend on at least one element from 

each row of at least k distinct rows. 

PROOF: We prove the first part of the lemma by induction on the number of rows. 

When k = 1, + = 1, which is true from Observation (1). This gives the basis for 22 

the inductive proof. Now assume that the result is true for k-1 rows; Le., assume that 

with the values of all elements in k-1 rows computations of at most - + - elements can be completed. In addition to the values of all the elements in the k-1 rows, 

let the values of all the elements in the kth row be known. For the kth row there is a 

(k-1)2 (k-1) 

2 2 
10 

unique diagonal element, whose value can be computed using only the information 

from that row (Observation (1)). From Observation (2), using the information from 

each of the previous k-I rows and the kth row, computations at one unique new element can be completed. Thus, by knowing the appropriate values at all the elements 

in the kth row, computations of at most k extra elements can be completed; i.e., the 

values of all elements in the k rows can be used to complete computations of at most 

Rk 

2 2 22 (k-1)2 +- (k-l) + k = - + - elements. This completes the inductive proof of the first 

part of the lemma. 

To prove the second part, assume that there is a set of + elements such that 2 

all the computations at these elements depend on values from less that k rows. Since 

the computations at any element depends on its own initial value, it is obvious that for 

the above assertion to hold, these + k elements must lie on less than k rows. 22 

Assume that k - E rows hold these elements and that information from k - e' rows is 

sufficient to complete the computations at these elements, for E 2 E' 2 1. From the first 

part of the lemma, all the information from any k - E' rows can be used to complete 

computations of at most - (k-')2 + - elements. This implies that there are at least 2 2 

&' 1 i.(k - - + -) > 0 elements in k - E rows where the computations cannot be completed 22 

using the information in the k - E' rows - a contradiction. Therefore information from 

k2 k 

22 at least k rows is needed in order to complete the computations at - + - elements. 
11 

Next the above described scheme with a total data traffic of O(m2 dp) in factoring 

the m x m matrix using p proccssors is shown to be optimal in its communication 

requirements. 

THEOREM 2 : Assuming that the computational load is to be uniformly distributed 

among p processors, the data traflc involved in computing the Cholesky factor of an m 

x m dense matrix is R(m2 GI. 

PROOF: Consider the computations corresponding to the elements in the set 

S = (ai, I i,j > - ). In Figure 3 the gray area FGH denotes this set of elements. The m 

2 

m3 

6 total computational work in factoring the m x m matrix is - + 0(m3) and that in factorm3 ing the part corresponding to area FGH is - + o(m3. This implies that out of p pro- 12 

Gj H 

Figure 3. 
12 

cessors, cl.p processors must be assigned to complete the factorization of area FGH, for 

m3 m3 some constant 0 < c1 < 1, to achieve a uniform load distribution of c2.- + o(-) work 

P P 

per processor. 

Let n be some partitioning of elements, that partitions the elements in the set S in 

cl.p subsets such that each subset has the same amount of computational work. It can 

be shown that for any such partitioning n there are constants c', c" > 0, such that the 

- - 

number of elements in any subset are in the interval [;a$, c"-$]. Ignoring the lower 

m3 order terms each subset defined by the partitioning n has c2-- amount of work. Now 

P 

the element u,,,~ in the set S has the maximum amount of computational work associated with it, which is proportional to m since u,,,~ is modified m times. Therefore, the 

number of elements in any subset must be greater than or equal to c'.d, for some 

P 

m constant c'. Similarly, for any element a. ,, , i.e., an element in the column - + 1, Sy+l) 2 

the computational work is the minimum over the set S. This minimum computational 

m work is proportional to - Hence the number of elements in any subset defined by n 2' 

must be less than or equal to c".E for some constant c". 

2 

P 

m2 ' Let ci.-, c I ci I c", be the number of elements in the ith partition defined by n. 

P 

d From Lemma 2, computations at any ci.- elements depend on information from at 

P 

.m rows where C; is a constant. Now for any element aij E S, j 

m Hence the values of at least elements from each of the ci.2 m is greater than - 2' + 
13 

rows is required in completing the computation of any subset of elements defined by 

C: m2 

2 + n. Thus the data traffic associated with each partition is at least -.- . Since there 

are cl.p such partitions, the data traffic involved is at least cl*s-m24j. 2 

Therefore under the condition of uniform load distribution, the data traffic 

involved in computing the Cholesky factor of an m x m dense matrix using p processors is Q(m2 $). 

0 

It should be noted that, in an order of magnitude sense, the computational work 

associated with each block in the SPOCC factorization scheme presented earlier is the 

same as that of a partition under the uniform load distribution scheme. Thus the data 

traffic associated with that scheme is optimal in an order of magnitude sense. 

4. Communication Requirements of a SPOCC Factorization Scheme for Sparse 

Matrices 

In this section the total data traffic of a SFQCC factorization scheme in factoring 

a sparse, symmetric, positive definite matrix A corresponding to a regular n-vertex 2-D 

grid graph is analyzed. A 9-point stencil, unless otherwise stated, is used. In the following discussion L denotes the lower triangular Cholesky factor of the mamx A. 

First, the worst case communication requirement complexity result of O(nT), which is 

independent of the number of processors used, is obtained and then a scheme that 

3 

1 +a 

reduces the communication requirement complexity to O(n 2), when nu, as 1, 
14 

processors are used, is presented. The result is shown to be optimal under the condition of uniform load distribution. 

Clearly, the communication requirement is the worst when the use of local 

memory is not allowed. Now consider the computations involved in computing a 

i-1 

nonzero element lij E L. Recall that the expression aij - lik.ljh is computed first, folk=l 

lowed by a single division by ljj. Thus, for each multiplication there is one subtraction 

operation, at most one division and three memory references and a constant overhead 

such as index computation. Let the values of all the elements of the lower triangular 

part of matrix A and those of L be stored in the shared memory. Any number of processors are allowed to participate in computing a nonzero element of the factor provided that no single operation is performed by more than one processor. The following theorem gives the bound on the total data traffic under these assumptions. 

Although the result is obvious, it is useful because it is independent of the number of 

processors used and gives the worst case bound on the data traffic even for the models 

of computation that are more restrictive. 

TIIEOREM 3 : The worst case data traflc associated with factoring the matrix A is 

PROOF: For the assumed model of computation, each multiplication operation requires 

at most a constant number of memory references. Thus, the total data traffic is 

I cons1 . number of multiplication operations. 

From [l], the number of multiplication operations associated with factoring matrix A is 

3 3 

o(nZ). Hence, the total worst case data traffic is O(nZ). 
15 

Returning to the model of computation used to get bounds for the dense matrix 

case, the above worst case bound on the data traffic is improved when the number of 

processors is no where a < 1. This is possible because now the local memory is used 

to store the data that is needed in more than one computations. 

First, bounds are obtained on the data traffic associated with computing the elements of the Cholesky factor along the columns corresponding to a generic "+" shaped 

separator in the nested dissection ordering of a 2-D regular grid graph. Let the ordering scheme be such that the vertices on vertical part of a "+" separator are ordered 

after those on the horizontal part. Let each of the two horizontal parts of the "+" 

separator be referred to as the horizontal sub-separator and the single vertical part be 

referred to as the vertical sub-separator. 

Let 6 = ( k I k ij and liL f 0, li,k E L ); Le., ~ is the set of all columns of the factor 

L to the left of the column j+l such that the elements in row i of these columns are 

k 

nonzero. Let q$ = ~ ; i.e., is the set of all the columns to left of column j+l Fi 

such that on each of these columns there is a nonzero element in at least one of the i 

through k rows of the factor. Let r represent any m-vertex sub-separator. It is assumed 

that all the vertices in any sub-separator are ordered consecutively. Let low(r) and 

high(r) be the indices of the lowest and the highest ordered vertices, respectively, on 

the sub-separator r. Note that high(r) - low(r) + 1 = m. The following lemma establishes some basic sub-separator related properties that are useful in analyzing the communication requirement s. 
16 

LEMMA 3 : Let r be any m-vertex sub-separator. (i) Corresponding to the vertices of 

I' tlicrc is u dense m x m tricin~ukir diagonal block in the Cholesky frictor. (ii) In the 

,fuctor L, the columns low(T) through high(r) contain at most four of-diagonal rectangular blocks with nonzero elements. Each of these blocks is of size at most (cl.mtcz) x m 

where cI I 2 and c2 I 3 are integer constants. Any nonzero element in these columns is 

in one of these four blocks or in the diagonal triangular block and no where else. 

PROOF: The first part of the lemma is obvious. In Figure 4, the sub-separator r 

separates the vertices in regions R, and R2. Since the vertices in these two regions are 

ordered ahead of those of r, the fill-in due to the elimination of vertices in regions R1 

and R, ensures a dense m x m mangular diagonal block bounded by columns low and 

high(0 as shown in Figure 5. 

I IR1 I 

r2 

r4 

Figure 4: Sub-separator r with four surrounding sub-separators. 
To prove the second part of the lemma, consider Figure 4 again. In that figure 

the thickness of the lines qualitatively indicate the separator levels in the nested dissection ordering. Let rl, r,, r3, and r, be the four partial sub-separators that surround the 

sub-separator r. Because of the nature of the nested dissection ordering the vertices of 

r are "connected" to only those higher ordered vertices that lie on rl, r,, r,, and r, 

and to no other vertices.? Thus, all the nonzeros on columns low(T) through high(r) in 

rows below high(r) are confined to only the rows corresponding to the vertices on rl, 

r,, r3, and r,. Furthermore, each vertex in r is "connected" to every vertex on these 

four partial sub-separators and hence, the four rectangular blocks are dense. This is 

shown schematically in Figure 5. It can be verified that if r is a horizontal m-vertex 

n 

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

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

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

Figure 5: Off-diagonal blocks with nonzeros corresponding to sub-separator r. 

t Verrcx u is said to be "connected to vcrtcx v. if there exists a path [u. ulr 9. .... uk, VI of largth me or more in Ihc 

grid graph such that indcx(y) < min(indcx(u).index(v)). for 1 5 I 5 k; Le. if fij # 0. fij EL where i = indCx(u). j = indCx(v). 
18 

sub-separator, then the surrounding box of vertices is of dimension (2~3) x (mi-2). 

Therefore there are two rectangular off-diagonal dense blocks of dimension at most 

(h+3) x m and the other two of dimension at most (mi-2) x m. Similarly, if r is a 

vertical m-vertex sub-separator, there are four off-diagonal rectangular blocks of 

dimension at most (m+2) x m in the factor. If r is not surrounded on all four sides 

then some these blocks will be missing. 

0 

Using the above stated sub-separator properties, the following lemma gives the 

data dependencies of the nonzero elements in any of the five blocks specified in 

Lemma 3. The lemma shows that the number of nonzero elements in any row i of the 

factor L is less than c-m where c is an integer constant and m is the size of the subseparator to which the vertex corresponding to row i belongs. It is then shown that for 

any row i, the computations at all the elements lij E L, low(r) sj 5 higho, for some mvertex sub-separator r require a total of less than cam nonzero elements from that row. 

Note that this count is independent of the sub-separator to which the vertex 

corresponding to row i belongs. Thus, the computations at all the elements in a row of 

any of the five blocks specified in Lemma 3, require only cam elements from that row 

irrespective of the relative location of the off-diagonal blocks in the factor. 

LEMMA 4 Let r be any m-vertex sub-separator. The nonzero elements from row i, 

i 2 low(q, required in completing the computations of all elements lij E L such that 

low(n I j s high(r), are those elements in row i on the columns in the set given by, 

ip[low(f-)k.gh(r) qYghm . For all i 2 low(r), the number of such useful nonzero elements kink r) 
19 

in row i is I c.m for some constant c. 

PROOF: Any nonzero element lij E L, i 2 low(T) and low0 I j I higho, is in one of the 

five blocks specified in Lemma 3. Hence to prove the result of this lemma only the 

rows that intersect one of these blocks need to be considered. The result for 

low(T) I i 5 high(r) is first proved followed by that for i > highu. 

When bw(Q I i S high(T), n q?") = qyMI, since, 

qPgh(0 c qlow(&,high(r) high I7 . By definition, the set qrghm contains all the columns that have a 

nonzero element in row i. Clearly, the nonzero elements from the row i required in 

completing the computations at all the elements lii E L, low(T) I~S higho, are on 

columns in the set q,Aigh(n. 

To measure the size of the set qrgqn note that it is bounded by the number of 

vertices ordered ahead of the vertex i and which are "connected" to vertex i. Using the 

recursive nature of the nested dissection ordering it can be verified that in the case of 

constant degree grid graphs and when fow(T) I i s high(r), the size of the set qr#hm is 

bounded by cm, where c is a constant dependent on the degree of the graph and on 

whether r is a horizontal or vertical sub-separator. If r is a horizontal m-vertex subseparator, then for a 5-point stencil, c = 7 and for a 9-point stencil, c = 11. When r is 

a vertical m-vertex sub-separator, the values of c are 5 and 7, respectively. This completes the proof when low(Q I i I high(r). 

The case where i>high(r) is considered next. As shown above, IIqfghm 11 

depends on the size of the sub-separator to which the vertex i belongs and hence when 
20 

i > high(T), I1 q?xh(n 11 can be much greater than O(m) where m is size of r.? However, 

when the computation of only those elements in row i that lie on columns fow(T) 

through high(r) :ire of concern, each of these computations consists of ii product of ;i 

nonzero element in row i and a nonzero element in one of the rows fowo through 

high(r) in the column highQ or in some other column to the left of it. Thus, for these 

computations only the columns that have a nonzero element in row i and in row j 

where /ow(T) s j I high(r), are of interest. Now the set e$pfixMn consists of all the 

columns that have a nonzero element in at least one of rows fow(T) through high(r). 

Similarly, qpxh(n consists of all the columns that have a nonzero element in row i. 

Clearly, the set i$zRhgh(r) n qYxh0 consists of all the columns which contain all the 

pairs of nonzero elements that can be usefully employed in completing the computations at all the elements I,,, low(T) I j I higho. Thus, the nonzero elements from row i, 

i > high(r), required in completing computations at all the elements fi, e L such that 

low(T) I j I high(r), are those elements in the row i on the columns in the set given by, 

hixh r) WO 

qiow(h.lu8h(r) n qi - 

To get a bound on the size of qkmfighn n qygh(D, consider the m-vertex horizontal sub-separator r shown in Figure 4. It is surrounded by sub-separators rl, r,, 

r,, and r,. Suppose that low(rI) 5 i I high(rl). Now, the set qE$.&ith(T) n consists of columns corresponding to vertices on r or corresponding to those vertices 

ordered ahead of them which are "connected" to at least one vertex in r and to the 

vertex corresponding to row i. Using the recursive ordering of the nested dissection 

t 11 (I 11 dcnotcs the cardinality or he number of elancnts in sct U. 
21 

scheme it can be shown that the number of such vertices is less than 7m. Thus, 

The same bound is obtained 

when hw(r4) I i I high(r4). If h(r9 5 i I high(r2) or low(r3) I i I high(r3) then it can be 

verified that, II fit$RhgMn n qrgh(n 11 I 3m. If r is vertical sub-separator the two 

11 qhlgm iO4n,,,i8,,(r) n qrgMn II I 7m, for low(rl) 5 i I high(rl). 

5 

2 bounds are 5m and -m respectively. 

The scheme for completing the computations in all five dense blocks associated 

with an m-vertex sub-separator r using p processors is now described. If p = 1, that 

processor accesses the necessary values and completes all the computations. If p > 1, 

the factorization corresponding to the m x m triangular diagonal block is first completed using all the processors and then the factorization corresponding to the four 

off-diagonal blocks. For the first part, as stated in the previous section for the dense 

matrices, the m x m dense diagonal block is partitioned into - ‘- - square blocks and 22 

m m 31 r diagonal triangular blocks each of size - x - where p = - + - and each of these 

r r 2 2’ 

p partitions is assigned to a unique processor. Each processor completes the computations corresponding to its partition by accessing the required data from the shared 

memory. For the purpose of factoring, the off-diagonal blocks are treated as if they 

were adjacent, and the resultant rectangular block is partitioned into p sub-blocks each 

c.m m of size - x - where c I 6 for a horizontal sub-separator and c I 4 for a vertical 4p dp’ 

sub-separator. Again each partition is assigned to a unique processor. 
22 

Using the results from the previous lemma, a bound is obtained on the data traffic 

associated with the computations performed as outlined in the above scheme. 

LEMMA 5 : The total data traffic associated with completing the computations, using 

p processors, in all the dense blocks within the columns low(T) through higho, for 

some m-vertex horizontal sub-separator r, is bounded by (53 -+- 11fi)m2-J7;. It is 

bounded by (28 + 8Qrn2-$ for an m-vertex vertical sub-separator. 

PROOF: Let r be an m-vertex horizontal sub-separator. Assume that this sub-separator 

is surrounded on all four sides by vertices ordered higher than any vertex on the subseparator. Such a sub-separator has the worst case communication requirements among 

all the m-vertex sub-separators. 

First, compute the data traffic associated with factoring the triangular diagonal 

12'2 

22 block using p = - + - processors. Each of the sub-blocks requires nonzero elements 

from at most 3 rows out of the m rows in the range low(T) through high(r) of the factor. No other information is needed. From Lemma 4, each of these rows has at most 

11m nonzeros. Thus the communication requirement of each partition is at most 

m 11m.2- and the total communication requirement of the triangular block is, bounded 

r 

42P 

by I142m2.$. 

Now consider the data traffic associated with the off-diagonal blocks. Each parti6m m 6m tion is of size - x -. Thus, each partition requires nonzero elements from - 

rows which are below the row high(T) in the factor. From Lemma 4, each of these 

rows has at most 7m nonzeros that are useful in completing the computations in any of 

$6 J7; 
23 

thc partitions. Each partition also requires information from - m rows from the region 

.Ir, 

/OW(~ through high(r). Each of these rows has at most Ilm nonzeros. Thus, the communication requirement of each partition is at most 7m.% + 11m- = - s3m2 and the 6 6 6 

total communication requirement of completing the computations at the off-diagonal 

blocks using p processors is less than or equal to 53mv. 

Adding the communication costs corresponding to the diagonal and the offdiagonal blocks we get the total data traffic associated with r to be less than or equal 

to (53 + 1 Ifi)d.+i 

A similar analysis can be used to compute the data traffic when r is an nt-vertex 

vertical sub-separator and can be shown to be bounded by (28 + 8a)d.G. 

0 

Applying the results from the above lemma, a bound is obtained on the total data 

traffic in factoring the sparse matrix A associated with the 2-D grid graph. First a 

description of the overall load distribution scheme among the processors and some 

notation are given. For the sake of simplicity, assume that the n-vertex 2-D grid graph 

- 1 - I 

is a n2 x n2 square grid and that there are nu processors, as 1. The vertices of the 

grid are ordered using the nested dissection ordering scheme such that the vertices in 

each square sub-grid are consecutively ordered. In Figure 1, such an ordering for a 7 

x 7 grid is shown. This particular ordering is chosen only for the sake of describing 

the scheme. The results presented here are applicable to other nested dissection order- 

-- la 1 a -- 

ing schemes as well. This ordering results in no sub-grids each of size n2 x n2 2. 
24 

The computations at all the nonzero elements corresponding to the vertices of each 

such sub-grid is assigned to a unique processor. After completing the work, two processors from adjacent sub-grids combine to complete the work corresponding to the 

n2 2-vertex horizontal sub-separator that separates the two sub-grids. This is followed 

-- In 

-- la 

by the work on the 2n -vertex vertical sub-separator. Four processors work on each 

such vertical sub-separator. This process is continued until the factorization of the 

entire matrix is completed. Whenever more than one processor is working on any 

sub-separator, the scheme presented earlier for a generic m-vertex sub-separator is 

applied. 

Let Th (m. p, k) represent the data traffic using p processors in completing the computations at all the nonzero elements lij E L in the columns corresponding to an mvertex horizontal sub-separator that is surrounded by higher ordered vertices on k sides. 

Let T, (m. p. k) represent the same for an m-vertex vertical sub-separator. From Lemma 

5, Th (m. p, 4) S (53 + lld2)m2.$ and 't, (m. p. 4) 5 (28 + 8fi)d-$. kt 'tr (m. p, k) 

represent the total data traffic, using p processors, in completing the computations 

corresponding to all the sub-separators within an m-vertex sub-grid that is surrounded 

by higher ordered vertices on k sides. 

The following theorem gives an upper bound on the total data traffic in factoring 

the matrix A associated with an n-vertex 2-D regular grid graph using no processors 

with the scheduling scheme as described above. 

THEOREM 4 : The total data traffic in factoring the n x n sparse matrix A using n" 

* 

a 1 +- 2 

a I +- 

processors is O(n z); i.e., T~ (n, n", 0) = O(n ). 
25 

- 1 - 1 1 

PROOF: On a n2 x n2 regular grid there is an nT-vertex vertical sub-separator and 

I 

1- 

2 two -n -vertex horizontal sub-separators (ignoring the additive constant -1). The 

vertical sub-separator is not surrounded by any vertices that are ordered after the vertices on the vertical sub-separator. Each of the two horizontal sub-separators are surrounded on one side. These three sub-separators subdivide the n-vertex grid graph into 

1 1 1- four sub-grids of size ln2 x -n2, each surrounded on two sides by higher ordered 

vertices. Thus, the total data traffic in factoring the corresponding matrix A is given 

2 2 

by * 

1-1 11 

22 4 '4 

- 1 1 

rR (n. na, 0) = T~ (n2, na, 0) + 2.5h ( -n2, -na, 1) + 4.5, ( --n -na, 2). 

A recursive expansion of the above expression contains data traffic terms for vertical sub-separators of different sizes that are surrounded on zero sides, two sides, three 

sides (in two different ways), and on all four sides by higher ordered vertices. It also 

contains data traffic expressions for horizontal sub-separators of different sizes surrounded in five different ways. To keep the analysis simple, it is assumed that all the 

I 1 

1- 

2 2 four sub-grids of size -n2 x 'nT are surrounded on all four sides. This 

simplification results in a conservative expression for the data traffic, but affects only 

the constant terms in the bound. Thus, 

Now, 
26 

Froin Lemma 5, it follows that, 

(1 1 +- 11 1 

4 '4 16 T# ( -n -nu, 4) I - (134 + 85&).n '. 

An analysis similar to that given in Lemma 5 yields 

I+% '. - 1 

T, (n2, na. 0) < 8% 

and 

a 1 I+- 

2 8 

1 1-1 

Th ( Tn 2, -na, 1) < -(22 + 25@n 2. 

Thus, we get, 

a 

2. 

1 +- 1 

2 7, (n. na, 0) 5 -(78 + 71@n 

In the following theorem the communication bound of the above presented 

scheme is shown to be optimal in an order of magnitude sense, by giving a lower 

bound proof. 

THEOREM 5 Under the condition of uniform load distribution, the data trafic in 

a I+- 

factoring the n x n sparse matrix A, using na processors, a < 1, is Q(n 2). 

PROOF: For a regular 2-D grid graph with n vertices, the separator size for nested 

- 1 

dissection ordering is n2 171. From Lemma 3, it follows that the factor L has an 

- 1 - 1 

n2 x n2 dense triangular diagonal block incorporated in it. From Theorem 2, the data 

traffic involved in completing the computations associated with the elements of this 

dense triangular block, under the condition of uniform load distribution using na pro- 
27 

I+" 

cessors, is Q(n 

ing the factorization of this dense block, the result follows. 

0 

2). Since the factorization of A cannot be completed without completFrom Theorems 4 and 5, it is clear that the load assignment scheme described 

here for factoring the n x n sparse matrix using na processors is optimal in an order of 

magnitude sense. Note that when nu, a > 1, processors are used, the data traffic bound 

given in Theorem 3 holds. 

5. Conclusions 

In this paper load distribution schemes for multiprocessor systems with local and 

shared memory are presented for factoring dense and sparse symmetric, positive 

definite matrices. The communication requirements of these schemes are analyzed and 

are shown to be asymptotically optimal for systems corresponding to regular 2-D problems where the vertices of the grid graph are ordered using the nested dissection ordering methods. For these schemes it is shown that the total data traffic in factoring an 

2+Q 

n x n dense, symmetric, positive definite matrix is O(n 2, when no processors are 

used. The total data traffic in factoring the n x n sparse, symmetric, positive definite 

I+Q 

matrix corresponding to a 2-D grid problem, is shown to be O(n 2, when no, a s I, 

3 

processors are used. When nu, a > 1, processors are used the total data traffic is O(nT). 

Under the condition of uniform load distribution, these results are shown to be optimal 

in itn order of magnitude sense. 
28 

In $1, George et al. have given a load assignment scheme Adr factoring the 

matrix A that has a total data traffic of O(nl+") when O(na) processors are used. In this 

assignment scheme, the computational work corresponding to an entire colunin of the 

matrix A is performed by no more than one processor. In other words the work associated with any column is treated as indivisible. Under the condition of uniform load 

distribution and a column-level indivisibility, in [SI and [8] it is shown that the total 

data traffic of that scheme is asymptotically optimal. The scheme presented in this 

I +Q 

paper reduces the communication requirement to O(n 2, by removing the constraint 

of column-level indivisibility. Here the indivisible work unit is the computation 

corresponding to a nonzero element in the factor. Note that this does not necessarily 

enforce the use of a finer grain computational model. However, the SPOCC factorization scheme allows an efficient use of up to O(n) processors before the total data traffic 

3 1 

reaches the maximum value of O(nT), whereas in the former scheme only up to O(nT) 

processors may be efficiently used before the total data traffic reaches the maximum. 

Another notable aspect of the scheme presented in this paper is that the reduction 

in the communication requirements is brought about by improving the utilization of the 

data accessed from shared memory by each processor. Consider the factorization of an 

m x m dense matrix. Let the data utilization of a data element accessed by a processor be defined as the number of computations in which that element is used by that 

processor divided by m. Since an element in the factor is needed in at most m computations, the maximum utilization of any data accessed is one. Let the aggregate data 

utilization for a processor be defined as the average utilization of the individual data 
29 

IIelcmcnts accessed by that processor. In the SPOCC factorization of an m x m dense 

2m2 matrix, each processor accesses at most - elements from the shared memory and G 

m each element is used in at least - computations. Thus, the utilization of each data G 

-1 accessed is - and so is the aggregate utilization of all the data accesses. On the 

other hand, with the column-level work assignment scheme, each processor accesses 

O(m2) elements from the shared memory. Of these, only O(2) elements have an utili- 

.iT 

P 

1 

P 

zation of one and the data utilization for the remaining elements is - which gives an 

1 

P aggregate data utilization of approximately -. Similar improvements in data utilizations are obtained in factoring a sparse matrix. 

It should be noted that the square shape of the submatrix partitions produce the 

best possible aggregate utilizations. For the algorithm considered here, the data dependencies are such that rectangular and square partitions give rise to high data utilizations. Since the square partitions have the minimum perimeter for a given area, the 

number of data elements accessed (which is proportional to the perimeter of the partition) for a given work load (which is proportional to the area enclosed), is also a 

minimum for the square partitions. 

An effect of the improvement in the aggregate utilization of data and the resulting 

reduction in the communication requirements is the segregation of the accesses to the 

shared data. Since the total data traffic in factoring an m xm dense matrix using p 

m2 

G processors is O(d. G), on an average each processor accesses only O(-) data. Note 
30 

V. K. Naik and M. L. Patrick. Communication requirements ofparallel sparse Cholesky faclorizatim. in prepration. 

that the total shared data is O(m2). Thus, on an average each element in the shared 

memory is accessed by O(+) processors. The column-level assignment scheme, however, has a total data traffic of O(d. p) and thus, on an average each processor 

accesses o(&) data or on an average each element in the shared memory is accessed 

by ~(p) processors. An obvious implication of this observation is that for the scheme 

presented here, not only is the total data traffic reduced but also the requests at individual shared addresses. This can have considerable impact on the performance of the 

systems with a large number of processors. 

In the factorization scheme presented here, the work assigned to each processor is 

balanced in an order of magnitude sense which is sufficient for carrying out the 

analysis. For practical purposes it is necessary to consider the effect of the constants 

involved on the load imbalance. Schemes for reducing the load imbalance caused by 

the constant terms in the column-level assignment scheme and in the SPOCC scheme 

are presented elsewhere.2 Some implementation aspects are also considered there. 

Finally, the scheme developed here for factoring a sparse matrix is applicable to a 

wider class of problems. The the communication requirements for these problems are 

found to be optimal as 

In summary, we have presented schemes that reduce the data traffic involved in 

factoring both dense and sparse, symmetric, positive definite matrices on multiprocessor systems with a two level memory hierarchy. These reductions are made possible 

by improving the utilization of each data element accessed. 
31 

ACKNOWLEDGMENTS 

The authors would like to thank Bob Voigt for his support and encouragement. 

REFERENCES 

(1) J. A. George, Nested dissection of a regular finite element mesh, SIAM J. Numer. 

Anal., vol. 10, pp. 345-363, 1973. 

(2) J. A. George, M. T. Heath, J. W. H. Liu, and E. Ng, Sparse Cholesky factorization on a local-memory multiprocessor, Technical report ORNL/TM-9962, Oak 

Ridge National Laboratory, Oak Ridge, Tenn, 1986. 

(3) J. A. George and J. W. H. Liu, Computer Solution of Large Sparse Positive 

Definite System, Prentice-Hall, Inc., Englewood Cliff, NJ, 1981. 

(4) J. A. George, J. W. H. Liu, and E. Ng, Communication reduction in parallel 

sparse Cholesb factorization on a hypercube, in Hypercube Multiprocessors 

1987, ed. M. T. Heath, SIAM Publication, 1987. 

(5) J. A. George, J. W. H. Liu, and E. Ng, Communication results for parallel sparse 

Cholesky factorization on u hypercube, Submitted to Parallel Computing, 1988. 

(6) G. H. Golub and C. F. Van Loan, Manix Computations, The Johns Hopkins 

University Press, Baltimore, MD, 1983. 

. 

(7) R. J. Lipton, D. J. Rose, and R. E. Tarjan, Generalized nested dissection, SIAM J. 

Numer. Anal., vol. 16, pp. 346-358, 1979. 

(8) V. K. Naik and M. L. Patrick, Analysis of communication requirements of sparse 

Cholesky factorization with nested dissection ordering, Roc. of the third SIAM 
32 

Conference on Parallel Processing for Scientific Computing, Los Angles, Dec. 1- 

4, 1987. 

(9) E. M. Reingold, J. Nievergeld, and N. Deo, Combinatorial Algorithms: Theory 

and Practice, Pren tice-Hall, Inc., Englewood Cliff, NJ, 1977. 

, 
1. Report No. 

NASA CK-181626 

ICASE Report No. 88-14 

2. Government Accession No. 3. Recipient's Catalog No. 

10. Work Unit No. 

505-90-2 1-0 1 

4 Title and Subtitle 

DATA TRAFFIC REDUCTION SCHEMES FOR SPARSE 

CHOLESKY FACTORIZATIONS 

5. Report Date 

February 1988 

6. Performing Organization Code 

Langley Research Center 

Hampton, Vi4 23665-5225 

- 

7. Author(s) 

Vijay K. Naik and Merrell L. Patrick 

4 14. Sponsoring Egency Code 

8. Performing Organization Report No. 

88-1 4 

I 

15. Supplementary Notes 

Langley Technical Monitor: Submitted to Int. Supercomputing 

Richard W. Rarnwell Conference 

9 Performing Organization Name and Address 

Institute for Computer Applications in Science 

Mail Stop 132C, NASA Langley Research Center 

Hampton, VA 23665-5225 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 

and Engineering 

Final Report 

11. Contract or Grant No. 

NAS 1-1 8107 

13. Type of Report and Period Covered 

Contractor Report 

%. Abstract 

17. Key Words (Suggested by AuthorM) 

communication requirements, shared 

and local memory, multiprocessor 

systems, sparse Cholesky factorizations 

Load distribution schemes are presented which minimize the total data traffic in 

the Cholesky factorization of dense and sp.me, symmetric, positive definite matrices 

on multiprocessor systems with local and shared memory. The total data traffic in factoring an n x n sparse, symmetric, positive definite mamx representing an n-vertex 

regular 2-D grid graph using na, a 5 1, processors is shown to be O(n *). It is O(nT), 

when no, a 2 1, processors are used. Under the conditions of uniform load distribution 

these results are shown to be asymptotically optimal. The schemes allow efficient use 

of up to O(n) processors before the total data traffic reaches the maximum value of 

O(nz). The partitioning employed within the scheme, allows a better utilization of the 

data accessed from shared memory than those of previously published methods. 

3 1 +o 

3 

18. Distribution Statement 

59 - Mathematical and Computer 

62 - Computer Systems 

Unclassified - unlimited 

Sciences (General) 

19. Security Classif. lof this report) 20. Security Classif. (of this page) 21. No. of pages 

Uncl ass i f i ed Unclassified 34 

22. Price 

A0 3 

NASA FORM 162S OCT 86