NASA Contractor Report 178016

leASE REPORT NO. 85-48

NASA-CR-178016

19860007475

leASE

MULTIPROCESSOR SPARSE L/O DECOMPOSITION

WITH CONTROLLED FILL-IN

... - .. .- "

,

,

.'.

,':

.... ". ,,'

...

Gita Alaghband

Harry F. Jordan

Contract No. NAS1-17070. Grant No. AFOSR 85-1089

October 1985

INSTITUTE FOR CO~WUTER APPLICATIONS IN SCIENCE AND ENGINEERING

NASA Langley Research Center, Hampton, Virginia 23665

Operated by the Universities Space Research Association

NI\SI\

National Aeronautics and

Space Administration

Langley Research center

Hampton. Virginia 23665

J/\••i\lI 6() ,-J ? '1':3" . i(~ '..'

.JJ,r;c,L.U I\ES[hf<Cd CEtHER

LJ CRARY, r~p.s/\

! !M"f':iON, VIRGINIA

m

Multlprocessor Sparse L/U Decomposition 4

with Controlled Fill-in

Gita Alaghban(I

and'

]larry F. Jordan

') l)epartment of Electrical and (.( mputer Engineering

University of Colorado

Abstract

l)urin_ I_/1; decoml)osition of a sparse matrix, it is possible to l)t.rf()rm ('(mil)Utati_)n on many diagonal elements simultane()uslv. Pivots that can l)e processed in l)ara]le] are related l)v a compatibililv relation and are gr(,ul,(,d in a

('ompatil)le set. The collection of all maximal compatibles yields difft, rent

maximum sized sets of pivots that can be processed in parallel." Generation of

the ma×ima] compatil)les is based on the information obtained from an

incompalil)le table. This table provides information about pairs of incompatible pivots. In this paper, generation of the maximal compatibles of pivot

elements for a class of small sparse matrices is studied first. The algorithm

involves a binary tree search and has a complexity exponential in the order of

the malrix. Different strategies for seh, ction of a set of compatil)h, pivots

based on the _larkowitz criterion are investigated. The competing issues ¢)f

parallelisn_ and fill-in generation are studie(l and results are l)rovided. A

technique for obtaining an ordered compatible set directly from the ordered

incoml)atil)le table is given. This technique generates a set of coml)atibh,

pivots with thel)rol)erty of generating few fills. A new hueristic algorithm is

then proposed that combines the idea of an ordered compatible set with a

limited binary tree search to generate several sets of compatible pivots in

• linear time. Finally, an elimination set to reduce the matrix is selected.

Parameters are suggested to obtain a balance between parallelism and fill-ins.

, Results of applying the proposed algorithms on several large aI)plication

matrices are presented and analyzed.

¢Research was supported in part by NASA Contract No. NASI-17070 and by the Air Force Office of

Scientific Research under Grant No. AFOSR 85-1089 while the authors were in residence at ICASE, NASA

Langly Research Center, Hampton, VA 23665.

:i.

Introduction

Solution of a linear svstem of equations is required in many al)l)li('ati(m

programs. One such area is tile VI,SI circuit simulation l)r(_grams, h;ver.v

computer-aided circuit analysis program includes a routine that solves a system of sparse linear equations. If implicit integration is used, at every tirol.

step one must solve a s.vstem of nonlinear e<l.uati+ms(usually by Newlon ileralion). At every iteration a system of linear equations mus! be s<dv+,(I.

Del)en<ling on the integration m+,:h<_(l,the number of times that a sparse system of linear equations needs t<) be s<)l\'e<!mav be large. If it is lmssibl+, 1+,

reduce tile solution time for the sparse s\'st¢,m, the total circuit analysis tim,

would l)e significantly reduced. One method f,,r solving such a system is the

factorization of the matrix into lower an+l Ul)l>crtriangular matrices fc)lh)wc(I

by forward and back suhslitutions.

One ]>romising area for advances in solution technique is lhe use of l)ara]-

lel COml)uters and parallel algorithms. Our l)revious work on parallelizing the

MA28 [1] sparse matrix package for the IIE;l) [2] multil)rocessor sugg+,sts that

sufli<.ient parallelism is not ol)tainable in sparse I,/U decompositi<,n withoul

processing multil)le l+iv+)tsin paralh, l [3]. lh_rallel pivoting st rat egi+,s h:l','e

been investigated bv ('alahan [.l] and n,orc rt,ccntlv by Wing and llaung [5].

[61.Jess and Kees [_:]and Peters [g]. Al[h(mgh the number ofol,erathms i,ossibh, in parallel may be large in a very sparse system, exploitation of all the

available parallelism may significantly in<.rease the generation of fill-ins (zero

element of the matrix becoming nonzero as a result of elimination). Sin<'e

fill-in increases the total COml>utation work, it is important to keep the

number generated under control. The purpose of this work is to stud,," spars+,

L/I_; decomposition on a multiprocessor by means of an algorithm which

exploits parallel pivots and keeps fill-in low. The class of sparse systems

guiding the study will be those arising from the simulation of VLSI <.ircuil.,

using a program such as SPICE [9].

\Ving and llaung in [5] represent the triangulation process by a directed

graph where the vertices represent a divide or update operation (operations

required for performing the triqngulation}, and the edges determine the precedence relation of the operations to be executed. By assigning level numbers

to the directed gral>h, thev identify all operations on the same level to be

done in parallel. They usea weighted <'ombinati+m offll-in cost an+1 depth of

computation in a heuristic to determine a nearly optimal pivot sequence.

\\'hile \Ving and ]laung identify all the operations that can be done in parallel. `'ve will identify all pivots that can be pr+>cessed in parallel at each stel).

" An issue that has not been discussed in tile literature is that in a sparse

matrix there are usually different sets of possible pivot candidates for each

, stel), and the sizes of these sets may well vary. It seems imi)ortant to study

these possibilities and the effect of parallel l)ivoting on apl)lieation matrices.

Algorithms identifying parallel pivot candidates are complex, so it, will be of

value t.o come up with such algorfl.hms only if the amount of parallelism in

circuit domain matrices is large enough to justify the computation required

to identify it.
In this l)al)er, x_t. assume a sharetl-meni,,r.v, NIINII) model for our 1,ar:llh,!

c()ml)utalion, ill which the total memory address space is accessible uniformly

to all parallel units (I)rocesses or individual pr(wessors). This COml)Ul:llicmal

model sl,ould provide s.vnchronization mechanisms to allow multiple memory

updates. Ifmvltil)le updates are aimed at the same memory cell, the penally

paid is a shorl dela.v in access time. Based on ll, is computational model, the

fiirsl h:llf of this paper is devoted to stud,,' the amo.nt of parallelisnl lha!

exists in application matrices. This is carried out by producing all l)ossible

sets of pivot (.andidalcs which cnn be processed in parallel at each step for a

number of small ma!rices. Observations are then made on different slra!egit,s

for choosing one of the sets produced at each step. and hence the gt.neralion

of fill-ins nnd possible parallel pivoting steps. The complete and del'liled

analysis of this study leads us into the second half of the paper, where we

describe a fast heuristic algorithm to produce a set of acceptable parallel

l)ivo! ean(ii(lales for reducing the matrix at (inch step. Issues inv,)lved in

balancing paralh'] work and fill-in generation are discussed and verified

through simula!ed results.

Parallel Pivot Candidates

The Iri:lngulnti()n method used here as men!ioned above will I)e sparse

l,/l" (lt,(.()ml)ositi()n. l"or simpli(.ily, we only consider lhe diagonal clemt,nls of

the ma!rix as pivot (.andi(lat e.'s. Note !hat I)i\'o!ing usually, refers !(_ unsvm-.

me!ric l)ermulali()ns of the matrix for swal)l)ing an off-diagonal ma!rix element with a diagonal clement. In this l)apcr, we art" only considering symmetric permutalions of the matrix. Even though we are not pivoting in the

al)ove sense, the le!ins pivot and pivoting are used throughout lhe pal)er to

refer to lhe diag()nal ehqncnt used to reduce the matrix at a given step anti a

svmme! ric l)crmu! alton respect ivelv.

In a sparse matrix, two pivots ai; and a.. can be processed in parallel if

a.. and a.. are both zero. In other words, _uring elimination, row j is not U

involved _1_1the elimination process taking place for pivot aii, and row i is not

involved in the process for a... This statement can only be true if wel)rovide

.L/ . ,

correct synchronizations for simultaneous update during the elimina!i,)n with

l)arallel pivot candidates:

1. During elimination, when processing pivots all,ayi,.., in parallel, it is l)()ssib]e thai an element of a nonpivot row needs to be updated bv all or

some of the parallel processes handling pivots i,j,.., for the current step.

In order for each process to obtain a coml)letely updated value, as a

result of a previous update, the update operation must be done asynchronously by" parallel processes. On the other hand, the order in which

parallel processes update an element is of no importance (except for

round off errors).

2. During elimination, when processing Divots a..,a..,.., in parallel, it is pos- * I! ]

sible that a fill-in is generated in position (m,n)J{ It is also possible that

more than one process tries t,o generate a fill-in in the same position

(re,n). The position (re,n) for the fill-in must be created once by one

process only, and other processes will update its value as in 1.
If two pivots all and aj). can be processed in parallel, and if aii and %.k

