Minicomputers and Large Scale Computations

2 Conjugate Gradient M e t h o d s for Solving

Downloaded by EAST CAROLINA UNIV on April 11, 2018 | https://pubs.acs.org Publication Date: June 1, 1977 | doi: 10.1021/bk-1977-0057.ch002

Algebraic Eigenproblems J. C. NASH Research Division, Economics Branch, Agriculture Canada, OttawaK1A0C5,Canada S. G. NASH Mathematics Department, University of Alberta, Edmonton, Alberta, Canada The problem which is addressed in this paper is that of find ing eigensolutions (e,x) which satisfy Ax=e Β

x

(1)

where A and Β are real symmetric matrices and Β is positive defin ite. The particular case where Β is the identity is also of inter est. This study will focus on methods using the conjugate gradient algorithm because this is very frugal of storage when used to solve linear equations [1] or function minimization problems [2]. The algorithm can be used in either fashion to solve the eigen problem (1): 1) By solution of the equations (A - kB) y = Β

x

i

(2)

i

with x = y / || y || (3) i+1

i

i

where the double vertical lines mean "the norm of", it is possible to find the eigenvector x (or y) having eigenvalue e closest to the shift k. This is the process of inverse iteration [3,4]· 2) By minimization of the Rayleigh quotient R

T

=

x

A

x

/

T

x

B

x

(4)

T

with respect to x, where the denotes transposition, the eigen vector x corresponding to the most negative eigenvalue e is found [5]· R takes on the value e at its minimum. Note that if Β is the identity, minimization of 1

R

=

T

x

(A

-

2

k1)

x

/

T

x

x

24 Lykos; Minicomputers and Large Scale Computations ACS Symposium Series; American Chemical Society: Washington, DC, 1977.

(5)

2.

NASH

AND

25

Conjugate Gradient Methods

NASH

g i v e s the e i g e n v e c t o r c o r r e s p o n d i n g to the e i g e n v a l u e w h i c h i s c l o s e s t t o t h e s h i f t k. A l t e r n a t i v e l y , r o o t - s h i f t i n g o r o r t h o g o n a l i z a t i o n a s d i s c u s s e d by S h a v i t t et^ a]_. 16] may be u s e d t o f i n d some o f t h e h i g h e r e i g e n s o l u t i o n s .


