o

NASA Contractor Report

ICASE Report No. 90--48

182072

ICASE

THE USE OF LANCZO'S METHOD TO SOLVE FHE

LARGE GENERALIZED SYMMETRIC EIGENVALUE

PROBLEM IN PARALLEL

(_'_AqA-L_-_)'_L;7=) ]'_- UCL G r LA\!C7U'5 M_[,H0 r_

Mark T. Jones

Merrell L. Patrick

IN P^k _ ..... eLL[ _inol R_:porL

t.g_L n,t3

q3,/ol

Contract No. NAS1-18605

July 1990

Institute for Computer Applications in Science and Engineering

NASA Langley Research Center

Hampton, Virginia 23665-5225

Operated by the Universities Space Research Association

National Aeronautics and

Space Administration

I_ngley R_rch C_nter

Hamplon, Virginia 23665-5225

The Use of Lanczos' Method to Solve the

Large Generalized Symmetric Eigenvalue

Problem in Parallel

Mark T. Jones*and Merrell L. Patrick *t

Abstract

The generalized eigenvalue problem, Kx = ._Mx, is of significant practical importance, especially in structural engineering where it

arises as the vibration and buckling problems. New software, LANZ,

based on Lanczos' method has been developed for solving these problems and runs on SUN 3, SUN 4, Convex C-220, Cray 2, and Cray

Y-MP systems.

Preliminary results of using the Force to obtain a multiprocessor

implementation of LANZ on MIMD paraUel/vector systems are reported here. A parallel execution time model of LANZ is defined

and used to predict the performance of LANZ as well as examine

hypothetical modifications to LANZ. The results of using dynamic

shifting to improve parallelism are presented. Finally, the results of

assigning a group of processors to separate shifts and finding all the

desired eigenvalues using LANZ in parallel are reported.

*Department of Computer Science, Duke University, Durham, NC 27706

t

This research was supported by the National Aeronautics and Space Administration under

NASA contract No. NAS1-18605 and AFOSR grant No. 88-0117 while the authors were

in residence at ICASE. Additional support was provided by NASA grant No. NAG-I-466.

1 Introduction

The generalized eigenvalue problem, Ks: = )_Mz, is of significant practical importance, especially in structural engineering where it arises as the vibration and

buckling problems. New software, LANZ, based on Lanczos' method has been

developed for solving these problems and been reported on in [1] [2]. LANZ

uses a technique called dynamic shifting to improve the efficiency and reliability

of the Lanczos algorithm [3]. Improved methods for solving symmetric indefinite linear systems and for finding eigenvalues of the tridiagonal matrices that

arise when using Lanczos' method have been developed for use in LANZ [4].

Loop unrolling techniques were used to obtain improved performance on vector

processing machines [5]. An improved version of Parlett and Scott's selective

orthogonalization algorithm [6] was used to maintain orthogonality of the Lanczos vectors. Implementations of LANZ on a Convex C-220 were used to study

the performance of LANZ and compare it with a subspace iteration code used

by structural engineers. In all cases tested LANZ had superior performance to

the subspace iteration code.

In this work, the Force is used to obtain a multiprocessor implementation

of LANZ on MIMD parallel/vector systems. Strengths and weaknesses of the

Force as a parallel programming language are examined in Section 2. A parallel execution time model of LANZ is defined and used to compare the actual

parallel performance of LANZ against the performance predicted by the model

in Section 3. It is also used to predict the effect of changes to LANZ without the actual implementation of the changes, e.g., it is used to compare the

performance of vectorized versus non-vectori_ed versions of LANZ on a Cray

Y-MP. Backward triangular matrix solution is shown to be a bottleneck in the

parallel version of LANZ and the tradeoffs involved in designing algorithms

for a parallel/vector machine are discussed in Section 4. The results of using

dynamic shifting to reduce the number of backward solves by increasing the

number of factorisations are reported in Subsection 5.1. Finally, the results of

assigning groups of processors to work in parallel on separate shifts are reported

in Subsection 5.2.

2 Parallel Implementation

LANZ consists of over 15,000 lines of FORTRAN source code with four major

computations: 1) factorization, 2) triangular matrix solution, 3) matrix-vector

multiplication, and 4) computation of the eigenvectors. Observation of the

runtime profile shown in Figure 1 of LANZ on one processor of the Cray YMP reveals that over 97 percent of the execution time is spent on these four

computations. 1

i For this profile, LANZ took fifteen steps to find the five lowest elgenvalues.
Computation Time Percentage

(seconds)

Factorization 13.98 67.9

Matrix Multiplication 3.09 15.0

Triangular Matrix Solution 2.73 13.3

Computation of Eigenvectors 0.27 1.4

Other 0.49 2.4

Figure 1: Execution profile of LANZ on a single processor of the Cray Y-MP

The parallelization of LANZ focuses on these four computations. Two constraints directed the selection of a parallel language to be used in carrying out the