can also 1)e processed in parallel, then ai;,aii, and akk cannot necessarily be

• processed ill parallel. Tile relation between parallel pivot candidates is

retie.rive and symmetric, but not transitive, and is thus a compatibility relation. Two pivots related in this way will siml)ly I)e said to I)e comi)z_lillh, in

what folh,ws. A consequence of the nonlransitivity (,f the coml,alil,ilily relalion is that it classifies tile elements of a set into nt,ndisjoint sul)sels, so that

all meml)ers ofa sul)set are compatible. These subsets arecalh, d compatil_ilit)" classes. Thus, in order to come up with all l)ossil)le sets of pivots that ('an

1)e l)ro('essed in parallel and are of maximum size, we need to find all maximal

COml)atibh,s. A maximal compatible is a eonlpatible that is not included in

any larger compatil)le.

To clarify the discussion, we define a boolean matrix B for each sparse

matrix A, such that:

b..= 1 qf a..: IJ€:O

bU..= 0 otherwise

'1

where b.. and a.. denote elements of B and A respectively. U tj

Several apllroaches for constructing the set of maximal compatil)h,s exist,

and the,v art, all based on construction of an incoml)atible table [10]. The

incoml)alible table specifies pairs of incoml)atible elements. Assume pivots

are taken from the diagonal elements of tile sparse matrix and are numbered

1 through n corresl)onding to diagonal elements of rows 1 through n. Now we

could repr(,sent the incoml)atible table as a table consisting of {n-l)(.(_lumns,

where each column i has (n-i) elements. Columns of the table correspond to

pivot elements of the matrix. Column one of the tal)le, (.orresl)(mding to

pivot numl)er one, is set to the bit vector resulting from oring row and

column one of the matrix B and keeping the last (n-i) elements. The same

process is repeated for pivot 2 (column 2 of tile table), for the submatrix

ol)tained from the original matrix with row and column one eliminated. For

every column of the table that is completely constructed, the corresponding

row/column of the matrix is eliminated. The process is repeated for all pivots

in order. ]t is important to note that the incompatible table is constructed for

a given ordering of tile sparse matrix. Thus, there are n! different incompati-

. ble tables for n! possible diagonal orderings of an n by n sparse matrix. In

what follows, we represent the incompatible table as an array of dimension n,

say implbi(n), with elements of the array being sets of at most, n elements

• each. Each set corresponds to a column of tile table. As an illustrative example. the incmnpatible table for the matrix A1 of Fig. 1.1.a is given in Fig.

1.1.b.

The maximal compatibles are found by combining the pivot-pairs from

the ineoml)atible table into larger groups with coml)atible elmrlents. Several

s.vsteinat ic approach es for extract ing th e maximal compatibles h ave been su ggested, and they all use all exhaustive search routine. The one approach that

seems to be more suitable for programming on a digital computer is one that

assumes initially that all pivot candidates can be grouped into one set. Then
4

1 2 3 4 5 6 7

1 X "

2 X X

3 x x x x

4 x

5 x x x

6 x

7 x x x

Matrix AI

Fig. l.l.a

2

3

4 x

5 x x

6 x

x x ]

1

7

1 2 3 4 5 6

Incompatible Table

Fig. 1.1.b

tile information from the incompatible table is used for contradictions and

splitting tile groups where necessary. This procedure involves searching a

binary tree. Initially. it is assumed that all pivots are coml)atible. They are

grmlped ill one set consisting of all pivot elements. This set will I)e at tire

ro¢_t of a binary tree, level zero. Next, the set of pivots incompatible with

pivot number one, ol)tained from the illcompatible table, is used to split the

set at the root into a left and aright set., constituting level one. The left set

consists of all elements of its parent, set at level zero, except those incoml_atihie with pivot one. Tile right set. consists of the same elements as the starting set. (parent set,), except pivot one itself. At. the next step, the incompatible information for pivot number two, is used t.o break each set at level one

into a left, and right set, for level 2. Furthermore, since the matrix is sparse,

some of the sets at, a given level will not split into smaller sets for some pivoting elements, but they may still consist of ineolnpatible elements and will

split for some later pivots. Consequently, tire binary tree corresponding to

this search will not always bea dense tree. This process is repeated until no
more splitting of the sets is possible. The I(,af s[,ts are tl, en checked and

every set included in a larger leaf set is eliminated. Tile. remaining st,lseonstit uteallpossiblemaximalcompatibles. Note that the length ofapath from

the root 1o a leaf could be at most n.

q,

The above process is shown for the example matrix of Fig. 1.1 in Fig i.2.

lnilially, pivols number 1 through 7 are grouped together as the starting set.

Column one of the incompatible table indicates that pivot 5 is incompatible

with pivot one. Thus the starting set, is splil into two sets (1,2,3,4,6,7) and

(2,3,,t,5,6,7). At the next level, these two sets are broken into four sets, each

using the incompatibility information for pivot number lwo from the table.

This process is <'ontinued until no more. splits are possible. At the end, the

[1,°-,3,4,5,6,71

[1,2,3,4,6,71 [2,3,4,5,6,71

[1 [1,3,4,6,71 12,5,61 I3,4,5,6,71

a d

/

[1.3.4,7] 11,4,6,71 [3,4,71 [-t,5,6,71

b c

14,5,61 [4,6,71

e

Binary Tree Search to Obtain

the Set of Maximal Conlpatibles

Fig 1.2
extra sels (3,4,7) and (4,6,7) which are included in tile maximal sets (b) and

(c) respectively, are eliminated. The remaining five sets are 1he maxilnal

compatibles.

A high level description of the above procedure is given below:

procedure MAXCOMP(sset,i)

Assuml)tions:

- pivot candidates ,are numbered from 1 to n.

- initially sset consists of all pivots in the matrix and

i is |he first, pivot.

while i<n do

begin

(*split sset into left and right sets*)

lset = sset- imptbl[i]

rset = sset- [i]

if {ls(,t not a coml)atible set) then

maxcomp(lset,i+l)

if(rset not a compatible set) then

maxcomp(rset,i+l)

end

In tile above procedure, many branches do no! need to be continued to tile

completion of ILL(,search, since they are included in other subtrees. Moreover, as will be descril)ed later, we only need Io produce conlpatible sets of

maximunl size. Thus. there art, many branches in this tree that couhl be

trimmed to limi! the amount of search. Even including these features, thi.,€

algorithm has exl)onential complexity, and only serves to obtain inf_,rmatiml

about sparse matri('e,,.

To slu(lv lhe isstles (li._cusse(! earlier, a ]'.\.";('AI, program was writ lea to

perform symbolic L/l: decoml)osition on a sparse matrix. Our objecli\e was

to study the effects of parallel pivoting so the program performs the de('ompost!ion lo tile last parallel step and does not continue if parallel pi\€_! candida!es are not available. The structure of the program is outlined below:

p rogra m I'I\'OTS1'71"

- l_ead in input ma!rix and constrv, ct matrix structure.

-('onstrucI all maximal compatibles.

-if parallel pivo!ing ix not I)ossil)le go to stop

-Pick a set of c'Oml)atib]e pivots to be processed

in paralh, l.

-Permute the matrix according 4o the parallel pivots for this step.

-reduce the matrix and insert the resullant fill-ins.

-Repeat.

-Stop.
An alysis Performed

In general, in matrices arising from circuits there, are man)' different sets

of compatible pivots of equal maximum size. Depending on how a set is

chosen to reduce the matrix at, each step, we obtain a different behavior in

generation of fill-in elements, and as a result, different possibilities for continuingparal]e] pivoting in the next steps. The issues of generation of fill-ins

and parallelism in pivoting have been studied. We used different strategies

• to select a set of COml)atible pivots and then obtained statistical information

from some circuil malrices generated from the SPICE circuit simulalion program.

The Markowitz criterion [11] is well known for minimizing the generation

of fill-ins in sparse matrices in sequentihl pr()gramming. It is based on the

fact that at sic 1) k, the maximum numl)er of fill-ins generated by ch(_osing at]

as pivol is (ri- l)(c;- 1) . ltere ri- 1 is the. number of nonzero elements oth(,r

than ai/ in lhe i-[h row of the reduced matrix, and c3.-1 is the number of

nonzero elemenls otht, r than a0. in column j of the reduced matrix. Markowitz selects: as pivot element at step k, the element which minimizes

(ri-1)(c/-1). Thel)roduct(r i-])(c:-l)d ois the Markowitz number of eh,ment

a... In what follows, we use the Markowllz idea as a basis for tho seh,ction of 13 •

a COml)atible pivol set.

In our first analysis we coral)are two different strategies for choosing a set i

, of compalil)le l)i\'()ts among all maximal coml)atil)les. In both (,ase_ we (.onl sider ()nlv the st,ls (_f nlaxinlunl size. The first strategy (called :klarl_ox\ilz

i sum) ('hooses lh:ll s(,I :among all sets of maxinlum size in which the sum (,f th,.

Mark()witz numbers ,,f all its elements is minimum. The prol)lem here is lh_lt

some of the l)iV_)Isi_l the set chosen for re(lu_.ing the matrix may g(,n_,r:ll(,

fill-ins in lhe same l),_siiions, and thus we overt,slim;ate the Markowitz count

for a l)urely sequenii:ll case. As an altcrnali\'e, a second strategy is employed

(called Ored Markowitz). ]tere, using the boolean matrix]3 corresl)ondingto