The a p p r o a c h e s a b o v e s u g g e s t ways o f f i n d i n g one e i g e n v a l u e and v e c t o r a t a t i m e . No d e f l a t i o n o r o r t h o g o n a l i z a t i o n t e c h n i q u e i s needed f o r p r o b l e m s w i t h d i s t i n c t e i g e n v a l u e s . O n l y a few v e c t o r s o f s t o r a g e a r e r e q u i r e d and t h e p r o g r a m s we have w r i t t e n run i n v e r y s m a l l memories and c o u l d be u s e d on v e r y s m a l l m a c h ines s u c h as the programmable desk top c a l c u l a t o r s . The a i m o f t h e work r e p o r t e d h e r e was t o f i n d a r e l i a b l e y e t s i m p l e method f o r s o l v i n g t h e s y m m e t r i c m a t r i x e i g e n p r o b l e m when t h e m a t r i c e s i n v o l v e d would n o t n e c e s s a r i l y f i t in main memory in t h e c o m p u t i n g d e v i c e . In t h i s r e g a r d , t h e methods o f N e s b e t [7] and S h a v i t t [6J had p r o v e d u n r e l i a b l e , w h i l e t h a t o f D a v i d s o n [8j was t o o compl i c a t e d f o r t h e t a r g e t m a c h i n e s . Conjugate

Gradients

B e a l e [9] has g i v e n a s h o r t but l u c i d d e r i v a t i o n o f t h e c o n j u g a t e g r a d i e n t s (eg) f a m i l y o f a l g o r i t h m s . The f u n d a m e n t a l idea i s to g e n e r a t e a set o f l i n e a r l y independent s e a r c h d i r e c t i o n s by means o f a r e c u r r e n c e r e l a t i o n s h i p . T h a t i s , g i v e n a f u n c t i o n S(x) and i t s g r a d i e n t £ ( x _ ) t o g e t h e r w i t h t h e j ' t h s e a r c h d i r e c t i o n t_., t h e n e x t s e a c h d i r e c t i o n i s d e f i n e d a s

ij+i

=

ij - £

α

6

where α i s some p a r a m e t e r , t h e f o r m u l a f o r w h i c h d e t e r m i n e s t h e m e t h o d . Once s e a c h d i r e c t i o n s e x i s t , p r o v i s i o n o f a mechanism f o r a c t u a l l y performing a l i n e seach completes the s p e c i f i c a t i o n of an a l g o r i t h m e x c e p t f o r s t a r t o r r e s t a r t c o n d i t i o n s . W i t h i n t h i s work we c h o s e a l w a y s t o s e t

ί ο - - a.

(7)

and i f t h e a l g o r i t h m had n o t c o n v e r g e d in η s t e p s , where η i s t h e o r d e r o f t h e p r o b l e m , t o r e s t a r t in t h e same way f r o m t h e c u r r e n t p o i n t o r v e c t o r .x. M i n i m i z a t i o n of the p a r t i c u l a r f u n c t i o n S(x)

= x

T

C x

+ 2 x

T

w +

(any

where C i s a g i v e n p o s i t i v e d e f i n i t e s y m m e t r i c a given v e c t o r s o l v e s the l i n e a r equations C χ = -

scalar) m a t r i x and w

is

(9)

w

M o r e o v e r , in t h i s c a s e i t can be shown t h a t c o n v e r g e in no more t h a n η s t e p s , w h i l e t h e

(8)

the a l g o r i t h m l i n e search

Lykos; Minicomputers and Large Scale Computations ACS Symposium Series; American Chemical Society: Washington, DC, 1977.

should

26

MINICOMPUTERS

simplifies

to

the computation o f h. J

=

g

T

g /

A N DLARGE

the

step

t! "~J

C t. -J

SCALE

COMPUTATIONS

length

(10)


S i m i l a r s i m p l i f i c a t i o n s e x i s t f o r t h e R a y l e i g h q u o t i e n t (k). T h e s e are d e a l t w i t h below. B o t h t h e s p e c i a l i z e d and g e n e r a 1 - p u r p o s e eg m i n i m i z a t i o n a l g o r i t h m s are n o t o r i o u s l y s e n s i t i v e to implementation d e t a i l s w h i c h b e d e v i l any c o m p a r i s o n s . T h e s e d e t a i l s a r e d i s c u s s e d a t g r e a t e r l e n g t h i n [10] where t h e a l g o r i t h m s a r e p r e s e n t e d i n s t e p - d e s c r i p t i o n f o r m . C h a p t e r 19 o f [10] i s p a r t l y f o u n d e d on the p r e s e n t s t u d y . Direct

Rayleigh

Quotient

Minimization

As a r e s u l t o f o t h e r r e s e a r c h , a c o n j u g a t e g r a d i e n t s program f o r g e n e r a l f u n c t i o n m i n i m i z a t i o n was a v a i l a b l e [11] and o u r i n i t i a l hope was t o a p p l y t h i s d i r e c t l y t o t h e m i n i m i z a t i o n o f (4) o r (5). The m a j o r d i f f i c u l t y t h i s r e v e a l s i s t h a t b o t h t h e s e f u n c t i o n s a r e homogeneous o f d e g r e e z e r o , i m p l y i n g t h a t t h e s c a l e o r n o r m a l i z a t i o n o f x_ i s a r b i t r a r y . T h i s c a u s e s t h e m a t r i x o f s e c o n d p a r t i a l d e r i v a t i v e s , o r H e s s i a n , o f t h e f u n c t i o n R(x) to be o n l y p o s i t i v e s e m i d e f ï n i t e , v i o l a t i n g some t o t h e a s s u m p t i o n s w h i c h u n d e r l y t h e eg a l g o r i t h m . B r a d b u r y and F l e t c h e r [5] a d o p t e d a s t r a t e g y w h i c h r e s t r i c t s x_ t o I i e on a c o n v e x s u r f a c e , t h u s r e m o v i n g t h e t r o u b l e s o m d e g r e e o f f r e e d o m a t t h e e x p e n s e o f some e x t r a work w i t h i n t h e p r o g r a m . One c o u l d a l s o f i x one o f t h e e l e m e n t s o f x_ and a l l o w a l l t h e o t h e r s t o v a r y , e x c e p t t h a t t h e c h o s e n component s h o u l d p e r h a p s be v e r y s m a l l compared t o t h e r e s t , r e s u l t i n g i n v e r y l a r g e a l t e r a t i o n s in t h e s e a t e a c h s t e p o f t h e a l g o r i t h m . R a t h e r t h a n employ s p e c i f i c c o n s t r a i n t s , we t r i e d v a r i o u s p e n a l t y f u n c t i o n s t o impose a n o r m a l i z a t i o n on x. T h e s e a r e summarized in T a b l e I. U n f o r t u n a t e l y we must r e p o r t t h a t none o f t h e s e t e c h n i q u e s c a n be c o n s i d e r e d g e n e r a l l y u s e f u l due t o t h e s l o w c o n v e r g e n c e t o t h e minimum o f t h e f u n c t i o n . Our t e s t s were p r i m a r i l y made u s i n g t h e b i h a r m o n i c m a t r i x d e f i n e d by Ruhe [12] w i t h Β t h e i d e n t i t y . P r i m a r i l y the general m i n i m i z a t i o n is i n e f f i c i e n t because t h e l i n e s e a r c h i s i n e x a c t . F o r t h e R a y l e i g h q u o t i e n t w i t h o u t any p e n a l t y f u n c t i o n t o m a i n t a i n n o r m a l i z a t i o n t h i s can be c o r r e c t e d s i n c e t h e g r a d i e n t o f (k) i s

£ = 2 (A χ -

R Β χ ) / x

T

Β χ

(11)

so t h a t t h e one d i m e n s i o n a l m i n i m i z a t i o n o f R(x_ + h t_) w i t h r e s p e c t t o t h e s t e p l e n g t h h i s a c c o m p l i s h e d by t h e s o l u t i o n o f a q u a d r a t i c e q u a t i o n . T h i s i s s t r a i g h t f o r w a r d , t h o u g h n o t t o be u n d e r e s t i m a t e d ( s e e A c t o n [13]). A p a r t from t h e s e a r c h , the R a y l e i g h q u o t i e n t m i n i m i z a t i o n can be o r g a n i z e d so t h e r e c u r r e n c e (6) g e n e r a t e s s e a r c h d i r e c t i o n s


2.

NASH

AND

27


NASH

which a r e conjugate w i t h r e s p e c t ot the l o c a l Hessian m a t r i x . T h i s a v o i d s s e a r c h d i r e c t i o n s which c a u s e l a r g e growth in the e l e m e n t s o f x. The a p p r o a c h i s due t o G e r a d i n [ 1 h]. H o w e v e r , he d o e s n o t s p e c i f y t h e d e t a i l s o f h i s i m p l e m e n t a t i o n . T h a t w h i c h we have a d o p t e d i s g i v e n in f u l l by Nash [ 1 0 ] .


Inverse

Iteration

by C o n j u g a t e

Gradients

E x c e p t by t h e u s e o f o r t h o g o n a 1 i z a t i o n o r p r o j e c t i o n t e c h niques R a y l e i g h q u o t i e n t m i n i m i z a t i o n cannot f i n d o t h e r than the e x t r e m e e i g e n s o l ut i o n s o f ( 1 ) . H o w e v e r , i n v e r s e i t e r a t i o n can f i n d an e i g e n s o l u t i o n w i t h e i g e n v a l u e c l o s e s t t o any p r e s c r i b e d s h i f t k p r o v i d e d a s t a r t i n g v e c t o r n o t t o t a l l y d e f i c i e n t in t h e r e q u i r e d e i g e n v e c t o r i s g i v e n , ( i n some c a s e s t h i s i s a n o n t r i v i a l p r o v i s i o n . F u r t h e r m o r e , f o r m u l t i p l e e i g e n v a l u e s we c a n o n l y f i n d one e i g e n v e c t o r f r o m t h e r e l e v a n t s u b s p a c e u n l e s s some o r t h o g o n a l i z a t i o n scheme i s e m p l o y e d . ) C o m p a r i n g e q u a t i o n s ( 2 ) and ( 9 ) shows method we s h o u l d make t h e identification C =

(A

-

kB)

;

w = -

Β χ.

;

that

to

u s e t h e eg

χ = y.

(12)

However, C was s u p p o s e d t o be p o s i t i v e d e f i n i t e , w h i c h i s now i m p o s s i b l e i f k i s g r e a t e r t h a n t h e s m a l l e s t (most n e g a t i v e ) e i g e n v a l u e o f o u r p r o b l e m ( 1 ) . Then one may t r y t o s o l v e t h e l e a s t squares problem (A

-

kB) (A T

-

kB)

y.

=

(A

-

kB)

T

Β χ.

(13)

w h i c h d o e s have a p o s i t i v e s e m i d e f i n i t e c o e f f i c i e n t m a t r i x but w h i c h i n c r e a s e s t h e amount o f work and w o r s e n s t h e n u m e r i c a l c o n d i t i o n of the equation system. Following Ruhe and W i b e r g [k] we have o p t e d t o i g n o r e t h e r e q u i r e m e n t t h a t C be p o s i t i v e d e f i n i t e , t h e r e b y d i s c a r d i n g t h e c o n v e r g e n c e r e s u l t s w h i c h have been p r o v e n f o r t h e eg a l g o r i t h m a p p l i e d t o l i n e a r e q u a t i o n s . To c o m p e n s a t e f o r t h i s , we must c h e c k t h a t t h e s t e p l e n g t h computed in ( 1 0 ) i s n o t t o o l a r g e , s i n c e t h i s w o u l d i m p l y t h a t t h e s h i f t k i s t o o c l o s e t o an e i g e n v a l u e . T h i s p o s e s t h e d i l e m m a : 1 ) For r a p i d convergence o f i n v e r s e i t e r a t i o n i t i s d e s i r a b l e have k a s c l o s e t o an e i g e n v a l u e a s p o s s i b l e . 2 ) F o r t h e eg a l g o r i t h m t o work w i t h o u t o v e r f l o w , k s h o u l d n o t be t o o c l o s e t o an e i g e n v a l u e . to

