c 

NASA Contractor Rcpt 178422 

ICASE REPORT NO. 87-75 

ICASE 

PARALLEL PIVOTING COMBINED WITH 

PARALLEL REDUCTION 

GI ta Alaghband 

fNASA-CR-178422) PARALLEL PIVOTIUG CGUlYBD ma-14662 

WITH PARALLEL REDUCTION Final Report 

41 P unclas 

WASA) 

CSCL U9B 

G3/61 0118199 

wnEraCt8 NO. NAS1-17070 and NASI-18107 ' 

December 1987 

INSTITUTE FOR COMPUTER APPLICATIONS IN SCIBMCE MID ENGINEERING 

NASA Langley Research Center, Hampton, Virginia 23665 

Operated by the Universities Space Research Assoclatim 

.. , i. r . i .: , 
Parallel Pivoting Combined with Parallel Reduction 

and Fill-in Control 

Gita Alaghband 

University of Colorado at Denver, 

Department of Electrical Engineering and Computer Science 

1100 14th Street (Campus Box 104) 

Denver, Colorado 80202- 

Parallel algorithms for triaogularization 01 large, sparse, and unsymmetric 

matrices are presented. The method combines the parallel reduction with a new 

parallel pivoting technique, control over generation of fill-ins and check for numerical stability, all done in parallel with the work being distributed over the active 

processes. The parallel pivoting technique uses the compatibility relation between 

pivots to identify parallel pivot candidates and uses the Markowitz number of 

pivots to minimize fill-in. This technique is not a preordering of the sparse matrix 

and is applied dynamically as the decomposition proceeds. 

Research was sPpportcd by the National Aeronautics and Space AdminisaatiOn under NASA 

Contracts No. NAS1-17070 and NAS1-18107 and by the Air Foxce Office of Scientific 

Research under Cantnct APOSR 85-0189 while the ubor was in residence at the Institute for 

Computer Applicatiau in Science and Engindng (ICASE), NASA Langley Research Center, 

Hampton, VA 23665. 
1. Introduction In this paper we present multiprocessor algorithms for solving 

large systems of linear equations where the coefficient matrix is sparse and unsymmetric. VLSI circuit simulation, structural analysis, partial differential equations, 

and chemical analysis are few examples of applications requiring the solution of 

such systems of equations. 

The algorithms described in the paper are designed for a shared-memory, 

MIMD model for parallel computation, in which the total memory address space is 

accessible uniformly to all parallel units. This computational model provides sync hronization mechanisms to allow multiple updates. If multiple updates are aimed 

at the same memory cell, the penalty paid is a short delay in access time. 

Given is a system of linear equations: 

Az=b (1.1) 

where the coefficient matrix, A, is large and sparse. This paper concentrates on a 

direct parallel solution method for solving (1.1) by factoring A into lower (L) and 

upper ( U) triangular matrices respectively. 

A = LU (1.2) 

The solution is then obtained by forward and back substitution steps: 

Ly = b (1.3) 

uz - y (1.4) 

To solve (l.l), a sparse matrix technique based on the following principles is used: 

a) Only the non-zero elements of A are stored. 
2 

b) 

c) 

Arithmetic operatiow arc. pcrfot:i!ed on non-zero elements only. 

During the decompopitiou fill-ius are generated, i.e. new nowzero e!emeuts 

3re created in thr pmxss of generdng zeros below the diagonai. The 

number of fill-ins is kept small. 

The three problems qt:itc.d above are all related. Even though only nonzeros need to be stored, fill-ius must he stored in the matrix structure. Therefore, 

minimization of fill-iu will result in minimization of the arithmetic operations and 

storage as well. One must find a permutation of the sparse matrix A to satisfy the 

above goal. The problem of fincling an optimum permutation to minimize fill-in is 

NP-complete [l]. and many heuristic algorithms have been developed to obtain 

near optimal solutions for this problem. Most of these heuristics find optimum 

permutations of the matrix which minimize fill-in in sequential solution process 

while they often minimize the amount of possible parallel work in parallel process. 

Therefore an ordering, or pivoting strategy to minimize some combination of fill-in 

and parallel execution time must be determined. The design of a heuristic algarithm which identifies a set of pivots to be processed in parallel while minimizing 

fill-ins is described in detail in [2], [3], and [4]. Other parallel pivoting strategies 

have also been suggested [SI, IS], [7], [8], [SI, [lo], [ll], [12], (131. In this paper we 

concentrate on parallel implementation of sparse LU decomposition procedure 

using the parallel pivoting technique described in 121, [3], and [4]. In this implementation pivots are tested for numerical stability as well as for sparsity. 

A brief description of the parallel pivoting algorithm is given in section 2. 

Section 3 describes the storage structure used in the implementation. In section 4 

the various parallel procedures to perform steps of triangularization are described 

and analyzed. In section 5 we represent actual performance results from the 

parallel implementation of the sparse LU decomposition on the HEP computer. 
3 

Finally, in section 6 some concluding remarks :we presenter!. 

2. Parallel Pivoting Algorithm The Triangulation of an nX tt ma!.iix 

A - [a,,] can be described By the following procedure. 

for K = I,?, ..., n- 1 and for each a,k#O 

Ojk - j>k 

Otk (2.1) 

For each pair a,k.kj # 0 

- Oik 'kj i>k, j>k (2.2) 

In (2.2) if o,,=O but a,k*atl fO, a fill-in is generated. It is obvious that if we have 

sufficient processors, the divide operations (2.1) for each column K can be done in 

parallel. Also, for each k the update operation (2.2) for all pairs 0,k'Okj fO can be 

done in parallel. Our experience in employing this approach has indicated that 

the sparsity of application matrices leaves parallel processes with little work to 

perform if only reduction for a single pivot is done in parallel [14]. During Sparse 

LI; decomposition it is possible to perform computation on many diagonal elements simultaneously. In parallel LU decomposition of general unsymmetric 

sparse matrices several key issues must be considered: 

a) Parallelism and fill-in are two competing issues and a balance between the 

two must be obtained. In other words minimizing fill-ins results in limited 

parallelism, and maximizing parallelism results in uncontrolled generation of 

fill-ins. 

b) A test for numerical stability of pivots must be made to ensure the accuracy 

of the solution process. 

c) In applications where the sparse linear system must be solved repeatedly, it 

must be possible to decompose structurally identical matrices using the 
4 

information produced for !.be Grst decomposition of such matrix. 

d) A storage structure snicabie Tor parallel processing must be determined. 

A heuristic algorithm has been designed in [2], [3], [4] which identifies parallel 