lhe sl)arse matrix under consideratiou, we ('()unt number of nonzeros in a vector that is the result of ORing rows of pivot candidates in the set and mullil)ly Ibis number by lhe number of non-zeros in a vector resulting from ORing

columns of the l)olenl ial pivots.

Coml)arison of the above strategies on our test cases shows thai the fir.,,t

method is almost ,_lways superior. Our results show that, in gen(,ral, i)v

minimizing the _larkowilz sum we always get fewer fill-ins gener:lle(! and

often more rows are reduced in parallel steps. This study has shown that the

. amount ,:)fparallelism in circuit matrices is quite high but that the generation

of fill-in terms is also quite high in most. cases when compared to the sequenlial runs on the same matrices. The number of potential pivots to be pro-

• cessed in parallel al each step seems to 1)eso high that we could l)rocess fewer

pivots in parallel in a step without limiting the parallel work considerably.

An experiment to study this l)ossibility is performed by picking the maximum

sized set with nlinimum Markowitz sum as was exl)lained above. This set is

then used to reduce the matrix, with the following analysis performed on the

set, of eompatil)le pivots.
-Discard the pivot with maximum Markowitz count and determine

number of fill-ins that would be generated as a result.

-Repeat the above procedure until no more pivots can be discarded

from tile set, either because tile set size is too small or because

all Markowitz sums are zero.

Although the above analysis of reducing tile size of the set of COml)atil)le

pivots was done for each step, tile actual elimination and fill for a slc 1) was

done using tile maximal compatible with lowest Markowitz sum. This

analysis is repealed at each parallel step and the results show that it is possible to decrease tile generation of fill-ins at this step significantly by reducing

[ the amount of parallel work slightly. In fact, discarding only onecoml)atible

j pivot results in a decrease of at least abo'ut one third in the number of fill-ins

that would be generated otherwise.

\\'e performed this analysis over all generated sets of compatible pivots

also. In this exl)eriment, we chose maximum sized set, with minimum ,Markowitz sum and used it for re(Io(.ing the matrix as described I)elow:

- for nil sets of maximum size do

find the i)ivol with maximum Markowitz count and

remove fr(,m the sol.

find Ill(, set ()f maximum size and minimum Nlarkowitz sum

and determine number (>ffll-ins that would be generated

from the prr)cessing <)fthis set.

- repeat the above l)rt)c't'ss.

Similar results were obtaine(i by applying the nbove two procedur(,s to ,mr

lest matrices. Ilence. we will use the first method for the next l)hase. "l'h:lt is.

tile next analysis is l)erformed on tile set of maximunl size and minimum Nlarkowitz sum.

Even though tile above experiment shows we can always generate fewer

fill-ins at a step I)v avoiding the maximum 1)ossible parallt, lism, it (lees not

indicate that this will not delay lhegeneration of fill-ins to later steps. In our

next exl)eriment, we choose the maximum sized set with minimum ._Iarkowitz

sum. l)ut this time we (lis(';_rd the pivot with maximum Markowitz cot,nt

from the set and use the rt,sultin_ set for elimination and fill generation. \Ve

will nlso repeat the pre\'ic)us analysis by reducing the set, size and determining

numl)er of resulting fill-ins. TI, is work confirms our previous result that by

discarding some of the paralh,I pivot candidates according to their high hlarkt)wilz count we decrease the total generation of fill-ins.

Results of Complete Analyses

A set. of circuits to be simulated by the SPIC,E circuit simulati()n l)rc)-

gram is available as a benchmark to test SPICE. \Ve used these circuits as

input, to SPICE and generated their corresponding matrices. These matrices

are used as test cases for analysis purposes. The frst circuit is a siml)le

differential pair and generates a 16 by 16 m:_trix with 57 nonzeros. Tile

matrices are of small sizes and the size range is between 12 by 12 to 24 by 24.

The complexity of our algorithm to generate all possible maximal sets of compatible pivots would not allow us to test larger matrices, but the generated
information produces valuable statistics about parallelism and circuit

matrices. An algorithm with t+,lerablecoml)lexity to produce a set of compatible pivots will irrvolve heuristics; therefore, it will not give total information about the matrix.

The results of comparison of Markowitz sum and Ored Markowilz strategies are summarized in Table 1.1 {tables are provided in appendix A at. lhe

end of this paper}. The first column describes the circuit, the order of the

matrix, and number of nonzeros. Tire second column indicates the l)aralle]

pivoting step. Columns 3 to 5 correspond to the Markowitz sum strategy

described earlier, and columns 6 through 8 correspond to ()red Markowitz.

The first column for each algorithm is the size of the maximum set of pivots

obtained at a step, tire second column is the minimum operation count

obtained for such a set, and the last column specifies the number of fill-ins

that are generated as the result, of processing the indicated set. ('.t_!umn g

indicates the total number of maximal comp:_iibles generated at each step.

The last two columns are information generated l)y the SI'I(..:E l)rogram

about the amoun! of fill-in generated and the l)er('eniage of the matrix whi(.h

is zero.

As can i)e seen from the table, in e`"ery case the second stra!egy r(,sulle<l

in equal or more fill-ins and equal or fewer l):,rall('l slel)s with fewer number

of r(_ws reduced. This indicates that the _l:lrkowitz sum is a I)eit(,r h(,urislic

for selecling the se! of pivots among mnn.v sets. This ('an be observed from

the 16 I)v 16 matrix of the differential pair cir('uil. In the firs! step, with sets

of size six, Markov,'ilz sum generated 6 fill-ins while Ored Markowilz generated 8. The pivot set chosen by the Markowitz sum generated fewer fill-ins

than the Ored _l:lrko',vitz algorithm, and, as can be seen, the Ored _larkowilz

resulted in twice as many fill-ins as the Marko`"'itz sum arrd fewer pivots `'vere

processed in parallel (1-1 for Ored Markowitz and 15 for M'arkov,'ilz sum}.

The same behavior resulted from the EC,I., COml_atible SCMITT trigger circuit

which produced an 18by 18 matrix. The nulnher of fill-ins at s!ep2_,fpnrallel triangulation is 10 for Ored _1arkowitz and onl.v 4 for the _larkoxvitz .,,urn