Ruhe and W i b e r g [k] n e v e r t h e l e s s u s e d e g - i n v e r s e iteration to r e f i n e e i g e n v e c t o r s , using s h i f t s v e r y c l o s e to the e i g e n v a l u e , a f t e r w a r d s r e f i n i n g t h e e i g e n v a l u e i t s e l f by means o f t h e R a y l e i g h q u o t i e n t . We have a p p l i e d t h e i r method t o t h e more d i f f i c u l t case where t h e s t a r t i n g v e c t o r may c o n t a i n o n l y a s m a l l component o f t h e d e s i r e d e i g e n v e c t o r and o n l y a c r u d e e s t i m a t e o f t h e e i g e n v a l u e may be a v a i l a b l e . In t h e c a s e where t h e s t e p l e n g t h ( 1 0 ) e x c e e d s some t o l e r a n c e ( t h e r e c i p r o c a l o f t h e s q u a r e r o o t o f t h e


28

MINICOMPUTERS

A N DLARGE

SCALE

COMPUTATIONS

m a c h i n e p r e c i s i o n ) o u r p r o g r a m h a l t s , s i n c e we have y e t t o be c o n v i n c e d t h a t a s u i t a b l e a u t o m a t i c s h i f t i n g p r o c e d u r e c a n be devi sed. As a c o n t r o l on t h e p r o g r e s s o f t h e i n v e r s e i t e r a t i o n we compute t h e r e s i d u a l s