parallelisation: 1) the language must require little rewriting of the non-parallel

source code, and 2) the resulting implementation should be transportable to

several shared-memory parallel architectures.

The Force is a set of extensions to FORTRAN that provides a shared memory model of parallel processing on MIMD shared memory architectures [7]. The

Force includes constructs for both fine- and coarse-grain parallelism. The Force

starts a process on each processor. Each Force process communicates through

shared variables and synchronizes using barriers and critical regions. Loop iterations are partitioned among Force members by prescheduling or self-scheduling

constructs. Subroutines must be declared as Forcesubs if Force constructs are to

be used within them. These "Forcesub"s are called using Forcecall statements

and all the processes must encounter this Forcecall statement and execute the

subroutine. The Force is implemented as a preprocessor to FORTRAN on the

following shared-memory computers: Cray Y-MP, Cray/2, Cray X-MP, Flex/32,

Alliant FX/Series, Convex C220, Sequent Balance, and the Encore Multimax.

The Force has been shown to be useful for implementing parallel linear algebra

algorithms [8].

The use of the Force to parallelize such a large program revealed a shortcoming in the language. The Force assumes that all the processors are simultaneously executing a segment of code unless a synchronization construct has

specified otherwise. This is a useful model when implementing a short algorithm such as a factorization subroutine [8]. However, in a large, practical code,

such as LANZ, which contains several different algorithms and a fair amount

of i/o, this is not the best model for computation. Most of the source code in

such a large application is executed very few times, if at all, and is inherently

sequential. The major portion of the execution time is spent in a few subroutines that execute factorization or matrix multiplication algorithms. In these

subroutines the Force model of parallel computation is the preferred one. Using

the Force required a large number of unnecessary synchronization constructs in
the inherentlysequentialcode. Although thisdid not have a largeimpact on

performancc, itwas a time-consuming and error-proneprogramming task. Also,

the Force does not allow a set of processorsto execute one segment of code in

parallelusing barriersto synchronize while another set of processorsexecutes a

differentsegment of code. The barriersin the Force assume that allor none of

the processorswillencounter each barrier.

The additionof a few capabilitiesto the Force language would alleviatethese

problems and provide added functionality.Ifan entiresubroutinecould be designated as sequential,and yet retainthe abilityto declareshared variables,then

the necd for synchronizationconstructsin inherentlysequentialcode would be

alleviated.The abilityto calla Force subroutine from a subroutine designated

as sequentialwould alsobe necessary.For example, the Force requiresthat the

main program in LANZ be declared a parallelroutine so that ithas the capabilityto declareshared variablesand callsubroutines that willbe executed in

parallel.However, nothing in the main program should be executed in parallel,

therefore,severalBarrier statements arc requiredto forcethe code to execute

sequentially.The abilityto declarethe main program to be sequentialwould

alleviatethe need for unnecessary Barrier statements and stillallow the main

program to declareshared variablesand callparallelsubroutines.The Forcecall

subroutine should be extended to allow the number of processesexecuting the

subroutine to be specified.This would allow setsof processorsto execute differentsegments of code in parallel.The barrierconstructwould alsohave to be

extended to accommodate situationsin which allthe processorsare not executing the same code. The utilityof thisaddition isillustratedin Subsection 5.2.

The addition of these constructswould be very useful,but may be difficultor

impossibleto implement as a preprocessor.

In addition,the Barrier synchronization construct provided by the Force

proved to be too expensive for use in the loops that occur in the factorization subroutine. In itsplace,a user-implemented construct which was much

lessexpensive was used. The Force barriercan not take advantage of the fact

that it residesin a loop, whereas, the user-implemented construct is specificallydesigned to execute in a loop.Conceptually,the synchronizationdesiredis

shown in Figure 2. Ifimplemented using the Force,line2 becomes a UBarrier"

statement and line4 becomes an "End barrier"statement. These statements

are expanded by the preprocessor and the code is shown in Figure 3, where

"LOCKON" and "LOCKOFF" are provided by the operating system and are

very expensive. The user-implemented constructuses a Force "Barrier" before

execution to initializevariablesand then uses the shared array "commun" to

synchroni_.c.The loop as implemented using the user-implemented construct

isshown in Figure 4, where "me" is a processor'snumber and "nprocs" isthe

number of processors. The user-lmplemented construct has the advantage of

not requiringa callto an operating system routine._ The construct is not a

2The i_tiallzation itep doel require such a call.
1)

2)

3)

4)

5)

6)

z)lo

fori= 1,ndo

execute code in parallel

wait forallprocessorsto reach thispoint

execute code with one processor

tellallprocessorsto resume executionat thispoint

execute code in parallel

continue

Figure 2: Loop with desiredsynchronization

C

C

I)

2)

3a)

3b)

3c)

3d)

3e)

3f)

3g)

4)

5a)

Sb)

sc)