with none being generated in the next steps. Ored Marko`'vitz generated -!

more fill-ins at step 3 a11d ,,','as not able to find an.v more parallel pivo! c:_ndidates, but the first strategy continued to do one more parallel step. Of

course, there are cases where both stra!egies produced similar or close resul!s.

as can be seen from lhe table. The table also indicates thai, in 1)arallel runs.

generation of fill-ins is much higher than in sequential runs of the SI"I('E pro-

• gram. At the same time it can be seen that the matrices generally do not

become dense rapidly, and parallel pivot candid:ales are available to almos!

the very last steps of the triangulation process.

The result of our next analysis is shown in Table 1.2. At, each step, a set

of maximum size and mininlum Markowitz sum is selected to reduce the

matrix. Fu_lhermore, from this set we repeatedly remove a pivot with maximum Markowitz count and compute the number of fill-ins that w_mld be

generated if this set were used toreduce the matrix. As can be seen from the

table, in every case it is possible to reduce number of fill-ins significantly 1)y

reducing the amoun! of parallelism slightly. For example, for the 16 bv 16
10

matrix with 57 nonzeros, we can see thai if we reduce the number of t.€)mpalibit pivots from 6 to -1 by removing the two pivots with highest M;irkowilz

count from the set, we can l)revent generation of more fill-ins. Also ill tilt'

last 24 1)3"2,-1matrix with 158 nonzeros, we can reduce the number of generated fill-ins by a factor of 2 (from 40 to 20), if we discard two pivots in step

one in the same fashion. This is a general result that can be observed from

the table for all cases and all parallel steps.

In tile next exl)eriment we confirm that it is l)(,ssible to reduce the total

generation of fill-ins, as opposed to just. at each step , I)2' using fewer than lhe

maximum aural)or of compatible pivots. In every (.as(, we have been al)le to

reduce the Iolal number of fill-ins by sonic fra(.lion (at ]east about ono third).

COml)are(I to the case where nlaximum l)aralh,]ism was eml)loyed . "]'heso

results are summarized in Table 1.3. ll(,re we chose to discard a l)i\'(_t fr(ml

the maximal ('t)ml)alibleset according to its highest Markowitz count. Ifih(,

maximal se! would not generate any fill-ins, because of a zero Markowitz sum,

we did not discard any pivots from the set. The Iota] number of fill-ins g(,n.

eraled for the first matrix (16 !)3' 16) is 2 which is one third of the amoun!

gen(,rate(l wilh our first experiment (6). This number was reduced fr(,m .10 t_,

26 for the case of the 24 l)2" 2-1 matrix with 151 nonzeros. In this (.:lse the

number or parallel st(,l)S was increased fmln 5 to 6, but the total number t,f

rt)ws thal could be rt,(lu(.ed in t]lt,S(,Slel)S rt,maiwt,d (.C)llStalll. In fact. in m,,st

cases, tile numl)er of parallel SlCl)S is ilL(.rt,ns(,(I. I)u! !he total number _,l"

l)ivots lhai could be l)ro('essed in lh(,se steps (]_('s not change mu_'h (no

change is greater than one addition or r(,(lu(.lic)n in the number of rc,du(.(,d

rows ).

Generation of Compatible Sets from the Incompatible Table

1_ is clear that in large sparse circuit matrices the number of l)ossib]e

pivots to be processe(! at each step will 1)emuch higher than our smal I exam1)lematrices, and lherefore, it will be possible to obtain enough paralh, I work

by jusl considering a sub-maximal set of compatible pivots at each step. The

algorithm described involves a co)replete binary tree search and has exl)onentia] complexity in the order, n. of the sl)arse matrix. In order to come up

with a good heuristic, we need to relax tILe requirement of finding the m,'lximalset of compatible pivots with minimum Markowitzsum. Asa conclusion

from the above analysis, we will have to reduce the size of tile set to decrease

the generation of fill-ins. Keel)tag these problems in mind, an acceptable set

would be one which has a large number of pivot candidates for parallel processinl_ and a low enough Markowitz sum. We now need to look for a procedure which tends to produce a number of compatible sets of reasonably

large size and low Markowitz sum. Having generated such sets, we can then

choose the best candidate among these coml)atible sets using the same criteria as before. In what follows, we will describe different issues which will

lead us to a good heuristic algorithm and a se! of parameters lo be used in

trading off between fill-in generation and the size of the set, of parallel pivo!

candidales.
11

So far, the information from the inconil)atible table has been used to

construct the maximal compatil)le sets of pivots in a complete binary tree

search algoritllm. A more careful analysis of tile incompatible tabh_ could

provide a set of compatible pivots without the need for searching the tree.

" As we know, this table gives information al)out tile incompatil)le pairs of

pivots. In other words, by looking at column i of tile table corresponding to

pivot i, we obtain all pivot numi)ers j>i where pivot j is incompatible with

pivot i, for a given ordering of the matrix. Note that, we are assuming pivots

are taken from the diagonal of the matrix and they are numbered 1 through n

corresponding to roars 1 through n of the matrix. Consequently, if column i

of the table is null, then the corresponding pivot number i is eoml)alible with

every pivot whose corresponding column 'lies to the right of column i. llence,

by scanning the incompatible table, we Call find a set of compatible pivots

whose corresponding columns in the table art, null. Clearly, pivots with such

a property are compatil)le and can begrouped ill a COml)alibleset. 17singlh,

representation of the incompatible table described earlier, the al,_x'e procedure can be formulated as:

scan implbl from right to left

for each colunin i of imptbl do

if( inlplbl; is empty) then

(*add the <.orresponding pivot to the set of compatibles*)

compset = compse! + It]

where compset is the set of compatible pivots whose corresponding c_,lulnnS

in the table are null. Now if there exists a pivot k such that the set of pivots

incoulpatible with it in column k of tile tabh,, is disjoint from tile set of

already constructed €.olnpatible pivots in compse! , then k is coral)at iblt, wit h

every pivot in conlpse! . Therefore, we can expand conlpsel I)3"adding k t_ it.

The above proced tire Call now be complet ely tit,scribed as:

scan imptbl from right to left,

for each column i of implbl do

begin

if ( imptbl i 1'7 compsel it: erupt.v) then

(*add [i] to the set of COml_atibles_)

compset = compset + [i]

else

delete row i of imptbl

end

The compat.ible set, compset , produced by this procedure, will be referred to

as an ordered compatible set from now on, since it is obtained by inapo_ing a

" specific ordering on the diagonal elements of tile matrix to gel, the incompatible table. As an example, the incompatible table of matrixA2 in Fig. 2.1.a is

given in Fig. 2.1.b. The compatible set, corresponding to tlle null columns of

the table consists of pivots 10 and 11. This set consists of 2, 10, and 11 after

the above expansion.

As was explained previously, our strategy for selecting a compatible set,

among all possible compatible sets of equal maximum size was to select the
12

1 2 3 4 5 6 7 8 9 10 11

1 x

2 x

3 x x x "

4 x x

5 x x x .

6 x x

7 x x x x

8 x x x x

9 x x x x

10 x x x x x

11 x x x x x

Matrix A2

Fig. 2.1.a

o

3

4

5 x

6 x x

X

x x

9 x x x

10 x x x x x [

I

11 x x x x x x I ]

I

1 2 3 -1 .5 6 7 8 9 10

Incompatible Table

Null columns: (10,11)

Compset: (2,10,11)

Fig. 2.1.b

one with minimum ._larkowitz sum. That is, to select, the set in which tile

sum of Markowitz numbers of the pivots in its set is mil_imum. If we consider

the set of compatible pivots constructed above directly from the incompatible

table, we see that it consists of pivots 2, 10, and 11, which in turn have Markowitz numbers 0. 4. and 12. In general, we would like to have a compatible

set consisting of pivots with as lowMarkowitz numbers as possible. It is also

clear that pivots with low Markowitz numbers generally have fewer incompatibilities. Moreover. by looking at the incompatible table of Fig. 2.1.b, wesee

the compatible pivots 10, 11 are obtained from the right end portion of the
13

table. This is usually tile case. sineeasweconslrue! columns of Illt, inc(mlpatible fable, we arc. left with a smaller submalri× to work with. l"hu_,, afler

Cmnl)leting each column, we have fewer incompalibles left for the e(m_truelion of the next column. These observations lead us to use a differenl order-

• ing in which the first column of the incompatible table has the maximum

number of incoml)atibles and as we work our way to the right end of the

• table, the number ofincoml)atibles will decrease io the minimum. Such an

ordering implies the resulting incompatible {able will have more null c_:lumns

clustered at the right end. So the ordered compatible set that ca1: be constructed from the ordered fable will be of a larger size and smaller ,Nlarkowilz

sum than the results of lhe above procedure. As a result of these argumenls.

we sorl the pivots in order of decreasing Markowitz numbers. I,rsing this new

ordering, we can conslru(.t a new in('oml)alible table with the firsl (.olumn

corresponding to the pivot xviih highest Markowitz number and the last

column corresponding to the pivot with lowes! Mnrkowitz numl)er. As an

example, the Mnrkowitz numbers and the new ordering of the pivols arc.

shown in Fib. 2.2.a for matrix A2 of Fig. 2.1.a. The corresponding ordered

incompatible table is given in Fig. 2.2.b. It can beseen from Fig. 2.2.b tha!

the collection of pivots corresponding Io null columns of the fable gives a

compatible set of size -1and Markowitz sum .I consisting of pivots 1, 2.3. and

4. This is in comparison with set of size 2 and .klarkowi!z sum 16 gent'raled

: from lhe unordered :,,icompatil)ie table of Fig. 2.1.1). After exl)andin_ lhi.,,

set. we pr,dnce a compatible set of size .5:and ._larkowilz sum 16 con,_isting (,f

pivols I, 2, 3, 4. and 9.

Limited Binary Search Tree

In this section, we will combine the idea of an ordered compatil)le set

with the tree search algorithm described earlier to obtain a limited tree

Markowit z

Pivo! Order Nu tuber

1 0 9

2 0 11

3 2 8

4 2 6

- 5 2 1o

6 4 7

7 3 3

8 9 4

9 12 5

10 4 1

11 12 2

Fig. 2.2.a
14

I1

8

6

10 x x x

7 x x

3 x x

4 x x x

,5 x x

I x . x ]

2 x I

9 11 8 6 10 7 3 4 5 1

Ordered lneompatil)ie "]'able

Null eolumns: (1.2.3,.1)

(_ompset , . (1,_,,_3,.t ..q)

Fig. 2.2.1)

seareh al_rithnl which l)r¢,duees ,qn ae('eptab]e se! of <'Oml)alible l)i\'_Is for

reducing !he' malrix. Given a _el of all pivo! elemen!s, we can n_,w direell.v

1)r_(lueeaset ofeompatii)lepivols from lheordered ineompalil)letable. This

ordered eompalil)le set is obtained for the initial s!arling set. at ihe ro_l _f

lhe binary search free. A child set in the tree is a subset of il._ parent sel. In

thi: eonlexl, every set at any" given point in lhe tree has fewer pivols lhan

the reel set. Sueh a set couhl be considered as a slarling set itself. Provided

we could produee lhe eorreel ineompatil)le lable for ihis sel, we eould generate ils corresponding ordered compatible set direetly from the new table.

The ineompalibh, lal)le for a given starling sel. Si, is the original lal)le

wilh those r_ws an(! ('€_lumlls (._rresl)on(ling lo the pivols abs,..nl from 5;'.!

e]iminaled. If we le! 5' be the iniliai set of all pivot ('and!dales in !he sparse

malrix an(I S. be an arbitrary slartingset, in lhe tree then the l)roeedurelo

ob!nin lhe ordered eoml)alible set for S; , eompset i , from an ul)dale(l and

ordered in(-_,ml)alil)le fable can be rel)resented as:
15

1. compsel i = empty

2. le_s = 5"-S.

I

3. for j = n down to 1 do

" 4. begin

5. if(jc S;)then

6. begin

7. tempset = implbL- less

8. lempset tempsel N compset i

9. if ( tempset = empty ) then

lO. compsel i = compsel i + [j]

II. end

12. end

where less is the set of pivots absent from Si. Line .5 allows only those

columns of the incompatible table whose corresponding pivot j is in ._.lol)e !

lc.-,ted for lhe ('ompatibility relation. Set less is used in line 7 to eliminale

rows corresponding to the absent pivots in S.. compset, holds the (.urront set ! !

of cc)mpatible pivots. A check for a new pivot being compatil)h, wilh lhose

already in compset i is made in line 9.

It is now possil)h, to l)rodu(.e an ordered compatible set for any set at any

arbitrary point in the tree dire(.tlv fr()m the incompatible table. (;iven a

starting set, our method of producing an or(lered compatible set tends to gent,rate a large set of low .X.larkowitz sum. Thus, we can produce a number of

ordt, re(! ('ompatil)le sets for many starting sets at different points in the tree

and choose the best candidate among them 10 reduce the matrix. Thefollov,'-

ing theorem will eliminate of some of the redundant work.

Th e 0 r enl

.-tll ordered compatible sets derived from the starting sets in the binary

search tree u'ith lerei L-1 or less are included in the ordered compatible sets generated from the sels at lerel L of the tree. (i.e.. it is only necessary to generate

ordered compatible sets for starting sets at level L to cover those at level I(L.)

Proof

Let 5" be the initial starting set at the root of the binary tree consisting

of all pivots, Pl .... ,P . Let So" $1 be the left and right children of 5` . Let