r_

=

(A -

sB) χ

(14)

where s i s t h e c u r r e n t e s t i m a t e o f t h e e i g e n v a l u e . The sum o f s q u a r e s rjr_ p r o v i d e s a c o n v e n i e n t c o n v e r g e n c e t e s t . Unfortunately d i f f e r e n t n o r m a l i z a t i o n s a r i s e in t h e i m p l e m e n t a t i o n s o f t h i s a l g o r i t h m and t h a t o f G e r a d i n , so t h a t t h e c o n v e r g e n c e c r i t e r i a a r e n o t d i r e c t l y c o m p a r a b l e . We have s o f a r c h o s e n n o t t o p e r f o r m t h e e x t r a c o m p u t a t i o n n e e d e d t o make t h e c o n v e r g e n c e c r i t e r i a ident i c a l . Hybrid

Algorithms

H a v i n g now s e v e r a l p r o g r a m s a t o u r d i s p o s a l , i t i s t e m p t i n g to c o n s i d e r , s a y , combining a f u n c t i o n m i n i m i z a t i o n w i t h i n v e r s e i t e r a t i o n to r e f i n e the v e c t o r . E a r l y though q u i t e e x t e n s i v e t r i a l s with the general minimization a l g o r i t h m using f u n c t i o n D o f T a b l e I f o l l o w e d by i n v e r s e i t e r a t i o n a s i n t h e p r e v i o u s s e c t i o n showed t h a t t h i s c o m b i n a t i o n c o u l d be more e f f e c t i v e t h a n f u n c t i o n m i n i m i z a t i o n a l o n e . However, i f o n l y t h e extreme e i g e n s o l u t i o n s a r e d e s i r e d , t h e G e r a d i n a l g o r i t h m can in o u r e x p e r i e n c e s u p p l y an a c c u r a t e s o l u t i o n more e f f i c i e n t l y than any c o m b i n a t i o n m e t h o d , i n c l u d i n g G e r a d i n w i t h i n v e r s e iteration. D i s c u s s i o n and Examples T a b l e l i l i s t s t h e programs e i t h e r d e v e l o p e d o r used f o r c o m p a r i s o n p u r p o s e s w i t h i n t h i s s t u d y . We t e n t a t i v e l y recommend t h e a l g o r i t h m s w h i c h u n d e r l y t h o s e c a l l e d GER and INVIT. T h e r e a d e r i s a d v i s e d t h a t t h i s c a n be no more t h a n a t e n t a t i v e r e c o m m e n d a t i o n s i n c e work i s s t i l l g o i n g on i n t h i s a r e a . F o r i n s t a n c e , the a u t h o r s understand that P r o f . G.W.Stewart o f the U n i v e r s i t y o f M a r y l a n d has been d e v e l o p i n g a l g o r i t h m s f o r t h e sparse eigenproblem. His bibliography [16] p r o v i d e s a w e a l t h o f p o s s i b l e d i r e c t i o n s f o r r e s e a r c h . A t t h e t i m e o f w r i t i n g , we a r e s t i l l i n v e s t i g a t i n g t h e eg methods o f F r i e d [YJ]· T h e r e r e m a i n t o be s e t t l e d t h e q u e s t i o n s o f (1) c o n v e r g e n c e c r i t e r i a and (2) s t a r t i n g v a l u e s f o r t h e s h i f t a n d t h e i n i t i a l v e c t o r . F o r t h e program GER c o n v e r g e n c e i s assumed when t h e g r a d i e n t norm s q u a r e d £ 0