5d)

5e)

5_)

5g)

s)

7)10

assume that nurnhere is0, the lock Iockvarisoff

and the lock lockwaitison

do 10 i: 1, n

execute code in parallel

calllockon(lockvar)

if(numhere.lt.(nprocs- 1))then

numhere : numhere + 1

calllockoff(lockvar)

calllockon(lockwalt)

endif

if(numhere .cq.(nprocs-1)) then

execute code with one processor

endif

if (numhere.eq.0) then

call lockoff(lockvar)

else

numhere-- numhere- I

calllockoff(lockwait)

endif

execute code in parallel

continue

Figure 3: Loop with Force synchronization

4
Oa) commun(me) = 0

Ob) Barrier

0c) End barrier

i) dolOi= l,n

2) execute code in parallel

3a) if (me.eq.1) then

Sb) j= 2

3c) 20 continue

3d) do 30 j = j, nprocs

Se) if(commun(j).ne.i) goto20

3f) 30 continue

4) execute code with one processor

sa) commun(1)= i

5b) else

5c) 40 continue

5d) commun(mc) = i

5e) if (commun(1).ne.i) then

5f) goto 40

5g) endif

5h) endif

6) execute code in parallel

7) 10 continue

Figure 4: Loop with user-implemented synchronization

general substitute for the "Barrier-End barrier" construct of the Force because

one initialization step must take place for every user-implemented construct.

This initialization step requires the equivalent of a Force UBarrier." The userimplemented construct is superior when used in a loop because the cost of the

initialization step is amortized over the number of iterations in the loop. Extensive experiments have shown that the user-implemented construct is superior

to the Force "Barrier" and to the synchronization constructs available in the

Cray autotasking package on both the Cray Y-MP and Cray-2 when executing

in dedicated mode[9].

5
3 Parallel Execution Model of LANZ

A parameterized parallel execution model of the computations in LANZ has

been constructed to allow the prediction of performance on parallel computers of varying characteristics, as well as for problems of varying size and type.

This model is based on the parallel version of LANZ outlined in the previous

section. One application of the model is the comparison of the actual parallel

performance of LANZ against the performance predicted by the model. This

ensures that the parallel implementation of LANZ is performing as expected.

Two other applications of the model are: 1) prediction of performance on different architectures, and 2) prediction of the effect of changes to LANZ without

the actual implementation of the changes. Given the parameters listed in Figure 5, the cost of execution on p processors can be estimated using the following

model:

T(V)= ,_jtj(p) + ,_.t1,(v)+ ,_.t.(p) + ,,,t_._(v) + -.t.o(v) + to(v). (1)

Because of its complexity, the model is split into the following submodels:

G8

tj(p) = _[s.Cp)+ ("'(" 2- 1)_. +_ _(_2 1)(c_ + c(1)+ (:)

j=l

¢18

c.=((o- _,,1,(_.- 1111+

i=1

gIs(P) = n[fs(p)+ asCas-1) 'as)]2 CCm+ C(1)'4-c"(--_ (3)

t_.(p) = _s[f.(p )-4- asCas2- 1)(c"_ + c(1) + as(p- 1)c. + (4)

_p(-_, as)]

2n t_,,,Cp) : --c.(1=(_,)+ f, Cp)+ c,.(2, P P p-- 1) (5)

1"1,s _ 1¢i.

t,_Cp) = ---_c,_:(p, 1) + (tmra + c¢,(n, 1) + e,jm(n)) (6)

_°(p) = [2c.,_(,_)+ 4c..(,,.,1)+ s_,_(,_,1)]+ (_1

,_i.(-_, _.)(t._(p) + 2_.(., 1)+ c_._(_)+ _..(., 1)1.

The operation, extended saxpy of size j, is used in the submodels and is defined

8.8

J

= • - _ _,,, (s)

i=l
where z is a vector, a_ is a scalar quantity, and each y_ is a separate vector. The

extended saxpy operation takes full advantage of the vector processors on the

Cray Y-MP.

The cost of a single factorization, s t/(p), is given Equation 2. The extended

saxpy operations in the second line of the equation dominate the computational

costs. To parallelize the computation, at each of the i7" steps a processor is

given _ extended saxpys to compute. Because this part of the computation P

parallelizes well and dominates the execution cost, good parallel speedup is

expected in the factorisation algorithm.

The models for triangular matrix solution, t/_(p) and _b_(P), are given in

Equations 3 and 5. 4 These computations do not parallelise nearly as well as

factorization. Two reasons for this poor parallelisation are: 1) the ratio of

computation to synchronization is much lower than for factorisation, and 2) as

the number of processors, p, increases, the efficiency of the vector operations

in this computation deteriorates. The reason for this deterioration in vector

performance is now given.

The dominant cost for forward triangular matrix solution is the extended

saxpy operation. This operation has been parallelized by splitting the vectors

into p pieces, thereby decreasing the vector lengths as p increases. Thus, as the

parallelism increases, the vector efficiency decreases.

The case is somewhat more complicated for backward triangular matrix solution because the dominant operations are j inner products, s This analysis

is based on the assumption that the inner products have been parallelised by

splitting the vectors into p pieces as was done for forward triangular matrix solution. If this is done, then the same analysis that was used in forward triangular

matrix solution holds.

Another approach is possible, however, if the size of the matrix-matrlx operations, j, can be arbitrarily specified. 6 The alternative approach is to assign _,

where j is evenly divisible by p, inner products to each processor to compute m

parallel. Thus, the vector lengths are unchanged as p increases. However, this

approach has two drawbacks: 1) if j is not evenly divisible by p, then poor load

balancing occurs, and 2) a processor on the Y-MP can compute several inner