• compset be the ordered coml)atible set obtained directly from the incompatible table for the set S . Similarly, let compset oand compset I be the ordered

compatible sets corresponding to SOand S 1 respectively.

A pivot Py can split a set S iff:

P.€ S and

{s_t of incompatibles with P. } f') S ;e empty.
1B

Assume P. splits S int_; SOand S l ; then"

So = S - {set of incompatibles with Pj } and

s,=s- 1 .1.

There are two cases to consider:

i. P. not in compsel

The table corresponding to S 1 consists of the same null columns and

compatible pivots as in compset so:

compsel = compset I.

ii. P. €compse!

"l'tlen we must have:

imptbIp 1'7compsel = empty

1 .

sinceP, is compatlblewith all pivots in compset . In this case, compsel o

obtaineJd from So is equal to compsel . \Ve know P. is in the set So and

that the incoml)atible table for So is the same as thg table for the l)arent

set S with those rows and columns corresl)onding to incompatil)]es of Pi

eliminated. Thus, all the compatible information which resulted in production of compset is transferred from the l)arent set S to So and consequently:

compsel = compsel o.

The above argument proves that, at level l, one of the sets S o or ,5"! will

produce the same ordered compatible set, asproduced by its parent set. This

proof holds for an)" two children of a set. In other words, at any point in the

tree, an ordered compatible set corresponding to a parent set is reproduced

by one of its children.

Induction on level verifies that generating the ordered compatible sets

for every set from the root. through level L of the tree does not produce any

more information than producing the ordered ('Oml)atible sets for every set. at

level L only.

As a consequence of the theorem, we generate all the sets at a given level

in the binary tree, and for each set, we produce an ordered compatible set.

from the ordered incompatible table. Among the generated compatible sets

we choose the set. of largest size and lowest Markowitz sum to reduce the

matrix and call the the resulting set the elimination set.

If we note that we split, each set at each level of the tree for a given pivot,

according to its incompatibility information, then generation of the starting

sets at different levels could be done in various ways:

i. We could split the starting sets using the original pivot ordering given by

the input sparse matrix. This would generate completely random

results.
17

it. Tile same ordering used to order the incoml)atibI e table could b(. used to

split tile sets. This left to right ordering does not seem to agr(,e with our

low Markowitz sum requirement. At, each split {level in the tree}, we

include one of the pivots, say Pi ' with highest number of incoml)atibh,s

• {highest Markowitz number) m the left subtree. This inclusion also

means we take a large number of pivots incoml)atible with Pi out of the

sets in the left subtree. These l)ivots that are incompatible with Pi have

lov,'er !\larkov,'itz numbers than Pi and could themselves 1.)eCOml)atil)le

with some other elements in the set. As a result, this ordering will produce a left set corrsiderably smaller in size than the resulting right set.

Moreover, the left set contains pivots of high Markowitz number which

would produce man)" fills if used to'reduce the matrix. Therefore, some

of the large compatible sets with small Markowitz sums cannot be generated from one of the sets in the left, subtree unless we search very dt,ep

in the tree. In this case, the desired compatible sets would be in one of

the right subtrees.

iii. A third alternative would be to split the sets with pivots in increasing

order of their ;\larkowitz numbers. Of course, in this ease, the incompatibility information of the pivots used to split a starting set is taken from

the right end of the incompatibletab]e. Thus the completeincompatibi].

it)" information for a pivot i is obtained by concatenating the row and

column i of t he table. This l)rocess, seems togivea bett:er balance to the

binary tree for the first few levels used to generate the starting sets

required in our algorithm. Furth(,rmore, it has the proI)erty that does

not ignore pivots of low Nlarkowitz numbers.

The high level description of this algorithm is given below:

Program Parallel Pivot ing

calculate/\larkowitz numbers of pivots in the

remaining unreduced matrix.

- SORT pivots in decreasing order of Markowitz numbers

- produce all startillg sets at level I.q_EVEL taking the

pivots to split the sets from the root to tTI_EVEL in

order of increasing *larkowitz numbers.

- for each set, at U'I_EVEL pr(,ciuce an ordered compatible set, from

the ul)dated ordered incoml)atible table.

- among the ordered compatible sets generated above choose the

maximum sized set with minimuna Markowitz sum (Elimination set,).

Here, ULEVEL is a preset, level number indicating the depth of the tree to be

searched. The algorithm is no longer exl)onential in time. An efficient imple-

• mentation of the required sort and set operations are important factors in

efficient execution of the algorithm. The set ol)erations used in tlreconstruction of the incompatible table are of order 1 (adding an element to the set. or

a test for membership). The incompatible table can therefore be constructed

in time nz, where nz is the number of nonzero elements of the matrix. Ceneration of an ordered compatible from tire incompatible table requires scanning n sets corresponding to the columns of the table, and performing intersection and difference operations on the sets. These operations are of order n
18

with a conslant factor equal to the inverse of tile numl)er of bits l)er ('ompulerword. The set ol)eralions are usually implemented in machinel:lnguage

or micro cod(, and thus have a small time factor. They could be considered to

have a constant time (rather than order of n) coml)ared to tile time taken to

execute a high level language statement. Production of all starting sets for

level ULEVEL takes a constant time. (-;eneration of an ordered compatible

for each starling set at, ULEVEL takes a constant times n as explained above.

For reasonable values of ULEVEL, all ordered coml)atible sets can be derived

in parallel for different starting sets. In the next section we will see that good

results are obtained for small, constant values of ULEVEL compared to n.

The complexity of the algorithm is bounded al)ove by the sorting; algorithm.

Thus, employing an efficient parallel sort would improve the performance of

the new algorithm.

Balance between Parallelism and Fill-in Generation

Even though the above procedure tends to 1)roduce large sets of low Markowitz sums, we still could optimize the generation of fill-ins by considering a

subset of the elimination set. Thai is, there could still be some room fortrading off I)etween parallelism and fill generation. To accomplish this task, we

need to come up with parameters to, control the number of pivots to be pro-

('essed in parallel and tile number of fills to be generated. One such l)nrameter could be the size of the set ofcompatil)iepiw_ts. By allowing a percentage

of the set Io be discarded, we can control the the number of COml)atible

pivots to a degree that does not limit our parallel work by too much. For

clarity, this parameter is called the shrinkage parameter and is used as a

lower limit to shrink the elimination set by a percentage of its size. A

different parameter could be an upper limit on the size of the elimination set.

This limit would allow just enou_ih work to keep our parallel processes busy.

Of course shrinking of the elimination set must not be done arbitrarily by

throwing pivots out oftheset. In general, wewould like to shrink our set by

discarding pivotstha! would cause generation of man)" fills. Such pivotstend

to have high Markowilz numbers. We already have pivots ordered according

to their Markowitz numbers. Wecould use this ordering to scan pivots with

highest Markowitz number in the elimination set and test against a threshold

value. If pivots with Markowitz numl)ers greater than a threshold exist and

if our shrinkage parameter allows, they are discarded from the elimination

set. Use of a threshold value will allowusnot toshrink aset, lhal consists of

all good pivot candidates of reasonably low Markowitz numbers. I"o serve

this purpose, the threshold value should be set in comparison with low and

high Markowitz numbers of the pivoting elements in the matrix. Again the

ordered Markowitz numbers of pivots can be used to set such a threshold

value conveniently. One way is to specify a fraction of candidates to be discarded from the elimination set. Consequently, we set the threshold to the

Markowitz number of a pivot, in a specific position in the list of pivot elements of the unreduced matrix (ordered by decreasing Markowitz numbers).

Any pivot above this point in the ordered list is considered to have a high

Markowitz number and therefore is a candidate for being discarded from the

set, and any pivot below this point is considered acceptable.
19

Pivots in tile elimination set are scanned in o,"der of their highest Markowitz number. If a pivot with ,\larkowitz nlJml_er greater than Ill(, threshold exists and if the set is not already of minimum size, it is discarded from

• the set.. Tile process is repeated until either no more pivots of large Markowitznumbersareleft, in the set or the set, cannot be further shrunk. In the

next section we present the result of different strategies and various parame-

" t.ers discussed here for a number of test matrices.