C)

x (A-kl) x

D)

x7(A-kl) x/x x +

T

2

2

+

(x x-l) T

T

Least s q u a r e s . T h i s i s not a Rayleigh quotient. Difficulties in f i n d i n g c e r t a i n s o l u t i o n s unless starting vector "good".

2

(x x-l) T

Extremely

2

poorly E)

x x/(x (A-kl) x)

Table

T

T

II:

to s o l v e the

Comments

Funct ion A)

29


NASH

2

+

(x x-l) T

Programs a p p e a r i n g

2

s l o w when

eigenvalues

separated.

An u n s u c c e s s f u l a t t e m p t t o separate close eigenvalues. Slow c o n v e r g e n c e .

in t h e

study.

GER

The G e r a d i n a l g o r i t h m

INVIT

Inverse i t e r a t i o n using conjugate g r a d i e n t s s o l u t i o n of the l i n e a r equations (2).

TQL

A BASIC v e r s i o n o f TQL1 i n [15] f o r t h e e i g e n v a l u e s o n l y o f a symmetric t r i d i a g o n a l m a t r i x .

SSE

Eigenvalues only using

a s implemented

o f a symmetric

Sturm s e q u e n c e s

tridiagonal

in [ 1 0 ] . for

matrix

[3].

NES

Nesbet's algorithm

MOR

Method o f o p t i m a l

CGT(B)

General conjugate g r a d i e n t s penalized Rayleigh quotient

[J]. relaxation

[6]. function minimizer using (B) o f T a b l e I w i t h ζ = 1.


30

MINICOMPUTERS

(This

is

lacking

not

normalized.)

a component

in

However,

the

in

A N DL A R G E

SCALE

some c a s e s

direction

of

the

this

desired

e s p e c i a l l y in c a s e s where t h e r e a r e s y m m e t r i e s (other than t h a t about the p r i n c i p a l d i a g o n a l )

in or

COMPUTATIONS

vector

is

eigenvector,

the the

matrices matrix

elements are i n t e g e r s . In s u c h c a s e s we have u s e d a p s e u d o random number g e n e r a t o r t o p r o d u c e e l e m e n t s in t h e interval if

