NASA Contractor Report 181863 

ICASE Report No. 89-40 

ICASE 

DATA TRAFFIC REDUCTION SCHEMES FOR 

CHOLESKY FACTORIZATION ON ASYNCHRONOUS 

MULTIPROCESSOR SYSTEMS 

Vijay K. Naik 

273 49 1109- 

June 1989 

Institute for Computer Applications in Science and Engineering 

NASA Langley Research Center 

Hampton, Virginia 23665-5225 

Operated by the Universities Space Research Association 

National Aeronautics and 

Space Administration 

Langley Research Center 

Hampton, Virginia 23665-5225 
Recently, ICASE has begun differentiating between reports with a mathematical or applied science theme and reports whose main emphasis is some aspect of 

computer science by producing the computer science reports with a yellow cover. 

The blue cover reports will now emphasize mathematical research. In all other 

aspects the reports will remain the same; in particular, they will continue to be 

submitted to the appropriate journals or conferences for formal publication. 

i 
I 

Data Traffic Reduct ion Schemes for C holesky Factorization 

on Asynchronous Multiprocessor Systems* 

Vijay K. Naik 

IBM Rescarcli, T.J. Watson Research Center 

P.O. Box 704 

Yorkton-n Hcights, N.Y. 12598, U.S.A. 

Merrell L. Patrick 

ICASE, NASA Langley Research Center 

Ha.mpton, \rA 23665, U.S.A. 

and 

Computer Science Department, Duke University 

Durham, N.C. 27706, U.S.A. 

Abstract: Communication requirements of Cholesky factorization of dense and sparse 

symmetric, positive definite matrices are anaIyzed. The communication requirement is 

characterized by the data traffic generated on multiprocessor systems with local and shared 

memory. Lower bound proofs are given to shorn that when the load is uniformly distributed 

the data traffic associated with factoring an n x n dense matrix using nu, a 5 2, processors 

is fl(n2+a/z). For an n x n sparse matrices representing a fi x fi regular grid graph 

the data traffic is shown to be n(n1+a/2), a 5 1. 

Partitioning schemes that are variations of block assignment scheme are described and 

it is shown that the data traffic generated by these schemes are asymptotically optimal. 

The schemes allow efficient use of up to O(n2) processors in the dense case and up to O(n) 

processors in the sparse case before the total data traffic reaches the maximum value of 

0(n3) and O(n3l2), respectively. It is shown that the block based partitioning schemes 

allow a better utilization of the data accessed fiom shared memory and thus reduce the 

data traffic than those based on column-wise wrap around assignment schemes. 

'This research was supported by the National Aeronautics and Space Administration under NASA Contract Nos. NAS1-18107, NAS1-18605, and by the Air Force Office of Scientific Research under Contract 

No, AFOSR 88-0117 while the second author was in residence at the Institute for Computer Applications 

in Science and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23665. 

iii 

PRECEDING PAGE BLANK NOT FILMED 
1. Introduction 

Consider the problem of solving a system of linear equa.tions 

Ax=b 

where A is an n x n symmetric, positive definite coefficient matrix, x is an n x 1 vector of 

variables, and b is an n x 1 vector of constants. Applying the Cholesky decomposition to 

A yields 

A = LL~ 

where L is lower triangiilar with positive diagonal elements [lo]. From this factorization, 

the solution to the system of equations is obtained by soIving the triangular systems 

Ly = b 

and 

1' L x=y 

Recently, efforts have been reported for efficiently parallelizing the various steps in computing the solution of the dense and sparse systems. Most of this work has concentrated on 

developing algorithms that extract as much parallelism as possible on specific architectures 

[14,18,19,1,5,3,9]. The main emphasis there is on distributing the computational load as 

evenly among the processors as possible and little attention is paid towards the data traffic 

complexity. 

In this paper we are interested in the parallel Cholesky decomposition schemes with minimum data traffic for factoring dense and sparse symmetric, positive definite matrices. The 

model of computation assumed for this purpose is that of a multiprocessor system with two 

level memory hierarchy such that each processor has local memory and all processors have 

access to a common shared memory. Accessing any nonzero element in the shared memory 

is assumed to generate a unit data traffic. No data traffic is generated in accessing the 

local memory. The total number of shared memory accesses from the beginning to the end 

of the algorithm is defined as the communication requirement or total data trafic of that 

algorithm implemented on the multiprocessor system. 

In [17] the communic'ation requirement of the Gaussian elimination algorithm implemented 

on threr different architectures is analyzed. For a bus architecture where a data element 

may be broadcast to all the processors in one step and counts as one unit data traffic 

independent of the number of processors receiving the data, the data traffic complexity is 