Analysis of the Results

The complexity of tile binary tree search algorithm to obtain maximal

COml)atible sets was such that it could not be run to completion for a 38 by

38 matrix. To verify the validity of otir heuristic program, we performed

every analysis described in this section on the small test matrices of Table

1.1. Recall that the new algorithm produces a number of starting sets for a

given level [UIA';VEL) of the binary search tree. For each starting set, an

ordered coml)aiil)le set is produced. Among the generated ordered compatible sets, the set with maximunl size and minimum Markowitz sum is seh:cted

as the eliminatiol, set at that parallel step. Two alternative ordering;s for

generation of siarting sets at ULEVEL were discussed earlier. For simplicity,

we call the algorithm to reduce a sp,arse matrix 17y compatible pivots using

the decreasing order of Markowitz numl)ers for starting set splitting,

DCOMP. Similarly, the algorithm which uses the increasing order of Markowitz numbers is called ICOMP.

Detailed information produced by I)COMP and ICOMP are presented for

three sparse matrices in Table 2.1. Column one of the table gives a description of tile sparse matrix undcrconsideration. Cohlmn 2, specifies the parallel step. Columns 3, 4, and 5 give the number of compatible pivots in the

elimination set, its ]_larkowitz sum and number of fill-ins generated at each

step for program DCOMP. Similar information is summarized in the next

three columns for program ICOMP. The information presented here is for

ULEVEL=4. The first two matrices have been completely analyzed in the

previous section and are presented here to show the validity of our proposed

algorithms. It is interesting to see that, for the first matrix, DCOMP produced exactly the same results as the complete tree search program. On the

other hand, ICOMP produced different results. Even though ICOMP produces a smaller conll)atible set in the first step. it finds larger sets in the next

steps and reduces the same number of rows (i.e., 21) in five parallel steps.

ICOMP generates 22 fills, almost half tile number produced by DCOMP (40)

or even the complete binary tree search algorithm (40). The same behavior is

observed from the second 24 by 24 matrix. The third matrix is obtained from

the circuit of an 8-bit full adder and is a 144 by 144 matrix with 616

nonzeros. Note that both a]gorithmsproduced an elimination set of 72 pivots

in the first, step, and so, half of the matrix can be reduced in parallel in one

step. In this case the advantage of ICOMP over DCOMP is not significant.

To see how variation of depth will affect the resulting compatible sets,

we ran both programs for values of ULI'_VEL between 2 and 5, for a number

of matrices. These results are summarized in Table 2.2. Again the first

column describes the matrix. Thc second colunan specifies ULEVEI,.
20

Colunlns 3 to 6 give related infornlation for tl,e I)(:OMP program, ilere the

first and the third columns specify number of l)arallel sleps taken Io reduce

the matrix and number of rows reduced in those steps; respectively. The

second column is tile average parallel work at each step and is obtained l)y

dividing total numl)er of rows reduced in parallel by the number of parallel

steps. The fourth column gives the total number of fill-ins generated by

parallel reduction. The next four columns of the table provide similar information for the IC()MP program. The last. two matrices of the fable are produced from the SPAR l)rogram, which is a structural analysis program [12].

These twomalriceshavea pecu]iarblock stru('lure. Our initial objective was

tostudy sparse m::Iri(.es arising from SI'ICE. "]'hese matrices ordinarily have

a random sl)arsity structure, but at the'same time, the limited connectivity

between nodes of the input, circuit results in a limited number of nonzeros per

row/column. The SPAR matrices will provide some insight into the behavior

of our heuristic algorithms for a wider (.lass of matrices.

It is clear from the table that, in almost, every case, IC:OMP produces

better results both in terms of number of rows reduced in parallel and

number of fill-ins generated. As was expected, D('O,kl]' finds elimination sets

of lower Markowitz sums as we search deeper in the tree. This is o])served

from the first 2-1 I)v 2.1 matrix and from the last SI'..\R generated 505 by 505

matrix. Ill the first matrix, ]COMI ) ])r_)(ll,ced 18 fills, reducing 21 rows in 5

parallel steps, while DCONIP generate(! more than twice the numl)er of fills

and reduced 20 rows ill 5steps. The number of fills decreases forDCOMPas

ULEVEI, is increased, while I¢_',OMP takes the opposite direction. This also

shows that reasonal)lv acceptable compatible sets, both in terms of size and

Markowitz sum. are generated for small values of ULEVEL and it is not

necessary to search very deep in the tree. The above observations hold for

every matrix presented in the table, except the 78 by 78 matrix produced by

the SPAR program. This matrix does not have characteristics typical of

SPICE generated matrices; but, as we will see in our next analysis, acceptable

results are produced for this matrix as well. Note that there are cases fortlle

DCO.MP program in which a higher average l)arallelism is indicated in the

table than for I('O._IP. In those situations, it is often tile case that fewer

rows have 1)een reduced by DCOhlP than by l(:Ohlt ).

The remaining analyses are performed on the ICOMP 1)rogram only,

since it produces better results. In what follows, a value of 4 is used for

ULEVI'_L. The next step is tostudy tile effects of varying the param(,ters proposed earlier to obtain a balance between generation of fill-ins and the

amount of parallel work. Results are summarized in Figures 2.3 to 2.6 for

four of the matrices of Table 2.'2. In these graphs, four different symbols are

used to represent four different values of the threshold parameter. Recall

that the threshold is set to the Markowitz number of a specific pivot in the

ordered list, of pivot candidates. On the gral)hs, the threshold value is given

as a fraction of the pivoting elements in the reniaining unreduced matrix,

ordered in order of decreasing Markowitz numl)ers. For example when the

threshold is 1/3. tile Markowitz number of the pivot residing in the 1/3 point

of the ordered list of pivot candidates in the unreduced matrix is obtained.

Any pivot in the elimination set with Markowitz number greater than this
21

0 1/10

• z: X 1,'3

o 4" 1,,2

,,,..1

° *"a.• 2/3

_,_- _ _ ®

_- "1" -'1"

! ! !

o _, , 30. 4O,

5HRI;_KAGEPARAMETER

Fill-ins vs. Shrinkage for l)iffcrent Thresholds.tlLE\:EL= 4

144 by 144 Matrix of 8-bit Full ,.\dder

Fig. 2.3

17. I ! It

0 _ _

SHRINKAGEPARAMETER

Fill-ins vs. Shrinkage for Different. ;l'hresholcls,ULEVEL= 4

24 by 24 Matrix of MOS Amplifier AG/DG

Fig. '2.4
22

, _ ® 0 0 0 Izle

X 1/3

+ 1/2

• 2/3

_. X X

_. + +

*

143. @

O 15. 311. 4O.

SHRINKAGEPARAMETER

Fill-insv.ShrinkageforDifferentThresholds,ULEVE1,=4

78 by 78 Matrix of SPAR Program

Fig. 2.5

55e_Lq

_ X

!_'77. •

i i i

SHRINI{AGEPARAMETER

Fill-ins vs. Shrinkage for Different, Thresholds,LlLEVEL= 4

g05 by 505 Matrix of SPAR Program

Fig. 2.6
23

value is a candidate to be discarded from tile set.

The graphs present information about the number Of generaled fill-ins

versus the shrinkage parameter. In each case, the analysis is performed for

threshold values of 1/10, 1/3, 1/2, and 2/3. For ever)' value of the threshold,

• shrinkage parameter values of 0, 5, 30, and 40 percent are considered, llere

we have limited our sets to be at, most of size 25. Thus any set which consists

. of more than 25 compatible pivols is reduced to size 25 by discarding pivots

of highest Markowitz numbers. From tile information presented in the

graphs of Fig. 2.3-2.6, it is apparent that we have been able to reduce the

generation of fill-ins in ever)" case. For the 2.1 by 2.1 matrix (Fig. 2.-11, the

number of fills is reduced by 1.1_, over the range of threshold values considered. This number is higher for tlie rest of our test cases. \Ve have

obtained an overall reduction of 30._, 36_, and 50,%in the number of fill-ins

produced by ICO.MP for the 8-bit full adder matrix and the two SPAR.

matrices, respectively. It is important to note that the number of parallel

steps laken to reduce the malrices and the number of rows reduced in those

steps did nol change considerably with a change in the above parameters. In

each case, tht, change was not more than one in the number of steps or lwo in

the numher of n, duced rows. Thus we have b_,_,nabh, to reduce the, _eneration of fill-ins considerably l)y givinB,up an insi_llificanl amount of the existing parallel capabiiily. In other words, by eml)loyin _ lhe above paramolers, a

better balance belween number of COml_alible l_i\'ols generated at dilh,reni

steps is achi_'ved.

A characteristic of these a!_plication sparse malrices is that many pivots

have equal Mark_witz numbers. When pivots are ordered, those wilh equal

_larkowilz numbers are clustered together. Therefore, when the value of the

shrinkage l)aralllt'ler is small, the threshold has no effect. This can be seen

for <lifferen! values of the threshold parameter for the shrinkage l>arameler

e<l+:alto 5_'_.ill :_l_.x"of the gral>hs of Fig. 2.3 lhrc>ugh Fig. 2.6. As this reduction percentage is increased, tilt, threshohl parzlmeler plays a more eIl't,ctive

role. In tilt, S-hi! full adder malrix, shrinking lhe set size by 5_ accounls for