( - 0 . 5 , 0 . 5 )

the

vector

of

ones

seemed

inappropriate.


The a l g o r i t h m s have been t e s t e d i n BASIC i n a 6 0 0 0 b y t e p a r t i t i o n o f a Data G e n e r a l NOVA o p e r a t i n g i n 2 3 b i t binary f l o a t i n g p o i n t a r i t h m e t i c . T h i s c o r r e s p o n d s q u i t e c l o s e l y t o some o f t h e s m a l l d e s k t o p c o m p u t e r s in s i z e o f memory. Some o f t h e programs IBM 3 7 0 / 1 6 8 digits). 150

GER, INVIT, and CGT( ) - were a l s o t e s t e d on an a n d / o r Amdahl V6 i n s i n g l e p r e c i s i o n ( s i x hexadecimal

GER and l i n e s of

INVIT a r e q u i t e s h o r t , BASIC d e p e n d i n g on t h e

being a p p r o x i m a t e l y 1 0 0 to code f o r initialization,

n o r m a l i z a t i o n and p r i n t i n g o f r e s u l t s . GER was i m p l e m e n t e d and t e s t e d o n c e in l e s s t h a n two h o u r s f r o m a s t e p - d e s c r i p t i o n r e c i p e . T h a t i s , l e s s t h a n two h o u r s e l a p s e d f r o m t h e moment t h e f i r s t l i n e o f BASIC was w r i t t e n u n t i l t h e c o m p u t e r p r i n t e d c o r r e c t r e s u l t s f o r a t e s t p r o b l e m . INVIT i s s l i g h t l y l o n g e r in l i n e s o f program, for vectors, took 4 9 6 2 programs

2022

but r e q u i r e s o n l y 1 9 1 8 b y t e s o f s t o r a g e compared t o GER. M o r e o v e r , i t n e e d s o n l y 6 i n s t e a d o f 7 w o r k i n g so t h a t the o r d e r 1 0 0 problem r e p o r t e d in T a b l e III b y t e s t o run w i t h GER but o n l y 4 4 2 6 w i t h INVIT. B o t h r e q u i r e s u b - p r o g r a m s t o compute t h e r e s u l t s o f t h e

mult i p i i c a t ions ν By way o f d e f i n e d by

=

A x ;

illustration,

Α..

=

£

=

Β x_

consider

- 2 + l/[(n+l)

2

+

the

i ] 2

(15)

tridiagonal

for

matrix

i=l,2,...,n. (16)

A

i,i

+

1

=

A

i

+

l,i

"

1

f

o

r

i - 1 . 2 , - . n - l .

w h i c h we c a l l t h e F r o b e r g m a t r i x . Β w i l l be s e t t o t h e i d e n t i t y . On t h e NOVA we computed t h e s m a l l e s t e i g e n s o l u t i o n o f t h i s m a t r i x f o r o r d e r s η = 4 , 1 0 , 5 0 , and 1 0 0 . The r e s u l t s , w i t h c o m p a r a b l e values f o r o t h e r programs, r e s u l t s were f o u n d on t h e

a r e g i v e n in T a b l e 1 1 1 . V e r y s i m i l a r IBM 3 7 0 . T h e y i l l u s t r a t e what we have

o b s e r v e d i n a l l c a s e s t o d a t e , t h a t i s , t h a t when o n l y t h e s m a l l e s t (or l a r g e s t ) e i g e n v a l u e i s d e s i r e d t h e p r o g r a m GER i s t h e most e f f i c i e n t y e t p r o d u c e s a c c u r a t e s o l u t i o n s . U s i n g t h i s program w i t h t h e m a t r i x (A to o b t a i n

a solution

-

k

I)

(17)

2

with eigenvalue

closest

to

k causes a


very

2.

NASH

31


A N D NASH

s u b s t a n t i a l slow-down in t h e r a t e o f c o n v e r g e n c e . T h i s i s h a r d l y s u r p r i s i n g s i n c e we have s q u a r e d t h e c o n d i t i o n number o f t h e m a t r i x a s w e l l a s d o u b l i n g t h e work i n e a c h m a t r i x - v e c t o r product s t e p o f t h e p r o g r a m . T h e g e n e r a l f u n c t i o n m i n i m i z e r CGT(B) p e r f o r m s r e a s o n a b l y w e l l on t h e s e m a t r i c e s but u s e s r o u g h l y f i v e t i m e s a s much e f f o r t . Inverse i t e r a t i o n (INVIT) i s a l m o s t a s t e d i o u s . However, t h e r a t e o f c o n v e r g e n c e h e r e i s g o v e r n e d by t h e rat io


ρ

=

(e

}

-

k)/(e

2

-

(18)

k)