products simultaneously more efficiently than it can one or two inner products,

therefore if £ is small, the operations will be inefficient.

p

A model, t,_m(p), for the cost of the multiplication of a sparse matrix 7 times

$Thls submodel is constructed under the assumption that a variable banded factorlzation

algorithm utilizing matrix-matrlx operations is being used. The advantages of matrlx-matrlx

operations are discussed in [2].

4The assumption is again made that rrmtrlx-matrix operations will be used.

sj is the size of the blocks in the matrlx-matrlx operations.

eThe block size can be arbitrarily specified if a symmetric positive definite matrix is being

factored, but not if a symraetric indefinite matrix is being factored. This is discussed in more

detail in [2].

7The assumption is made that M, if the vibration problem is considered, and Ka, if the

buckling problem is considered, are sparse.
a vector is given in Equation 5. The computation is parallelised b_ assigning

columns of the matrix to each processor. Each processor uses saxpy operations

to compute the contribution of each of its columns to the result vector. Each

processor stores the sum of the contributions of its columns in its own vector.

After all the processors have finished computing their own vectors, each computes a portion of the result vector from the contributions of each processor

using saxpy operations. Because only one synchronization step is required, a

speedup of almost p on p processors can be expected.

This algorithm is superior to the obvious algorithm of assigning rows of the

matrix to each processors and letting each processor compute a single element of

the result vector. Although the obvious algorithm requires no synchronization, it

does require the use of vector inner products which are far slower than the saxpy

operations used in the previous algorithm. For large p the cost of communicating

the partial results would become too expensive and the obvious algorithm would

be superior.

The model for the cost of computing an eigenvector, t,c(P), is given in Equation 6. The major computations used to compute an eigenvector in LANZ are:

1) yi = Qjs_, which is a full matrix multiplication, and 2) the normalization of

y_ to ensure that y/TMlli : 1. The full matrix multiplication can be partitioned

in a fashion similar to sparse matrix multiplication with similar results to be

expected. The normalization requires a sparse matrix multiplication, a vector inner product, and a vector division. Because the dominant computations,

matrix multiplications, parallelize well, good speedup can be expected.

The model for the other computations in LANZ, to(p), includes the cost

of the n-length vector operations. These operations have not been parallelized

in LANZ. The cost of reorthogonalisation is reflected in the second group of

terms in Equation 8. Because the number of reorthogonalizations varies so

widely between problems, this model approximates the cost based on the number

of Lanczos steps and the number of eigenvalues computed. Because sparse

matrix multiplication is the dominant computation in this portion of _o(p),

good speedup can be expected from it.

One application for this model is to ensure that the parallel implementation

of LANZ is performing as expected. For a comparison of the model against

the implementation, LANZ was run on a medium size eigenproblem, r_-12054,

where the ten lowest eigenvalues were sought on an eight processor Cray Y-MP.

An examination of Figure 6 reveals that the implementation times are very close

to the predicted times from the model.

4 Analysis of the Parallel Implementation

When the performance data from LANZ in Figure 6 is transformed into the

speedup curve in Figure 7, a plateau in the speedup curve can be observed as p

exceeds four; there are several reasons for this early plateau.

8
Parameter Description

n order of the eigensystem

p number of processors

13 average bandwidth of the linear system

a average number of non-zeros per row of the linear system

as average block size during the factorization

f,(i) cost of synchronization given i processors

c_p(i, j) cost of j simultaneous/-length vector inner products

c0==(_) cost of sparse saxpy operation with _ non-zeroes

e,,, cost of single multiplication

e_m(i) cost of/-length vector multiplication

e_ cost of single addition

cs®(i, j) cost of/-length j size extended saxpy operation

n, number of Lanczos steps

n_ number of eigenvalues

n/ number of factorizations

t/(p) time for 1 factorization on p processors

i¢_(p) time for 1 forward triangular matrix solution on p processors

tbs(p) time for 1 backward triangular matrix solution on p processors

tram(p) time for 1 sparse matrix multiplication on p processors

tee(p) time for 1 eigenvector calculation on p processors