m_sl of l he re<lucli_n in till-il_:., and afler that, the changes are not

significant. ]:t,r lhe lwo matrices generated I)v tilt, SPAR program, results

are more evenly (lislributed over tile changes in the above parameters (Figures 2.5 and 2.6). Different values of these l)arameters do not start It) affect

the resulls until we allow a large fraction of pivots to be discarded in the 2-1

by 24 matrix of Fig. 2.4.

• Recall tha! \re have restricted the compatil)le sets to be at most of size

25. For the matrices thai generate considerably larger sets, the threshold

parameter does not play an imt)ortant role i)ecause an upper limit is imposed

" on the size of tile elimination set. To observe the effects of tile threshohl

parameter, the upper limit is ignored. Resulls for the same values of the

threshold and shrinkage parameter for the 8-bit full adder and the 505 by .505

SPAR matrix are presented in Fig. 2.7 and Fig. 2.8. The other two matrices

do not generate as large sets and, therefore, are not aff_,cted by the removal of

the upper limit parameter. A threshold value of 1/10 does not cause any

reduction in fill-ins in either case. This suggests that the generated
24

:_. 0 0 0

!

0 1/10

X 1/3 •

"F 1/2

_-_,. _ • 2/3

_0" X ×

-_. + +

2tO. I

I I !

0 S. ,30. 40.

SHRINKAGE PARAMETER

Fill-ins vs. Shrinkage for Different ThreshoI(Is,ULEVI_,L = 4

no upper Limit, 144 by 144 Matrix of 8-bit Full Adder

Fig. 2.7

_ _ 0 0

I

I'4

h

,__.. X X

q427. +

4172. •

3gtT. •

3._e I I !

D 5. 30" 40.

BHRINKAGE PARAMETER

Fill-ins vs. Shrinkage for Different Thresholds.l.rLEVEL= 4

no tipper Limit, 505 by ,505 Matrix of SPAR Program

Fig. 2.8
25

COml)atible sets consist of pivoting elements of lower Markowitz numbers

than the highest one tenth. These graphs show more variation with the

threshold parameter, even t.hou_oh a higher number of fill-ins is pr¢_duced.

Without exception, the fewest fill-ins are produced for tile largest, values of

the threshold and the shrinkage parameter.

The number of fill-ins produced by ]('OMP compares reasonably with

sequential runs on the same matrices. For exan_l_le, the sequential run on the

" 8-bit full adder matrix produced 166 fills, and ICOMI' produced 196, which is

an increase of about 18%. In general, as expected, the number of fill-ins produced by ICOMP is higher than the sequential results, but the difference is

not. great.

Coneluslon

Solution of sparse systems of equations is essential in many application

programs. Often such a system has to l)e solved repeatedly. In this paper we

veritied that in sparse matrices arising from eh, ctronic circuits it, is possible to

do computations on many diagonal elements simultaneously. A <'(m_l)lete

analysis of some test matrices, done by generating all maximal COml)alible

sets of pivot elements, indicated the existence of many compatibh, l)ivols in

these matrices. \Ve have shown our l,esl matrices do iiot become full during

the (lecoml)osition. Furthermore, it was shown thai many parallel (.omputalion steps are possible, and during these steps, the matrix is often reduced

coml)letel.v. The (.ompeting issues of l>arallel pivoting and till-in generation

have i)een studied, and we verified through examples that it is possible to

reduce the l)roduction of fill-ins by removing some of the parallel pi\'oT e:indidates from the elimination set on the basis of high Markowitz numbers. A

heuristic algorilhm was then proposed to l)rodu('e large (.onapatible sels of low

hlarkowitz sums I)\" a combination of an ordered partial tree search strategy

and generation of ordered conlpatible st,ls. l)ifferenl orderings to produce lhe

orderc(l compatil)le sets were suggested, and their advantages and disadvantages were discussed and verified lhrough the simulated results. Anuml)erof

parameters to provide a balance l)etween generation of fill-ins and the

amount of p::ral]e] work were suggested, and their effects were determined in

the simulated results.

The incompatible table required by the algorithm can be constructed in

time nz (number of nonzero elements of the matrix). Production of starting

sets fora given ULEVEL takes a constant time. For ULEVEL small and (.onslant compared to n, generation of ordered coml>atibles from starting sets is

" of order n set intersection and difference ol>erations. Assuming efficient

implementation of the set operations is available, the heuristic algorithm has

- a complexity l)ounded above by the sorting algorithm required in the program. Thus, employing an efficient, parallelsort 1)rogram would iml)rovethe

total performance of the new algorithm. Nevertheless, our results show that

many compatible pivots are produced for parallel reduction of the sparse

matrices, and the process can be repeated until the matrix is ahnost completely reduced. In cases where the matrix is not completely reduced, the

remaining submatrix is of such a small size that parallel operations have little

effect. Significant reduction in generation of fill-ins is obtained by varying
26

t,he proposed l)arameters. Moreover, as the result of these parameters, a

better balance between the number of coml_atible pivots generated at

different steps was achieved, while the reduction in parallel work proved to

be insignificant,.
27

REFERENCES

[1] Duff, I. S., "MA28: A Set of Fortran Subroutines for Sparse Linear E(lualions," AERE Reporl R.8730, London (1977).

[2] Jordan, 11. F., "Parallelizing a Sparse Matrix Package," Iieporl CSDG

8I-1, Computer System Design Group, Eleclrical and Coml)uter

Engineering Department, University. of Colorado, Boulder (1981).

[3] Alaghl)and, G. and Jordan, H. F., "Parallelizing a Sparse Matrix Package," Report CSDG 83-3, Computer System Design Group, Electrical and

Coml)uter Engineering Departnlent, University of Colorado, Boulder

(June 198.3}.

[-1] Calahan, D. A., "Parallol Solution of Sparse Simultaneous Linear Equations," t'roc. 11th Annual Allerlon (:oaf. Circuits and System Theory, pp.

729-735 (Oct. 1973).

[5] iluang, J. \\'. and Wing, O., "Ol),timal Parallel Triangulation of a Sparse

.Matrix," IEEE Trans. on Circuits and Systems, v. CAS-26, no. 9, pp.

726-732 (Sept. 1976).

[(i] Wing, O. and lluang, .!. W., "A Computation Model of Parallel Soluiion

of Linear Equations," IEEE Trans. on Compnlcrs, v. C-29, no. 9, pp.

632-638 (July 1980).

[7] Jess, J. A. G. and Kees, It. G. lkl., "A Data Structure for Parallel I,/!"

Decomposition," IEEE Trans. on Computers, v. C.-31, no. 3, pp. 231-239

(March 1982).

[8] Peters, F.J., "Parallel Pivoting Algorithm for Sparse Symmetric

Matrices," Parallel Compuling, v.1, North-llolland, pp. 99-110 (19S4).

[9] Nagel, L. W., "SPICE2: A Computcr Program to Simulate Semiconductor

Circuits," Memorandum Et?L-31S:20, Electronics Research Laboratory,

University of C,alif., Bei'keley, C.A 9-1720 (1975).

[10] Lewin, I)., Logical Design of Switchin9 Circuits, American Elsevier Publishing, New York, (1974).

[11] Markowitz, H. M., "The Elimination Form of the Inverse and It's Appli-

• cation to Linear Programming," Management Sci., v. 3, pp. 255-269

(April 1957).

. [12] "SPAR," NASA CR 158970-1, Engineering Information Systems Inc., San

Jose, CA (Dec. 1978).
28

APPENDIX A
29

Parallel Pivoting Strategies

Matrix Step !Vlarkowitz Sum Ored Marl:owitz No.of SPICE

" NO.of OP NO.of NO.of OP NO.of Sets Fill-in %

Pivots Count Fill-in Pivots Count Fill-in Sparse

Differential 1 0 23 0 6 48 8 50

Pair

2 4 4 0 4 1,5 4 14

I

16 by 16 3 3 5 0 2 0 0 4

Nz 57

4 2 2 0 2 2 0 I

total 4 15 6 14 12 0 80.27

Cascaded 1 5 12 6 5 24 6 24

RTL

Inverter 2 -1 8 4 4 9 4 4

12 I)y 12 3 2 2 0 2 1 0 l

Nz 34

total 3 11 10 11 10 0 79.88

ECL 1 8 39 16 8 80 16 122

Colll ]):it ible

SCIIMITT 2 4 18 4 4 30 10 16

Trigger

3 2 8 0 2 16 4 6

18 by 18

Nz 66 4 2 5 0 2

• total 4 16 20 14 30 6 80.05

TABLE 1.1
3O

Parallel Pivoting Strategies

Matrix Step MarkowitzSum OredMarkowitz INo.of SPICE

NO.of OP NO.of NO.of OP NO.of Sets Fill-in %

Pivots Count Fill-in Pivots Count Fill-in Sparse .

MOS 1 4 0 0 .! 0 0 20

Memory

('ell 2 4 0 0 4 0 0 6

13 by 13 3 3 6 2 3 4 2 3

Nz 47

total 3 11 2 11 2 0 76.02

MOS 1 8 169 36 8 156 36 173

Aml)lifier