pivot candidates and allows the matrix to be reduced for multiple pivots simultaneously while it minimizes fill-ins. It is a dynamic algorithm which can be 

:Ipplied at any point in the decomposition phase and does not require a preordering of the input matrix. It allows pivots to be tested for numerical stability. 

Therefore at any point during the reduction, if numerically unstable pivots are 

encountered, unsymmetric permutations can be performed. The algorithm can 

then be applied to the remaining unreduced submatrix. This technique also allows 

structurally identical matrices to be decomposed using the information generated 

during the decomposition of the first matrix. In subsequent decompositions a test 

for numerical stability should be made. If the test is not satisfied, an off-diagonal 

permutation can be made and the parallel pivoting algorithm can be applied anew 

to the unreduced matrix only. 

Here we will concentrate on parallel implementation of this algorithm and 

will not go into a detailed description of its design. A complete and detailed 

description and analysis is available in [2], (31, [4]. In what follows a brief description of the algorithm and the required steps is given. The procedure to implement 

each step is presented in detail in section 4. 

Pivots that can be processed in parallel are related by a compatibility 

relation and are grouped in a compatible. In other words pivots f,,, f,,, fu are 

compatible and can be processed in parallel if and only if elements 

ai,, aji, ad, ak,, (lit, ak; are all zero. The collection of all maximal compatibles 

[15], [16] yields different maximum sized sets of pivots that can be processed in 

parallel. Several methods for generating maximal compatibles exist and they are 
5 

all based on the constructiou of an inplicntion (incompatible) table. The incompatible table gives information ihut pr.irs of incompatible pivots. Produetior of 

all maximal compatibles iar-oive! siLinnry tree search and is exponectial in thc 

order of the, mat,ri.u. This problem is solvcd by a technique which generates an 

"ordered incompalible table" haset1 on the Markowitz number [l;] of the pivot. candidates. 

The Markowitz crittriou is a heuristic for minimizing fill-ins in sparse 

matrices in sequential programming. It is based on the fact that at each step k, 

the maximum number of Ell-ins generated by choosing alj as pivot is 

(rl- l)(c,- l), where (rl- 1) and (r,- 1) are the number of nonzeros other than alj 

in row i and column i of the reduced matrix. Markowitz selects as pivot element 

at step k, the element which minimizes (II- l)(c,- l), which is called the Markowitz number of element aij. 

An "ordered Compatible" can then be produced directly from the 

ordered incompatible table without the need to search the tree. The resulting set 

of compatible pivots has the property of generating few fills. The heuristic algorithm combines the idea of an ordered compatible with a limited binary tree 

search to generate several sets of compatible pivots in linear time. An 

elimination set to reduce the matrix is generated and selected on the basis of a 

minimum Markowitz sum number (sum of the Markowitz number of pivots in a 

compatible). Several parameters are introduced to trade off parallelism for fill-in 

which can be controlled by the program. In summary the algorithm requires the 

following steps: 

1. An incompatible table is constructed by scanning the sparse matrix. 

2. Pivots are ordered according to their Markowitz numbers. 
6 

3. A limitcd binary tree search produces several starting ~ats at a given 

level (I;LEVEI,) of tbe tree. 

4. .An ordered ronrpatible is generated for each sti?rtiug set at ‘JLLL’EL 

from the corresponding ordered incompatible table. 

5. The ordered ronrpatible of maximum size and minimum Markowitz sum 

is selected as the elimination set to reduce the matrix. 

6. A set of program parameters can be applied to the resulting 

elimination set to furt,her minimize fill-in. 

3. Storage Structure The basic global data structure used in the parallel LU 

decomposition program is described below. Each element of the matrix structure 

consists of five fields: the real numerical value, the row index, the column index, a 

pointer to the next element in the row, and a pointer to the next element in the 

column. The incompatible table is represented by an array of dimension n, order 

of the matrix, with elements of the array imptbl being sets of n elements each. 

Each set corresponds to a column of the table. Column i of the table, irnptbli, 

holds the incompatible information for pivot i of the matrix. Note that the parallel pivoting algorithm considers only the diagonal elements as pivot candidates. 

Unsymmetric permutations are possible in between parallel pivoting steps. 

cornpsi holds the resulting elimination set. -. 
Type Definition: 

ptr= rnotpac; 

matpac * record 

pointer ~ype to 3 matrix element. 

wal : real; red value. 

row : integer; row index. 

column : integer; column index. 

ne : ptr; pointer to next element in row. 

nr : ptr; pointer to next element in column. 

end; 

roef=(r,c); row aod column list. 

sefs=set of l..n; set tjpe 

Variables: 

A: array( rocl,l..n)o/ ptr; matrix structure. 

no/r,nojc: array( l..n)o/ integer; number of nonteros 

imptbl: array( l..n)oj sets; 

compst: sets; elimination set. 

in row and column. 

incompatible table. 

4. Parallel LU Decomposition In order to write efficient parallel programs 

one must consider the underlying parallel architecture to which the program is to 

be applied. In an MIMD environment parallelism must be applied at the highest 

possible level in the program in order to effectively exploit the underlying parallel 

hardware. In our design and implementation we have used the idea of universal 

parallelism due to Jordan (181, [IS] which is based on writing parallel programs 

assuming that all the parallelism needed by the programmer exists throughout the 

program execution. A set of parallel programming constructs known as "the 

Force" implemented for several shared-memory MIMD computers [ 181, [IS] are 

used in the implementation of the algorithms presented in this paper. 

A high level block diagram of the program is given below. The entire 

LU Dceompor~tion program is executed by NPROC processes. These processes 

can be created by a driver routine. The parallel routines are specified by a Forcccall followed by a brief description of their function on each box. Therefore the 

body of each Force subroutine is executed by NPROC processes in parallel. After 
8 

the program is completely exerlited, the parallel procesaes are joined in the driver. 

The flowchart consists of two major loops, parallel pivoting loop and single pivoting loop. The parallel pivoting loop is executed as long as the program 

can find compatibles of more than one pivot, otherwise the single pivoting loop is 

executed. During parallel pivoting steps only diagonal elements are considered as 

pivots and unsymmetric permutations are not permitted. In single pivoting steps 

unsymmetric permutations are allowed and hence any matrix element can be considered 3s pivot. 

In the remainder of this section we describe the parallel algorithms involved in the 

LU Decompoaition program by stepping through the flowchart in the given order. 

4.1. Parallel Sort A sorting routine is required to sort the pivots in decreasing 

order of Markowitz number. The ordered list of pivots is used at several points in 

the parallel pivoting algorithm: in the construction of an 

ordered incompatible table, in the construction of a partial binary tree search, and 

finally it is used to trade off parallelism for fill-in by discarding a fraction of compatible pivots with Markowitz number higher than a given threshold value in the 

ordered list. 

The sorting algorithm used for this purpose is the Batcher sort [20], [21]. 

Batcher's sorting scheme is somewhat like Shell's sort but, the comparisons are 

made in a novel way so that no propagation of exchanges is necessary. The 

amount of bookkeeping needed to control the sequence of comparisons is rather 

large. All comparisons/exchanges specified by a given iteration can be done simultaneously. In the procedure below processes are prescheduled over a range of 

indices and they perform the comparison/exchange operations in parallel. As can 

be seen from the algorithm at each iteration we must compute the range of 

, 
d 

I ni t i:r I i 7~ 

I>nr: icr 

P:lrpv== true 

End Bnrrirr 

I / 

.) 

F’nrrllrl pivot.ing loop 

N - 

SiiigIr pivoting loop Pnrpiv? 

Rack solve 

find next pivot: 

diagonal pivoting 

first, complrte 

I Forcecall S II’PS 1.’ 

Forcecall .S”.’fRO I!’ 

permute pardlcl pivot rows 

to put uirrtris in pivot order 

+ 

find next pivot: 

diagonal pivoting 

first, complrte 

I Forcecall S II’PS 1.’ 

+ Fo rccc a I1 r c drc ce 

I 

’ I 

Forcecall SORT i 

sort pivots in decreasing 

order of Alarkowitz no. , 

? 

Fo rc cc a I 1 itip a ir 

construct the incompatiI)le 

table 

I: o rc ec a I 1 Co my 3 et 

produce elimintion set 

U I 

Parpivm false 

stability and discard 

unstable pivots 

for fill-in 

I Forcecall SIVPSY save pcrrnat ation information 
10 

nonadjacent pairs for compariwo. 'Thlu is 3 rather large overhead but is performed in parallel by prorewe\ 

of sequential code is ceecicd to 2.dlust r.iic rug? of indices for the uext iteration. 

In bet weer. tbe iterations however, a large qcction 

c 

Procedure Batcher awl 

Global l,p,q,r,d 

Barrier 

t= [lognl 

p" 2'-1 

q= 2"', roo, d=p; 

End barrier 

while (p > 0) do 

begin 

Presched DO ' i= I, n - d 

(com pute correct index ) 

q- (i- l)/p 

j= p'q +r +i 

if(i 5 (n- d)) then 

End Presched DO 

Barrier 

if( p f q ) then 

compare and exchange; 

d-q-p 

9' gR 

1-p 

P' lP/21 

else 

endif 

End barrier 

endwhile 

It has been shown that with enough parallel operations, sorting is completed in 112 [log n1 ( [log n1+ 1) steps. The sequential work and the small critical 

section used in implementation of the barrier construct will dominate the parallel 

work unless n is very large. 

The semantics for Barrier construct are such that all procelwcs pause when they reach the Barrier. 

After all have arrived, one process executes the section of code enclosed by Barrier-End barrier pair. After 

the sinKly executed code section is complete. 311 processes will resume execution after the End barrier. 

* Preached DO loop causes the body of the loop enclosed between it and the matching End Preached 

DO to be executed in parallel for different values of i. Instance8 of the loop body must be independent for 

different values of i. 
11 

4.2. Parallel Inyair The incompnfible table is constructed in thig routine. 

Each columo of this table corresponds to a pivot or the matrix and c*,,ntaim t,he 

list of pivots incompatible wi:,h the pivot under consideration. This informaiuon is 

used in the ronatruct,ion of the partial binary tree search described in the next section. Aasume pivots are numbered 1 to n Corresponding to diagonal elements of 

row9 1 tbrorigh n of the matrix ordered with decreasing hlarkowitz number. 

Columo i or the incompatible table corresponds to pivot number I of the matrix. 

Each column of the table can be constructed independently by a parallel process. 

Parallel processes are prescheduled over a loop of indices (i) corresponding to 

diagonal pivots of the matrix. Each process, say i, scans the row-column pair 

corresponding to pivot i. If a nonzero element alj or a,i is encountered a mark for 

pivot j is entered in row j of column i of the incompatible table, indicating pivot 

i is incompatible with pivot i. No process synchronization is required since each 

process is responsible for scanning row-column pairs of different diagonal elements 

and updating the corresponding columns of the table. Figure 4.1.b shows the 

ordered incompatible table for the sparse matrix of Figure 4.1.a. The algorithm 

csn be described as: 

Procedure Incompatible table 

Global imptbl[ 1 ..n] of set 1 ..n; 

Global n,nrem ; 

Presched DO i= nrem,n 

scan row i for any nonzero 0,) 

if fJ not in impfbl( i) then 

add f, to imptbl( i) 

scan col i for any nonzero a,, 

if P, not in irnptb4 i) then 

add P, to imptbl( i) 

End Presched DO 

The construction of the incompatible table requires scanning NZ nonzeros 

of the matrix. As can be seen from the procedure the only set operations required 

are addition of a new element to a set and a test for membership. These 
12 

Markowitz 

Number Pivot 

1 0 

2 0 

3 2 

4 2 

5 2 

6 4 

7 3 

8 9 

9 12 

10 4 

11 12 

8 

9 

11 

Order 

9 

11 

8 

6 

10 

7 

3 

4 

5 

1 

2 

12 3 4 5 0 7 8 9 10 11 

X 

.y 

X X X 

X X 

X X X 

X X 

X X X X 

xx X X 

X X xx X 

X X xx 

X X xxx 

Slatrix A1 

Pivots Ordered with Markowitz Number 

Figure 4.1 a 
13 

11 

8 

0 

10 

3 

4 

5 

1 

2 

- 

I 

9 11 8 6 10 7 3 4 5 1 

Figure 4.1.b Ordered Incompatible Table 

operations are O( I), therefore the incompatible table can be constructed in 

O( NZ/NPROC) time with NPROC parallel processes. 

4.3. Parallel Compset The procedure that produces the ordered compatibles 

has two major parts. The first part, generates several starting sets at a given level 

(C'LEVEL) of the binary search tree. The second part produces an 

ordered compatible for each of the starting sets from the incompatible table. 

The binary tree search is a systematic approach for extracting the maximal compatibles. Initially, it is assumed that all pivots are compatible. They are grouped 

in one set consisting of all pivot (diagonal) elements. This set, S, will be at the 

root of a binary tree, level zero. Next, the set of pivots incompatible with the 

pivot of minimum Markowitz number, PJ, obtained from the incompatible table, 

imptbl, is used to split S into a left S, and a right S2 set, constituting level one. 

S, consists of all elements of its parent S except those incompatible with Pi. S, 

consists of the same elements as S except PJ itself. At each level of the binary 

tree sets are produced by splitting the parent set into left and right sets, taking 

pivots in increasing order of Markowitz number from the ordered list of pivots to 

split the sets. This process continues until we have produced all starting sets, 
14 

,~~'/~"~",-' through .S(''-E"E[* - I at. level ULEVEL,. The p;rrtiaI binary tree 

search for the example matrix of Figure 4.1 and for ULEVELs3 is shown in Figure 4.2. ~1s We can see eigot starting sets are produced for this level of the tree. 

Note that set 5 and 6 are the same as their parent set simply because the parent 

set could not have been split for pivot number 10. Different ordering3 of pivots for 

splitting the nodes of the binary tree are considered in [2], [3], and [4]. 

In the second part of this procedure an ordered compatible is generated for 

each of the starting sets. This is done by scanning the incompatible table 

corresponding to each starting set in decreasing order of Markowitz number of 

pivots in the starting set. 

The incompatible table for a given starting set, SI, is the original table with those 

rows and columns corresponding to pivots absent from SI eliminated. 

For each starting set, SI, its corresponding incompatible table is scanned. Any 

pivot P, whose corresponding column in the incompatible table, imptblp,, is null is 

added to the ordered compatible, compset,. In addition any pivot P, for which 

imptblp n compset,= empty is also added to compset, since compset, does not contain any pivots incompatible with P,. Finally the ordered compatible of maximum 

size and minimum hlarkowitz sum is selected as the elimination set to reduce the 

matrix. The ordered compatiblea corresponding to the eight starting sets above 

are given in Figure 4.3. Any of these sets can be selected to reduce the matrix in 

parallel. Among these ordered compatibfes compset5, compset6, and compset8 are 

of maximum size (5). The set with minimum Markowitz sum will tend to generate 

fewer fill-ins. Therefore compsetd or cornp8et6, which ever is produced first, will 

be selected as the elimination set. 

1 
15 

[1,3.4,5,6.7,8,9,10,11] 

Figure 4.2 Partial Binary Tree Search 
16 

comnpset, = [2.3,7,10], Markowitz sum z- I) 

cornpaet- = [;',3.5.';;, Xtsrkowicz sum = 7 

compact3 = i2,3,9], Markowitt sum = 14 

cornpael, = [2.3,4,3!, Mwkowitz sum = A 

rompel, ~1.2,3,7,10), Markowitz sum = 0 

compact, = [1,2,3,i,lO!, Markowitz sum 9 

compret, (1,2,3,1(3], Markowitt sum = 6 

cumpaeta = [1,2,3:4,9], Markowitz sum - 11 

The Ordered Computibles and Their Markowitz Sum Number 

Figure 4.3 

To produce the starting sets at C'LEVEL, processes are assigned to the 

nodes of the part.ial binary tree from the root to level ULEVEL-1. A process cannot start. to split a set until the set is produced by the parent process. To accomplish this, a lock is assigned to each set from the root to ULEVEL-1. The lock is 

initialized to false except for the root set. As soon as a process has completed generation of a child set, it sets the lock for the child set to true allowing the next 

process t,o proceed. This is done by selfscheduling processes over the work. Production of starting sets as described above is embodied in the first self-scheduled 

loop in the algorithm below. 
17 

Procedure Compact 

Global lock( nset); 

Global imptbl( 1,n) set of I ..n; 

Global S(1..2X met); 

Global compsf( 1,neel) set of l..n; 

Local less set of I..n; 

Local tempset set of l..n ; 

nsd :number of sets from root to CU,EVEL- 1. 

Selfsched DO is 1, (2(mE'EL-') - 1) 

wait until lock( i) true; 

t.ake the next pivot, PI, with Lowest 

Markowitz number to split set, : 

produce left set. set lock(2X i) to true; 

produce right set, set lock(2x i+ 1) to true; 

End Selfsched DO 

Barrier 

End barrier 

For each starting set, S,, produce an ordered compatible, compacti 

presched DO is 21'LEtEL-l 2Uf-EmL-1 

7 

compset, = empty 

for j= n down to 1 do 

le88 = s- SI 

begin 

if ( P,B SI )then 

begin 

ternpaet = irnptbl, - less 

tempset = tempset n rompaet, 

if ( tempel- entpty ) then 

cornpeet, = compset, + [PI] 

end 

end 

find a local maximum. 

End Presched DO 

Critical ' max 

find a global maximum 

End critical 

Generation of ordered compatible8 is done by prescheduling processes over 

the sets at ULEVEL. Each process is responsible for keeping an updated copy of 

the ordered compatible of maximum size and minimum Markowitz sum it pros A procew tder the next unassigned value of i se soon aa it is free. This tends to even the work lord 

' Mutual exclusion is accomplished by critical sections, begun by Critical statement md ended by 

over processes when the execution time of the loop can vary significantly for diflcrent i values. 

End critical. 
18 

I 

diicr4. In order to do this, proccwcs c-xeciite 3 section of code to obtain 3 local 

maximum. After processes have completed the execut,ion of the prescheduled loop 

body, they execute a critical section to obtain a global maximum (eliminatton set). 

Production of K starting sets for a given C'LEVEL takes a constant time. 

For IJLEVEL small and constant c.ompared to n, generation of 

ordered compatibler from starting sets is of order n set intersection and diEerence 

operations. Assuming efficient implementation of the set operations is available, 

O( setop), the heuristic algorithm has a complexity of O(K.n-setop), where setop 

can be assumed to be constant. Employing NPROC processes will reduce the execution time of the second prescheduled loop by l/NPROC. Of course, the complete execution time cannot be improved by l/NPROC because of the synchronization code used in waiting for locks to become true in the barrier code and in the 

critical section to find 3 global maximum. As LrLEVEL is increased the number of 

parallel processes that. can be effectively used increases but at the same time the 

complexity of the algorithm will increase. It is important to choose ULEVEL such 

that the amount of parallelism provided by the underlying hardware is effectively 

exploited to speed up the execution time and not to add to its size. 

4.4. Numerical Stability and Trade off Parameters A test for numerical 

stability is done by prescheduling processes over the parallel pivot candidates in 

the elimination set. Each process searches its pivot column for the maximum 

entry, Vmaz. A pivot is numerically stable if 

pivot tolerance < lpivot valuel and 

u x I Vmaz I < !pivot value I 

pivot tolerance and u are user defined values which define the desired accuracy. If 

a pivot does not satisfy the test, it is discarded from the elimination set. When no 
19 

more parallel pivobs exist, i c. the single pivot.ing loop, if unstable pivots arc still 

present an accept.able pivot is obtained by a complete pivoting strategy. Thus an 

unrylnmetric permutat ion is performed to put the matrix in the new pivot order. 

Note that for matrices in which diagonal pivoting becomes impossible at an early 

stage during the decomposition, it is possible to switch to single pivoting code for 

a few steps. During single pivoting steps unsymmetric permutations are possible 

and will change the matrix structure. So parallel steps may become possible 

\ 

again. 

In previous papers we have shown that it is possible to minimize generation of fill-ins significantly by reducing the amount of parallel work slightly 

according to some criteria [2], [3], [4]. Trading off parallelism for fill-in is done 

according to the size of the elimination set and a number of parameters: 

1. Shrinkage parameter: By allowing a small percentage of the elimination set 

t,o be discarded we can cont.rol the number of compatible pivots to a degree 

that does not limit our parallel work by too much. 

2. Upper limit parameter: This limit would allow just enough parallel work to 

keep our parallel processes busy. 

3. Threshold parameter: In shrinking the size of an elimination set only pivots 

with hlarkowitz number higher than a threshold value in the ordered list of 

pivots may be discarded. Pivots with low Markowitz numbers do not tend to 

generate many fills and need not be discarded. 

If trade off is possible then pivots are discarded from the elimination ret 

asynchronously by parallel processes. Of course it is not necessary to use a very 

tight synchronization. It is possible to calculate the number of pivots with the 

highest Markowitz number to shrink the elimination set and to let parallel 

i processes to discard these pivots without synchronization. This approach is more 
20 

parallel but less flexible on thr *ize of the resulfing elimination art. 

i 

4.5. Parallel SWPSV Once prrrsllel pivots are determined the macrix is permuted to the new pivot order. Thc permutation information is saved iu procedure 

SWPSV. The destination addre~ses of parallel pivot candidates are stored such 

that at each parollel step roas and celumus are swapped only once. At each step 

indices of pivots in the elimination et are checked against the upper left hand 

corner of the unreduced matrix. Let nrem be the index of the next row to be 

reduced in the remaining unreduced matrix. Let npiv be the number of pivots in 

the elimination set. If i is the index of 3 pivot in the elimination aet, then if 

i < nrem + npiv - 1 

there is no need to permute row and column i. Therefore pivot i can be marked to 

avoid unnecessary permutationq. This is accomplished in a prescheduled do loop 

over the pivots in the elimination Jet. Having identified the necessary permutations, the required information is stored in permutation vectors in a self-scheduled 

loop. Each process obtains the index of a row to be swapped and updates the 

corresponding entries in the permutation vectors. This is a simple routine and is 

parallelized over the compatible pivots in the elimination set( npiv). The order of 

the sequential routine is O( npiv). since it only involves exchanges of entries in the 

permutation vectors for parallel pivot candidates. 

4.6. Parallel SWPROW The actual permutation of rows and columns is performed by routines S WPRO IC' and S WPCOL. Parallelism could be most effective 

if rows of the parallel pivots are completely permuted first followed by column 

permutations. A single step row-column permutation would involve many changes 

in the row and column lists of the matrix structure and would require tremendous 
21 

amount of synchronizing (*ode 'Thus SIYPHOIY permutes all the pivot row9 ii; 

parallel, and S WPCOL pcrinuteu all pivot columns. Barricr 3yactro2izzt iou IS 

necessary in between the r3lls to the two routines. These routines are syrumet,ric 

in the function they perform, so only SU-PROW is described. 

The rows nnd columns of the pivots in the cfirniriotion set are to be permuted w;t.h others in the remaining matrix such that the matrix is in the pivot 

order with any ordering of elements within an eliminationset. !n permuting parallel pivot rows the next row pointer and the row index fields must be updated. So 

after all parallel pivot. rows are swapped with their destination rows, each column 

of the matrix will be in increasing order of row indices obtained from the permutation vector. This suggests that we can sort the columns of the matrix according to 

the new ordering given by the permutation vector. Of course not every matrix 

column has to be sorted. Only columns having a nonzero element in any of the 

rows involved in permutation must be sorted. Therefore by constructing a bit 

vector which is the result of the union of the boolean vectors corresponding to the 

permuting rows (parallel pivot rows and their destination rows), we gather the 

indices of columns bo be sorted. Construction of the boolean vectors for parallel 

pivot rows and their destination rows is done in parallel be prescheduling 

processes over these rows. The union operation is then performed sequentially. 

Note that the union could be done in parallel using a parallel tree sum computation method. This would involve much storage for the intermediate results but 

would speed up the operation. 

Next, every column having bit position in the resulting bit vector set must be 

sorted. Each column consists of very few nonzero elements due to the sparsity of 

the matrix and hence a simple bubble sort can be used efficiently to sort the 

columns. The sorting of columns are independent operations and can be done 

simultaneously by parallel processes. This is done by self-scheduling the processes 
22 

over the work. Each proces c~sccules R small critical section to obtain the index 

of the next column to be sor:cd. The algorithm description follows: 

Procedure SN'PRO It' 

Global brow: array11 .2X r~piw] of acfs; 

bit vectors of rows to be permuted. 

Global colindez: aetJ; tit vector of column indices 

to be sorted. 

Presr hed DO ia I, 2 X npiv 

(initialize the boolean vectoru) 

brouj i)= 0 

End Presched DO 

Barrier 

End barrier 

Presched DO i= 1, 'LX npicr 

obtain index, j, of the row to be permuted. 

scan row j and for each nonzero ut& 

add k to brou(i) 

End Presched DO 

Barrier 

colindez=brow( 1)U brow(2)u . . * U brow(2x npiu) 

End barrier 

Critical neitcol 

get R local column index. j, from rolindez to be sorted. 

End critical 

if (j is a valid index) then 

sortnest: 

bubble sort column j using the information 

from the permutat ion vector. 

go to sort next 

endif 

The sort is O( n$), where nz is an average number of nonzeros per row or 

column. The sort must be done for all columns in the rolindez vector. This 

number is usually a multiple of n:, say Knz, and in the worst case could be n, the 

order of the matrix. It also involves the set operations union and the next element 

from colindez which is implementation dependent. The next element operation is 

performed within the loop and can be done in O(1). Thus on the average the 

number of operations in this routine is of order of: 

o( ~d) 
23 

If K.nz parallel processes exist,, SWPRO W can be done in time O[nt)?. 

4.7. Parallel reduce The numerical decomposition and insertion of fill-ins for 

the elimination set is parallelized by self-scheduling processes as follows: 

The reduction of each row for a pivot in the elimination aet is performed 

by a parallel process. If no more rows are left to be reduced for this pivot, the 

reduction process for the next row can be started by parallel processes looking for 

more work. The logic to do this is contained within a critical section of code. A 

process obtains a local pointer to the next row having a nonzero in the pivot 

column it is processing, advances the global pointer and exits the critical section. 

If no more rows are left to be updated for the pivot under consideration, the process advances a global pointer to the next pivot in the elimination set and obtains 

the next row pointer for the new pivot in the same manner and exits the critical 

section. Thus processes work in parallel over rows of a single pivot first and over 

the parallel pivot candidates in the elimination aet next. 

Each process is responsible for dividing the nonzero element in the pivot 

column by the pivot and subtracting a multiple of the pivot row from the row, j, 

being updated. The process must also check for a possible fill-in and insert it if 

necessary. The search for a possible fill-in must be done atomically so that parallel processes do not try to insert the same element in the same position with 

different values. Means must be provided to allow only one process to insert a 

fill-in alk and others to update the element after it is inserted. This can be done 

by locking the row j and column k of the matrix such that it can only be searched 

by one process at a time. The locking of a row and column is done by a critical 

section on elements of two asynchronous arrays, one for rows and one for columns. 

Any two shared arrays, for example nojr and nojc, that are not used throughout 
24 

t t11s process can be risen far this purpose. Of course only one element can be 

Inserted in the matrix 5tructrire 31 :\ time since the insertion causes changes in the 

pointer st riicture for rows a~d columns other than j and k alone. It is important 

to note that the probability of searching for the same element by more than one 

process at the s3me time is very low The synchronization described above is 

necessary for correct solution and does not increase the execution time by much. 

'The updating of an element must also be done atomically by processes. Again the 

probability of more than one process trying to simultaneously update the same 

element is low. This synchronization can be done by simple scoreboarding which 

must be available in machines matching our computational model. 
25 

Procedure reduce 

Barrier 

ii - index of the first pivcit iii tho ordered ma~rix. 

End harrier 

Critical nett 

c getwork: 

(get next parallel work! 

isii locd pivot index 

(get a local row index to Le reduced) 

j- nettrow in pivot columc ii 

if(i not valid) and (more pivots j then 

begin 

i= ii+ I ; advance to hhe next pivot and update 

ii= i 

goto local; get the wxt. row. 

local: 

, the glob31 pivot pointer. 

endif 

update global nertrow pointer information. 

End critical 

if( j valid) then 

begin 

uJI = uJl/uII; 

scan pivot row i and for each n,k: 

(check for a possible fill-in) 

(lock row i and column k 1 

Critical nofr(j) 

Critical nofc(k) 

divide by pivot 

if( QJk. not in matrix) then 

Critical insert 

insert the element 

End critical 

End critical 

End critical 

atomically update u,k 

ujk= - 

goto getwork 

endif 

Reduction of the matrix for a single pivot requires a complete scan over 

the matrix which can be done in time NZ. For npiv pivots the time would be proportional to npiv-NZ. NZ changes as fill-ins are encountered or a3 the matrix gets 

smaller due to the reduction. DUB [.I] reports that over a wide variety of matrices, 

the number of arithmetic operat ions performed has been empirically observed to 

be about 7*/4n where T is the number of nonzeros in the decomposed form. The 

value of T is generally not known a priori. Experience has shown that a value 
26 

5/',?N% is a *atisfactory e5tiruat.e for T although an esbinate of ordcr nlog n is more 

realivtic for problems arising from PDEs in two dimensions. 

5. Implementation Results The program IS implemented for the HEP (Heterogeneous Element Processor) pipelined shared-memory computer built by Deneicor, 

Inc. [E]. The resuits prcscnt.ed here are on a single PEM (Process Execution 

Module). Thc execution pipeline on the IiEP has eight steps. In the HEP the 

degree of simultaneous execution is limited by the length of the various pipelines 

and may be characterized by an average pipeline length. Thus on a single PEM a 

program may be expected to speed up by no more than 7.5 to 9.5 over single 

stream execution. The LU Decomposition program has also been simulated on a 

Vax 11/780 and tested on many application matrices arising from electronic circuits and structural analysis producing successful results [2], [3), 141. Here we 

represent the timing results of running the program over a 141 by 144 matrix from 

the circuit of an 8-bit full adder and employing different values for trade off 

parameters. Figure 5.1 represents the execution time of LU Decompoaition program for different numbers of processes from 1 to 25. The result is for the case 

when maximum parallelism is used. For NPROC=l, the matrix is completely 

reduced in 10 parallel steps. The number of compatible pivots at each step is 72, 

25, 16, 11, 6, 5, 3, 2, 2, and 1 respectively. Note that in the first step half of the 

matrix is reduced in parallel. The execution time decreases with an increase in the 

number of processes up to NPROC-11. In fact there is UNPROC reduction in 

execution time for small values of NPROC as new processes make efficient use of 

the execution pipeline. This decrease in execution time bottoms out as the pipeline becomes full. The slope of the linearly rising tail of the curve indicates the 

length of time spent in critical sections in various points in the program. A com- 
27 

EXECLITIBN TIHE 

I SECBNOS I 

3. SGO 

3.000 

2.500 

2.000 

I. 500 

I .ooo 

0.500 

0. 

Y 

Y 

X 

X 

Figure 5.1 Execution Time vs. Number of Processes 

No Trade off 

144X 144, NZ= 616,8-Bit Full Adder 
28 

plete model for analysis of parallel programs can be found in (2.31. Defining the 

speed up to be: 

T(1) 

T( YPROC) Sa 

where T( 1) is the time to execute the program with one process and T( "VPROC') is 

the same time using 'VPROC' processes. Then a speed up of 4.82 is obtained for 11 

processes. Note that t,his is not speed up measured with respect to the best 

sequential algorithm, but only gives insight to the parallelism in this program. 

For a small number of processes, execution time versus number NPROC 

of processes can be represented as: 

C? T(NPR0C) = C,+ NPROC 

where G, represents the sequential portion of the work and C, the parallel portion. A simple least squares fit to determine C, and C, is applied to a linear portion of the execution time versus NPROC curve to estimate the degree of parallelism. This analysis shows that the code is 87% parallel. Figure 5.2 shows the execution time versus XPROC for individual routines. As can be seen there is a 

sharp increase in the slope of the sort curve for large NPROC which indicates 

parallel processev spend more time in the critical section than doing useful parallel 

work. Of course a reason for this behavior is the small value of N, order of the 

matrix. As N increases the slope becomes smaller. A speed up of 3.3 for 8 

processes is obtained for the sort. The degree of parallelism for this routine is 

58%. 

For the Compact routine speed up is 5.7 for 10 processes, and the code is 93% 

parallel. Here the value of ULEVEL is 4, so there are 10 starting sets for which an 

ordered compotibfe must be produced in parallel. Therefore every time number of 

active processes divides 10 evenly a sharp decrease in execution time is observed. 

Speed ups for reduce and SWPROW are 4.4 for 11 processes and 0.5 for 12 
29 

EXEillTlBN TlnE 

I ZECBNOS I 

1.190 r 

1 .ooo 

0.900 

0.800 

0.700 

0.600 

0.500 

0.400 

0.300 

0.200 

0.100 

0. 

CBHPSET x 

S0RT 4 

SWPRBW 0 

REDUCE 

Figure 5.2 Execution Time vs. Number of Processes 

for Individual Routines, No Trade off 

144X 144, NZ=616,8-Bit Full Adder 
30 

. .- 

procrqws rcrpectively. Tlw drgn 1. of p:ir:~lIvliqin for t,tlcsc two rc,iitint.+ WI- 8 j";, 

and 93% respectively. Our nuiricricai results about number of availabie paraiic! 

pivots indicate the existence of L~YQ parsilel operations snd it is clear that by 

moving to comp1ltei.s with more parallel units or PEMs in the cast? of HEP a 

higher degree of parallelisoi can be achieved. 

The program is noncletermirlistic when executed in parallel. In many cases 

there are several sets of cquai maximum size and minimum Markowitz sum. 

Depending on the number of processes and their relative speeds, one of the candidate ordered compatibles will be selected as the elimination set. Thus diflerent 

results are produced. The number of fill-ins generated for single stream execution 

is 280. For different values of NPROC this number is in the range of 89% to 

103% of the fill-ins produced sequentially. 

Table 5.1 shows the result of running the program on the same matrix when trade 

off parameters are used. The values of parameters for this case are given below: 

Threshold 1/3 

Shrinkage Parameter 30% 

Upper Limit 25 

CrL E \'EL 'I 

Table 5.1 

Routine NPROC Speed up 

LU Decompoae 9 5.81 

Bateher 7 3.5 

Com pa et 18 6.1 

SWPROW 14 0.54 

reduce 9 4.15 
The higher speed up indicates that by employing the above paramcters a better 

balance between number of compatible pivots generated at different steps is 

achieved. A reduction of 23% in 811 iu obtained as the result of the above parameter variations which compares reasonably with results from the best sequential 

program (100). The fill-in can further be decreased by assigning different values to 

trade off parameters. 

The results of running the program on a 505 by 505 matrix produced from 

SPAR, a structural analysis program [24] is given in Figure 5.3 and 5.4. In the 

NEP, for every new process a local stack area is allocated. This area depends on 

the amount of local storage and some other system parameters. Due to the limited 

size of the available memory the program could only be run on this matrix for up 

to 7 processes. The trade off parameters for this run have the following values: 

Threshold 213 

Shrinkage Parameter 40% 

Upper Limit 00 

UL E VEL 4 

Using a least squares fit the degree of parallelism for each routine is summarized in 

Table 5.2. 

Table 5.2 

_- 

Routine Degree of Parallelism 

LU Decompose 95% 

Botcher 93.37% 

Compset 91.45% 

SWPROW 92.22% 

reduce 85.44% 
60.000 

SO. 000 

40.000 

30.000 

20.000 

1o.coo 

0. 

32 

EXECUTION TlY 

I SECBROSl 

70. OOC CBMPSET X 

SBRT 

SWPRBW 0 

REOUCE I 

Figure 5.3 Execution Time vs. Number of Processes 

for Individual Routines 

ULEVEL=4, Shrinkage=40%, Upper Limit=80, Threshold- 2/3 

505 x 505, NZ-5889, from SPAR Program 
33 

EXECLlTl'dN TIHE 

. NPRK 

Figure 5.4 Execution Time vs. Number of Processes 

Comparison of Parallel Pivoting and Reduction 

ULEVEL - 4, Shrinkage-40%, Upper Limit-60, Threshold- 2/3 

505X 505, NZ- 5880, from SPAR Program 
34 

The 93% parallelism from t,he sort indicates that, for large values of 'C', parallel 

processes spend more time perfarming parallel operations than in the critical section. 

Figure 5.4 compares the total execution time spent in routines to Bnd the parallel 

pivoting candidates and execution time of the rest of the program for LU dccomposition. As can be verified from thib figure the time spent to Bnd the parallel 

pivots is much less than the time to perform the decomposition. This difference 

increases as the size of the matrix becomes larger, verifying the advantage of 

parallel pivoting. 

6. Conclusion A set of parallel algorithms for performing LU decomposition of 

general unsymmetric sparse matrices for shared-memory MIMD computers has 

been presented. The sparse LU decomposition technique employs a parallel pivoting strategy to solve the problem of having enough parallelism in sparse matrices. 

The main features of the heuristic algorithm can be summarized as follows: 

-It can ident,ify a good set of parallel pivots in linear time. 

-It is 3 stepwise algorithm and can be applied to any submatrix of the original matrix. Thus it is not a preordering of the sparse matrix and is applied 

dynamically as the decomposition proceeds. 

-Pivots can be tested for numerical stability and unsymmetric permutations can be performed if necessary. 

-Trade off between parallelism and fill-in is possible under several program 

controlled parameters. 

-The information produced by the algorithm can be stored to decompose 

structurally identical matrices. 
35 

We have presented the parallel reduction combined wit.h parallel pivoting 

tecl~uic~ue, control over the gcoc-ration of fills and check for numerical stability, all 

done in parallel with work being diqtribiited over the active processes. The program verifies that it is actually possible to do parallel pivoting in sparse matrices 

on multiprocessors and take advantage of the existing parallelism in the problem 

and in the hardware. The timing aualysis of the routines indicate that every routine has been eflectively parallelized. The small slope in the execution time versus 

number of processes of L Ir Decornpoaition program which represents the amount 

of synchronization overhead verifies the effectiveness of parallelization and 

machine utilization. 

Acknowledgment The author would like to express her sincere appreciation to 

Dr. H. F. Jordan for his support and constructive comments and advice. 

, 
36 

Refer en ces 

(11 M. Yannakakis "Comput ins the Minimum Fill-in is NP-complete," SIAM J. 

.4lg. DISC. .Wath. 2. pp. 77-79. 1981. 

[2] C. Alaghband ",\fulfrproce~aor Spnrae LCr Decomposition with Controlled 

Fill-in." Ph.D. thesis, Gepartment of Electrical and Computer Engineering, 

1:niversity of Colorado, Boulder, May 1980. 

(31 G. Alaghband. H. F. Jordan "Spar3e Gaussian Elimination with Controlled 

Fill-in on a Shared hlemory Xlultiprocessor," ECSE Technical Report 86-1-5, 

Electrical and Computcr Engineering Department, LJniversity of Colorado, 

Boulder, November 1986. 

[4] C. Alaghband and H. F. Jordan "Multiprocessor Sparse L/U Decomposition 

with Controlled fill-in," IC.4SE Report No. 85-48, NASA Langley Research 

Center, Hampton, Virginia 23065, 1985. 

[5] D. A. Calahan, "Parallel Solution of Sparse Simultaneous Linear Equations," 

Proc. 11-th Annual .4llerton Con/. Circuits and System Theory, pp. 729-735, 

Oct. 1973. 

[6] J. W. Huang and 0. lying "Optimal Parallel Triangulation of a Sparse 

Matrix," IEEE Trans. on Circuits and Systems, vol. CAS-20, No. 9, pp. 726- 

732, Sept. 1970. 

[7] 0. Wing and J. W. Huang "A computation Model of Parallel Solution of 

Linear Equations," IEEE Trans. on Computers, vol. (2-29, No. 9, pp. 032-038, 

July 1980. 
37 

, 

[8] Y. F. Zhou "Optim31 Farnilel Triangulation of a Sparse Matrix- A Graphical 

Approach," IEEE 19HI Syrup. un Circuits and Syatcms. 

[O] K. Nlskajim,r ".4 Graph 'I'hroretical .4pproach to Parallel Triangullt.ion of a 

Sparse Asymmetric hiatrix," Proceedings of 1984 Con/. orb Informution Science and Systcnis. 

(10; J. A. C;. Jess and H. G. hl. Kees "A Data Structure for Paralle; LL; Decornposition," IEEE Trona. on Computers, vol. C-31, no. 3, pp. 231-239, March 

1982. 

[l I] F. J. Peters "Parallel Pivoting Algorithms for Sparse Symmetric Matrices," 

Parallel Computing 1, pp. 99-1 10, 1984. 

[ 121 M. R. Leuze "Parallel Triangularization of Substructured Finite Element 

Problems," IC.4SE Report no. 84-47, Sept. 1984. 

[13] I. S. Duff "Parallel Implementation of Multifrontal Scheme," Argonne 

National Laboratory, hfat bematics and Computer Science Division, Technical 

Memorandum no. 49, hlarch 1985. 

[l4] G. Alaghband and H. F. Jordan "Parallelizing a Sparse Matrix Package," 

Report CSDC 83-3, Computer System Design Group, Electrical and Computer Engineering Department, University of Colorado, Boulder, June 1983. 

[IS] Z. Kohavi "Switching and Finite Automata Theory," Computer Science Series, 

Second Edition, McGraw Hill Book Company, 1978. 

[IS] D. Lewin "Logical Design of Switching Circuit#," American Elsevier Publishing, New York, 1974. 

[17] H. M. Markowitz "The Elimination Form of the Inverse and its Application 

to Linear Programming," Management Science, 3, pp. 255-269, 1957. 
38 

[18] H. F. Jordan "Parallel Computation with the Force," ICASE Report no. 85- 

45, NASA Langley Research Center, Hampton, VA, October 1985. 

, 

I191 H. F. Jordan "Structuriug Parallel Algorithms in an MIMD, Shared Memory 

Environment," Parallel Computing, May 1986. 

[20] K. E, Bat.cher Proc. AFlPS Spring joint Computer ConIerence, 82, pp. 307- 

314, 1968. 

[Zl] D. E. Knuth "The -4rt of Computer Programming, Sorting and Searching," 

vol. 3, Addison-Wesley, 1973. 

(22) J. S. Kowalik The HEP Supercomputer and Its Applications, Ed., MIT Press, 

1985. 

[23] H. F. Jordan "Interpreting Parallel Processor Performance Measurements," 

Report CSDG 85-1, Computer System Design Group, Electrical and Computer Engineering Department, University of Colorado, Boulder, November 

1985. 

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

Jose, CA, Dec. 1978. 

. 
Report Documentation Page 

19. Security Classif. (of this report) 

Unc la s si f i ed 

1. Repon No. I 2. Government Accession No. 

21. No. of pages 22. Price 

A03 

20. Security Cld. (of thii pa I 

Unclassified'@ 40 

NASA CR-178422 

ICASE Report No. 87-75 

4. Titk and Subtitle 

PARALLEL PIVOTING COMBINED WITH PARALLEL 

REDUCTION 

9. Performing Organization Name and Address 

Institute for Computer Applications in Science 

Mail Stop 132C, NASA Langley Research Center 

Hampton, VA 23665-5225 

wring Agency Neme and Address 

gtional Aeronautics and Space Administration 

Langley Research Center 

Hampton, VA 23665-5225 

and Engineering 

12. S 

3. Recipient's Catalog No. 

5. Report Date 

December 1987 

6. Performing Organization Code 

8. Performing Organization Report No. 

87-75 

10. Work Unit No. 

505-90-21-01 

11. Contract or Grant No. 

NAS1-17070, NASl-18107 

13. Type of Ropon ad Period Covered 

Contractor Report 

15. Suppkmenta 

Langley ?eN2Eical Monitor: 

Richard W. Barnwell Journal 

Submitted to Parallel Computing 

I Report 

16. Abstract 

Parallel algorithms for triangularization of large, sparse, and unsymmetric 

matrices are presented. The method combines the parallel reduction with a new 

parallel pivoting technique, control over generations of fill-ins and check for 

numerical stability, all done in parallel with the work being distributed over the 

active processes. The parallel technique uses the compatibility relation between 

pivots to identify parallel pivot candidates and uses the Markowitz number of 

pivots to minimize fill-in. This technique is not a preordering of the sparse 

matrix and is applied dynamically as the decomposition proceeds. 

17. Ke Words ISuggestd by AuthorIs)) 

mu1 tlprocessor , Gaussian elimination, 

parallel pivoting 

18. Distribution Statement 

61 - Computer Programming and 

64 - Numerical Analysis 

Unclassified - unlimited 

Sof tware 

I 1 

NASA FORM 1- OCT 88