to(p) time for the other calculations in LANZ

Figure 5: Parameters required by parameterized model
14

12

8

Time in

Seconds

I0

6

Model

Implementation

I I I t I I I I

1 2 3 4 5 6 7 8

Number of Processors

Figure 6: Comparison of the implementation vs. the model on a Cray Y-MP

10
5

Speectu_

4

3

2

1 2 3

4 5 6

Numt3er of Processors

7 8

Figure 7: Speedup curve of LANZ on a Cray Y-MP

11
First, the parallel implementation of LANZ did not sacrifice the vector performance of the code in order to obtainbetterparallelspeedup. Ifthlshad been

done, then the speedup of the parallelimplementation on p processorsover the

parallelversion on one processorwould be excellent;however, the uniprocessor vectorimplementation would be fasterthan the parallelimplementation on

four processors.The multiplesegmented computational unitsin each processor

of the Cray Y-MP are,in essence,usurping much of the parallelisminherent

in LANZ. The plateau in the speedup curves for LANZ would occur laterif

vectorlzationwas not used. To illustratethispoint,parameters for the model

were gathered from the Cray Y-MP with vectorizationand without vectorization.The performance of the vectorizedand unvectorizedversionsis plottedin

Figure 8. From thisfigureitisclearthat the plateauin speedup can be delayed

significantly.This would be foolish,however, because the vectorizedversion is

always fasterthan the unvectorized version. The ultimate goal should be to

minimize execution time,rather than to maximize parallelefficiency.

Second, a comparison ofthe sequentialexecutionprofilein Figure 1 with the

parallelexecutionprofilein Figure 9 of the same computation on four processors

of the Cray Y-MP revealsthat the triangularmatrix solutionsarethe bottleneck

in the parallelversion of LANZ. The percentage of execution time required

for triangularmatrix solveand the sequentialcomputations increasedfrom 16

percent on a singleprocessor to 38 percent on four processors. This is not

unexpected, because it was stated in the analysisof the implementation in

the previous section that triangularmatrix solutiondoes not parallclizewell.

Two reasons for thispoor parallelizationcan be observed in the submodels for

forward and back triangularmatrix solution:i) the low ratioof computation

to synchronization,and 2) the decreasingvectorlengthas p increases,resulting

in a degradation of vectorperformance.

Third, the type of problem being run willbe a significantfactorin determining the speedup. Ifthe problem requiresa largenumber of Lanczos steps

and only one factorizationto converge to the sought aftereigenvalues,then the

triangularmatrix solutionswilldominate the runtime and cause poor speedup.

However, if the eigenvalues in a particular frequency range are being sought,

at leasttwo factorizationswillbe necessaryand, possibly,only a few Lanczos

stepswillbe required,resultingin excellentspeedup. Therefore,the distributionof the eigenvaluesas well as which eigenvaluesare being sought willaffect

the speedup.

5 Improving Parallel Performance

Two approaches for improving the parallelperformance of LANZ willbe described. These approaches improve performance by either: 1) reducing the

number of Lanczos stepstaken and thereforereducing the number of triangular

matrix solutions,or 2) allowing more than one triangularmatrix solutionto

12
70 -

6O

5O

40

Time in

Seconds

30

2O

10

.... l I 1 i , I I I

1 2 3 4 5 6 7

Number of Processors

I

8

Figure8:Estimatedperformanceofvectorizedvs.unvectorizedversions

13
4

3

Speedup

2

with dynamic

shif_

- without dynamic

shifting

I I i I

1 2 3 4

Number of Processors

Figure 11: Speedup curves for LANZ with and without dynamic shifting

16
t of the lowest eigenvalues are required, or 2) all the eigenvalues in a particular

range are required.

If all the eigenvalues in a range are required, then the selection of shifts is

straightforward. A shift on each endpoint of the range is always required to

compute the number of eigenvalues in the range. Additional shifts can then be

selected in the interior of the range, if more than two groups of processors are

available. Because the eigenvalues closest to a shift are found first, the groups

working on the endpoints will calculate some eigenpairs outside the specified

range. Initially, two assumptions will be made: 1) the eigenvalues are linearly

distributed and, 2) each group will compute an equal number of eigenvalues

around its shift. Given these assumptions, the initial shifts, _k, will be chosen

such that each group will be assigned an equal amount of work. Given the

endpoints of an interval, a and b, and the number of groups, rn, the equations

for computing the rn/_'s are (if rn > 1):

_, = ,, (9)

/_k a + (k b - a

= - 1) _2-y_1 lot 1 <k <r_

/-trn = b.

If rn -----1, then the shifts can be chosen in one of two ways: 1) two shifts will be

chosen, a and b, or 2) three shifts will be chosen, a, b and the midpoint. The

second method will result in fewer Lanczos steps than the first method, but will

require an extra factorization. Because the cost of factorization and the cost of

a Lanczos step depends on the problem, the selection of a shift selection method

