ML19305E704

Modified Subspace Computation for Eigenvalues & Eigenvectors.

ML19305E704
Person / Time
Site:	Sequoyah
Issue date:	06/30/1978
From:	Jesse Rollins, Renee Taylor PMB SYSTEMS ENGINEERING, INC.
To:
Shared Package
ML19305E692	List: ML19305E691 ML19305E694 ML19305E704
References
IEB-79-14, NUDOCS 8005200293
Download: ML19305E704 (18)

Text

.

a05260.d93 i

..h I

i MODIFIED SUBSPACE C0tIPUTATI0ti .

FOR EIGEllVALUES AfiD EIGEtiVECTORS l

l l

l by -

Robert L. Taylor and James M. Rollins -

1 l

l Consulting Report 78-1 -

.==. . . . _

j . .

I Prepared for .

PMB Systems Engineering Inc.

500 Sansome Street San Francisco, California 94111 y

i i: ,

?;;. -

.=

p l

[~,dhlab_____..__.___..__. . April 1978 (!bd. June 1978) . _ _ _ . , _ . .

' ++- - .

., .u. , , , , l

,w * .  !

1. . - = .'=*b.,.,..,

. - . 2: - : - : -- - .

l Table of Contents

~

=

. b:~ . Pace Ta bl e o f Co n t e n ts . . . . . . . . . . . . . . . . . . . . . . . . . i .

1. Introduction . . . . . ....................1

2. Isolation of Subgroups by Bisection . . . . . . . . . . . . . 4

3. Subspace Iteration on a Subgroup . . . . . . . . . . . . . . . 6

. 4. Selection of Start Vectors . . . . . . . . . . . . . . . . . . 8

5. Modification of Convergence Test . . . . . . . . . . . . . . . 10

6. Removal of " Ghost" Vectors . . . . . . . . . . . . . . . . . . 11 References ...........................12 Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 e

e e

e

(=  :

e 0

=*

Y

. 'W

1. F

$ 9

+ -.

, - 4 ... .

4

_. __ _ , .__...: -. .7 g

.g , ,

a . ' ~ ~ . . - . -~.. , , - * *  % . . . . . . . . ,

X "~ '_ t - u j -

... -p.

' fg '

g m g e "EO'1 . g e a y l we* -

amm- ene .m eme. *e *===ame e . - . ' _ . . ' N, -. . , . s. .

_ ... --,- ~ ~ - - - -

. , _ , . . . -. . . . . . . - - - . . - - ~ - - - - - - - - - - . . . - . . - - - - - - - .. -. . - - -

1. Introduction IL The calculation of eigenpairs for large structural systems often consumes a large amount of , total computer time. Much effort has been

~

expended on procedures to reliably extract selected eige'npairs at a reasonable cost.

One popular procedure is the subspace method as given by Wilson and Bathe [1] and implemented in the SAPIV structural. analysis system [2]. '

The general eigenproblem in the SAPIV program is defined by Av = MyA (1) where A and M are the symetric stiffness and " lumped" mass matrices (i.e., M_ is diagonal), y_ is the set of eigenvectors which are mass ort'honormal, and A_ is the array (diagonal) of eigenvalues. The diagonal arrays are stored internally in vectors.

kn a typical problem the stiffness and mass arrays are o

~

~

size, and it is required to find the lowest set of eigenvalues. Often the number of eigenvalues 4 required is much smaller than the total number ..

and it is in this regard that the subspace method is an extremely viable algorithm for extracting the required eigenvalues and their associated vect.or.

The basic subspace algorithm can be sumarized as follows:

(1) T Perform an LD_L decomposition o A. 7 For the k th (2) iteration, set '

-- }i .,

. ,,, d; ;.f~-

and solve the inverse iteration steps - '

u (jij _ z_ '

Ax_N+I = N-\ . . - _ _ - . .

l. k+1 ._;-

~

for E

- l2' , 'S

- - - ,. . . . . . , - / -- _

(3) Construct the subspaces is = k y . _A* (_x+1)T _A _x_k+1 and ,

M* = k k ~

(x . +1 )'T M x +1 .

(4) Solve the reduced general eigenproblem -

. _A* R = h* R 3*

48 k

. (5) Compute y_k+1, x +1 R , then repeat (2) to (5) until " converged."

In use the set of vectors y_ consists of the p-vectors corresponding to the values of the p-lowest eigenvectors desired, plus some extra vectors which are of ten called " guard vectors." Let q be the size of the total number of vectors retained. Further, let n be the size of A and H, and m the half-band of A_. ,

Thus, the operation counts for each step in the iteration process are:

2 (1) fnm operationstoconstructL_0[ decomposition.

~

k nq operations to get y and2mnqoperationstoget[*I (foniard solution and backsubstitution).

2 (3) fnq + nq operations to get M*, and2 fnq + 2nmq operations to get A*. -

(4) Solution of the general eigenproblem is:

C(i)i q', where C(i) and i are iteration, parameters in the

, double Jacobi algorithm and depend on the threshhold and number of iterations to get convergence. c. ;

a 2

(5) nq multiplications to get y,k+1 , [;

a .

,.y.

es=- - Usually convergence will be achieved in a' few iterations; however, as the number of vectors increases, the number of iterations will also increase.

.< S I

j W g gie WN OFM@M *~

• • w- .==-4 m.. . o wm ep =

. . 3 This is primarily due to the fact that convergence of the subspace is dominated by the ratio of the eigenvalue just outside the subspace (A qq in Fig. 1) to the largest eigenvalue desired (A i p in Fig.1). As I the cize of the subspace inc~r eases, this ratio reduces (unless an 1

' extremely large set of guard vectors is retained). Another related i

I disadvantage is that the lower values converge before-the-larger-values, thus much wasted effort is expended to get the larger values.

~

It is clear from the operation counts that situations for which q 2 m become inefficient in the sense that subspace iterations are  !

expensive relative to the triangular decompositions. Furthermore, as q increases the cost of the Jacobi solution goes up extensively - the number of iterations to diagonalize A* and !!* are' of order 10-20 for some large problems, C(i) is a constant depending on the threshhold

,.. algorithm used in the Jacobi routine. If q > 0.1 n, the .lacobi iteration itself can exceed a factorization cost. '

. From the above discu[sion, it is clear that as the number of sought eigenvalues increases' the sub: pace becomes prohibita51y more and more costly and some modifications are required. In thi!. report we discuss one such modification to the algorithm provided in sap'IV. The modifica-tion consists of dividing the . subspace into small g oups and then using l a shifted subspace algorithm to get the eigenvalues and vectors in each small group. A binary search roatine using the Storm sequence property I

.is used to determine the bounds on the groups. Onte the groups are isolated, the subspace is used with shifts to determine"the eigenvalues w.

, and vectors in each small group. In this regard, special; considerations ph=

are necessary to select the Sitial subspace subgroup vectors in order to n jt . .

achieve convergence. Some ncw difficulties are also encountered in that

n. -

i

m 1

" ghost" vectors occasionally appear in the interval. A strategy is

([= employed to detect and remove these vectors. Finally, consideration of the I/O is necessary since, in general, many more iterations are necessary to Det all the vectors. We discuss the aspects of our modi-fications in more detail in the following sections. -

2. Isolation of Subcroups b Bisection - --- -

- - - - - -~

In order to establish the subgroups of size s r (i.e., we limit the number of vectors in each subspace to be equal to or less than r), a binary bisection algorithm is employed. The value of r is determined as r = max (4,i[!42BAND

.2 _. .

+ 1) , -

l r = min (r, core space limitation) ,

.-... where MBAflD = m the semi (half) band of A. The chosen value is to balance costs of a subspace iteration and the triangular decomposition. For very 1arge (n) problems, it may b'e necessary to further reduce r to utilize I/O more efficiently or because of core space limitations.

t s

' Once the value of r is determined, it is necessary to know how many eigenvalues are required and estimate the largest eigenvalues wanted (upper bound estimate) or a cut-off frequency (i.e., a frequency'above I which no eigenvalues are wanted). If no cut-off (upper bound) frequency

( is specified, the program will use the Gershgorin disk's center value

( h;$ ) '

S = CUTOFF = 2 max 4 (Hii ). 2' i = 1. n  ? ...

' Only values with nonzero M$$ are used. (fi.B. this simple algorithm does not 5.5.E

~

sufficeifaconsistentmassisused).

• l i, ,

--s 0

l .,

~

5 Starting from the cut-off value, a Sturm check is made using 'the

!?

triangular decomposition of the shifted stiffness. Accordingly, we compute

=

C A - SM =

EDET ,

and count the number of sign changes in D (thisgivestheSturmcheck).

We should note that C is not positive definite and that we do not use pivots or double precision in the decomposition. Consequently, there is '

a remote possibility of the decomposition failing (none has appeared to occur at this. time). From the Sturm property, we now know how many values are below S.

We record the shift and the number of values by storing in AT(tiVAL) =

S .

where ilVAL is number of values less than S. This records the necessary data to recover the shift and Sturm property (AT is in comon /EM/ and currently is dimensioned 500; it can be increased to a much larger value withoutincreasingcoredema'nds).

We next halve the interval which has the largest number 4 of eigenvalues; thus, in the second step we set S-(i.e., SHIFT)toS/2andrepeattheprocess.

We continue bisection until

. we isoi.ste all the subgroups, i.e., see Fig. 2 [3]. Each subgroup has less than or exactly r-values between recorded shifts. It is possible  !

that two adjacent subgroups combined have less than r-values. Consequently, l we regroup these into a single gre :: recording a final group configuration,

.which is then output to the line printer. At this time we are prepared i l

to start the subspace computations for each subgroup in hurn. We begin with the lower values and work to the upper ones. .

5 a-5 .

=-

6

. 3. Subspace Iteration on a Subaroup r'

After isolating each subgroup, the subspace alscrithm is used to determine the eigenvalues and vectors in the subgroup. The algorithm given above is modified in two basic ways: -(a)weuseonly ~ d 2r vectors in each subspace - r values serve as guard vectors and (b) we shift to the mid-interval as evaluated by the bisection bounds (see Fig.3). We can summarize the steps for a typical subgroup as: -

~

(1) Compute 't.,D,[ from A - og = C, where a is the shift to the mid-interval.

(2) For the'kth iteration set, k k y,g and solve (A_ - aM ) [+I = [ ,

k where now 'v are the subspace vectors for the current subgroup; i.e. , 4 x n.

(3) Constrgct the subspaces '

k+1

.- .P " (Ek+1)T A,1 and -

i,4 (/+1)T g y+1 which.are now 4 x 4 arrays (symmet'ric). - '

. (4) Solve the reduced general eigenproblem:

h (e + =r- ) n - c a r k '

(5) ComputeU+I=x+1 gandrepeat($)to,(5)untilA*+A, (i.e., values in interval converge). -

(6) Repeat (1)to(5)forallsubgroups. <-

a Some further coments are necessary on this algorithm. These deal I

~

i with how to select " good" start vectors, which will then control conver-1 . m= -

t r gence, how to deal with any " adverse" vectors which appear (we call these l

l

, ) -

7

" ghost" vectors), how to utilize core and I/O effectively, and last, some f coments on computed A* and 5*. He begin with the last two, namely, stor-age I/O and computing 5* and B*.

We do not want to compu'te X* as shown above. The reason is that we will need two copies of an array of size i in core and a block of A or we will have to make a large number of I/O operations with A; both are cbjec-tionable. If we look at the original implementation (i.e., see Fig. 4a),

this is precisely what was done - two copies VL and VR and a block of A.

If q is large, the size of a block (HEQB) is controlled by the subspace computation of A_* and B_*. We now have small size subgroups and will store all of 1 in core plus a block of ,A_(_C_) and all the mass. In this way we

, will reduce I/O to just A_, being transferred by blocks once f6mard and once backward. To get all of 1in core restricts the s' ubspace size further

.= for large problems where I/O would begin to dominate costs. Thus, the I

modified core utilization for subspace computation of A* and B* is as shown in Fig. 4b. We note particularly that v_ appears on left and right k

side of equations 4 - no second copy being necessary. This implies that i ,

y k, and all are stored in the same array. This is a substantial change from that originally given in Fig. 4a. .

~ In order to use this storage scheme, we must look further at how X*

canbecomputed;5*canbecomputedsimplys,inceHisdiagonal. Using the definitions for C and A*, we note that 9

~

k* ( ) E( ) g. , -

W

~. .

_ .... - ~ . . . . ~ . .-.

k+1)T T k g +1) -E ,

.. g . - .

Now let us look at the solution for x ,

S

.a

\ . .

, . . . . . . _ E . w w- a .-e = ^- * " ' ' ~ ~ * - ~ ' ' ' " ' ~~

. 8 C, ([d ) = [ = [ D_ [ T ( C)

('

which we compute by a foniard solution

_Lh*

, - ~ '

3_ .

and the backsubstitution 2

.L_ (_x_k+1)=D-I c

We note now that z_ can be used to compute X* as follows. Substitute for

[x I the D-I z_ and obtain

\

zT pD)g2

~ _

A* =

z is stored in y also. -

.r - This is identical (with exception of reciprocal diagonal n.itrix) to the algorithm fo'r getting 5_*. It can be computed as we compute the forward <

, i solution to the equations by blocks since D-I is diagonal (othemise we .

would need more 1 than one block at a time) and all z is in core. The comple-tion of backsubstitution gives x.k+1 and we can then compute 5* again easily since M is~ diagonal. These revisions are incorporated into the rewritten subroutine REDBAK. Core revision is necessary in SAP 4, MODES, SSPACEB,

-t

• i where NEQB is determined together with r.

4. Selection of Start Vectors ,

] The selection of ~v for each subgroup will affect (s.ignificantly) the subspace computations in each subgroup.

If a v is sele 3ted which.has a

, very poor content of the eigenvectors in the subgroup, convergence may

! _ ' .... never occur in a reasonable number of iterations (e.g., s.ay 16). Several

[8

9 types of vectors have been tried.

These include random vectors and selected unit vectors similar to that originally used. The vectors were mass orthogonalized in some calculations using a Gram-Schmidt process (G-S).

Random vectors did not perform reliably - convergence not being achieved or eigenvalues were missed. The later subgroups' especially suffered.

Since very few vectors are used, the unit vectors were often not rich .

enough" to get all the eigenvalues. The strategy employed was to use '

.the largest values of

• "ii

'i (Ag j - aMjj) ,

as the unit vectors (note they are orthogonal). In the first subgroup H_

is the first J_. -

In performing the calculations, the vectors of the guard values which were above the p'revious s'ubspace subgroup were used also.

~

When fewer than r values were used, there were often cases in which no_ values were above.

By using r values we usually have at least one being broupht forward.- This definitely improves the behavior - especially if the unit vectors are mass I

orthogonalized to these vectors. The start vectors currently used consist of the vectors above the previous subgroup, with the remaining vectors being made up from all the elements of . .

1 v

N '

_ ii .

(Agg - aM$$)

~

c..w g

divided equally into the y" as shown in Fig. 5. These vectors are all made mass orthonormal by the G-S process.

We note that they remain mass ortho-

- (E== normal (but not stiffness) for all subsequent iterations; thus, the only

10 loss in orthogonality which can ccme in is that between subgroups.

( The start vectors in each subgroup are also orthogonalized to the previous converged vectors to ensure that they are Ritz vectors. In theory this should preclude'the possibility of ever having " ghost" vectors.

In practice, however, we usually compute the eigenpairs to less than

- - ~ - - - - - -

-~

. machine precision. Thus, it is possible-for-the-vectors-to-still-contain . -

small portions of all the eigenvectors. In setting the start vectors ,,

~

these "small errors" may eventually result in start vectors which have .

significant components of eigenvectors below the subgroup and these lead to " ghost" vectors., A " ghost" vector is merely a vector which-is a com-- '

~

bination of vectors whose eigenvalues are both above and below the subgroup.

The combination gives a value in the subspace (not an eigenvalue) which corrupts the result. The effect of these vectors can be minimized by

,_- computing eigenvalues to a precision of at least half the machine word length. In the.next two sections, we discuss the procedure adopted to detect and remove the " ghost" vector from the calculations. The setting

of the start vectors is performed in SI:;VEC; ORTH 0G.is used to make them N

mass orthogonal and NORM makes them mass orthonormal.

5. Modification of Convergence Test To reliably detect good and bad vectors in a sLbgroup, we have modified the convergence test in EIGSOL. We define t'he to)erance as follows: ,

~

RT01.V

= min j

[(ABSD(i)(DL(j)-D(i))\q, j y., ,

=

This ensures that " ghost" vectors moving through the group do not affect a

, determination of good vectors. In practice we then find that " ghost" vec-

~~ 3g-.

tors produce large RTOLY compared with good vectors.

3 , .

- ~

6. F.e= val of " Ghost" Vectors ,

In some subgroups we have noted that more eigenvalues exist (between the boundaries defined by bi,section) after a few iterations than there should be in the interval. If iteration continues, these' tend to converge to an eigenvalue already in the interval and remain. In some instances, they may eventually " drift" to the next value and thei'next till they pass out of the interval. We call these " ghost" vectors. If we did not know -

better from the Sturm property during bisection, we might think they were.

eigenvalues - they are not! They must be removed and the sooner the better. More or .less by. trial and error, we have found that by the end of four subspace iterations in a subgroup we c$n detect the presence of bad copies quite reliably. If full print out is used, this appears at the end of each iteration as the values of L and NF. The value of L is l'= the number of values currently between the boundaries t ' the subgroup, and NF is the numberi which shoul,d be in the interval.

If at the end of four iterations we detect ghost vectors, we delete the vector in the interval which has the largest change in the eigenvalue from the previous iteration.

If after two more iterations we still have a ghost, another deletion is made. This will reduce the number of guard vectors but, to date, has , l worked well. The routine SWEEP is called from EIGSOL to delete a vector.

Y

.g S

p . *-

I

.?

~ - :Y_

,~

.'* , , . .-?- , - ~ ~

. . .A ~:.- ~ .

~

e

,. w *h. g a. A ...gpq vm -

s '49*'

eAM*SD#N%T*****"*"~*****"#E

e - - -

12

([~ _P,eferences

1. K. J. Bathe and E. L. Oilson,tiumerical Methods in Finite Element Analysis, Prentice-Hall, Englewood Cliffs, llew Jersey,1976.

2. E. L. Wilson, K. J. Bathe, and F. E. Peterson, "SAPIV - A Structural Analysis Program for Static and Dynamic Response of Linear Systems,"

EERC Report 73-11, University of California, Berkeley,1973 (revised 1974).

3. A. Jennings and T.J.A. Agar, " Hybrid Sturm Sequence and Simultaneous i -

Iteration Method," Symp. on Appl. of Comp. Meth in Engr., University of Southern California, Los Angeles,1977.

P

.~

J

• q -

e

. . 1 t -

1 i

{ -a

! . j l .

ie W ..

\'

~

• fw \

l l l l . . .=. . .-

l J. . gy. . ,

. .- ..c. _ .

._ ~ .. ._

. - - -. - ,. - . _ u . ... ..

13 p- .

I l

l i

. i

.A .. .

1 l

Eigenvalues to " Guard" .

-< A-  : 4 -

V be computed Values

}  :  :  :  :  : y --.

A A A p-1Ap Ap+1 1 2

• *

• A A q Aq+1 FIG. 1. SUBSPACEEIGEilVAlbES

.=-

l 's Eigenvalues to

~

4 be computed

' J' >=

Subgroup 1 Subgroup L -

y-  ;

A A w '

2 x <

1 p  ;

.-y -

FIG. 2. SUBGROUPS FOR SUBSPACE

. s . .

T.: ..: .

$. 'T

...,..;em,.. ..

~ , . , . .. -.. . _ . , -

_s . ,

e .

s -

=~i. . ~ . - - . . . ~ U ~. $. . "' -

._ . f.'iCEi .

.~:  ;. . . - .

l D b .

%e e

S .m 4 -

mag e-g 6.

. - %;..,53 c --

)

I.... l 1

. l

. 1 i

l l

Shift a

-)

\

Upper Limit From Bisection n

. . . - _ w_ sj .

, Lower Limit From Typical Sobgroup k Guard Values Bisection p Jg .

~ -

( - -

A

j k -

Ak +r - A

(~~~ .

. S e

b

.4

. 8

' FIG. 3. TYPICAL SUBGP.0UP IN SUBSPACE

'i e

.6 1, .

Y .

e f s k..~; ~-

..g.... . .

g . ...e .++.m. ..

, , , , RE ,

mw e.

f

_. O * * '== "

I -- - - -

15 (G

.-_ 9_ ,

_ 4 -

XM A VL == VR

- \ .. _ _ , ,_ ..

(a) Original " Blocking" Scheme .

5 4

~ .

~

. l

""~

I XH (b) Modified " Blocking" Scheme Information In-Core ,

}:.

Information in Backing Store '3 e

@ FIG. 4. ARRAYS I!! HIGH SPEED CORE FOR StfBSPACE COMPUTATIO l

16 7, -

' .C

. m m L .M

• =

L O 2 l's o e .

M U h

',/

u c 3>

C2 &

• r=

L '

g .. - .

O c) E 4 Cs L ~ ~~

O '  % %O

• O L CJ C .

> c:s > 0 c C .C M M c OM L m L3 L O

• J C) 4

• C m  % U n O L CJ "U .O.

r O >

CJ M . h.

CL > O. M E CJ D *O L L O CJ 't3

(.;: D L L- mM a

w :o m L -

' /[. _ '

o o

'2 #o

- *1 Ek 1 hk l l

. 1

'l *

- tion-Zero Region Usino - -

ii (Ag j - aligg) d e

, FIG. 5. START VECTORS FOR A SUBGIt0UP

(:';  ;

. s 5

(===.

~:

N e n

,} -