shown to be n(n,'). For a nearest neighbor ring network, where each transmission of a data 

element across a link of the ring coants as one unit data traffic, the data traffic complexity 

is shown to be n(n2 sp), where p is the number of processors on the ring. The data traffic 

complexity for a nearest neighbor grid network is shown to be R(n2&). In all the cases 

it is assumed that no element is computed in more than one processor; Le., recomputation 

1 
is not permitted. By using a different proof technique than that given in [17], it is shown 

here for the assumed model of compiitation that the data traffic complexity of the Cholesky 

factorization scheme is n(n2fi). The proof for the lower bound holds even if recomputation 

is allowed provided each processor is assigned at lcast n3/6p amount of work. Although we 

do not prove it, our resiilt holds for other variations of Gaussian elimination as well. 

In [4], a parallel sparse factorization scheme is given for local memory multiprocessor systems. This scheme has a total data traffic of O(nl+" log, n,) using na processors. This resiilt 

is improved to O(n'+") in [7]. In this paper we prescnt a factorization scheme that has a 

total data traffic of O(n'+"/2). Our main results and the organization of the paper is as 

follows. 

In the following section, the data dependencies involved in the Cholesky factorization of a 

dense matrix are discussed and a parallel assignment sclicme is presented. It is shown that 

the data traffic associated with that scheme is O(n,' - Ji;) when an n x n dense matrix is 

factored using p processors. By giving a proof on the lower bound for the data traffic, in 

Section 2.4 it is shown that under the condition of uniform load distribution the computation 

time and data traffic complexities of the assignment scheme are asymptotically optimal. In 

Section 3., the case of factoring sparse, symmetric, positive definite matrices is considered. 

The sparse matrices considered here are restricted to only those matrices that represent the 

graphs arising in finite difference and finite element applications. In Section 3.5, a block 

based parallel factoring scheme for sparse matrices is presented. The data traffic in factoring 

an n x n sparse matrix corresponding to a 2-dimensional regular grid graph is shown to be 

O(n Jlr). In Section 3.6 a lower bound on the data traffic in factoring the sparse matrix 

is shown to be n(n Ji;). These results can be extended to other graphs that satisfy an 

f(n)-separator theorem [13]. Preliminary versions of the results given here appear in [15] 

and [16]. 

For the sake of clarity, in the following discussion the dense ma.trix is assumed to be of size 

m x m and the spa.rse ma.trix of size n x n. 

2. Parallel factoring of dense symmetric, positive definite 

matrices 

The basic algebraic scheme considered here for factoring an m x rn symmetric, positive 

definite matrix A is the column version of the Cholesky decomposition method [lo]. An 

outline of this algorithm is given next and the data dependencies are discussed. Following 

that a partitioning scheme with optimal data traffic is presented. 

In the following discussion, values in row i refer to the values of the elements on and to the 

left of the diagonal. Simila.rly, a;,, (u*,j) represents all the elements in IOU' i (in column j) 

I 

2 
of t,he lower tria.ngular part of the ma.trix under considera.tion. 

2.1 The Cholesky factorization 

for j = 1 until m, do 

begin 

.. Initialize lij = ai$, E = 1, - - * rm 

for k = 1 until j - 1 do 

for i = j until m do 

1. * 1. . - 1. . - 1 i,3 - r,j i,k 3,k ; 

for k = j + 1 until m do 

Ikj = ik,jflj,j 

end 

In the above algorithm for cla.rity, the values of lij are shown separately from those of ui,j. 

In practice lij ma.y overwrite aij. 

Clearly, in the Cholesky factorization scheme outlined above, computing the elements in 

a column j of L requires values of the elements in columns 1 through j - 1 of L and the 

values of the off-diagonal elements in j are used for computations of elements in columns 

j + 1 throiigh m of L. Specifically, computing an element lij in L requires all the values 

from the set, 

Aij = {lj,jl I 1 5 j' 5 j} U {li,jl I 1 5 j' < j} U {di,j}. 

Moreover, the steps of the innermost loop, where a product of two elements of L is subtracted from ui,j, may be performed in any order. Once I;,; is computed, it is used in the 

compiitation of every element in the set, 

Ai,j = {~i,jl I j < j' 5 i} u {~jt,i I i 5 j' 5 m}. 

Again the element lij may be used in any order in the computations of the elements in the 

set Ai,j. 

3 

I 

I 
2.2 A partitioning scheme for Cholesky factorization 

Without loss of gcnerality, assume that p = (r2 + r)/2 where t is an integer. The lower 

triangular part of the matrix A is divided into p partitions by taking r vertical and r 

horizontal sections each of size a, where a = m/t. All except r of the resulting p partitions 

are square blocks of size a x a. The remaining t partitions which Lie on the diagonal of the 

matrix are a x 8 triangular blocks. Each of the partitions is assigned to a single processor. 

Initially, each processor reads the data for its partition fiom shared memory into its local 

memory. The computations proceed in parallel according to the column version Cholesky 

algorithm as fullow. The t processors in charge of the partitions containing the left most 

B x B 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. This is continued until the entire 

factor is computed and written into the shared memory. This partitioning and factorization 

scheme for dense symmetric, positive definite matrices is referred to as the block orienled 

column Cholesky-factorization scheme or simply as the BLOCC scheme. 

(i-l)*mtl i.. 

Figure 1: The data tra.ffic associated with block I 

2.3 Data traffic complexity of the BLOCC scheme 

First consider the data traffic associated with computing the elements in a generic square 

block I in the factor L shown in Figure 1. In that figure, the darkened area represents the 

4 
data elements that me reqiiired for the computations in block I. Let block I be bounded by 

the columns (i - 1)s + 1 and i. s, and by the rows (j - 1).9 + 1 and j s where 1 5 i 5 t and 

1 5 j 5 T. The following lemma provides bounds for the communica.tion cost of block I. 

Lemma 1 A total data trafic 0f(2i-1/2).s2+s/2 is nece.9sary and suficientfor computing 

the elements in the aquare block I; it is the same for all square blocks bounded by the columns 

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

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

PTOO~: See [15]. I 

Using these results, a. bound on the total data tra.ffic is obtained, as shown next. 

Theorem 1 The total data traffic associated with the BLOCC scheme for factoring an 

m x m dense symmetric, positive definite matrix using p processors is O(m2fi). 

Proof: The total data tra.ffic associated with all the blocks bounded by columns (i - 1)s + 1 

and i 8, 1 5 i 5 T, is given by, 

(T - i) x (data traffic associated with a square block) 

+ (data traffic associated with a triangular block bounded by the given columns). 

From Lemma. 1 and the fact that there a.re t such column partitions, we get the total data 

tra.ffic involved in factoring the 

r 

i=i 

Ignoring the lower order terms, 

m x m dense matrix using the BLOCC scheme as: 

the total data traffic is 

r 

= (2t * i - 2i2) s2 

i=l 

Since p = (t2 + ~)/2 and s = m/r, the total data traffic is O(m2&). 

2.4 A lower bound on the data traffic complexity 

1 

Theorem 1 gives an upper bound on the data._traffic associated with factoring the m x m 

matrix using the BLOCC assignment scheme. In the next theorem a lower bound on the 

5 
data traffic in factoring the dense symmetric, positive definite matrix is established. Before 

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

of the factor. Consider the computations at any element a,j as shown in Figure 2. To 

compute the corresponding element I,j, values at all the elements in row' j and the values 

of elements in column 1 through j of row i are needed. Thus, if aij is an off diagonal 

element then 2j valiics are needed for the computations and if it is a diagonal element (i.e., 

i = j) then j values are required. There are three observations to make regarding these 

computations as follo\r-s: 

i. The values at all the elements in any row i are needed to complete the computations 

corresponding to the diagonal element a;,;; no other values are needed for computing l,,i. 

Moreover, no other element in the factor can be computed by knowing only the values in 

row i. 

ii. If i and j are any two rows such that i > j, then the values of the elements in these two 

rows are used to complete computations at exactly one off-diagonal element aij. Values 

from no other row are needed to complete the computations at that element. 

iii. For the computations at a subset of elements spread over kr rows and A, columns, values 

from at least max(k,, k,) rows are needed. 

Figur 2: Data dependencies for the computations at ai,j and 

ments in subset S 

t eleThe last observation follows from the fact that the computations require values fiom h, 

rows as well as fiom the k, rows that correspond to the k, columns. However, some of the 

6 
I 

t 

k, rows and the k, row mag overlap. In addition, note that if j' is the leftmost of the k, 

columns then at least all the values in columns 1 through j' on max(k,, k,) rows are required 

in the computation of the elements in the subset under consideration (sce Figure 2). 

The above observations and thc result established in the following lemma are used to get a 

lower bound on the data traffic in computing the Cholesky factor. 

Lemma 2 Let 1.T' be the amount of computational work which is to be distributed uniformly 

among p processors and let a be any constant less than one. For any subset of this computation consisting of Ti'/2 amount of work, there are at least (1 - a) .p/(2 - a) processors 

each assigned a W/2p or more work from that subset. 

Proof: The work is uniformly distributed among p processors and hence each processor 

is assigned lV/p amount of work. Now let S be a subset consisting of W/Z amount of 

computational work. All p processors may be assigned some portion of work from S. Let 

w; be the computational work from S assigned to processor pi, where 0 <_ wi 5 W/p. 

Therefore. 

M' 

i=l 

and W/2p is the average work from S performed by each processor. Thus, there is at least 

one processor that is assigned W/2p or more work fiom S. Let a be a constant less than 

one and suppose that x processors are assigned at least a - lY/2p amount of work from S. 

Each of the x processors may be assigned at most lV/p work from S. Now there are p - x 

processors that compute less than a . W/2p amount of work from S. Therefore, 

Solving the inequality for x we get, 

1-a 

2-a a!?--- - P. 

Thus, there are at least (1 - a) p/(2 - ct) processors ea.ch computing a - W/2p or more 

a.mount of work from S. I 

In the following theorem a bound on the da.ta traffic associated with computing the Cholesky 

factor of an mxm, matrix is esta.blished. Note tha.t the result holds under stronger conditions 

than required by the model of computa.tion assumed here. 

Theorem 2 Let A be a dense, m x m symmetric positive definite matrix that resides in the 

common memory. If the computational work is uniformly distributed among p processor.?, 

then the data traffic involved in computing the Cholesky factor ofA is n(m2.Jp3. Forp 2 6, 

the data traffic is R(m2 &) even if the initial values of matrix A are in the processor local 

memory before the computation begins. 

7 
&. m 

Figure 3: Data. tra.ffic associated with factoring of elements in region 

FGH 

Proof: Suppose that matrix A is initially stored in the common memory. The initial value 

of each element is fetched at least once by the processors for computing the factor. Thus, 

at least m2/2 a.mount of data tra.fFic is associated with the Cholesky factorization. Hence 

when the number of processors is a constant, there is nothing to prove. 

In the following it is proved that even if the matrix A is distributed initially among the 

processors according to their work assignment, the data traffic is n(m2 - Jp3 when the 

number of processors is greater than 16(a + 4/a - 4)/3 for any constant a less than one. 

Since there exists an a less than one such that 16(a+4/a -4)/3 is less than six, the proven 

result holds for all p 2 6. Consider the computations corresponding to the elements in the 

set S = {aij I i, j > rn/2}. In Figure 3 the region FGH denotes this set of elements. The 

total computational work in factoring the m x m matrix is mS/6 + m2/2 + m/3 and that 

corresponding to the elements in set S is m3/12 + m2/4 + m/6. Thus, the amount of work 

associated with S is exactly half of the total work. If a is a constant less than one, then 

from Lemma 2, there are at least (1 - a) -p/(2 - a) processors each computing at least 

a m3/12p amount of work from the region FGH. 

Let 11 be the set of processors each with at least a.m3/12p amount of work from the 

region FGH and let pi E II. Now the computational work associa.ted with any element in 

FGH is at most m (work corresponding to element urn,,,, ). Therefore at least a - m2/12p 

elements in the region FGH are assigned to processor pi. Let z be the number ol rows 

on which the elements assigned to processor pi lie. This implies that there are at least 

[(a - m2)/(12 -p - z)1 columns on which the elements assigned to pi lie. Therefore, from the 

8 
i 

I 

Observation (iii) above, data from at least max(z, [(a - m2)/(12 -p. z)]) number of rows in 

the region DFGE arc required for completing the computations from region FGH assigned 

to processor p,. Without loss of generality, assume that the quantity 12 .p. z divides a - m2 

evenly. Now, the quantity max(z, (a. m2)/(12 -p - z)) is minimum when z = 6. m/m. 

Thus, the computations in processor pi require at least all the values of the elements on 

the fie m/a rows in region DFGE. In region DFGE each row has m/2 elements and 

thus the values of at least & - m /( 4 fi) elements from the region DFGE are needed 

in processor p, for performing computations in the region FGH. 

Now processor p, mag also be assigned some work from the region DFGE in addition to that 

in FGH. Hence to complete the proof, it is necessary to show that of the fi. m2/(4 - &) 

elements needed by processor pi at least c m2/& elements are not available locally, where 

c is a constant less than one. In that case the data traffic associated with processor pi is 

at least c - m2/& and since there are at least (1 - a) .p/(2 - a) such processors, the total 

data traffic in computing the Choleskp factor of an m x m dense matrix is n(m2 - fi). We 

complete the proof by showing in the following that p, accesses at least c - m2/fi non-local 

elements from region DFGE for completing the computations in the region FGH. 

Processor p, is assigned at least a. m3/12p amount of work from the region FGH. Since 

each processor is assigned m3/6p amount of work (the uniform load distribution condition), 

p, performs at most (2 - a). m3/12p amount of work in the region DFGE. The data traffic 

associated with processor p; in completing the work in the region FGH is a minimum 

when all the elements from region DFGE assigned to pi lie on the &- m/a rows. 

Furthermore, to reduce the data traffic, as many elements on these rows as possible should 

be assigned to processor pi. Now the computational work corresponding to any element 

a;j is j; that is the work associated with an element on the leftmost column of the matrix 

is the smallest and it increases for elements on any row from left to right. Therefore the 

data traffic associated with p; is a minimum when it is also assigned the computational 

work corresponding to the elements in the leftmost columns on the chosen IORS of region 

DFGE. Let E be the number of the leftmost columns on which the the elements from region 

DFGE that are assigned to p, lie. The shaded region shown in Figure 3 corresponds to the 

elements which minimize the data traffic for processor pi. Since processor p, performs at 

most (2 - a) - m3/12p amount of work in DFGE, the condition on k is given by, 

fi.m k (2 - a) e m3 -. xi 12P mF I 

;=1 

i.e., 

4(2-a) m2 1 

It can be verified that there is a constant greater than one, such that if p is greater 

than 1G(a - 4 + 4/a)/3, then the right hand side of the above inequality is at most 

m/2P for all values of m. This gives a bound on k. Therefore work corresponding to 

at most (fi. m2)/(4,0 - &) elements in the region DFGE may be assigned to processor 

pi which will minimize its data traffic for the computation in the region FGH. Hence of 

9 
the (6. m2)/(2 - a) elements needed by processor pi for completing the computation 

in the rcgion FGH, at lcast (1 - I//?). fi. m2/(4,&) elements are not available locally. 

Thus, if thc number of processors, p, is greatcr than 1G(a - 4 + 4/(r)/3 for any a less than 

one, then the data traffic associated with proccssor p, is at least c m2/& for some constant c less than one. Since there are at least (1 - a).p/(2 - a) such processors, the result 

follows. 1 

2.5 Remarks on the BLOCC Scheme 

Assuming that each step of the innermost loop in the Cholesky decomposition costs one 

computa tional time unit and ignoring the costs associated with other steps, the scquential 

computation time for factoring the rn x m matrix A is m3/6 + O(m2). The BLOCC scheme 

described above has a computation time of m3/2p + O(m2/p), where p is the number of 

processors uscd. As shown in Theorem 1 the associated data traffic is less than fi-m2.&/3. 

Thus, the time and the data traffic complexities of the BLOCC scheme are optimum in an 

order of magnitude sense. However, the computational load in the BLOCC assignment 

scheme is not perfectly balanced. The processors that compute elements in the partitions 

that are towards the left side of the matrix L finish computation earlier than those that are 

on the right. This balance mag be improved in several different ways, but at the cost of 

increasing the data traffic. In one such scheme the columns of the matrix are assigned to 

each processor in a wrap around fashion; that is, columns i,p + i, . . . , m - p + i are assigned 

to processor i. All the elements on any column of L are computed by a single processor. 

Let this assignment scheme be referred to as the wrap around assignment scheme. In this 

scheme the computation is distributed more evenly among the processors than that in the 

BLOCC scheme. The computation time is reduced to m3/6p + O(m2), provided m is at 

least p(p + 3)/2. However the data traffic associated with this scheme is m2 - p/2, which 

is suboptimal. In [ll] and [2] the wrap around assignment scheme is recommended as 

a preferred method for computing the factor on a multiprocessor system because of its 

good load balancing properties. Their analysis does not take into account the cost of the 

associated data traffic, which must be taken into accvunt for reducing the overall execution 

time. 

3. Parallel factoring of sparse, symmetric positive definite 

matrices 

In this section a partitioning and assignment scheme is presented that computes the factors 

of an n x n matrix, associated with a 2-d regular grid graph, using na, (Y 5 1, processors 

with a total data traffic of O(n1t"/2). Then it is shown that the data traffic in factoring the 

matrix is fl(n'+'/2) when the load is distributed uniformly among na processors, (Y _< 1. 

10 
I 

It is also shown that in any scheme tha.t requires nu processors, a 2 1, the data. traffic is 

0 ( n, /2 ) . 

For fa.ctoring a sparse ma.trix efficiently proper ordering of the matrix is essential. Ordering 

of the matrix to be fmtored also determines the data. dependencies a.nd hence the data traffic 

a.ssocia.ted with any partitioning and a.ssignment scheme. For ma.trices associated with 

regular grid gra.phs, nested dissection is a. well known ordering scheme [6]. In the following 

a. few ba.sics of this ordering scheme a.re briefly described a.nd some notation is introduced 

that is necessary for the a.nalysis presented later. In the following it is assumed that the 

rea.der is familiar with the elementa.ry concepts underlying the nested dissection algorithms, 

a.nd the terms such as elimimtion order and the fill associated with the elimination process. 

It is also a.ssumed tha.t the reader is fa.milia.r with the basic gra.ph theory concepts rela.ted to 

matrix representations of systems of equations, in pa.rticnlar, the notion of vertices, edges, 

sepa.ra.tors, subgraphs of a. gra.ph, a.nd the correspondence between the vertices and the 

rows and columns of the ma.trix, between the edges a.nd a.nd the non-zero elements, a.nd the 

a.dded edges a.nd the fill-in during the factorization of the matrix. For details see [12] and 

[6] and the references therein. 

3.1 Nested dissection method as applied to 2-d grid graphs 

A nested dissection method may be viewed as a divide-and-conquer algorithm on an undirected graph. It relies on finding a small set of vertices, called the separator set, in the graph 

such that the removal of these vertices divides the graph approximately in half. Informally, 

the nested dissection method orders the vertices of the graphs as €ollo~~s. The vertices in 

the separator set are ordered last. Then the vertices in the subgraphs obtained irom the 

original graph by removing the separator are ordered recursively. In [12] a nested dissection 

algorithm is given for ordering the vertices of any graph G such that G and all subgraphs 

of G satisfy a fi-separator theorem. The ordering produced by this algorithm guarantees 

a O(n1ogn) fill and O(n3l2) sequential operation count for a system corresponding to an 

n-vertex graph G. In [8] a nested dissection algorithm is given for ordering the vertices of 

a graph G that has a &-separator decomposition.’ For a detailed treatment of the nested 

dissection methods and for the relevant practical applications see [6]. 

For the sake of simplicity and clarity, here only the systems corresponding to fi x fi 

regular grid graphs are analyzed. However, the techniques developed for analyzing data 

traffic complexities are applicable to other systems where thc nested dissection method can 

be used to give a “good” ordering. In the following, the nested dissection method used 

for ordering the vertices in a ,h x fi regular grid graph is briefly described. In the 

discussion, the grid graph is sometimes simply referred to as the grid and a subgraph of the 

grid graphs is referred to as a subgrid. For the rest of the discussion, it is also assumed that 

> 0 if G has a 

6-separator C and every connected component of G’- C has a 6-scparator decomposition. 

‘A graph G is said to have a h-separator dccompoaition for constants a < 1 and 

11 
1: I1 

.I 

Figure 4: A 7 x 7 grid with nested dissection ordering 

the vertices of the grid are connected according to a 9-point stencil, unless otherwise stated. 

Let V be the set of vertices of a 6 x fi regular grid graph. Without loss of generality, 

assume that fi = 2' - 1 for some integer 1. Let SO be the set of 2' - 1 vertices on a vertical 

mesh line, the removal of which partitions V into two subgrids, 171 and 1% such that the 

vertices of both the subgrids are arranged in a (2' - 1) x (2I-l - 1) mesh. The vertices of 

SO are numbered from n - 6 + 1 to n in any order. Suppose that VI is the left subgrid 

and Vz is the right subgrid. Let S1 be the set of vertices on a horizontal mesh line that 

dirides 1'1 into two equal parts each containing (2I-l - 1)2 vertices that are arranged along 

a (2I-l - 1) x (2'-' - 1) square mesh. Similarly, let S2 be the set of vertices from 1.5 which, 

when removed, produce two equal halves from 15. Both SI and S2 contain 2I-l- 1 vertices. 

Let the vertices in S1 be numbered from n - 2fi + 2 to n - 3&/2 + 1/2 and those in S2 

he numbered from n - 3&/2 + 3/2 to n - fi. Thus, the removal from 1' of the vertices 

in the set SO U SI US, partitions V into four (6 - 1)/2 x (6 - 1)/2 subgrids. The 

separator set So U S1 U S, is referred to as the "+"-separator for the grid corresponding to 

1'. The middle vertical part of the "+"-separator is referred to as the vertical sub-separator 

and each of the two horizontal halves of the "+"-separator is referred to as the horizontal 

sub-separator. All the vertices of the four subgrids are numbered by recursively identifying 

and ordering the vertices on the "+"-separators of the subgrids induced by the vertices 

ordered so far. The recursion stops when a subgrid has only one vertex on it. For any 

"+"-separator, there is a vertical sub-separa tor and two horizon tal sub-separators. With 

the above described ordering scheme, for any given "+"-separator, the vertices on the two 

12 
horizontal sub-separators are given numbers that are smaller than those assigned to the 

vertices on the corresponding vertical sub-scparator. Thus, we say that the vertices on 

a horizontal sub-scparator are ordered ahead of the vertices on the corresponding vertical 

sub-separator or that the vcrtices on the vertical sub-scparator are ordered after those on 

the horizontal sub-separators. An example of ordering the vertices in a 7 x 7 grid is shown 

in Figure 4. Observe that the grid is recursively partitioned into four subgrids by a set of 

vertices that form a "+"-separator. 

To label the subgrids and the separators of the grid graph, we use the notation given in 

[7]. Each subgraph and the separator that induces the subgraph are given a level number 

depending on the recursion level of the nested dissection on which the subgraph is ordered. 

Under this scheme the original grid is called a level-0 (sub)grid. The four subgrids of siie 

(fi - 1)/2 x (fi - 1)/2 are the level-1 subgrids. The "+"-separator that partitions the 

level-0 grid into the four level-1 subgrids is called the level-1 "+"-separator or simply as the 

level-1 separator. Thus, if n is equal to (2' - 1)2, there are 1 levels of subgrids numbered 0 

through 1 - 1 and 1 - 1 levels of separators, numbered 1 through 1 - 1. 

In the following it is a.ssumed tha.t the matrix to be factored is ordered using the nested 

dissection scheme and thai the symbolic factorization step is already completed. 

3.2 Cholesky factorization scheme revisited 

Consider the Cholesky fa.ctorization scheme described in Section 2.1 for factoring a. sparse 

symmetric, positive definite matrix. 

Clearly, the main difference between factoring a sparse and a dense matrix using the 

Cholesky factorization scheme is that in the former case there is no need to modify column j by all columns to the left of it. Specifically, column j is modified only by columns 

k €or which lj,k # 0. Moreover, i€ column k modifies column j, only the nonzero elements 

of column k need to be fetched. Exactly which elements are needed is formalized later. In 

Figure 5(a), the zero-nonzero structure of L, corresponding to the vertices of the separators 

on the first two levels, is shown schematically. The shaded areas represent the nonzeros. 

The corresponding grid is shown in Figure 5(b). It is clear from the figure that only certain 

values from certain columns are needed for computing an element of the factor. 

Another important difference is tha.t, because of the ordering applied, several columns may 

be computed simultaneously. As stated earlier, column i and row j of the matrix corresponds 

to a vertex vi in the elimination graph and the factoring of the matrix corresponds to the 

elimination of the vertices. Thus, all the vertices on the level 1-1 subgrids may be eliminated 

simultaneously followed by those on the level 1 - 2 and so on. This observation is useful in 

extracting parallelism in the factorization step. 

13 
' ' '. 

I, 

Figure 5: Structure of L 

3.3 The worst case data traffic complexity 

In this section a bound on the worst case data traffic complexity for factoring the matrix A 

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

memory is not allowed. Thus, an upper bound on the worst case data traffic is obtained 

by assuming that the values of all the elements of the lower triangular part of matrix 

A and those of L, as well as any intermediate results, are stored in the shared memory. 

Suppose also that any number of processors are allowed to participate in computing a 

nonzero element of the factor provided that no computation is repeated. Consider the 

computations associated with a nonzero element l;j E L. Recall that in computing li,j, first 

aij -xi=: li,k-lj,k is evaluated and then the resulting value is divided by lj,j. 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. Therefore, in the worst case, 

each multiplication operation in the Cholesky factorization is associated with a constant 

amount of data traffic. The following theorem gives a bound on the worst case total data 

traffic. In the proof of the theorem, the result given in Theorem 8.1.8 of (61 is assumed. 

That theorem states that the number of operations required to factor a matrix associated 

with an n-vertex 2-D grid ordered by nested dissection is given by 829n5l2/84 + O(n.log n). 

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

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

of computation that are more restrictive. 

14 
Theorem 3 The worst case data traffic associated with factoring the rnatriz A is O(n3l2). 

Proof: Associated with each multiplication operation in the factori7ation there are at most 

a constant number of memory refcrences. Suppose that k memory references are involved 

per multiplication. Thus, the total data traffic is 

5 k . number of multiplication operations. 

Now, the number of multiplication operations associated with factoring matrix A is O(n3/') 

[6]. Hencc, the worst case total data traffic is O(n3/'). I 

Note that the above theorem is applicable to all the graphs for which a fi-separator 

theorem holds. 

3.4 Data dependencies for the sparse Cholesky factorization 

The worst case bound on the data traffic established in Theorem 3 can be improved for the 

model of architecture assumed in the case of the dense matrices. In that model, no element 

is fetched more than once from the shared memory and hence the values of the elements used 

in more than one operation are stored in the local memory associated with the processor. 

To maximize the potential of such a model, it is necessary to clearly understand the data 

dependencies involved. The vertices of the grid are ordered using the recursive nested 

dissection scheme. Hence it is sufficient to investigate the data dependencies involved in 

computing the elements of L in the columns corresponding to the vertices in a generic 

"+"-separator. This is accomplished in the next two lemmas. 

Let 7: = {klk 5 j a.nd li,k # O,/i,k E L}; i.e., 7: is the set of a.U columns of the factor L to 

the lcft of the column j + 1 such that the elements in row i of these columns are nonzero. 

Let fj!,k = u:=; 7;; i.e., q;,k is the set of all the columns to left of column j + 1 such tha.t on 

ea.ch of these columns there is a. nonzero element in at 1ea.st one of the i through IC rows of 

the factor. Let I' represent a.ny m-vertex sub-separator. It is assumed that all the vertices 

in any sub-sepa.rator are ordered consecutively. Let low(I') and high(I') be the indices of 

the lowest and the highest ordered vertices, respectively, on the sub-sepa.rator I'. Note that 

high(I') - low(I') + 1 = rn. In Figure 6, a. sub-sepa.ra.tor I' is shown. This sub-separa.tor 

separa.t,es the vertices in regions R1 and Rz. The diagonal and off-diagonal non-zero blocks 

a.ssocia.tcd with this sub-sepa.ra.tor a.re shown in Figure 7. 

The following lemma establishes some basic sub-separator related properties that are useful 

in analyzing the communication requirements. 

Lemma 3 Let I' be my m.-verter sub-separator. (i) Corresponding to the vertices ofr there 

is a dense m x m triangular diagonal block in the Cholesky factor. (ii) In the factor L, the 

15 
Figure 6: Sub-separator I' with four surrounding sub-separators 

columns low(r) through high(I') contain at most four of-diagonal rectangular blocks with 

nonzero elements. Each of these blocks is of size at most (c1 . m + c~) xm where c1 5 2 

and c2 5 3 are positive integer constants. Any nonzero element in these columns is either 

in one of these four blocks or in the diagonal triangular block. 

PTOO~: The first part of the lemma is obvious. In Figure 6, the sub-separator I' separates 

the vertices in regions R1 and Rz. Since the vertices in these two regions are ordered ahead 

of those of I', the fill due to the elimination of vertices in regions Rl and Rz ensures a dense 

rn x m triangular diagonal block bounded by columns low(r) and high(I') as shown in 

Figure 7. 

To prove the second part of the lemma, again consider Figure 6. In that figure, the thickness 

of the lines qualitatively indicates the separaior levels in the nested dissection ordering. Let 

rl, I'2, r3, and r4 be the four pa.rtia1 sub-separators that surround the sub-separator I'. 

Because of the nature of the nested dissection ordering, the rertices of I' are "connected" to 

only those higher ordered vertices that lie on I'll I'Z, I'3, and I'4 and to no other vertices.' 

Thus, all the noneeros on columns low(I') through high(I') in rows below high(l?) are 

confined to only the rows corresponding to the vertices on I'l, I'2, I'3, and r4. Furthermore, 

is said to be "connected" to vertex u if there exists a path [a9u~ta?,...,a~,wJ of length one 

or more in the grid grrph such that indez(r,) < min~ndez(r),indet(u)), for 1 < t 5 I.; in ouch a cane, 

I,,, E L is A non-zero, where i = indez(n), j = indez(u). 

'Vertex 

16 
each vertex in I’ 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 7. It 

can be verified that if r is a horizontal m-vertex sub-separator, then the surrounding box of 

vertices is of dimension (2m + 3) x (m + 2). Therefore there are two rectangular off-diagonal 

dense blocks of dimension at most (2m + 3) x m and the other two of dimension at most 

(rn + 2) x m. Similarly, if I’ 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 of these blocks will be missing. I 

Figure 7: Off-diagonal blocks with nonzeros corresponding io 

sub-separator I’ 

From the above lemma it is clear that, in computing the nonzero elements in the columns 

corresponding to the vertices on a sub-separator, only the data dependencies of the elements 

in the four rectangular blocks and the diagonal triangular block need be considered. This 

is accomplished in the following lemma where a bound is derived on the amount of data 

required in computing the nonzero elements lying on a given row and on one of the fire 

blocks. 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 sub-separator 

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, lotu(I’) 5 i 5 high(r), for some rn-vertex 

sub-separator I’, require a total of less than c - m nonzero elements from that row. Note 

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

17 
row i belongs. Thus, the computations a.t all the elements in a row of any of the five blocks 

specified in Lrmrna 3 require only c m. elements from that row, irrespective of the relative 

location of the off-diagonal blocks in the factor. 

Lemma 4 Let I' be any rn-vertez sub-separator. The nonzero elements from row i, i 2 

lo~u(I'), required in completing the computations of all elements 1;j E L ~vch that low(I') 5 

j 5 h,igh,(I'), are thoae elements in row i on the columns in the set given by, T,I,~~~(~),~~~~(~) n vi 

For all i greater than or equal to low(I'), Ilfj,ow(r),high(r) n qhigh(I')II ia at most c-m for some 

constant c. 

-high(r) high(r) . 

high(r) 

Proof: Any nonzero element 1;j E L, i 2 low(I') and low(r) 5 J' 5 h,igh(I'), is in one of the 

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

that int,ersect one of these blocks need to be considered. The result for low(r) 5 i 5 high(I') 

is proved first followed by that for i > high@'). 

high(1') high(r) high(r) high(r) high(r) 

conta.ins all the columns that have a nonzero element in row 

 hen lo4r) 5 i 5 high(r), q,ow(r),high(r) fl vi = vi , since, vi %w(r),high(r)' 

By definition, the set vi high(r) 

i. Clearly, the nonzero elements from the row i required in completing the computations at 

all the elements 1;j E L, low(I') 5 j 5 high(I'), are on columns in the set vi high(r) 

high(r) To measure the size of the set vi , 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 low(I') 5 i 5 high(I'), the size of the set vi is bounded by c - m, 

where c is a constant dependent both on the degree of the graph and on whether I' is a 

horizontal or vertical sub-separator. If I' is a horizontal m-vertex sub-separator then, for 

a 5-point stencil, c is equal to 7 and, for a 9-point stencil, c is equal to 11. When I' is a 

vertical rn-vertex sub-separator, the values of c are 5 and 7, respectively. This completes 

the proof when low(I') 5 i 5 high(I'). 

high(r) 

high(r) The case where i > h,igh.(I') is considered next. As shown above, 117; 11 depends on the 

size of the sub-sepa.ra.tor to which the vertex i belongs and hence, when i > high(I'), 

Ilq;high(r))l ca.n be much greater thaa O(m) where m is size of I'. However, when the 

compiita.tion of only those elements in row i that lie on columns low(I') through h,i.qh(I') a.re 

of concern, ea.ch of these computaiions consists of a. product of a. nonzero element in row i 

and a nonzero element in one of the rows low(I') through high(I') in the column high(I?) or 

in some other column to t.he left of it. Thus, for these computa.tions, only the columns that 

have a nonzero element in row i a.nd in row j, where lo?u(I') 5 j 5 high(I'), a.re of interest. 

The set v,ow(r.),high(r) consists of all columns tha.t ha.ve a. nonzero element in a.t 1ea.st one 

of the rows low(l?) through high(I'). Similady, vi high(r) consists of all the columns that 

ha.ve a. nonzero element in row i. Clearly, the set q,o,,,(r),high(r) -high(I') n v;high(r) consists of all the 

-high(T) 

18 
i 

columns which conta.in all the pa.irs of nonzero elements that must be used in completing the 

computations a.t all the elements l;j, lot@) 5 j 5 high,(I'). Thus, the nonzero elcrnents 

from row i, i > high(I'), required in completing computations a.t all the elements 1;j E L 

such that low(I') 5 j 5 high(I'), a.re those elements in the LOW i on the columns in the set 

given by vlow(r),high(r) n vi 

-high(r) high(r) 

-high(T) high(r) To get a. bound on the size of vlOu, r),high(r) n vi , consider the m-vertex horizontal 

sub-sepa.rat,or I' shown in Figure 6. t is surrounded by sub-sepa.rators I'll I'2, I'3, a.nd r4. 

Suppose high(r) n q!igh(r) consists of columns that lom(I'1) 5 i 5 high(rl). The set fjlow(I.),high(r) 

corresponding to vertices on I' or corresponding to those vert.ices ordered ahead of them 

which a.re "connected" to a.t. least one vertex in a.nd to the vertex corresponding to row i. 

Using the recursive ordering of the nested dissection scheme it can be shown tha.t the number 

of such vert,ices is less t,han 7m. Thus, ~~ijlow(r),high(r) n vi 11 5 7m, for low(I'1) 5 i 5 

high(rl). The sa.me bound is obtained when low(I'4) 5 i 5 high(I'4). If low(I'2) 5 i 5 

high(rz) or lo~(I'3) 5 i 5 high(r3) then it can be verified tha.t, llfj~o,,,(r),high(r) high(r) n 1; high(r) 11 - 

I 

f 

high(r) high(r) 

3m. If I' is vertical sub-separator the two bounds are 5m and 5m/2 respectively. 

3.5 A partitioning scheme with minimum data traffic 

In this section a partitioning scheme for computing the factor of the sparse matrix A is 

described. Suppose that an n x n matrix is to be factored using n" processors, a 5 1. 

The vertices of the fi x fi grid graph corresponding to this matrix are ordered using 

the nested dissection method described earlier. Assuming a = (2' - 1)2, the ordering results 

in 1 levels of subgraphs and 1 - 1 levels of "+"-separators. If the original fi x fi grid 

is considered to be on level 0, then on level i there are 22i level-i subgraphs each of size 

(2'-i - 1) x (2'-i - 1). Without loss of generality, assume that a - 1 is an integer. Thus, in 

the partitioning scheme described here, all the vertices on a level-al subgraph are assigned 

to the same processor. In that scheme, initially each processor independently computes the 

elements in the factor corresponding to a (2('-'')' - 1) x (2('-")" - 1) subgraph which are 

separated from one another by the level-aI separators. Once the elements in the columns 

corresponding to the vertices on the level-(/ - 1) through level-al separators are computed 

locaJly, a processor pi combines with three other processors to compute the elements on the 

columns of L corresponding to the vertices on the level-(a/ - 1) "+"-separator. The two 

horizontal sub-separa tors are computed by two processors and the vertical sub-separator of 

that level is computed by all four processors. The next lower level "+"-separator is computed 

in parallel by sixteen processors from the four neighboring groups. This is continued until all 

the vertices are eliminated. On each level of computation each group of processors computes 

the elements of the factor independent of the other groups. The elimination of the vertices 

on the vertical sub-separator of level-1 ifi computed in parallel by all processors. This 

corresponds to factoring a fi x 6 dense matrix. The computations corresponding to the 

level-i separator, i < a-l, are performed as follows. The computations corresponding to the 

19 
vertices on level-(al - k) "+"-sepa.ra.tor, 1 5 k < a. I, are completed by p = 22'k processors 

working in parallel. Using all t.he a.va.ila.ble processors, the factorization corresponding to 

thc m. x m tria.ngn1a.r diagonal block is first completed. Then the processors are used to 

compute the elements corresponding to the four off-dia.gona1 blocks. For the first part, the 

BLOCC factorization scheme described for the dense ma.trices is used. The m x rn dense 

diagonal block is partitioned into r2/2 - ~/2 square blocks and r dia.gona.1 triangular blocks 

ea.ch of size m/T x m/T whcre p = r2/2 + ~/2, and each of these p prtitions is assigned to a 

unique processor. Ea.ch processor completes the compiita.tions corresponding to its partition 

by accessing the required data. from the shared memory. For the purpose of faa.ctoring, the 

off-diagonal blocks a.re treated a.s if they were a.djacent, and the resulta.nt rectangular block 

is pa.rtitioned into p sub-blocks each of size c m/fi x m f&, where c 5 6 for a horizontal 

sub-separator and c 5 4 for a vertica1 sub-separa.tor. Again each pa.rtition is assigned to a 

separate processor. This process is repea.ted on the next lower level "+"-separator. Thus, 

in the assignment scheme described hcre, each processor is assigned a new subblock on ea.ch 

level and the size of the subblock assigned to a processor va.ries fiom one level to the next. 

Let this pa.rtitioning scheme be referred to a.s the sparse block oriented column Choleaky 

factoriza.tion scheme or simply as the sparse BLOCC scheme. Note that the underlying 

numeric algorithm is the column oriented ChoIesky factoriza.tion. 

Data traffic associated with an m-vertex sub-separator 

The sparse BLOCC scheme, described above, may be considered as a sequence of steps, each 

step corresponding to the elimination of vertices on the "+"-separators of some level. Initially, a single processor computes all the non-zero elements corresponding to a "+"-separator 

in the factor. As the computation proceeds, more than one processor work together to compute the elements corresponding to a "+"-separator. On any such step, first the non-zero 

elements in the columns corresponding to the horizontal sub-separators are computed and 

then those in the columns corresponding to the vertical part are computed. Here we analyze 

the data traffic associated with any one step, on which p processors combine together to 

compute the elements corresponding to a sub-separator. 

By Lemma 3, for any sub-separator I? there are at most five non-zero blocks in the columns 

corresponding to the vertices on I'. The number of non-zero blocks is five when I' is enclosed 

within a rectangular box formed by the sub-separa tors with vertices that are ordered after 

those on I' (see Figure 6). The following lemma gives a bound on the data traffic associated 

with computing the elements in the coliimns corresponding to such sub-separators. Not all 

sub-separators are enclosed by such rectangular boxes. In such cases there are less elements 

to be computed and consequently there is less data traffic. For the sake of simplicity of the 

analysis, it is assumed that no element of the factor nwded in the computation of the five 

non-zero blocks is initially in the local memory of any of the p processors. Thus, the data 

traffic given below is a conservative estimate. 

20 
Lemma 5 Let I' be uny m,-vertex sub-separator and p be the number of processors available 

for computing the ekmenta of the factor in all the non-zero blocks within the columns lo?iJ(I') 

through high(I'). If I' is an m-vertex horizontal sub-seprator, then the associated data 

traffic is at most (53 + ll&)m2 - fi. If it i.9 a vertical aub-separator, then the data trufic 

is at most (28 t 8fi)m2 fi. 

Proof: Let I' be an m,-vcrtcx horizontal sub-separator that is enclosed completely within 

a rectangular box formed by the sub-separators whose vertices are eliminated after the 

rcrtices of I'. Such a sub-separator has the worst case communication requirements among 

all thr m-ver tex sub-scpara tors. 

First, consider the data traffic associated with computing the elements of the factor in the 

triangular diagonal block using p = r2/2 + r/2 processors. Each of the sub-blocks requires 

nonzero elements from at most 2m/r rows out of the m rows in the range loiu(I') through 

high(I') of the factor. From the proof of Lemma 4, 

each of these rows has at most llm nonzeros. Thus the communication requirement of 

each partition is at most llm - 2m/& and the total communication requirement of the 

triangular block is bounded above by 11&m2 . &. 

No other information is nccdrd. 

Now consider the data traffic associated with the off-diagonal blocks. Each partition is 

of size 6mlfi x m/&. Thus, each partition requires nonzero elements from 6m/& rows 

which are below the row high(I') in the factor. From the proof of Lemma 4, each of these 

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

partitions. Each partition also requires information from rn/@ rows from the region lou!(I') 

through high(I'). Each of these rows has at most llm nonzeros. Thus, the communication 

requirement of each partition is at most 7m .6m/& t llm . m/& = 53m2/& and the 

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

using p processors is less than or equal to 53m2J?i. 

Adding the communica.tion costs corresponding to the diagonal and the off-diagonal blocks 

we get the total daia traffic associated with I' to be less than or equal to (53+11&)rn2.fi. 

A similar analysis can be used to compute the data traffic when I' is an m-vertex vertical 

sub-separator and can be shown to be bounded a.bove by (28 + 8fi)m2. fi. I 

The total data traffic of the sparse BLOCC scheme 

Applying the resiilts from the above lemma, a bound is obtained on the total data traffic of 

the sparse BLOCC scheme. First some notation is introduced. Let q(m,p, k) represent the 

data traffic using p processors in completing the computations at all the nonzero elements 

1;j E L in the columns corresponding to an m-vertex horizontal sub-separator that is 

21 
surrounded by higher ordered vertices on k sides. Let r,(m,p, k) represent the same for an 

m-vertex vertical sub-separator. From Lemma 5, q,(m,p,4) is at most (53 + 1lh)rn' .J?i 

and ~,,(m,p,4) is at most (28+8&)m2-fi. Let ~~(m',p,k) represent the total data traffic, 

using p processors, in completing the computations Corresponding to all the sub-separa tors 

within an m,'-vertex sub-grid that is surrounded by higher ordered vertices on k sides. Note 

that the quantities Th and 7, represent the data traffic corresponding to the vertices on a 

horizontal and a vcrtical sub-separator, respectively, whereas rg represents the data traffic 

corresponding to the vertices on an entire sub-grid. 

The following theorem gives an upper bound on the total dat,a traffic in factoring the 

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

scheduling scheme as described above. 

Theorem 4 The total data traflc in factoring the n x n sparse matrix A, using n" 

processors, is ~(n'+"/'); i.e., ~~(n, nu, 0) = ~(n'+"/'). 

Proof: On an n1l2 x 71'1' regular grid there is an n'/'-vertex vertical sub-separator 

and two n1/'/2-vertex horizontal sub-separators (ignoring the additive constant -1). The 

vertical sub-separa tor 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 

by such vertices only on one side. These three sub-separators subdivide the n-vertex grid 

graph into four sub-grids of size n'/'/2 x n1/'/2, each surrounded on two sides by higher 

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

given by, 

1 1 11 

2 2 44 7g(~lnp,o) = ~,,(~1/2,?&Q,0) + 27h(-?&'/2,-nQ,1) +47g(-n,-n0,2). 

A recursive expansion of the above expression contains data traffic terms for vertical subseparators 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-separa tors of different sizes surrounded in five different ways. 

To keep the analysis simple, it is assumed that all the four sub-grids of size n1I2/2 x n'j2/2 

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, 

1 1 11 

2 2 44 Tg(n,n",O) 5 Tv(n1/2,.",0) + 2Th(-n'/', -nOJ) + 4.r,(-n, -nQ,4). 

Now, 

From Lemma 5, it follows tha.t, 

11 

4 '4 16 Tg( -n -n", 4) 5 

22 
An analysis simi1a.r to that given in Lemma. 5 yields 

and 

Thus, we get 

3.6 A lower bound on the data traffic complexity 

In the following theorem, the communica.tion bound of the spa.rse BLOCC scheme is shown, 

by giving a lower bound proof, to be optimal in an order of ma.gnitude sense. 

Theorem 5 Under the condition of uniform load distribution, the data trafFc in factoring 

the n x n sparse matrix A, using no processors, (Y 5 1, is fl(n1+a/2). 

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

ordering is n1l2 [12]. From Lemma 3, it follows that the factor Z has an n1l2 x nl/' dense 

triangular diagonal block incorporated in it. F'rom 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 no processors, is fl(nlta/'). Since 

the factorization of A cannot be completed without completing the factorization of this 

dense block, the result follows. I 

From Theorem 4 and Theorem 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 ma, (Y > 1, processors are used, the data traffic bound 

given in Theorem 3 holds. 

4. Concluding remarks 

In this paper we have analyzed the data dependencies in the Choleskg factorization of 

dense and sparse symmetric, positive definite matrices. The model of computation assumes 

23 
a multiprocessor system with a memory hierarchy. Based on this analysis it is shown that 

iindcr the condition of uniform load distribution the data traffic associated with factoring 

an n x n dense matrix is fl(nZta/') when nu, cr 5 2, processors are used. The same 

is shown to be ,(TI'+~/'), cr 5 1, for factoring an n x n sparse matrix representing a 

2-dimensional regular grid graph where the vertices are ordered using the nested dissection 

ordering methods. Block based partitioning schemes are presented that asymptotically 

achieve these bounds on the data traffic. 

The sequential computat,ion time for fa.ctoring the n x n sparse matrix A is 829n3/'/84 + 

O(n1og n,) [6]. As sta.ted for the dense matrix ca.se, the assumption here is that the computation cost of each step of the innermost loop is one and costs involved in the other steps are 

ignored. Under the sa.me assumption, it ca.n be shown that the computa.tion time for the 

sparse BLOCC scheme is a.t most 283n3/'-"/4 if nQ processors are used. In [7], a parallel 

scheme for factoring the ma.trix A on a multiprocessor system is given that is a.na.logous 

to the wrap around a.ssignment scheme described in the Section 2.5 for dense matrices. 

This scheme has the property of distributing the work evenly a.mong the processors. The 

computa.tion time to fa.ctor the sparse matrix A on nu processors with the wra.p around 

scheme is at most 197n3/2-a/4. However, the data tra.ffic associa.ted with that scheme is 

less tha.n or equal to 183n1+"/4. Note that the difference in the computation time with 

the BLOCC scheme a.nd with the wrap around a.ssignment scheme is less than a. factor of 

two. The BLOCC scheme is able to compute the fa.ctor efficiently in the case of the sparse 

ma.trices because the processors are now assigned blocks in a. wra.p around fashion which 

tends to distribute the 1oa.d evenly. On the other hand, the da.ta traffic associated with the 

BLOCC scheme is an order of ma.gnitude less than tha.t for the wra.p a.round assignment 

scheme. Moreover, in the former scheme, as many as n processors ma.y be used before the 

total data tra.ffic reaches the ma.ximum value of O(n31'), whereas in the Mer scheme only 

up to nl/' processors ma.y be used efficiently. The implications of the reduced data traffic 

on the performa.nce are as follows. 

The sparse BLOCC scheme reduces the communica.tion requirement to O(n'+"/') by removing the constra.int of column-level indivisibility. Here the indivisible work unit is the 

computa.tion corresponding to a nonzero element in the factor. The reduction in the communica.tion requirements is brought a.bout by improving the utiliza.tion of the da,ta. accessed 

from sha.red memory by each processor. Consider the fa.ctorization of a.n rn x rn dense 

ma.trix. Let the data utilization of a data element accessed by a processor be defined as 

the number of computations in which tha.t element is used by tha.t processor divided by 

m. Since an element in the fa,ctor is needed in a.t most m computa.tions, the maximum 

utilization of a.ny da.ta. a.ccessed is one. Let. the aggregate data utilization for a processor be 

defined a.s the a.verage utiliza.tion of the individual data. elements a.ccessed by tha.t processor. In the BLOCC scheme applied to a.n rn x m. dense matrix, ea.ch processor a.ccesses 

a.t most 2m,2/Jli elements from the sha.red memory and ea.ch element is used in a.t least 

m/fi computations. Thiis, the utilization of ea.ch da.ta accessed is at least l/fi and so is 

the aggregate utiliza.tion of all the data accesses. On the other ha.nd, with the column-level 

work assignment scheme, each processor accesses O(rn') elements from the shared memory. 

Of these, only O(m/p) elements ha.ve a. utilization of one a.nd the da.ta utilization for the 

24 
remaining elements is l/p which gives an aggregate data utilization of approximately l/p. 

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 utili7ations. 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 x m dense matrix using p processors is 

O(m2 - &), on an average each processor accesses only O(m2/&) data. Note 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(m2 -p) and thus, on an average each processor accesses O(m2) data or on 

an average each element in the shared memory is accessed by O(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. 

As a final remark, note that the data traffic analysis for the sparse BLOCC scheme exploits 

the fact that the underlying graph satisfies a &-separator theorem. Thus, similar schemes 

may be developed for any class of graphs satisfying an f(n)-separator theorem [13]. In 

such cases the data dependencies, the fill, and the computation time depend on f(n). In 

[12] the fill and the bounds on the sequential computation time for various values of f(n) 

are listed. Here we state the bounds on the corresponding data traffic when the systems 

are computed using n" processors. The data traffic of factoring a matrix corresponding to 

an n-vertex 3-dimensional regular grid using nu processors is O(R~/~+"/~). For that case 

the computation time is O(n2-"). In general, the total data traffic for computing a factor 

of a matrix corresponding to a &dimensional grid is O(n2a+a/2), where c = 1 - l/d. The 

computation cost is O(n3a-"). For an n-vertex finite element graph* with no element having 

more than k boundary vertices, the total data traffic in factoring the associated matrix is 

O(k2 - n1+"I2). The computation time is 0(L3 - n3j2-"). All these quantities are optimal in 

an order of magnitude sense. 

*A finife elcmrnf graph is any graph formed from a planar embedding of a planar graph by adding all 

possible diagonals to each face; i.e., there is a cliqtie corresponding to each face of the embedded planar 

graph. 

25 
References 

[l] P. Charrier and J. Roman. Study of the Parallelism Induced by a Nested Di~ection 

Method and of its Implementation on a Message Passing Multiprocessor Computer. 

Technical Report 1-8722, Universite de Bordeaux I, Talence, France, 1987. 

[2] J. A. George, hi. T. Heath, and J. W. H. Liu. 

a shared-memory multiprocessor. 

1986. 

Parallel cholcskg iactorization on 

Linear Algebra and Its Applications, 77:165-187, 

[3] J. A. George, R.I. T. Heath, J. W. H. Liu, and E. Ng. SoZutaon of Sparse Positive 

Definite Systems on a Shared-Memory Multiprocessor. Technical Report OItNL/TM10260, Oak Ridge Na.tiona1 Laboratory, Oak Ridge, Tenn., 1987. 

[4] J. A. George, hi. 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. 

[5] J. A. George, hl. T. Heath, E. Ng, and J. W. H. Liu. Symbolic cholesky factorimtion 

on a local-memory multiprocessor. Parallel Computing, 5:85-95, 1987. 

[GI J. A. George a.nd J. W. H. Liu. Computer Solution of Large Sparse Positive Definite 

Systems. Prentice-Hall, Inc., Englewood Cliff, NJ, 1981. 

[7] J. A. George, J. 117. H. Liu, a.nd E. Ng. Communication reduction in pardel sparse 

cholesky factorization on a hypercube. In M. T. Heath, editor, Hypercube Multiproces.~ors 1987, pages 576-586, SIAM Publication, 1987. 

[8] J. R. Gilbert and R. E. Tarjan. 

Numerical Mathematics, 50:377-404, 1987. 

The a.nalysis of a nested dissection algorithm. 

[9] J. R. Gilbert and E. Zmijewski. A Parallel Graph Partitioning Algorithm for a 

Message-Pas.~ing Multiprocessor. Technical Report TR 87-803, Department of Computer Science, Cornel1 University, Ithaca., NY, 1987. 

[lo] G. H. Golub and C. F. Van L0a.n. Matraz Computations. The Johns Hopkins 

University Press, Baltimore, MD, 1983. 

[ll] hi. T. Heath. Parallel Cholesky Factorization an Me.lrsage-Passing Multiprocesaor 

Environments. Technical Report ORNL-6150, Oak Ridge National Laboratory, Oak 

Ridge, Tenn., 1985. 

[12] R. J. Lipton, D. J. Rose, a.nd R. E. Tarjan. Generalized nested dissection. SIAM 

Journal of Numerical Analysis, 16:346-358, 1979. 

26 
[13] R. J. Lipton and R. E. Tarjan. A separator theorem for planar graphs. SIAM 

Journal of Applied Mathematics, 36:177-189, 1979. 

[14] J. W. H. Liu. Coinputa.tiona1 models and task scheduling for parallel sparse cholesky 

factorization. Pa.ralle1 Computing, 3:327-342, 1986. 

[15] V. K. Naik. On the Computation and Communication Radeoffs and their Impact 

on the Performance of Asynchronous Multiprocessor Sy.qtem.9. PhD thesis, Computer 

Science Department, Duke University, Durham, NC., 1988. 

[16] V. K. Naik and AI. L. Patrick. Analysis of communication requirements of sparse 

cholesky factorization with nested dissection ordering. In G. M. Rodrigue, editor, 

Parallel Processing for Scientific Computation, chapter 2, pages 9-14, SIAAI Publication, 1989. 

[17] Y. %ad. Commiinica.tion complexity of the gaussia.n elimination algorithm on multiprocessors. Linear Algebra and Its Applications, 77:315-340, 1986. 

[18] P. H. M'orley and R. Schreiber. 

a.rray. 

Environments: Parallel, Vector, and Systolic, 1986. 

Nested dissection on a mesh-connected processor 

In A. M'ouk, editor, Proceedings of the ARO Workshop on New Computing 

[19] E. Zmijewski a.nd J. R. Gilbert. A Parallel Algorithm for Large Sparse Symbolic 

and Numeric Cholesky Factorization on a Multiprocessor. Technical Report 86-733, 

Depa.rtment of Computer Science, Cornell University, Ithaca, NY, 1986. 

27 
Report Documentation Page 

1. Report No. 2. Government Accession No. 

NASA CR-181863 

ICASE ReDort No. 89-40 

3. Recipient's Catalog No. 

4. Tide and Subtitle 

Data Traffic Reduction Schemes for Cholesky 

Factorization on Asynchronous Multiprocessor 

Systems 

7. Authorlsl 

Vijay K. Naik 

Merrell L. Patrick 

Langley Technical Monitor: 

Richard W. Barnwell 

5. Report Date 

June 1989 

6. Performing Organization Code 

8. Performing Organization Report No. 

89-40 

10. Work Unit No. 

Final Report 

9. Performing Organization Name and Address 

Institute for Computer Applications in Science 

Mail Stop 132C, NASA Langley Research Center 

Hamp ton, VA 23665-5225 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 

Langley Research Center 

Hampton, VA 23665-5225 

and Engineering 

Proc. of 1989 Am Internationa 

Conf. on Supercomputing 

505-90-21-01 

11. Contract or Grant No. 

NAS 1- 1 8 10 7 

NAS1-18605 

13. Type of Report and Period Covered 

Contractor Report 

14. Sponsoring hgancy Code 

16. Abstracl I 

Commiinica tion rcqiiircmen ts of Choleskp factorization of dense a.nd sparse/ 

17. Key Words (Suggested by AuthorW) 

Chdles ky factorization , communication 

requirements, asynchronous multiprocessor systesm. local and shared 

memory 

18. Dirtribution Statement 

59 - Mathematical and Computer Science: 

Unclassified - Unlimited 

NASA-Langley, 196s 

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

Uncl as s i fed Un c las s i fie d 31 

22. Price 

A0 3