is problem dependent.

In this section, the focus is on the range case, however, a brief description

of the algorithm for selecting shifts when the lowest t eigenvalues are sought is

now given. Because only a star_ing point, 0, for the search is available, a group

of processors will start with a shift of 0. As estimated eigenvalues are generated

by the first group, more shifts can be selected based on the estimated eigenvalues. The shifts would be generated as soon as enough estimated eigenvalues

are available to give a good indication of where shifts will be useful. Other

groups of processors will begin work on these newly generated shifts. A tradeoff

exists between putting the idle groups of processors to work as early as possible

and making sure that the newly generated shifts are in areas of the eigenvalue

spectrum where useful work is done. In addition, it would be desirable to have

as many processors as possible work on the first factorization and then let a

subset of these processors work on the computation of the Lanczos steps at this

first shift. This algorithm is not as efficient as the algorithm for the range case,

because not all the processors can begin working at the same time. For the

group approach to be effective, a fairly large number of eigenvalues must be

sought, otherwise, one group could quickly find the desired eigenvalues, s.

Sin this case, the normal parallel I, ANZ program is efficient becatrse few Lanczos steps

17
Each group computes eigenpairs until:

a) it is finished with the sub-interval between its shift and the shift

immediately to its right and it is finished with the sub-interval

between its shift and the shift immediately to its left, or

b) the number of allowable steps is exceeded

c) the storage capacity of the group is exceeded

A group is finished with a sub-interval if:

d) all the eigenpairs in that sub-interval have been found, or

e) the group finds a converged eigenpair that is on the other side

of an eigenpair found by another c=rrently working group.

Figure 12: Algorithm for group approach

The algorithm for a group is given in Figure 12. Condition e allows the

groups to negotiate the portion of the sub-interval that each will compute. This

has a significant advantage over assigning half of the sub-interval to each group

when the eigenvalues in the sub-interval are concentrated near one of the shifts.

If each group is assigned half of the sub-interval, then the groups will take a

very unequal number of steps, resulting in poor load balancing.

The data that need to be communicated between the groups are: 1) the

shifts, 2) the inertia count at each shift, 3) the eigenpairs, and 4) the left and

right endpoints of the range of eigenvalues computed by each group. This data

is conceptually communicated via message-passing, however, shared variables

protected by critical sections are used in the Force implementation.

To make good choices for the size and number of groups, as well as to

predict the performance of the group approach, the model described in Section 3

will now be extended. The model for predicting the performance of the group

approach is

T(p) = max (:g,(1, p,, n,,, _,,)), (10)

i=l,ra

where the parameters are specified in Figure 13. Of course, the optimal number

of groups, rrh and processors per groups, i_, given p cannot be known = priori for

most problems. However, the model can be used to make general statements

regarding the selection of m and p_. From the model in Section 3 and the

analyses in Section 4, it can be concluded that the parallel efficiency of a group,

gi, will deteriorate as Pi and n,_ 9 increase. However, although the efficiency

are needed.

cLn general, em n. i increaaes, n. i increaaes.

18
Parameter Description

t_,(_, b,c, d) for group i, the Tv from equation 1 where,

a is the number of factorizations,

b is the number of processors,

c is the number of Lanczos steps,

d is the number of eigenvalues found,

rn the number of groups

pi the number of processors in group i

n0_ the number of Lanczos steps taken by group i

n_ the number of eigenpairs calculated by group i

Figure 13: Parameters forgroup model

decreases as Pi increases, the execution time does decrease) ° In general, as the

the number of groups, rn, increases, the number of Lanczos steps per group, n,_,

decreases, xl The tradeoff, with fixed p, between the number of groups and the

number of processors per group can be essentially characterized as trading off

decreased factorization time for more steps per group, e.g. as Pi increases, and

m decreases, the factorization time per group decreases, however, because the

group is now responsible for computing more eigenvalues, more Lanczos steps are

required. To illustrate this point more effectively, the time for factorization and

the number of Lanczos steps have been plotted against the number of processors

in group i. In the graph shown in Figure 1412, p is assumed to be held constant,

and, therefore, m is implicitly decreasing as p_ increases.

The group approach as described in this section has been implemented on

the Cray Y-MP using the Force parallel language. Due to constraints imposed

by the Force, however, each group consists of only one processor. Although

this constraint does not allow full use of the group approach, results from this

implementation, combined with the model and results from the parallel implementation in Section 2, allow very accurate assessments of the performance of

the group approach. To make these estimates, n,_ and n,, for group i in the

group implementation were used to calculate the time for group i given different

values of PiFirst, performance results from the implementation on the Cray Y-MP will

be given. The group approach using 1... 4 processors was run for three separate

situations on the same problem13: 1) all the eigenpairs in a range that contained

1°The decrease continues until a plateau is reached. After this plateau is attained, the

addition of more processors will only increase execution time.