where e^ i s t h e e i g e n v a l u e c l o s e s t t o k a n d e « t h e n e x t c l o s e s t . INVIT i n i t s m a j o r i t e r a t i o n c o n v e r g e s a s f a s t a s powers o f ρ tend to z e r o . F o r t u n a t e l y , as t h e e i g e n v e c t o r begins to dominate

Table Program

III :

Minimal

Order

GER INVIT TQL SSE CGT(B) NES MOR

4 4 4 4 4 4 4

GER INVIT TQL SSE CGT(B) NES MOR

10 10 10 10 10 10 10

GER INVIT TQL SSE CGT(B)

50 50 50 50 50

GER INVIT TQL SSE

100 100 100 100

eigensolution

Eigenvalue

0.350142

of

Rayleigh

Froberg's matrix. quotient

0.350144 0.350144

Matrix

produc

5 27

(7)

0.350145 0.350144 0.350144 0.350144 0.350505

7.44379E-2 7.44434E-2 7.44405E-2

3.48618E-3 3.49247E-3 3.48753E-3

8.86582E-4 9.64072E-4 8.89313E-4

36 10 >1000

7.44406E-2 7.44407E-2

1 1 43 (6)

7.44407E-2 f a i 1 ed 7.48813E-2

57 >1000

3.48733E-3 #.4873E-3

26

3.48748E-3

123

8.89216E-4

51 242 (4)

8.89202E-4

147

(5)

The f i g u r e i n p a r e n t h e s e s a f t e r t h e number o f m a t r i x p r o d u c t s i s t h e number o f i n v e r s e i t e r a t i o n s . In a l l c a s e s t h e s h i f t k=0.


32

MINICOMPUTERS AND LARGE SCALE COMPUTATIONS

x^ , the eg linear equations solution is usually accomplished in many fewer than η matrix-vector products. This has also been observed by Ruhe and Wiberg [4]. The number of inverse iterations, starting form a shift k=0, has been included in parentheses behind the number of matrix-vector products in Table III. Acknowledgement s


While this study used very minimal computing resources, we wish to acknowledge time provided on the Data General NOVA and IBM 370 (Datacrown Ltcl. ) at Agriculture Canada and the Amdahl V6 at the University of Alberta. Literature Cited [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17]

Hestenes M.R. and Stiefel E . , J. Res. Nat. Bur. Standards Section Β (1952) 49, 409-436. Fletcher R. and Reeves C.M., Computer Journal (1964) 7, 149-154. Wilkinson J.H., "The algebraic eigenvalueproblem ,Clarendon Press, Oxford, 1965. Ruhe A. and Wiberg T., BIT (1972) 12, 543-554. Bradbury W.W. and Fletcher R., Numer. Math. (1966) 9.,259-267. Shavitt I., Bender C.F., Pipano Α., and Hosteny R.P., J. Computational Physics (1973) 11, 90-108. Nesbet R.K., J. Chem. Phys. (19657) 43, 311-312. Davidson E.R., J. Computational Physics (1975) 17, 87-94. Beale E.M.L., in Lootsma F.A., "Numerical methods for nonlinear optimization", 39-44, Academic Press, London, 1972. Nash J.C., "Compact numerical methods: linear algebra and function minimization", To be published, probably in 1978. Nash J.C., "Function minimization with small computers", submitted to ACM Trans. Math. Software, 1976. Ruhe Α., in Collatz L., "Eigenwerte Probleme", 97-115, Birkhäuser Verlag, Basel, 1974. Acton F.S., "Numerical methods that work", 58-59, Harper & Row, New York, 1970. Geradin M., J. Sound. Vib. (1971) 19, 319-331. Bowdler H., Martin R.S., Reinsch C., and Wilkinson J.H., Numer. Math. (1968) 11, 293-306. Stewart G.W., in Bunch J.R. and Rose D.J., "Sparse matrix computations", 113-130, Academic Press, New York, 1976. Fried I. J . Sound. Vib. (1972) 20, 333-342. 11


Minicomputers and Large Scale Computations

Recommend Documents