A('/DC, 2 5 117 4 ,5 6-t 4 12

3 3 0 0 3 0 0 5

24 by 24

Nz 154 4 3 0 0 3 0 0 4

5 2 13 0 2 9 0 3

total 5 21 40 21 40 10 73.76

MOS 1 8 178 40 g 156 40 149

Amplifier

Transient 2 4 91 4 .1 72 4 14

3 3 0 0 3 0 0 5

24 by 24

Nz 158 4 3 0 0 3 0 0 4

5 2 25 2 2 16 4 3

total 5 20 46 20 48 22 71.20

TABI,I': 1.1

Continued
3

1

Fill-inStatistics

Matrix Step Set size OP.Count Fill-in ' 16 by 16 1 6 23 6

Nz575 142

.450

2440

3350

4220

12by1215 126

Nz34484

342

200

2484

342

200

3220

18 by 1818 39 16

Nz 667 3010

6 216

5 122

,180

24 184

392

200

3280

42 ,50

TABLE 1.2
32

Fill-ln Statlstles

Matrix Step Set size OP.Count Fill-in

13 by 13 1 4 0 0

Nz 47 2 4 0 0

3 3 6 2

2 2 0

24 by 24 1 8 169 36

Nz 154 7 120 24

6 84 18

5 59 10

4 34 6

3 18 2

2 9 0

2 5 117 4

4 68 0

3 3 0 0

4 3 0 0

5 2 13 0

2-1 by 24 1 8 178 40

Nz 158 7 129 28

6 93 20

5 68 14

4 43 8

3 18 4

9 2

2 ,1 91 4

3 5,5 4

2 3 0

3 3 0 0

4 3 0 0

5 2 2.5 2

TABLI: 1.2

Cont in ued
33

PARALLELISM vs. FILL-IN

Matrix Step Max Reduced MAX Ol ) Reduced Fill-in Mrkowitz SI'ICI';

Size Size Count OP.Count Sum

16 by 16 1 6 5 23 1.t 2

• Nz 57 2 5 .1 13 4 0

3 4 3 6 2 0

4 2 2 0

lotal ,1 14 . 2 6 0

12 by 12 1 5 4 12 8 4

Nz 34 2 4 3 _ 4 2

3 3 2 6 2 0

4 2 2 0

total 4 11 6 l0 0

18 by 18 I 8 7 39 30 10

Nz 66 2 6 .5 26 17 6

3 3 2 11 5 0

4 2 2 0

lotal 4 16 16 20 6

I

TA BI,I-; 1.3
34

PARALLELISM vs. FILL-IN

Matrix Step Max Reduced MAX OP Reduced Fill-in Markowitz SI)IC] .,

Size Size Count, OP.Count Sum .

13 by 13 l 4 4 0 0 0

Nz ,17 2 4 4 0 0 0

3 3 2 6 2 0

4 2 o 0

tot :tl 4 12 . 0 2 0

2.1 by 2.1 1 8 7 109 120 2,1

N z ! 5.1 2 5 4 93 57 2

3 3 3 0 0 0

4 3 3 0 0 0

5 2 2 0

O 2 2 0

lolal 0 21 26 -10 lO

2-t b.v 2-1 1 8 7 178 129 28

Nz 158 2 5 -1 1-10 91 4

3 3 3 0 0 0

4 3 3 0 0 0

5 3 2 22 13 4

6 2 2 0

lotai 6 21 36 -lO 22

TABLE 1.3

Continued
35

DCOMP ICOMP

Malrix Step No.of OP. No.of No.of OP. No.of

Pivots Coun! Fill-in Pivols Coun! Fill-in

MOS 1 8 169 36 4 0 0

:_ml)lifier 2 5 117 4 7 54 18

AC/DC 3 3 0 0 5 5.t 4

4 3 0 0 3 0 0

24 by 24 5 2 13 0 2 10 0

Nz 15-t

Total 5 21 ,10 21 22

MOS 1 8 178 40 4 0 0

Amplifier 2 4 119 10 7 61 17

Transient 3 3 0 0 5 60 8

,t 3 0 0 3 0 0

24 by 2-1 5 2 25 0 2 13 1

Nz 15,_

Tol:tl 5 20 50 21 20

8-Bi! 1 72 449 150 72 449 150

Full Adder 2 25 258 86 25 2.17 70

3 16 99 0 1(l 99 0

4 11 110 12 11 110 20

1.t.t b.v 1.1.1 5 6 75 22 6 75 22

Nz 016 0 5 4t3 8 5 40 8

7 3 29 -t 3 29 4

8 2 10 0 2 10 0

9 2 8 0 2 8 0

"l'ol a l 9 1-I2 282 142 280

Comparison of the Two Proposed Ordcrings

for ULEVI';I_= 4.

TABLE 2.1
36

DCOMP ICOMP

Avg. Rows No.of Ai, g. No.of

Matrix ULEVI'I, Steps Par. Fill-in Steps Par F_ows Fill-in

2,1 by 24 2 5 4 20 4.1 5 4.2 21 18

Nz 154 3 5 4 20 ,t4 5 4.2 21 22 "

4 5 4.2 21 40 5 4.2 21 22

5 5 4.2 21 .10 6 3.6 22 26

2.1 i)y 2-t 2 5 ,1 20 50 5 4 20 20

Nz 158 3 5 4 20 50 5 ,1.2 21 20

4 5 4 20 50 5 4.2 21 26

5 5 4 20 50 5 4.2 21 2[)

52 by 52 2 5 9 45 128 6 7.8 47 121

Nz 166 3 5 9 45 128 6 7.8 47 121

4 5 9 45 128 6 7.7 46 137

5 6 7.8 47 123 6 7.7 46 137

144 by 144 2 9 15.8 142 282 9 15.8 1-12 258

Nz 616 3 9 15.8 142 282 8 17.6 141 26.1

4 9 15.8 1,12 282 9 15.8 1.12 280

5 9 15.8 1.I'2 282 9 15.8 1.12 2_0

78 b.v 78 2 10 7.7 77 198 8 9.3 7.1 238

Nz "598 3 8 9.3 7-t 200 10 7.5 75 271

4 9 8.4 76 199 9 8.3 75 287

5 10 7.7 77 202 9 8.3 75 292

505 by 505 2 36 13.6 -188 5802 .10 12.3 -t93 5-132

Nz 5889 3 37 13.2 .190 5811 .11 12.1 .197 55(il

4 34 14.3 485 5802 41 12 489 5583

5 38 13 .195 5785 36 13.4 48.1 5695

Comp:trison of D('OMI' and I('OMP for Different levels

TABLI; 2.2

1. Rrport No. NASA CR-178016

ICASE Report No. 85-48

I2. Gov..-nment Acceuion No. 3. Recipie-nt'l Clotllog No.

4. T,lIr and Subt,lIr

Multiprocessor Sparse L/U Decomposition

with Controlled Fill-In

7. Authorh)

Gita Alaghband and Harry F. Jordan

5. Rrport Data

October 1985

6. Perform,ng Orpaniution Code

8. Prrforming Or9iJniution Rrport No.

85-48

...

1-------------------------------, 10. Work Unit No.

9. Performing Orpanization Namr and Address

Institute for Computer Applications in Science

and Engineering

Mail Stop l32C, NASA Langley Research Center

Hampton, VA 23665-5225

12. Sponsoring Agrncy Name and Address

National Aeronautics and Space Administration

Washington, D.C. 20546

11. Contract or Grant NQ.

NAS1-17070, AFOSR 85-1089

13. Type of Rrport and Period Covered

Contractor Report

14. Sponsoring Agrncy Code

505-31-83-01

15. Supplementary Notes

Langley Technical Monitor:

J. C. South Jr.

Final Report

16. Abstract

Submitted to IEEE Trans. on Computers

During L/U decomposition of a sparse matrix, it is possible to perform

computation on many diagonal clements simultaneously. Pivots that can be

processed in parallel are related by a compatibility relation and are grouped

in a compatible set. The collection of all maximal compatibles yields

different maximum sized sets of pivots that can be processed in parallel.

Generation of the maximal compatibles is based on the information obtained

from an incompatible table. This table provides information about pairs of

incompatible pivots. In this paper, generation of the maximal compatibles of

pivot elements for a class of small sparse matrices is studied first. The

algorithm involves a binary tree search and has a complexity exponential in

the order of the matrix. Different strategies for selection of a set of

compatible pivots based on the Harkowitz criterion are investigated. The

competing issues of parallelism and fill-in generation are studied and results

are provided. A technique for obtaining an ordered compatible set directly

from the ordered incompatible table is given. This technique generates a set

of compatible pivots with the property of generating few fills. A new

hueristic algorithm is then proposed that c.ombines the idea of an ordered

compatible set with a limited binary tree search to generate several sets of

compatible pivots in linear time. Finally, an elimination set to reduce the

matrix is selected. Parameters are suggested to obtain a balance between

parallelism and fill-ins. Results of applying the proposed algorithms on

several large application matrices are presented and analyzed.

•

17. Key Words (Suggested by Author(s))

multiprocessor

Gaussian elimination

parallel pivoting

18. Distribution Statement

61 - Computer Programming & Software

64 - Numerical Analysis

Unclassified - Unlimited

19. Security Oauif. (of this rrport)

Unclassified

20. Security Classif. (of this pa9&1

Unclassified

21. No. of Pages

38

22. Price

A03

N-30S For sale by the National Technical Information Service. Springfield. Virginia 22161