11However, this is very dependent on the distribution of the eigcnvaluea. Also, in general

_.m. nj. is not constant for different m in the same problem.

=L •

l_No units for the y-axis axe shown because the interest is only in the direction of the curves.

lathe order of the cigenproblem was 1824 and the average seml-bandwidth was 127. The

19
Factorization

Time

Number of

Lanczos Steps

Per Group

Number of Processors Per Group

F_gure 14: The tradeofl" be_.ween factorization time and the number of Lanczos

steps

2O
20 eigenpairs were sought, 2) all the eigenpairs in a range that contained 40

eigenpairs were sought, and 3) all the eigenpairs in a range that contained 60

eigenpairs were sought. In cases 1 and 2, good speedup was achieved. However,

a dropoff in speedup was seen in case 3 when m was increased from three to

four because the majority of the eigenvalues were concentrated in the far left

end of the range. This concentration of eigenvalues resulted in very poor load

balancing. Speedup curves for all three cases are plotted in Figure 15.14

To give the reader a better understanding of how the group algorithm operates, a short discussion of the division of work among the groups for a specific

example is now given. First, two aspects of the algorithm that should be noted

are: 1) because eigenvalues closest to a shift are found first, the groups working

on the outermost shifts will likely find eigenvalues outside the desired range,

and 2) two groups may compute the same eigenpalr. The division of the eigenvalue spectrum among the processors for case one is is shown in Figure 16. An

examination of the division of work using two groups in the figure reveals that

groups one and two computed several eigenvalues outside the desired range.

Also, group one computed one eigenvalue in the range of group two and group

two computed one eigenvalue in the range of group one, therefore, although

group one found 21 eigenvalues, only 20 eigenvalue are seen in its range.

To assess the performance of the group algorithm with different values of

p_ and m, results from the implementation in this section are combined with

results from the implementation in the previous section. The runtime of the

group algorithm is the runtime of the group that takes the most Lanczos steps.

If all the groups are of equal size, then the number of Lanczos steps that each

group takes is independent of p_. Let

r_m = max n,,. (11)

The implementation described in this section was run on a problem 16 and r_

was observed for m -- 1... 4. Then, the implementation described in the previous section was run on the same problem for rim, m = 1.-. 4, steps on p_,

i -- 1...4, processors. The execution time observed for each of these runs

would be the execution time of the group algorithm. The speedup curves for

m = 1... 4 groups are shown in Figure 17. From this graph, it can be concluded

that, in general, the group algorithm is superior to the algorithm described in

the previous section. It can also be concluded that, in general, the best performance of the group algorithm will not be attained with m = p.

size of the problem was restricted due to the limited stack size of the Cray Y-MP.

14_ re=l, then the group approach is equivalent to the parallel implementation described

in Section 2.

tSThe three cases are described in the prcvlous paragraph.

t6A vibration problem with n--12054 and an average bandwidth of 328. All the eigenpairs

in a specified range containing twenty eigcnpairs were sought.

21
4

3

Speedup

2

Case 2

Case 1

X I 1

1 2 3

Number of Groups

Case 3)

1

4

F_gure 15: The performance of the group implementation

22
Run 1 One Group Run 2

46 steps 23 steps

24 _.'s 13 _.'s

._ h,,(,, , , ,._ ,

t _2. _ t _ , r L|Ll: I i l J _ I

a /'\ b

! !

)

Two Groups

Group 1 Group 2

39 steps 38 steps

21 X's 16 _.'s /

_,,,,,,,,,,._,,,,,,, ,,_, ,,,, , , ,._ ,

L'"'"'"'-_' """ ''JL' '"' ' ' "_ '

\ a /\ b

t J

, I

)

Three Groups

Group 1 Group 2 Group 3

32 steps 30 steps 22 steps

14 ;L's 11 _'s 10 _.'s

,,,,, , , I I I ,, ,, ,,_ ,

a b

Four Groups

Group 1 Group 2 Group 3 Group 4

23 steps 22 steps 20 steps 18 steps

8 _.'s 10 _,'s 7 _.'s 7 _.'s

(,,,,_., h,,(.,, ,v,, ,h 4, .,, h,(, ,._ ,h

a b

a is the left end of the range

b is the right end of the range

X denotes a shift

23

F_.g'_re '.6: T_e _.iv_on oi" work ."or case one
1 2 -"

10

8

Speedup

6

4

2

D

I I i ! I I I I

2 4 6 8 10 12 14 16

Number of Processors

Figure 17: The estimated performance of the group approach

24
6 Concluding Remarks

The use of the Force parallel programming language on a large, practical code

revealed some shortcomings in the language and several remedies were proposed.

The parallelization of a Lanczos-based solver for the generalized eigenproblem

was discussed. An analytical performance model for this code was used to

explain the tradeoffs that were made when implementing the eigensolver on

a parallel-vector computer. The analytical model also was used to show that

a bottleneck in the parallel eigensolver is the forward and back substitution

algorithms. Two algorithms for improving the parallel performance were given.

The parallel eigensolver presented in this report is suitable for parallel computers with a moderate number of processors as well as for parallel-vector computers with a small number of processors. However, due to the bottleneck imposed by the forward and back substitution algorithms, this eigensolver is not

suitable for large scale parallelism. To exploit large-scale parallelism, the authors are seeking to eliminate this bottleneck by pursuing research into iterative

methods for solving (K - aM)_. = y.

References

[1] M. T. Jones and M. L. Patrick, "The Use of Lanczos's Method to Solve the

Large Generalized Symmetric Definite Eigenvalue Problem," Technical Report 89-67, Institute for Computer Applications in Science and Engineering

(ICASE), NASA Langley Research Center, Hampton, VA, 1989.

[2] M. T. Jones, The Use o] Lanczos' Method to Solve the Generalized Ei#enproblem. PhD thesis, Department of Computer Science, Duke University,

1990.

[3] R. G. Grimes, J. G. Lewis, and H. D. Simon, _The Implementation of a

Block Shifted and Inverted Lanczos Algorithm for Eigenvalue Problems in

Structural Engineering," ETA-TR-39, Boeing Computer Servies, Seattle,

WA, August, 1986.

[4] M. T. Jones and M. L. Patrick, "Bunch-Kaufman Factorization for Real

Symmetric Indefinite Banded Matrices," Technical Report 89-37, Institute

for Computer Applications in Science and Engineering (ICASE), NASA Langley Research Center, Hampton, VA, 1989.

[5] M. T. Jones and M. L. Patrick, "Factoring Symmetric Indefinite Matrices

on High-Performance Architectures," Technical Report 90-8, Institute for

Computer Applications in Science and Engineering (ICASE), NASA Langley

Research Center, Hampton, VA, 1990.

[6] B. N. Parlett and D. S. Scott, "Lanczos Algorithm With Selective Orthogonalization," Mathematics of Computation, vol. 33, pp. 217-238, 1979.

25
[7] H. Jordan, "The Force," Computer Systems Design Group, University of

Colorado, 1987.

[81 M. T. Jones, M. L. Patrick, and R. G. Vo_gt, "Language Comparison for Scientific Computing on MIMD Architectures," Report No. 89-6, Institute for

Computer Applications in Science and Engineering (ICASE), NASA Langley

Research Center, Hampton, VA, 1989.

[9] E. Poole, "Personal communication," 1990.

26


Report Documentation Page

1. Report No.

NASA CR-182072

ICASE Report No. 90-48

4. Title and Subtitle

2. Government Accession No.

THE USE OF LANCZO'S METHOD TO SOLVE THE LARGE

GENERALIZED SYMMETRIC EIGENVALUE PROBLEM IN

PARALLEL

7. Author(s)

Mark T. Jones

Merrell L. Patrick

9. Performing Organization Name and Address

Institute for Compuer Applications in Science

and Engineering

Mail Stop 132C, NASA Langley Research Center

Hampton, VA 23665-5225

12. Sponsoring Agency Name and Add_ss

National Aeronautics and Space Administration

Langley Research Center

Hampton, VA 23665-5225

15. Supplementary Notes

Langley Technical Monitor:

Richard W. Barnwell

3. Recipient's Catalog No.

5. Repo_ Date

July 1990

6. Performing Organization Code

8. Performing Organization RepoR No.

90-48

10. Work Unit No.

505-90-21-01

11. Contract or Grant No.

NASI-18605

13. Ty_ ofRe_and Period Covered

Contractor Report

14. Sponsoring ,_,gency Code

Final Report

16. Abstract

The generalized eigenvalue problem, KX = _Mx, is of significant practical

importance, especially in structural engineering where it arises as the vibration

and buckling problems. New software, LANZ, based on Lanczo's method has been developed for solving these problems and uns on SUN 3, SUN 4, Convex C-220, Cray 2, and

Cray Y-MP systems.

Preliminary results of using the Force to obtain a multiprocessor implementation of LANZ on MIMD parallel/vector systems are reproted here. A parallel

execution time model of LANZ is defined and used to predict the performance of LANZ

as well as examine hypothetical modifications to LANZ. The results of using dynamic shifting to improve parallelism are presented. Finally, the results of assigning

a group of processors to separate shifts and finding all the desired eigenvalues

using LANZ in parallel are reported.

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

generalized eigenvalue problem,

Lanczo's method, parallel

19. SecuriW Cla_if. (of this repot)

Unclassified

NASA FORM 1626 OCT 86

18. Distribution Statement

61 - Computer Programming and Software

64 - Numerical Analysis

Unclassified - Unlimited

20. SecuriW Cla_if, (of this page)

Unclassified

21. No. of pages

28

22. Price

A0 3

NASA-L_ngleY, 19eM}