Iterative Methods for Solving Large Sparse Systems of Linear Algebraic Equations

Computational Techniques for Differential Equations

623

J. Noye (Editor)

0 Elsevier Science Publishers B.V. (North-Holland), 1984

ITERATIVE METHODS FOR SOLVING LARGE SPARSE SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS

LEONARD COLGAN South A u s t r a l i a n I n s t i t u t e o f Technology, Adelaide, South A u s t r a l i a

I I

!

1

I

I 0 opt

I I

I I

2

-0

Len Colgan

624

CONTENTS

1.

INTRODUCTION.

2.

THEMODELPROBLEMS

....................................

625

...............................

625

..............................

627

MATRIX PRELIMINARIES 4. 4.1. 4.2. 4.3. 4.4.

BASIC ITERATIVE METHODS. . . . . . . . . . . . . . . . . . . . . . . . . . . . 628 Jacobi method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 630 Gauss-Seidel method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 632 634 Successive overrelaxation (SOR) method ...................... Symmetric successive overrelaxation (SSOR) method. . . . . . . . . . . . . . 638

5.

CHEBYSHEVACCELERATION . . . . . . . . . . . . . . . . . . . . . . . . . . .

6.

CONJUGATE GRADIENT ACCELERATION

7.

REDUCED SYSTEMS (RED-BLACK SYSTEMS)

8.

ALTERNATING-DIRECTION IMPLICIT (ADI) METHODS

9.

COMPARISON OF ITERATIVE METHODS

. . . . . . . . . . . . . . . . . . 645

. . . . . . . . . . . . . . . 649

. . . . . . . . . 651

. . . . . . . . . . . . . . . . . . .654

APPENDICES. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Appendix A: A list of symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . Appendix B . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Appendix C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Appendix D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Appendix E . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Appendix F . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . BIBLIOGRAPHY

640

.....................................

656 656 658 661 664 668 673 677

Iterative Methods

1.

625

INTRODUCTION Because o f t h e s u b s t a n t i a l i n c r e a s e i n speed and e f f i c i e n c y o f modern d i g i t a l computers, n u m e r i c a l methods i n v o l v i n g i t e r a t i v e processes have r e q a i n e d much p o p u l a r i t y . F o r example, t h e d i s c r e t i z a t i o n o f p a r t i a l d i f f e r e n t i a l e q u a t i o n s , u s i n g e i t h e r f i n i t e d i f f e r e n c e s o r f i n i t e elements, i n v a r i a b l y leads t o t h e problem of s o l v i n g a l a r g e system o f l i n e a r equations. I n p a r t i c u l a r , e l l i p t i c p a r t i a l d i f f e r e n t i a l equations i n two dimensions, o r more s o i n t h r e e dimensions, y i e l d l i n e a r systems which, because o f t h e v e r y l a r q e number o f equations w i t h v e r y few non-zero c o e f f i c i e n t s , sugqest an i t e r a t i v e method r a t h e r than a d i r e c t method. The a v a i l a b l e l i t e r a t u r e on i t e r a t i v e techniques i s vast, whi!e t h e number o f s u b s t a n t i a l l y d i f f e r e n t methods i s e q u a l l y vast. Hence a s u i t a b l e s e l e c t i o n has been made, t a k i n g i n t o account t h e s i m p l i c i t y o f d e s c r i p t i o n , ease of i m p l e m e n t a t i o n on a computer, e f f i c i e n c y i n t h e use o f s t o r a g e and time, a v a i l a b i l i t y o f s o f t w a r e , and most i m p o r t a n t l y , p e r s o n a l preference. Consequently, r a t h e r than d e t a i l i n g a succession o f methods and t h e i r many v a r i a n t s , a s m a l l r e p r e s e n t a t i v e s e l e c t i o n has been i n c l u d e d w i t h a l i s t o f r e f e r e n c e m a t e r i a l f o r those i n t e r e s t e d enouqh t o seek f u r t h e r in f o rma t i on.

A

b r i e f summary o f t h e c o n t e n t s o f t h i s a r t i c l e i s as f o l l o w s : The f i r s t sect i o n i n t r o d u c e s two model problems f o r subsequent d i s c u s s i o n , and t h i s i s f o l l o w e d by t h e necessary m a t r i x p r e l i m i n a r i e s . L i n e a r s t a t i o n a r y i t e r a t i v e methods of t h e f i r s t degree a r e then described, a l o n g w i t h t h e i r correspondinn convergence c r i t e r i a . The p a r t i c u l a r methods d e t a i l e d are t h e b a s i c Jacobi, Gauss-Seidel, Successive O v e r r e l a x a t i o n (SOR), and S y m m t r i c Successive Overr e l a x a t i o n (SSOR) methods. O f these, i t i s f a i r t o say t h a t o n l y Successive O v e r r e l a x a t i o n has any p r a c t i c a l c l a i m f o r d i r e c t implementation. However, t h e Jacobi and SSOR methods a r e i n c l u d e d because they can be a c c e l e r a t e d subs t a n t i a l l y u s i n g Chebyshev o r c o n j u g a t e g r a d i e n t a c c e l e r a t i o n procedures, which comprise t h e n e x t two s e c t i o n s . Then f o l l o w s a d e s c r i p t i o n o f how t h e s i z e o f t h e system can p o s s i b l y be reduced u s i n g a s o - c a l l e d " r e d - b l a c k " o r d e r i n n on r e c t a n q u l a r domains. A1 t e r n a t i n g - D i r e c t i o n I m p l i c i t (ADI) methods are r e p r e s e n t e d by t h e Peaceman-Rachford method. Throuahout these d i s c u s s i o n s , i t s h o u l d be understood t h a t v a r i a n t s o f these b a s i c methods can y i e l d f u r t h e r improvements. For example i t i s p o s s i b l e t o use b l o c k - i t e r a t i v e methods ( o r l i n e - i t e r a t i v e methods) i n c e r t a i n circumstances, and these can produce w o r t h w h i l e t i m e s a v i n g s when compared t o t h e fundamental p o i n t - i t e r a t i v e methods. These a r e d e s c r i b e d e x t e n s i v e l y i n many r e f e r e n c e s i n t h e b i b l i o n r a p h y . F i n a l l y , a b r i e f comparison o f s e l e c t e d i t e r a t i v e methods i s i n c l u d e d . Append i x A c o n t a i n s a l i s t o f t h e symbols used i n t h i s c h a p t e r , and subsequent Appendices l i s t p o s s i b l e FORTRAN s u b r o u t i n e s (and o u t p u t s ) t o implement a number o f t h e methods discussed.

A s p e c i a l acknowledgement must be made t o P r o f e s s o r David M. Younq o f t h e Center f o r Numerical A n a l y s i s , a t t h e U n i v e r s i t y o f Texas a t Austin. I n prep a r i n g a r e v i e w a r t i c l e such as t h i s , i t was necessary t o s e l e c t n o t a t i o n , d e f i n i t i o n s , p a r t i c u l a r methods and examples f r o m t h e l i s t o f r e f e r e n c e m a t e r i a l . P r o f e s s o r Young's i n f l u e n c e i n t h i s area i s e s p e c i a l l y s i q n i f i c a n t and much o f t h e f o l l o w i n g m a t e r i a l can be a t t r i b u t e d d i r e c t l y t o him. 2.

7HE RODEL PROBLEMS Consider t h e two-dimensional Poisson e q u a t i o n

a2u

7 +

ax

a2u ay

= -1, w i t h u(x,y)

=

0

on t h e boundary o f t h e u n i t square Osx
Len Colgan

626

Using a uniform g r i d of m s h - s i z e h , we can approximate the given e l l i p t i c p a r t i a l d i f f e r e n t i a l equation a t a p o i n t (x,y) by the usual f i n i t e d i f f e r e n c e formula U( X , Y - ~ ) + U ( X - ~ , Yx ), Y - )~+Uu ((x+h ,y)+u(x,y+h) =

-l.

“Model Problem 1” corresponds t o h = = 0.2. Then there would be 16 unknowns corresponding t o the function values a t the i n t e r n a l nodes. Let t h e s e be 16, where the natural ordering i s used. denoted by u, , i = l ,

...,

0 0 - 1 0 0 0 0 0 0 0 0 0 0 o”u,’ -1 4 - 1 o 0 - 1 o o o o o o o o o n u Z 0 - 1 4-1 o 0 - 1 o o o n o o o o o u3 0 0 - 1 4 0 0 0 - 1 0 0 0 0 0 0 0 0 u4 - 1 0 0 0 4 - 1 0 0 - 1 0 0 0 0 0 0 0 us

‘4-1

0-1 0 o 0-1 0 0 0-

0 0 0 0 0 0 0

0 0

0 0

.o

0 0 0

0 0 0

i.e.

A!

=

0 0 0 0 0 0

’1’

’0.04’

1 0.04 1 0.04 1 0.04 1 0.04 0 - 1 4 - 1 0 0 - 1 0 0 0 0 0 0 u6 1 0.04 o 0 - 1 4 - 1 o 0 - 1 o o o o n u7 1 0.04 1 0 0 - 1 4 0 0 0 - 1 0 0 0 0 US = h z 1 = 0 . 0 4 0 - 1 0 0 0 4 - 1 0 0 - 1 0 0 0 US 1 0.04 0 0 - 1 0 0 - 1 4 - 1 0 0 - 1 0 0 u10 1 0.04 0 0 0 - 1 0 0 - 1 4 - 1 0 0 - 1 0 u11 1 0.04 0 0 0 0 - 1 0 0 - 1 4 0 0 0 - 1 u12 1 0.04 0 0 0 0 0 - 1 0 0 0 4 - 1 0 0 u13 1 0.04 0 0 0 0 0 0 - 1 0 0 - 1 4 - 1 0 u14 1 0.04 0 0 0 0 0 0 0 - 1 0 0 - 1 4 - 1 u15 1 0.04 0 0 0 0 0 0 0 0 - 1 0 0 - 1 4,(u16,
b.

The c o e f f i c i e n t m a t r i x A i s s y m t r i c and p o s i t i v e d e f i n i t e , sparse, and has a band-width o f 9. "Model Problem 2" corresponds t o h =

$J = 0.025.

I n t h i s second model problem, we have a l a r g e sparse l i n e a r system Ay = b, where t h e m a t r i x A i s o f o r d e r 1521, and t h e v e c t o r y c o n t a i n s t h e unknown values u , i = 1 , 1521, numbered i n t h e n a t u r a l o r d e r i n g . A has a bandw i d t h of' 79, and t h e r e a r e a t most f i v e non-zero e n t r i e s i n each row o f A.

...,

A p a r t f r o m equations d e r i v e d f r o m p o i n t s a d j o i n i n g a boundary, we have t o s o l v e a system ( h = 1/40) o f t h e f o r m

I t e r a t i v e Methods

627

I n g e n e r a l , numerical s o l u t i o n s o f ( e l l i p t i c ) p a r t i a l d i f f e r e n t i a l e q u a t i o n s l e a d t o l i n e a r systems w i t h t h e c o e f f i c i e n t m a t r i x h a v i n g a s t r u c t u r e somewhat s i m i l a r t o t h e ones generated b y t h e model problems. However, f i n i t e element schemes o f t e n produce a m a t r i x w h i c h i s l e s s sparse and w i t h l e s s p a t t e r n . I n t h e s u b r o u t i n e s l i s t e d i n t h e Appendix, t h e m a t r i x A i s t o be i n p u t i n an unsymmetric sparse form, a l t h o u g h f o r t h e second problem i t w o u l d c l e a r l y be more economical t o use a s y m m t r i c sparse s t o r a g e scheme. For t h e f i r s t model problem, t h i s w o u l d r e q u i r e t h e d a t a t o be i n p u t as

A 14, -1, -1, -1, 4 , -1, -1, -1, 4 , -1, -1, -1, 4, -1, -1, 4 , -1, -1, -1, -1, 4 , -1, -1, -1, -1, 4, -1, -1, -1, -1, 4 , -1, -1, 4, -1, -1, -1, -1, 4, -1, -1, -1, -1, 4, -1, -1, -1, -1, 4 , -1, -1, 4 , -1, -1, -1, 4 , -1, -1, -1, 4, -1, -1, -1, 41 J A l l , 2 , 5 , 1, 2, 3, 6 , 2, 3 , 4 , 7, 3, 4 , 8, 1, 5 , 6 , 9, 2, 5 , 6, 7,10,3,6,7,8,11,4,7,8,12,5,9,10,13,6,9,10,11, 14, 7, 10, 11, 12, 15, 8, 11, 12, 16, 9 , 13, 14, 10, 13, 14, i 5 , 11, i 4 , 15, 1 6 , 12, 15, 161 ISTARTIl, 4, 8, 12, 15, 19, 24, 29, 33, 37,

42, 4 7 , 51, 5 4 , 58, 6 2 , 651

Althouqh t h e s o l u t i o n o f t h e second model problem, h =,,1/40, i s more t y p i c a l of “ l a r g e sparse systems o f l i n e a r a l g e b r a i c equations , t h e f i r s t problem, w i t h o n l y 16 v a r i a b l e s , has been i n t r o d u c e d t o enable some o u t p u t s t o be shown i n t h e Appendix.

3.

MATRIX PRELIMINARIES L e t B be a r e a l NxN m a t r i x , and ,x an N-dimensional r e a l v e c t o r . (a)

Then p ( B ) , t h e s e c t r a l r a d i u s o f B, i s t h e maximum o f t h e moduli o f t h e eigenvalues o f d e t m(B) and M(B) ( o r , i f unambiquous, merely m and M) r e p r e s e n t t h e s m a l l e s t and l a r g e s t eiqenvalues o f B, s o t h a t m(B)
(b)

For

(c)

We d e f i n e t h e 2-norm o f ,x by

Il_xllz

11 11,

m a t r i x norm

=

(x&=

(

i

P ( B ) 2 IlBll

x:)”

i=1

The 2-norm o f B i s g i v e n b y

1 1 ~ 1 =1 ~ I~(B~B)I+, and t h i s m a t r i x norm IIB_xlI2 IIBIIZ =

i5t-Q

-.

Clearly /IB_x/I2 (d)

i s s u b o r d i n a t e t o t h e v e c t o r norm i n t h a t

II,x112

5

/lBl/2//x2(j/2.Also, i f B i s svmmetric, t h e n / / B / / z= p(B).

For any n o n s i n g u l a r m a t r i x !J, we d e f i n e t h e Id-norm o f t h e v e c t o r ,x and of t h e m a t r i x B b y

Len Colaan

628

and

IIBllw

=

IIWBW-'112.

Hence, i f I4BW-l i s symmetric, t h e n l j B l l w = P ( B ) . (e)

B i s symmetric p o s i t i v e d e f i n i t e i f and o n l y i f B i s symmetric and z'Bp0 f o r a l l _x#_O. Then a l l t h e eigenvalues o f B a r e p o s i t i v e , and &here e x i s t s a unique symmetric p o s i t i v e d e f i n i t e m a t r i x , denoted by B , whose square i s equal t o B . Also, B i s symmetric p o s i t i v e d e f i n i t e i f and o n l y i f t h e r e e x i s t s a n o n - s i n g u l a r m a t r i x X such t h a t B = X'X.

F O

( f ) B i s i r r e d u c i b l e (indecomposable) i f i t cannot be p u t i n t o t h e f o r m ( G H), where F and H a r e square, by simultaneous row and column permutations.

(9)

A m a t r i x B i s s a i d t o have "Pro e r t A i f , by a p p r o p r i a t e simultaneous row and column p e r m u t a t i o n s , - d k d k i r i t t e n i n t h e f o r m

where D1 and D2 a r e square diagonal m a t r i c e s . (h)

B = Ibll I has weak diagonal dominance i f Ibl,

I

N 2

1 /bil 1

1 =1

(l
I #1

and f o r a t l e a s t one i, t h e i n e q u a l i t y i s s t r i c t . Then i f B i s a r e a l symmetric m a t r i x which i s i r r e d u c i b l e , b >O f o r N, and B has weak diaqonal dominance t h e n B i s symnekkic i=l, positive definite.

...,

The m a t r i c e s A i n t h e l i n e a r systems A! = ,b d e r i v e d f r o m t h e model problems a r e i r r e d u c i b l e , have " P r o p e r t y A", have weak diagonal dominance and a r e s y m t r i c positive definite.

4.

B A S I C ITERATIVE METHODS

Consider a l a r g e sparse l i n e a r system

A!

=

b,

(4.1)

where A i s a g i v e n r e a l , n o n s i n g u l a r NxN m a t r i x , and b i s a g i v e n r e a l Ncomponent vector. Such systems t h a t a r i s e f r o m d i s c r g t i z a t i o n s o f p a r t i a l d i f f e r e n t i a l equations y i e l d a m a t r i x A t h a t i s u s u a l l y symmetric p o s i t i v e d e f i n i t e , and so t h i s assumption s h a l l b e made. I n any case, such a c o n d i t i o n i s o f t e n necessary t o guarantee convergence o f t h e i t e r a t i v e methods considered.

A l i n e a r s t a t i o n a r y i t e r a t i v e method o f t h e f i r s t degree can be c o n s t r u c t e d I t h e o r e t i c a l l y ) i n t h e f o l l o w i n q manner. Suppose t h e r e e x i s t s a n o n s i n o u l a r m a t r i x Q which i n some sense can be imagined'as b e i n g an approximation t o A, b u t , u n l i k e A, i s r e a d i l y i n v e r t i b l e . More p a r t i c u l a r l y , i t i s p o s s i b l e t o s o l v e c o n v e n i e n t l y any l i n e a r system o f t h e f o r m Q_x = ,c. F o r example, Q may be

Iterative Methods

629

a diagonal, t r i - d i a g o n a l , or lower t r i a n g u l a r matrix, o r t h e product of two o r more such matrices. Q i s c a l l e d a s p l i t t i n g matrix. The system (4.1) i s c l e a r l y equivalent t o

Q_u = (Q-A)!

b

+

which leads t o the i t e r a t i v e process - (Q-A)_~'"' +

Q_U(nC1)

i s not usually found e x p l i c i t l y , t h i s can be written as

Although Q-'

-

=

where B = I

-

U(n+l)

k.

k,

f

and ,k

Q-'A

(4.2) =

Q-'b.

Because the i t e r a t i o n matrix B i s constant, this process i s c a l l e d s t a t i o n a r y . (Information concerning nonstationary methods, such as Richardson's method, can be found i n Forsythe and Wasow (1960); o r Young (1954), (1971).)

...

The basicmethod ( 4 . 2 ) is conver ent i f the sequence g"', $", _u('), converqes t o the t r u e solutiongy f o r a l l i n i t i a l vectors y c O ) , and t h i s l i m i t i s independent of the choice o f $'). A necessary and s u f f i c i e n t condition f o r convergence i s t h a t p ( B ) < l .

The e r r o r vector a f t e r n i t e r a t i o n s i s

5 ( n ) = u( n )

-lj.

The residual vector i s

c'"'

=

b -

Ag'"',

(from (4.1))

and t h e pseudo-residual vector i s

6'"'

=

-

E(n)

-

-

,,ln+l)

Since y = Bg + E( " + I )

-

+ ,k

=

=

=

k,

B g

-

,,(nl

yln)

, (from ( 4 . 2 ) ) , in t h i s case.

then from (4.2) we have (n)

and hence

$ ( O ) . (")

Consequently,

IIL I12

(0)

II

IIB" 1 1 2 ,

f o r any g")

+ 0.

112

The average r a t e of convergence f o r n i t e r a t i o n s i s defined by 1 R,(B) = - iT Tn I I B n ( ( z a I f Rn(B1) <

Rn(B2),

then BZ i s i t e r a t i v e l y f a s t e r f o r n i t e r a t i o n s than E l .

Then, t o reduce an i n i t i a l error by a prescribed f a c t o r , sav l o m 6 , i t would require approximately

Len Colgan

630

'n 10' i t e r a t i o n s .

"lp.7-

However, the average r a t e of convergence, Rn(B), i s r a r e l y available. The asymptotic r a t e of converaence, denoted by R(B), i s defined by R(B) = lim R ( B ) ,

n-

n

and i t can be shown t h a t R(B) =

- In

(4.3)

p(B).

(Since IIBn112 2 Cp(B)In f o r a l l n t l , then R ( B ) which IIBn112 < 1.)

2

Rn(B)

for any positive n f o r

This asymptotic r a t e of convergence R(B), f o r a convergent matrix B , i s certainly t h e simplest p r a c t i c a l way i n modern usage o f measuring t h e rapidity of convergence. Hence, a crude estimate o f th_e,number o f i t e r a t i o n s required t o reduce the norm of t h e e r r o r by a f a c t o r 10 i s given by (4.4)

However, this estimate can often y i e l d misleadingly low values f o r n. example, see Varga (1962). 1 4.1

Jacobi Method

[For

I t e r a t i v e Methods

(n+l)

UlS

(n+l)

+ u):

= %(u(,;)

+ u';:

u15

=

L(u::'

-t

u:!")

=

L(u\;)

+ ui;)

-t

631

0.04)

+ u'~:) + 0.04)

u):

-t

0.04)

Using an i n i t i a l e s t i m a t e _ u C o ) = Q, s u b r o u t i n e JAC (see Appendix B), a c o r r e s p o n d i n g c a l l i n g program, and a s u i t a b l e s t o p p i n g c r i t e r i o n (see l a t e r ) , t h e o u t p u t i n Appendix B i s obtained. The s p e c t r a l r a d i u s (SR i n t h e s u b r o u t i n e ) i s e s t i m a t e d , and s t o p p i n g i s i n v o k e d when

F o r t h e second model problem, t h e Jacobi method i s e q u i v a l e n t t o r e w r i t i n g (2.1) i n t h e i t e r a t i v e form

We can express t h e m a t r i x A as t h e m a t r i x sum A=D+L+U, where D i s t h e diagonal m a t r i x such t h a t d,l = a , 1 , and L and U a r e r e s p e c t i v e l y s t r i c t l y l o w e r and upper t r i a n q u l a r NxN m a t r i c e s , whose e n t r i e s a r e those r e s p e c t i v e l y below and above t h e main d i a s o n a l o f A. We can t h e n r e w r i t e (4.1) as DL.~ = -(L+U)!

+ ,b,

which i n t u r n y i e l d s t h e Jacobi i t e r a t i o n method

-

U(n+l)

= - 0 - l (L+U)y'"'

+ D-'_b.

Hence our general n o t a t i o n corresponds t o a s p l i t t i n g m a t r i x i t e r a t i o n matrix

BJ :-D-l(L+U)

I

=

-

D-' A, and

k,

Q,

:D, an

z D-'_b.

I t s h o u l d be p o i n t e d o u t t h a t i t i s p o s s i b l e t o n o r m a l i z e t h e i n i t i a l m a t r i x problem so t h a t t h e diagonal elements a r e u n i t y (see Varqa (1962)), and t h i s a v o i d s numerous d i v i s i o n o p e r a t i o n s d u r i n g t h e i t e r a t i v e process.

The Jacobi method w i l l converge p r o v i d i n g p(B,) < 1. However, t h i s i s n o t guaranteed, even though t h e m a t r i x A m i g h t be s y m m t r i c p o s i t i v e d e f i n i t e , because, a l t h o u g h i t can be shown t h a t t h e eiqenvalues o f B a r e r e a l a n d < 1, i t i s p o s s i b l e t h a t m(B,) < -1. (The Jacobi method converse if both A and 2D-A a r e symmetric p o s i t i v e d e f i n i t e . )

A u s e f u l s u f f i c i e n t c o n d i t i o n f o r convergence i s t h a t A i s i r r e d u c i b l e and has weak diaqonal dominance w i t h 4 , > 0 ( a s t r o n g e r c o n d i t i o n t h a n p o s i t i ve d e f i n i teness)

.

Another s u f f i c i e n t e x t r a c o n d i t i o n besides p o s i t i v e d e f i n i t e n e s s i s f o r t h e I n t h i s case t h e e i g e n v a l u e s o f B, a r e m a t r i x A t o have " P r o p e r t y A". hl f o r [A, I < 1. T h i s w o u l d i m p l y e i t h e r z e r o o r e x i s t as r e a l p a i r s t h a t p(B,) = M(B,) < 1.

*

Len CoZgan

632

Because o f i t s importance i n t h e i t e r a t i v e methods s t i l l t o be discussed, we denote XJ = M(BJ) t o be t h e l a r g e s t p o s i t i v e e i g e n v a l u e o f BJ. The Jacobi i t e r a t i o n m a t r i c e s f o r t h e model problems possess a l l these p r o p e r t i e s . I n f a c t , i t can be shown t h a t X, = cosnh i n each case. For h = 0.2,

For h = 0.025,

0.809017.

XJ

XJ = 0.996317.

Hence t h e appro2imate number o f i t e r a t i o n s r e q u i r e d t o reduce t h e e r r o r by a f a c t o r 10- i s n =

I n lo6 - I n cos ah'

For h = 0.2. t h i s i m o l i e s 66 i t e r a t i o n s . and i s v e r i f i e d i n ADDendix 8. I J i t h h = 1/40, t h i s y i e l d s n = 4475. [ i i t h g(") = ,O, i t a c t u a i l y takes 4479.1 We n o t e t h a t t h e converqence r a t e o f t h e Jacobi method appears t o b e very slow, I t can be shown f o r any g e n e r a l i z e d O i r i c h l e t problem t h a t t h e number o f i t e r a t i o n s r e q u i r e d t o achieve a s p e c i f i e d l e v e l o f converqence i s p r o p o r t i o n a l t o l / h 2 f o r a small mesh-size h. The f a c t t h a t t h e Jacobi method can be a c c e l e r a t e d , u s i n q Chebyshev a c c e l e r a t i o n , t o y i e l d a p r o p o r t i o n a l i t y o f l / h , i s adequate reason f o r i t s i n c l u s i o n here. The n e x t two b a s i c methods, Gauss-Sei del and Successive O v e r r e l a x a t i o n , a l t h o u g h i n d i v i d u a l l y f a s t e r t h a n Jacobi, cannot be a c c e l e r a t e d i n t h e same way.

NOTE:

Suppose o u r b a s i c Jacobi i t e r a t i o n was o f t h e f o r m

U(n+l)

=

i

(n) alul-k

(n) +

a2ul-l

+ a3 u(") + a4 u(") + c, f o r ai > 9, i=l,.., 1+1 1+k

4.

T h i s might correspond t o a general e l l i p t i c p a r t i a l d i f f e r e n t i a l equation where, say,

k = l - 1 h ' and u(x,y) i s g i v e n on t h e boundary o f t h e u n i t square. shown t h a t

+

XJ = 2(-

G) cosnh.

Our model problems correspond t o a, = 4.2

Then i t can be

4, i=1,..,4.

Gauss-Sei del Method For t h e f i r s t model problem, t h e Gauss-Seidel i t e r a t i v e method i s r e p r e s e n t e d by

d"+l)=

&(u:"'

+ '":it +

U:n+l)

=

g(,,\"+l)

,,in+')

=

L(~:"+') + u:")

Ul"+l

u:"")

)

+ + +

+ u\"' & ( J l n + l+ ) ~(6") +

= &(u:"+') =

+ 0.04) u:")

+ +

0.04) 0.04)

0.04)

u:"'

t 0.04)

I t e r a t i v e Methods

633

See Appendix C f o r a p o s s i b l e s u b r o u t i n e GS, w i t h an o u t p u t corresponding t o y(O) = 0, and a s t o p p i n q c r i t e r i o n s i m i l a r t o t h a t f o r t h e Jacobi method. I t r e q u i r e d 34 i t e r a t i o n s . For t h e second model problem, t h i s i s e q u i v a l e n t t o r e w r i t i n g (2.1) it e r a t i ve form

i n the

i.e. t h e l a t e s t e s t i m a t e s f o r a l l components of g a r e always used i n subsequent computations. I n m a t r i x form,

Hence t h e s p l i t t i n g m a t r i x i s BGS = - ( P L ) - ' U = I - (D+L)-'A,

*

Q,, E and

P L , the i t e r a t i o n matrix i s 3 (D+L)-l_b.

kGs

I t can be shown t h a t i f A i s symmetric p o s i t i v e d e f i n i t e , t h e n t h e Gaussconverge. (The eigenvalues need n o t be r e a l , and Seidel i t e r a t i o n complex eigenvalues t e n d t o produce an i r r e g u l a r convergence.) However, i f A i s " c o n s i s t e n t l y ordered" (see K i n c a i d and Young (1978); o r Varga (196?), f o r a p r e c i s e d e f i n i t i o n ) and has " P r o p e r t y A", t h e n t h e eigenvalues o f ,B, a r e r e l a t e d t o t h o s e o f B, i n t h e f o l l o w i n g way: t h e e i g e n v a l u e s o f ,B w i l l be r e a l , and f o r each nonzero p a i r 2 XI o f B,, t h e r e corresponds e i qenval ues 0, XIz f o r BGS

.

Hence, i n t h i s case,

and t h e Gauss-Seidel i t e r a t i o n converges a b o u t t w i c e as f a s t as t h e Jacobi i t e r a t i o n . I t w o u l d r e q u i r e about h a l f as many i t e r a t i o n s t o a c h i e v e t h e

Len CoZgan

634

accuracy s p e c i f i e d i n (4.4);

R(Ba,)

i.e.

= 2 R(B J

1. .

T h i s i m p l i e s about 33 i t e r a t i o n s f o r h = 0.2, and about 2240 i t e r a t i o n s f o r h = 1/40. (Using s u b r o u t i n e GS, t h e 2240 i t e r a t i o n s have been verified.) The number o f i t e r a t i o n s i s again, o f course, p r o p o r t i o n a l t o l / h 2 . The Gauss-Seidel method i s i n c l u d e d here m a i n l y t o i n t r o d u c e t h e n e x t procedure, Successive O v e r r e l a x a t i o n . B e f o r e d o i n q so, a b r i e f comment on I d e a l l y , we m i g h t wish t o guarantee a s t o p p i n g c r i t e r i o n i s warranted. that

say, w i t h an added s a f e g u a r d c o n d i t i o n taken n s nmx nmx *

, f o r some

prescribed

However, we cannot measure e i t h e r ( 1 ~ ' " ) ( ( 2 o r ( I L I ( ~ z p r e c i s e l y . Ue can f i n d t h e pseudo-residual b c n ) a t ea& sten, b u t t h i s can o n l y be used e f f e c t i v e l y i f we have a good e s t i m a t e f o r t h e s p e c t r a l r a d i u s o f t h e i t e r a t i o n m a t r i x . I t s h o u l d be p o i n t e d o u t t h a t a s t o p p i n g t e s t o f t h e form

not i m p l y

does 4.3

t h e accuracy s p e c i f i e d above.

Successive O v e r r e l a x a t i o n (SOR) Method For t h e second model problem, t h i s i s e q u i v a l e n t t o r e w r i t i n g (2.1) i n t h e i t e r a t i v e form n + l)

(,,

I

= W{4(Ui-,,

("+I) +

+

(n+') u,.l

+ u;:ie+

u:+":

h2)l

+ (l-o)ul("'.

The parameter o i s c a l l e d a r e l a x a t i o n parameter. The new e s t i m a t e u , ( ~ + ' ) can be regarded as a w e i g h t e d average o f a t e r m o b t a i n e d i n a manner s i m i l a r t o t h e Gauss-Seidel concept and t h e c u r r e n t e s t i m a t e d " ) . T h i s i t e r a t i o n i s l i n e a r , sta+,ionary and c o m p l e t e l y c o n s i s t e n t w i t h ' (4.1), i n t h e sense t h a t i f i t converges, i t gives t h e t r u e s o l u t i o n y. I n m a t r i x form,

-. ,,(n+l)

(E+wL)_u'"*'' _u("+"

-

= w ~ ~ - ' ( - ~ ~ ' " + ' ) IJ_u'"' = C-WU

= (DewL)-'C-wU

The s p l i t t i n g m a t r i x i s

+ (l-w)Dl_u'"'

+

+ (1-w)Dly'"' QSOR

BsoR :( D+oL)-lC-wth(l-w)D1 and

+ !)I

+ (i-w)_U'") o,b

+ (E+wL)-'w_b.

:$0 + L, t h e I t e r a t i o n m a t r i x i s =

I

-

( i D + L)-lA,

635

I t e r a t i v e Methods

kso,

(2+ L I - ' ~ .

F ( D + ~ L ) - '=~ ~1

For w = 1, the method i s merely Gauss-Seidel. I t can be shown tha t i f A i s symmetric positive d ef i n i t e and irreducible, then p(Bso,)
2 A,,

r

X2,

...,

k

XN1 ( X i

> 0 ) and N2 zero

Corresponding t o 2 A, 1 > 0, there are two eigenvalues fying the relationship lJ

+ w

-

1=

WAG,

i.1

o f BsoR

satis(4.3.1)

or equivalently,

For o > 1, t h i s usually produces complex eigenvalues N2 eigenvalues of BsoR ar e a l l equal t o (1-u).

i.1.

The remaining

I n p a r t i c u l a r , we know t h at w = 1 corresponds t o the Gauss-Seidel i t e r a t i o n . Then equation (4.3.1) confirms the e a r l i e r statement concerning the r a t e s of convergence of the Jacobi and Gauss-Seidel methods under the conditions being assumed; i . e . p(BG,) = C p ( B J ) 1 2 = Xi. There i s an optimal value of w, denoted by wept, t h a t minimises the spectral radius o f BSOR, depending on A,; i.e.

w

OPt

minimises - A max s A 5xJ J

i

12* (4.3.2)

2 p(BSoE) = w

-

1, i f

w t u,,~.

Len Colgan

636

T h i s graph shows how p(BS0,) v a r i e s w i t h choices o f w between 0 and 2, and how t h e r a t e o f convergence would improve d r a m a t i c a l l y as w approaches w f r o m below. OPt

, l e t p(B R ) = i. I f we were For convenience o f n o t a t i o n , when o = w t o know w, t , which we can c e r t a i n l y f i 8 8 ' v i a (4.3.57 i f we know A,, then t h e r a t e o? convergence o f t h e SOR method i s a b o u t 2d times t h e square r o o t o f t h e r a t e o f convergence o f t h e Jacobi method (and hence about t w i c e t h e square r o o t o f t h e r a t e o f convergence o f t h e Gauss-Seidel method). I n p a r t i c u l a r , t h e number o f i t e r a t i o n s r e q u i r e d t o o b t a i n t h e e r r o r r e d u c t i o n p r e v i o u s l y p r e s c r i b e d i s now p r o p o r t i o n a l t o l / h . For t h e model problems, A, = cos Hence, from (4.3.2),

w

= 1.8545 f o r h = 1/40.

which i s

I-.

opt

1~

h.

2 = ___ l+sinnh, which i s = 1.2596 f o r h

=

0.2,

and

The s p e c t r a l r a d i u s o f BsoR w i t h w = wa p t i s

0.2596 f o r h = 0.2,

and = 0.8545 f o r h = 1/40.

Using

we expect t o r e q u i r e 11 i t e r a t i o n s f o r h = 0.2 and 88 i t e r a t i o n s f o r h = 1/40. However, u s i n g t h e t r u e v a l u e o f wept as i n p u t s i n t o t h e s u b r o u t i n e SOR o f Appendix D, i t i s found t h a t 15 and 114 i t e r a t i o n s r e s p e c t i v e l y a r e a c t u a l l y r e q u i r e d . T h i s i s because t h e i t e r a t i o n m a t r i x Bso,, c o r r e s ponding t o Q ~ i s~n o,t a b l e t o be d i a g o n a l i s e d , b u t has a Jordan Canonical Form which leads t o a s m a l l e r r e d u c t i o n i n t h e e r r o r s t h a n p r e d i c t e d by t h e s p e c t r a l r a d i u s alone. For e f f e c t i v e use o f t h e SOR method, i t i s e s s e n t i a l t h a t we have an a c c u r a t e e s t i m a t e f o r A,, and thereby G p t . I f we o v e r e s t i m a t e Q ~ t h~ e decrease i n t h e a s y m p t o t i c r a t e o f convergence i s n o t t o o s e r i o u s . However, an underestimate of h p tproduces a c o n s i d e r a b l e r e d u c t i o n i n t h e

,

Iterative Methods

637

convergence r a t e . I f t h e m a t r i x A has “ P r o p e r t y A “ , i t i s i n t e r e s t i n g t h a t t h e v a l u e o f wept i s independent o f t h e o r d e r i n which t h e unknowns a r e a c t u a l l y determined f o r most p r a c t i c a l o r d e r i n g s . E s t i m a t e s f o r A, can sometimes be o b t a i n e d a p r i o r i , o r e l s e an a d a p t i v e procedure can be employed f o r t h e a u t o m a t i c d e t e r m i n a t i o n o f A, as t h e a c t u a l i t e r a t i v e procedure i s b e i n g c a r r i e d o u t . (See Young (1971) o r Hageman (1972).)

A p o s s i b l e a d a p t i v e procedure t h e n can e a s i l y be i n c o r p o r a t e d i n t o an SOR code u s i n g t h e f a c t t h a t g i v e n e i t h e r A, ( t h e s p e c t r a l r a d i u s o f B J ) o r ( t h e s p e c t r a l r a d i u s o f ,,B, w i t h w = Q ~ ~ a)l l , t h r e e q u a n t i t i e s A,, p and wept can be determined v i a r e l a t i o n (4.3.2) as w e l l as (4.3.3) which i s d e r i v a b l e from (4.3.1) and t h e f o r m u l a f o r

i=

p(Bso,).

The f o l l o w i n g a d a p t i v e procedure i s n o t r i g o r o u s i n t h a t i t r e q u i r e s f o r a l l values o f w used w h i l e we seek t h e optimum v a l u e . liowever, i n p r a c t i c e , i t has been q u i t e s u c c e s s f u l .

w < woqt

(i)

Choose a value w so t h a t w < wept i s c e r t a i n ; f o r example, we c o u l d s t a r t w i t h w = 1.

(ii)

Apply a few steps o f SOR, and t h e n use t h e r a t i o o f c o n s e c u t i v e pseudo-residuals

where

g n ) = u_(”cl)

-

d n ), as

an estimate f o r

i.

(iii)

Use t h i s e s t i m a t e f o r i, and t h e c u r r e n t w, t o c a l c u l a t e a corresponding e s t i m a t e f o r A, v i a (4.3.3)

(iv)

Use t h i s A, e s t i m a t e t o c a l c u l a t e an improved v a l u e f o r wept v i a (4.3.2). Then r e t u r n t o (ii).

P r o v i d e d we a r e s e n s i b l e i n d e c i d i n g how o f t e n t o a p p l y t h i s a d a p t i v e process, t h e v a l u e o f o,, can be o b t a i n e d f a i r l y e a s i l y w i t h o u t adding much c o s t t o t h e o v e r a l l Successive O v e r r e l a x a t i o n method. T h i s process t h e o r e t i c a l l y r e l i e s on t h e e s t i m a t e w approaching uOpt f r o m below. As i t does so, t h e convergence r a t e of t h e a c t u a l SOR method improves noticeably.

A

s t o p p i n g t e s t which i s v e r y n e a r l y e q u i v a l e n t t o (n)

llE 112 -___ <

say, i s as f o l l o w s :

A f t e r t h e nth i t e r a t i o n , d e f i n e [~(n),

if

-

I < R‘”’

< 1,

Then t h e SOR process t e r m i n a t e s when

Len CoZgan

638

T h i s t e s t i s found t o be s a t i s f i e d w i t h i n a c o u p l e o f i t e r a t i o n s a f t e r t h e d e s i r e d accuracy has a c t u a l l y been achieved. Appendix 0 shows t h e o u t p u t corresponding t o u s i n g a n i n i t i a l i n p u t of w = 1, and t h e n a d a p t i n g t o f i n d wept. I t takes 17 i t e r a t i o n s , which i s was known i n advance. ( O b v i o u s l y o n l y 2 more t h a n needed i f t h e t r u e W, a b e t t e r i n i t i a l e s t i m a t e f o r w would 6e s u p p l i e d , i f a v a i l a b l e . ) The s u b r o u t i n e SOR does n o t a l l o w t h e v a l u e o f w t o change t o o o f t e n , and i n t h i s case t h e r e i s o n l y one change f r o m t h e i n i t i a l w = 1 t o w = 1.28213 a f t e r 5 i t e r a t i o n s . Suppose we make i t even more d i f f i c u l t t o change t h e v a l u e o f w by a l t e r i n g t h e 1 1 5 t h l i n e o f t h e s u b r o u t i n e ( w i t h 0.01) t o read IF(ABS(

SR-PSNORM/PSPREV) .GT

.o . O ~ O ~ ) T H E N

i . e . change 0.01 t o 0.0001. Then t h e o n l y change i n w i s a f t e r 14 i t e r a t i o n s t o y i e l d a new w = 1.25967. A l t h o u g h t h i s i s a l m o s t e x a c t l y wept, i t takes a t o t a l o f 24 i t e r a t i o n s b e f o r e convergence i s o b t a i n e d . B a s i c a l l y , t h e more i t e r a t i o n s performed w i t h a g i v e n W, t h e b e t t e r w i l l be t h e n e x t e s t i m a t e f o r wept. However, we c a n n o t a f f o r d t o use a poor v a l u e o f w t o o l o n g . Obviously, a compromise i s r e q u i r e d . S i m i l a r l y , f o r t h e second model problem w i t h h = 1/40 and 1521 v a r i a b l e s , t h e a d a p t i v e SOR method ( w i t h 0.0001) r e q u i r e s 154 i t e r a t i o n s , compared t o 114 ifuopt was known. Using an i n i t i a l w = 1, t h e a d a p t i v e method changes w t o w = 1.8000 a f t e r 26 i t e r a t i o n s , t o w = 1.8500 a f t e r 34 i t e r a t i o n s , and t o w = 1.8549 a f t e r 130 i t e r a t i o n s , t h e r e b y c o n v e r g i n g i n 154 i t e r a t i o n s . (Compare w i t h t r u e wept = 1.8545.) 4.4

Symmetric Successive O v e r r e l a x a t i o n (SSOR) method Each i t e r a t i o n o f t h e SSOR method c a n be regarded as two d i s t i n c t o p e r a t i o n s . F i r s t l y , i t uses a f o r w a r d SOR i t e r a t i o n i n which t h e unknowns a r e computed s u c c e s s i v e l y i n t h e i r n a t u r a l o r d e r i n g , and t h e n a backward SOR i t e r a t i o n i n which t h e unknowns a r e computed i n t h e o p p o s i t e o r d e r . The same v a l u e o f t h e r e l a x a t i o n f a c t o r W, 0 < w < 2, is used i n b o t h p a r t s o f each SSOR i t e r a t i o n . L e t BBSOR be t h e i t e r a t i o n m a t r i x a s s o c i a t e d w i t h t h e backward sweep. Then each SSOR i t e r a t i o n , depending upon a r e l a x a t i o n f a c t o r W, c a n be w r i t t e n i n a m a t r i x f o r m such as

where and

BsoR BBSOR

= ( DtwL) = (O+wU)-'

-'C -wU+(

1- w ) D I

[ -wL+( 1-w)Dl.

Hence, t h e i t e r a t i o n m a t r i x f o r t h e Symmetric Successive O v e r r e l a x a t i o n method i s g i v e n by

I t e r a t i v e Methods

Bs ,OR

5

( D+wU)

-'C -wL+(

1-w) 0 1(D+wL )

-'C -mu+(

639

1 -w) 0 1.

I t can b e shown t h a t t h e c o r r e s p o n d i n g s p l i t t i n g m a t r i x i s

I n t e r e s t i n g l y , i f A i s a symmetric p o s i t i v e d e f i n i t e m a t r i x ( i . e . LT=U), and 0 < w < 2 , t h e n t h e s p l i t t i n g m a t r i x Q i s i t s e l f symmetric p o s i t i v e d e f i n i t e . T h i s f a c t can be used t o prove t h a t , i f 0 < w < 2, a l l t h e eigenvalues h. o f ,,B a r e r e a l and l i e i n t h e i n t e r v a l 0 5 X < 1. (See Young (1971).) The r a t e o f convergence o f t h i s SSOR method i s n o t p a r t i c u l a r l y s e n s i t i v e t o t h e p r e c i s e c h o i c e o f t h e r e l a x a t i o n f a c t o r W, and so t h e e x a c t optimum v a l u e o f w i s n o t necessary. A number o f good values f o r w have been expounded, depending upon t h e i n f o r m a t i o n a v a i l a b l e . I n any case, i t i s known t h a t t h e SSOR method converges s a t i s f a c t o r i1y prov i d e d

P(D-'LD-~u)

5

b.

T h i s p a r t i c u l a r c o n d i t i o n i s s a t i s f i e d f o r themodel problem, and f o r most examples o f g e n e r a l i s e d D i r i c h l e t problems p r o v i d e d t h e n a t u r a l o r d e r i n g o f t h e mesh p o i n t s i s employed. I f we know t h a t P(D-~LD-'U) 5 i j , and we also know A,, t h e l a r g e s t e i g e n v a l u e o f t h e c o r r e s p o n d i n g J a c o b i i t e r a t i o n m a t r i x , t h e n a good v a l u e f o r w would be W*

=

2

1+

4 2 m

(See Young (1971). )

Moreover, f o r t h i s v a l u e to*, i t can be shown t h a t

I n a more general s i t u a t i o n , we m i g h t be a b l e t o f i n d a n upper bound f o r

p(D-lLD-lU), say B ( L ~ ) ,and an upper bound f o r X J , say a good v a l u e f o r w would b e W*

1,.

I n t h i s case,

2

=

1+w i t h a c o r r e s p o n d i n g bound f o r p(BssoR).

As a t y p i c a l comparison between t h e e f f e c t i v e n e s s o f t h e SOR method u s i n g wept, and t h e SSOR method u s i n g w*, c o n s i d e r t h e model problems. Here h., = cosah, and we know 4 i %. Hence w*=--

2 1 + fi(1-cosnh)'

which i s v e r y c l o s e t o

2 1t 2 sin

ah

Len CoZgan

640

w

=

1

t

2 sinnh

Comparing these, we can see t h a t t h e SOR Method converges about t w i c e as q u i c k l y as t h e SSOR method. T h i s c o n d i t i o n i s t r u e f o r most l i n e a r systems d e r i v e d f r o m g e n e r a l i z e d D i r i c h l e t problems f o r which 5 J-h, and t h e number o f i t e r a t i o n s necessary i s a p p r o x i m a t e l y p r o p o r t i o n a l t o l / h , f o r small h. Besides needing t w i c e as many i t e r a t i o n s , t h e SSOR method r e q u i r e s t w i c e as much work p e r i t e r a t i o n . However, t h e s e d i s advantages can be overcome i n two ways. F i r s t l y , i t i s p o s s i b l e t o e x t r a p o l a t e t h e SSOR method, which c a n square t h e s p e c t r a l r a d i u s o f t h e i t e r a t i o n m a t r i x . B u t most i m p o r t a n t l y , s i n c e t h e eigenvalues X o f 8 a r e r e a l w i t h 0 IX < 1, we can a c c e l e r a t e t h e convergence o f t h e SSORSoR method by employing e i t h e r Chebyshev a c c e l e r a t i o n or c o n j u g a t e g r a d i e n t a c c e l e r a t i o n (see ahead).

B

A s u b s t a n t i a l improvement i n t h e convergence r a t e c a n be achieved i n t h i s manner. T h i s i s n o t p o s s i b l e w i t h t h e SOR method. S i m i l a r l y , we c a n a c c e l e r a t e t h e Jacobi method, b u t n o t Gauss-Seidel. 5.

CHEBYSHEV ACCELERATION (See Hageman and Young (1981); o r Varga (1962).)

K i n c a i d and Young (1978);

Consider a b a s i c s t a t i o n a r y i t e r a t i v e method o f t h e f o r m

-

$")

U("+l) = B

t

,k,

where t h e eigenvalues o f t h e i t e r a t i o n m a t r i x B a r e r e a l . F o r example, we m i g h t use t h e Jacobi i t e r a t i o n w i t h B, = I Symmetric Successive O v e r r e l a x a t i o n i t e r a t i o n w i t h BsSoR = (D+wU)-'C-wL+( l-w)DI(D+wL)-'C-d+(

-

D-'A,

or the

l-o)OI.

L e t m = m(B) and M =M(B) r e p r e s e n t r e s p e c t i v e l y t h e s m a l l e s t and l a r g e s t eigenvalues X o f B, so t h a t m

2

X

2

M, f o r a l l X.

The Chebyshev a c c e l e r a t i o n t e c h n i q u e does n o t assume t h a t t h e m a t r i x A o f t h e o r i g i n a l l a r g e sparse l i n e a r system i s n e c e s s a r i l y p o s i t i v e d e f i n i t e . I n s t e a d , i t m e r e l y r e q u i r e s t h a t t h e i t e r a t i o n m a t r i x B i s s i m i l a r t o a syrrmetric m a t r i x w i t h r e a l eigenvalues l e s s t h a n one. i.e.

M

< 1.

Even i f t h e b a s i c method d i v e r g e s because m < -1, t h i s a c c e l e r a t i o n process converges and i s s t i l l v e r y e f f e c t i v e . In a l l cases, i t produces a c o n s i d e r a b l e improvement i n t h e r a t e o f convergence. etrization
Fundamentally, we r e q u i r e t h e e x i s t e n c e o f a n o n s i n g u l a r s m a t r i x W such t h a t WBW-' i s symmetric w i t h eigenvalues X

I t e r a t i v e Methods

641

E q u i v a l e n t l y , we r e q u i r e W such t h a t w(I-B)w-~

= I - WBW-'

i s symmetric p o s i t i v e d e f i n i t e , For example, suppose A i s symmetric p o s i t i v e d e f i n i t e , and l e t W = 0'. b a s i c Jacobi method w i l l d i v e r g e i f m(BJ) < -1.) Then

(The

W (I - B ~ 1w-l = D'(D-'A)D-'' = CTC

where C = AHD-''. S i m i l a r l y , f o r t h e SSOR method, l e t

D-'(D+~U),

w=-

o

<

i 2.

Then W(I-BssoR)W-' i s symnetric p o s i t i v e d e f i n i t e . choose W = QrsoRor W = A' i n t h i s case.

I n f a c t , we c o u l d a l s o

Even though i t i s n o t e s s e n t i a l , i t i s c o n v e n i e n t t o assume t h a t A i s symmetric p o s i t i v e d e f i n i t e i n o r d e r t o a n a l y z e t h e e f f e c t i v e n e s s o f Chebyshev a c c e l e r a t i o n . I t i s a polynomial a c c e l e r a t i o n , and i s termed a s e m i - i t e r a t i v e method w i t h r e s p e c t t o t h e c o r r e s p o n d i n g b a s i c i t e r a t i v e method. An o u t l i n e o f t h e t h e o r y behind Chebyshev a c c e l e r a t i o n i s as f o l l o w s : Suppose we have a b a s i c i t e r a t i v e method t h a t generates a sequence o f v e c t o r s v ( ~ ) ; i . e . we have a n i n i t i a l v e c t o r y ( O ) , and t h e n B

v(l) = !(')

-

V(n+l)

=

f

y(n)

_k e t c , such t h a t

+ _k.

A t t h e nth i t e r a t i o n , we w i s h t o choose c o n s t a n t s anI , i=O,l,..,n define

and t h e n

i n such a way t h a t t h e sequence o f v e c t o r s u(") tends towards t h e t r u e was solution o f Ay = b i n some o p t i m a l f a s h i o n . I f t h e i n i t i a l v e c t o r y"' i n f a c t u , t h e n c l e a r l y v ( " ) = u f o r a l l n 2 0, and so we r e q u i r e i n t h i s case t h a t u_cn) = 9 f o r aT1 n 2 0.

u

n

Hence we have t h e added c o n d i t i o n t h a t n

Since y =

1

i = o

1 anl 1 = 0

an12, we have

= 1, f o r each n 2 0.

Len CoZgan

642

i.e.

2")= Pn(B) E ( ' ) , "

C

pn ( x ) =

where t h e polynomial Pn(x) i s given by

anlxi

, for n

2

O.

1 - 0

The o n l y r e s t r i c t i o n we have imposed so f a r on t h e polynomials Pn(x) i s t h a t P n ( l ) = 1. L e t W be t h e synunetrization m a t r i x corresponding t o t h e m a t r i x B. Since WBW-l i s symmetric, then WP,(B)W-l i s a l s o symmetric, and hence IIPn(B)Ilw = p(Pn(B)) (see M a t r i x P r e l i m i n a r y ( d ) ) . Thus,

imp1 ies

A t t h e nth i t e r a t i o n , our o b j e c t i v e i s t h e r e f o r e t o choose t h e constants a,, so as t o minimise p(Pn(B)).

..,

L e t XI = m, A z , of Pn(B) a r e P n ( h l ) ,

A = M < 1 be t h e eigenvalues o f B. The eigenvalues Hence, f o r each n, we wish t o minimise

.., Pn?XN).

Instead o f t h i s p r o b l m , consider a s i m i l a r m i n i m i s a t i o n problem which d e t e r mines t h e v i r t u a l s p e c t r a l r a d i u s p(Pn(B)); we minimise

given t h a t P n ( l ) = 1. The s o l u t i o n o f t h i s problem i s c l a s s i c a l , and can be expressed i n terms o f Chebyshev polynomials, T,(x). Because o f a l i n e a r transformation, and t h e c o n d i t i o n P n ( l ) = 1, we f i n d t h a t

C l e a r l y t h i s expression could be used t o generate t h e c o e f f i c i e n t s but t h i s would be cumbersome, and, besides, t h e necessary storage o f a l l t h e vectors y C o ) , y ( l ) , , y ( n ) would be unmanageable. T h i s Chebyshev a c c e l e r a t i o n process can be reduced t o a second degree i t e r a t i o n by making use o f t h e second order recurrence r e l a t i o n f o r Chebyshev polynomials; v i z .

...

T ~ ( x )= 1 Ti(x) = x T,+,(x)

= 2x T,(x)

-

TnVl(x),

n

t

1.

Using t h i s recurrence, and a c e r t a i n amount o f algebra, t h e Chebyshev

Iterative Methods

643

a c c e l e r a t i o n process can be expressed s i m p l y as t h e second-degree s e m i - i t e r a t i v e method : L e t m, M < 1 be t h e s m a l l e s t and l a r g e s t eigenvalues r e s p e c t i v e l y o f t h e general i t e r a t i o n m a t r i x B. D e f i n e c o n s t a n t s

2 Then, g i v e n

2-M-m

,for n

an a r b i t r a r y

u _ ( ~ + ' )=

*

w,+,{y(B~~'")+k)

+

2

0 let

( l - y ) ~ ' ~ )+ } ( l - ~ ~ + ~ ) { ~ - " ,

where t h e sequence o f a c c e l e r a t i o n parameters w

n

W1

= 1;

Wp

=

l-c,a'

n+l

=-

1-&0'w~

for n

i s defined by 2

2.

I f m and M a r e c o r r e c t l y chosen, t h e n Chebyshev a c c e l e r a t i o n converges cons i d e r a b l y f a s t e r t h a n t h e b a s i c method. One p o s s i b l e disadvantage i s t h a t i t r e q u i r e s s t o r a g e o f two i t e r a t i o n v e c t o r s a t each stage. A f i r s t - d e g r e e f o r m does e x i s t , b u t t h e a c c e l e r a t i o n parameters t h e n tend t o become l a r g e , p o s s i b l y making t h e i t e r a t i o n u n s t a b l e due t o r o u n d i n g e r r o r .

I t c a n be shown t h a t t h e sequence o f which i s g i v e n b y

W ,

parameters tends towards a l i m i t , ,w

From t h e d e r i v a t i o n o f t h e a c c e l e r a t i o n procedure, and knowing t h e maximum v a l u e o f t h e Chebyshev p o l y n o m i a l s o v e r t h e r e l e v a n t i n t e r v a l , i t c a n t h e n be shown t h a t t h e v i r t u a l s p e c t r a l r a d i u s i s

O f course, p ( P m ( B ) )

l+r" --

5; 2rn'2

2r""

f o r l a r g e n. Consider t h e second model

problem, where M = -m = cosah, h = 1/40, when we a p p l y Chebyshev a c c e l e r a t i o n t o t h e b a s i c Jacobi method. Then y = 1, u = M, and so r = - 1-s in r h l + s i nah'

To a c h i e v e t h e accuracy p r e s c r i b e d e a r l i e r , we r e q u i r e n

zw

185 i t e r a t i o n s ,

(Using s u b r o u t i n e JSI o f Appendix E, which c a n b e c o n f i r m e d by experiment. w i t h t h e t r u e values o f M, rn i n p u t , and a s u i t a b l e s t o p p i n g t e s t , 187 i t e r a t i o n s were r e q u i r e d . ) The number o f i t e r a t i o n s i s p r o p o r t i o n a l t o l / h f o r s m a l l h, and t h e r e a r e a p p r o x i m a t e l y 1% t o 2 t i m e s as many as f o r t h e Successive O v e r r e l a x a t i o n method u s i n g I,.I,,~.

Len Colyan

64 4

I f we were t o use Symmetric Successive O v e r r e l a x a t i o n , and a p p l y Chebyshev a c c e l e r a t i o n , a s i m i l a r a n a l y s i s shows t h a t t h e number o f i t e r a t i o n s i s now p r o p o r t i o n a l t o 1/& f o r small h, and i n t h e case o f t h e second model problem, o n l y 26 i t e r a t i o n s a r e r e q u i r e d . T h i s seems t o be v e r y few, c o n s i d e r i n g we are, i n e f f e c t , l o o k i n g f o r 6 s i g n i f i c a n t f i g u r e s f o r t h e 1521 v a r i a b l e s .

However, t h i s needs t o be balanced by t a k i n g i n t o account t h e amount o f work i n v o l v e d i n each i t e r a t i o n , a l o n g w i t h t h e t o t a l t i m e and s t o r a g e r e q u i r e m e n t s . I t should be noted t h a t t h e f o r n w l a e f o r t h e c o n s t a n t s y and u, as w e l l as f o r t h e a c c e l e r a t i o n parameters 4, i n v o l v e t h e s m a l l e s t and l a r g e s t eigenvalues m, M o f t h e b a s i c i t e r a t i o n m a t r i x B.

Also, a s t o p p i n g t e s t t h a t i s f r e q u e n t l y used i n o r d e r t o t e r m i n a t e t h e Chebyshev a c c e l e r a t i o n procedure a t about t h e c o r r e c t i t e r a t i o n i s o f t h e form

Hence, i f t h e values G f m, M a r e n o t known, e s t i m a t e s w i l l have t o be found. Suppose these a r e i, M r e s p e c t i v e l y . I t can be shown t h a t Chebyshev a c c e l e r a t i o n i s n o t r e a l l y s e n s i t i v e t o e r r o r s i n v o l v e d i n i, e s p e c i a l l y i f ii 5 m, b u t i t i s v e r y s e n s i t i v e t o e r r o r s i n fl. I n p a r t i c u l a r , as fl approaches M from below, t h e r a t e o f converges improves d r a m a t i c a l l y . A d a p t i v e procedures have been developed which generate a c c u r a t e e s t i m a t e s f o r f 4 while the actual a c c e l e r a t i o n process i s i n o p e r a t i o n . I n t h g case o f Jacobi w i t h Chebyshev a c c e l e r a t i o n , i t i s o f t e n s u f f i c i e n t t o l e t m = -R, w h i l s t i f SSOR i s t h e b a s i c method, fi i s u s u a l l y l e f t as z e r o throughout. O f course, i n t h i s l a t t e r case, a good v a l u e w* o f t h e o v e r r e l a x a t i o n parameter a l s o needs t o be found a d a p t i v e l y . D e t a i l s o f these v a r i o u s a d a p t i v e procedures can be found i n Hageman and Young (1981); o r Grimes, K i n c a i d , MacGregor and Young (1978). A p o s s i b l e s u b r o u t i n e JSI t o implement t h e Jacobi method w i t h Chebyshev a c c e l e r a t i o n , and which adapts t o f i n d M, i s l i s t e d i n Appendix E. I t r e q u i r e s an i n i t i a l e s t i m a t e f o r M ( c a l l e d SRJ i n t h e s u b r o u t i n e ) , which i s o f t e n zero, and t h e n adapts. There a r e i n b u i l t r e s t r i c t i o n s t o p r e v e n t t h e parameters b e i n g changed t o o r e g u l a r l y .

[ B r i e f l y , i f n-s i s t h e number o f i t e r a t i o n s performed s i n c e t h e e s t i m a t e o f M was l a s t changed a t i t e r a t i o n number s, t h e n d e f i n e , u s i n g some o f t h e n o t a t i o n i n t h e p r e v i o u s two mentioned r e f e r e n c e s , qT =

2r( n- 1 ____

/ 2

l+r'n d

Y

=

Zl/(n--s)

=

(x

t

)

r/x)/(l

+r)

and t h e n a new v a l u e f o r

fl i s g i v e n by

I t e r a t i v e Methods

+(Fl + where

lii =

m,

645

m + Y(2-8-i)),

fl a r e t h e previous e s t i m a t e s , with t h e option of l e t t i n g a new

-8.

(All the above v a r i a b l e s a r e , in e f f e c t , t h e most straightforward way of solving t h e "Chebyshev equation". ) 1 For t h e f i r s t model problem, i t r e q u i r e s 22 i t e r a t i o n s i f m , M a r e known, and 26 i t e r a t i o n s i f they a r e t o be found adaptively (see Appendix E ) . S t a r t b g with M = 0, m = -1, and t h e r e a f t e r l e t t i n g til = -8, t h e subroutine changesM t o 0.7624 a f t e r 2 i t e r a t i o n s , and then t o 0.8087 a f t e r 6 i t e r a t i o n s . (Compare with true M = 0.8090.) For t h e second model problev, i t r e q u i r e s 187 and 218 i t e r a t i o n s r e s p e c t i v e l y . For more information on Chebyshev a c c e l e r a t i o n , and i t s comparison w i t h t h e Successive Overrelaxation i t e r a t i v e method, one can r e f e r t o Golub (1959); or Golub and Varga (1961). 6.

CONJUGATE GRADIENT ACCELERATION Hestenes and S t i e f e l (1952), i n "Methods of conjugate g r a d i e n t s f o r solving l i n e a r systems", J . Res.Nat.Bur.Standards, introduced t h e c l a s s i c conjugate g r a d i e n t concept. I t i s a p a r t i c u l a r c a s e of e a r l i e r c o n j u g a t e - d i r e c t i o n methods, and can be derived by modifying the " s t e e p e s t descent" o p t i m i s a t i o n technique. Moreover, f o r an Nth order system A x = b, i t would converge t h e o r e t i c a l l y i n a t most N i t e r a t i o n s . However, a f t e r a c e r t a i n amount of research f l u r r y t h a t immediately followed, l i t t l e p r a c t i c a l implementation of this conjugate g r a d i e n t method was used u n t i l r e c e n t l y . Renewed i n t e r e s t has been sparked by broadening the scope of the o r i g i n a l method of Hestenes and S t i e f e l . Rather than regarding the conjugate gradient method a s an i s o l a t e d one, i t can be considered a s possibly t h e simplest c a s e of a c c e l e r a t i n g a b a s i c i t e r a t i v e method using "conjugate g r a d i e n t a c c e l e r a t i o n " . More s p e c i f i c a l l y , consider t h e b a s i c s t a t i o n a r y Richardson i t e r a t i v e technique, where t h e s p l i t t i n g matrix Q andthe symmetrization matrix W ( i f A i s symmetric p o s i t i v e d e f i n i t e ) a r e both t h e i d e n t i t y matrix I , and t h e i t e r a t i o n matrix i s I-A. Then t h e o r i g i n a l conjugate g r a d i e n t method can be regarded a s conj u g a t e g r a d i e n t a c c e l e r a t i o n applied t o this Richardson i t e r a t i o n . Moreover, a conjugate g r a d i e n t a c c e l e r a t i o n can be applied under the same c o n d i t i o n s a s required f o r Chebyshev a c c e l e r a t i o n , and i n v a r i a b l y i t converges in fewer i t e r a t i o n s . Also, the a c c e l e r a t i o n parameters do not involve eigenvalues, and hence can be generated automatically, although t h e l a r g e s t eigenvalue M of t h e i t e r a t i o n matrix may eventually be required f o r a s u i t a b l e stopping t e s t simil a r t o t h e ones outlined i n e a r l i e r methods. The f a c t t h a t the method theore t i c a l l y terminates i n N i t e r a t i o n s i s unimportant f o r l a r g e s p a r s e problems, because any p r a c t i c a l i t e r a t i v e technique n e c e s s a r i l y must converge in conside r a b l y fewer steps than this. For s i m p l i c i t y , we s h a l l assume t h a t A i s symnetric p o s i t i v e d e f i n i t e . Actually, a s i s the s i t u a t i o n with Chebyshev a c c e l e r a t i o n , i t i s s u f f i c i e n t merely t o r e q u i r e t h a t a symmetrization matrix W e x i s t s f o r the i t e r a t i o n matrix 8. e.g. Richardson (W=I), Jacobi (W=D"). SSOR ( W = l / J z T D-%(D+wU)), under s u i t a b l e assumptions. Let c(n) = b i s arbitrary.

-

$n+l)

dn)

A be the residual vector a f t e r n i t e r a t i o n s . Suppose Then t h e o r i g i n a l conjugate g r a d i e n t method i s of t h e form =

dn) + amp( n ) ,

IJ'"

Len CoZgan

646

in'

where the a r e suitably chosen direction vectors, and a r e generated iteratively by the relation

y)

=

(")

The parameters an =

t o )=

r(O) *

.

8, a r e determined by

01,.

(r(") * , -r ( " )

(P n , (

, given

(n-1) 'n-1.P

,Ap' n , )

, or

equivalently

(p(")

,rCn)

(t")

I

and

Hestenes and Stiefel suggest that the f i r s t expressions for % and & are simpler to compute, b u t the second values for a, and 8, tend t o yield better results. I n practice, the iterations would proceed thus: given hence

y(O)

u(~+')

$n-l),

T

compute r("), followed by, in order,

p_

(n)

, an and

The validity of the classic conjugate gradient method l i e s in the f a c t that the residuals a r e orthogonal, i.e. ( $ I ) ,c ( 1 ) ) = 0 for i # j ; and the direction vectors a r e "A" or "conjugate", i.e. (A"P('),A"E(')) i # j.

=

0 f o r i # j , or equivalently, (p(",Ap(")

- orthogonal", =

0 for

The property of orthogonal residuals has as an immediate consequence the theoretical termination property mentioned e a r l i e r . NOW, suppose we decide to eliminate the direction vectors p(") from formulae above. This yields a second-degree i t e r a t i v e process for u_ involving the residual vectors; ,,(n+l) *

-

un+1

{yn+l"

where (fn)

-

'n+1

(tn) ,A$n))'

n,

+

U_(n)

1t

p+",,

U_(n-l),

)

and the parameters wn a r e defined by the sequence w1

=

1

As can be seen, t h i s has a form somewhat similar to the Chebyshev acceleration process. In f a c t , t h i s can be regarded as a special case of the following more general conjugate gradient acceleration method.

Iterative Methods

647

Consider a general l i n e a r s t a t i o n a r y i t e r a t i v e method o f the f i r s t degree t o solve A! = 5 , o f the form

corresponding t o a s p l i t t i n g m a t r i x Q, where B = I - Q-'A. L e t W be a nons i n g u l a r symmetrization m a t r i x such t h a t W(1-B)W-' i s symnetric p o s i t i v e d e f i n i t e . L e t $(") = B u(") + k - y(") be the pseudo-residual vector.

-

Then the conjugate g r a d i e n t a c c e l e r a t i o n procedure as applied t o t h i s basic method i s given by:

u_('),

Given an a r b i t r a r y

-

U(n+l)

n

for

-

w n + l {Yn+l (BY

=

wn+l

'yn+1

0 let

2

(n)

(n)

+k)

+ ( l - y n + l ) y ( " ) 1 + (l-llJ+l)U_(n-l)

+ u'"'1

+

(l-wn+l)u_(n-'),

where the sequences of a c c e l e r a t i o n parameters yn and w, yn+l -

1 -

. -6'

n,

w1 =

1

6( n l TWTWB_G( e)

(6.1) a r e defined by

and

TWTW$'n'

1

I t can be shown t h a t t h i s procedure minimises t h e CWTW(I-B)IWnorm o f t h e e r r o r vector E ( ~ )when compared w i t h any general polynomial a c c e l e r a t i o n procedureapplied t o t h e same basic i t e r a t i v e method. I n p a r t i c u l a r , i t would o f t e n r e q u i r e fewer i t e r a t i o n s than Chebyshev acceleration. I t i s worthwhile t o p o i n t o u t a few r e l e v a n t p r o p e r t i e s o f t h i s conjugate g r a d i e n t a c c e l e r a t i o n process: (a)

I f we consider t h e basic Richardson i t e r a t i o n w i t h Q = I , B = I-A, k = b_, 6'"' = r C n )W, = I etc., i t can be seen immediately t h a t t h e conjugate Gradienf a c c e l e r a t i o n reduces t o the second-degree form o f the c l a s s i c conjugate g r a d i e n t method.

(b)

U n l i k e Chebyshev acceleration, t h e symnetri4ation m a t r i x W i s needed i n t h i s process. More p r e c i s e l y , we r e q u i r e W W , and upon i n s p e c t i n g W f o r the Jacobi and SSOR methods, t h i s i s u s u a l l y simple o r can be e l i m i n a t e d by a s u i t a b l e p r e c o n d i t i o n i n g ( s c a l i n g ) o f the o r i g i n a l l i n e a r system. I t can be shown t h a t t h e pseudo-residuals a r e W orthogonal i.e.

(c)

( w & ~ ' , w & ' ' )= 0, f o r i f j .

The a c c e l e r a t i o n parameters do n o t r e q u i r e knowledge o f t h e eigenvalues However, s i n c e a stopping t e s t i s o f the f o r m

m(B), M ( B ) .

Len CoZgan

648

an estimate f o r M(B) must be found, Adaptive procedures t o f i n d M(B) have been developed. Again, d e t a i l s can be found i n Hageman and Young (1981); o r Grimes, Kincaid, MacGregory and Young (1978). (d)

Since B&"' is needed i n the evaluation o f t h e parameter y n + l , i t i s b e n e f i c i a l i f we can avoid having t o determine By(") a l s o . This can be achieved by using the second form (6.1) f o r the d e f i n i t i o n of U ( n + l ) , and generating the pseudo-residual vectors $ ( " ) by means o f thewsecond degree r e c u r s i o n

-

&n+l)

-

-

(n)

w ~ tyn+lB$ + ~

+

(l-Yn+l

1~+, ,+~- l (

I$'

n-1) 9

which can be derived q u i t e simply. (e)

There appears t o be more work necessary i n the c a l c u l a t i o n o f t h e a c c e l e r a t i o n parameters on each i t e r a t i o n o f the conjugate g r a d i e n t a c c e l e r a t i o n than f o r t h e corresponding Chebyshev a c c e l e r a t i o n . A1 so, more storage i s required. Nevertheless, f o r a wide range o f problems, conjugate gradient a c c e l e r a t i o n uses l e s s computer time, even when the eigenvalues f o r Chebyshev a c c e l e r a t i o n a r e known i n advance. However, counter-examples t o t h i s can be constructed.

(f)

Suppose t h e m a t r i x A i s n o t necessarily symmetric, b u t AT + A i s symmetric p o s i t i v e d e f i n i t e . We can then use a s p l i t t i n g m a t r i x Q = J~(A~+A)

and so o b t a i n a basic i t e r a t i v e method

gu_cn+l)

= %(A~-A)!(") t

b,

where Q i s symmetric p o s i t i v e d e f i n i t e . For general problems, Q may not be e a s i l y i n v e r t i b l e , and consequently i t may r e q u i r e an i n t e r n a l i t e r a t i o n procedure w i t h i n each basic i t e r a t i o n before u _ ( ~ + ' ) can be computed Furthermore, the eigenvalues o f the i t e r a t i o n m a t r i x B = 4Q-'(AT-A) a r e purely imaginary. However, a eneralized conjugate g r a d i e n t a c c e l e r a t i o n rocedure has been developed, :orresponding t o t h i s basic i t e r a t i v e i r o c e s s . [See Concus and Golub (1976). 1 L e t _6(n) = Q-'I-;(AT-A)g'"' + 5?1 - u_ ( n ) be t h e pseudo-residual. method can be expressed i n t h e form

Then t h i s

where the a c c e l e r a t i o n parameters a r e given by W1

=

1

Appendix F contains a s i m p l i f i e d subroutine JCG t o implement Conjugate Gradient Acceleration a p p l i e d t o the b a s i c Jacobi i t e r a t i v e method. As w i t h JSI i n the previous section, i t can adapt t o f i n d M, t h e l a r g e s t eigenvalue o f B , which i s required i n the stopping t e s t . To do t h i s , Hageman and Young (1931) have shown t h a t i t i s convenient t o form ( t h e o r e t i c a l l y ) a symnetric t r i d i a g o n a l m a t r i x , o f order n, a f t e r n

I t e r a t i v e Methods

649

i t e r a t i o n s , r e p r e s e n t e d by

Then t h e l a r g e s t e i g e n v a l u e o f t h i s t r i d i a g o n a l i s t h e new e s t i m a t e f o r M. The s u b r o u t i n e NEWTON i s s u f f i c i e n t t o accomplish t h i s . Once M i s found t o a r e q u i r e d accuracy, no f u r t h e r e s t i m a t e s a r e sought. Because M i s o n l y r e q u i r e d i n t h e s t o p p i n g t e s t and n o t i n t h e a c t u a l i t e r a t i o n s , t h e number o f i t e r a t i o n s f o r convergence o f t h e JCG method i s independent o f t h i s . The o u t p u t f o r t h e f i r s t model problem i n d i c a t e s t h a t o n l y 3 i t e r a t i o n s were necessary, and a t t h e t i m e o f stopping, t h e a d a p t i v e technique had achieved M = 0.809017 ( p r e c i s e l y ! ) . For t h e second model problem, i t r e q u i r e d o n l y 63 i t e r a t i o n s , d u r i n g w h i c h M was found c o r r e c t l y on t h e 27th i t e r a t i o n , and no f u r t h e r e s t i m a t e s f o r M were sought a f t e r t h a t .

7.

REDUCED SYSTEMS (RED-BLACK SYSTEMS) Consider a s i m p l e r v e r s i o n o f t h e model problems where we s h a l l use h = 4. Then A i s o f o r d e r 9. Suppose we number t h e i n t e r n a l nodes i n t h e " r e d - b l a c k " o r d e r i n g , as f o l l o w s : 1.1)

(0.1

(1.0)

'

0 0 0 0

4

0 0 0 0 . 4 0 0 0 . 0 4 0 0 . 0 0 4 0 . 0 0 0 4 .

1 1 1 0

0

0 -1 0 0 -1 0 -1 -1 -1 -1 0 -1 0 -1 -1

................... .

-1 -1 0

-1 -1 0 0 - i -1 -1 -1 0 < O 0 - 1 - 1 -

0

0

-1

1

. . . .

4 0 0 0

0 4 0 0

0 0 4 0

0 0 0 4

I n g e n e r a l , a system Au_ = b_ i s s a i d t o be a " r e d - b l a c k system" i f i t can b e p a r t i t i o n e d so t h a t

where D1, D2 a r e d i a g o n a l m a t r i c e s o f o r d e r N1, N2 r e s p e c t i v e l y such t h a t N1, N 2 a r e t h e numbers o f " r e d " and " b l a c k " p o i n t s , and N = N 1 + N P . I f we assume A i s symmetric p o s i t i v e d e f i n i t e , then K = G . C l e a r l y , whenever f i n i t e d i f f e r e n c e d i s c r e t i z a t i o n s a r e used o v e r r e c t a n g u l a r domains, we can arrange f o r t h e l a r g e sparse l i n e a r systems t o be " r e d - b l a c k ' ' systems.

Len Colgan

650

Then A!

can be r e w r i t t e n as

= Di

21

[K

hi

+ G Fz =

+ D Z Y Z = !z,

which suggests an i t e r a t i o n o f t h e form: Given

gZo), p e r f o r m =

f D 1 $1"' I

-G

+

$2")

(7.1)

S i n c e D;,' 0;' a r e t r i v i a l l y found, t h e sequence o f u_:"' v e c t o r s can be computed a l o n g w i t h y $ " ) . E l i m i n a t i n g 45") f r o m t h e two equations, t h e above i t e r a t i o n process can t h e o r e t i c a l l y be regarded as $2 n+')

t&")

= D;'KD;'G

+ ( ~ ; ' b ~- D;'KD;'~~).

The s i z e o f t h i s system i s N z , which i s about h a l f t h e s i z e o f t h e o r i g i n a l system. Using t h e general n o t a t i o n f o r b a s i c i t e r a t i v e methods, t h i s Reduced System can b e expressed i n t h e f o r m $zn+')

=

BRSt2")

+ ,,k

,

krhere BRS E D;'KDY1G, k I n p r a c t i c e however, t h e z (D;'bz - D;'KD;'b,). f o r m (7.1) would be us;% f o r the"computations d u r i n g each i t e r a t i o n , and i t i s c l e a r l y possible t o store the v e c t o r s i n t h e same computer l o c a t i o n s as Moreover, t h e amount o f work i n v o l v e d i n each s t e p o f t h e Reduced the System i t e r a t i v e method i s about t h e same as f o r t h e b a s i c J a c o b i method.

.

As

d2")

-r

y2, t h e n

Suppose W = 0.:

W BW , -'

+

gl.

Then = D:(D;K ' D;G ' )D;" =

D;"KD;'

=

CTC where C

GD"; = D;"GO;",

s i n c e KT = G.

Since W BW , '= CTC i s symmetric, then llBRsIlw= p ( B R S ) , and t h e eigenvalues a r e r e a l and non-negative. A l s o , u s i n g t h i s " r e d - b l a c k " o r d e r i n g , t h e o f ,B b a s i c Jacobi i t e r a t i v e method can be w r i t t e n i n t h e f o r m

and so

Now,

I t e r a t i v e Methods

651

i . e . t h e y 2 values a f t e r two Jacobi i t e r a t i o n s a r e t h e same as a f t e r one Reduced System i t e r a t i o n . I t can be shown t h a t p(B,2) = p ( B i S ) = p ( B R S ) . When t h e m a t r i x A i s s y m n e t r i c p o s i t i v e d e f i n i t e and has " P r o p e r t y A", we know t h a t p ( B J ) < 1. Hence, under these c o n d i t i o n s , p(BR,) = [ p ( B J ) l 2 < 1. Because o f t h e s e p r o p e r t i e s , we can c l e a r l y a c c e l e r a t e t h i s Reduced System method u s i n g e i t h e r Chebyshev a c c e l e r a t i o n o r c o n j u g a t e g r a d i e n t a c c e l e r a t i o n . I n a l l cases, t h e number o f i t e r a t i o n s i s a b o u t one h a l f as many as needed i n t h e c o r r e s p o n d i n g Jacobi process t o achieve t h e same accuracy. The g a i n s i n t i m e and s t o r a g e requirements a r e o f a s i m i l a r o r d e r . For example, t h e Chebyshev a c c e l e r a t i o n procedure a p p l i e d t o t h e reduced system c o u l d b e c a r r i e d o u t i n t h e f o l l o w i n g way: Given ui')

,

. . . , uy) - ,

f i n d din) u s i n g Dlu_(l"' the expression D,&"'

= -G

dZn)+ b l ,

y l n ) + b2

= -K

and t h e pseudo-residual

&"'

b y means o f

- y (2n ) .

Then t h e Chebyshev a c c e l e r a t i o n i s o f t h e f o r m

e,

A" 2 9 , u

where c o n s t a n t s y = g i v e n r e c u r s i v e l y by

=

and t h e a c c e l e r a t i o n parameters on a r e

T h i s i s , o f course, e q u i v a l e n t t o t h e general form, because m = 0 and M = Xf, To o b t a i n t h e p r e s c r i b e d accuracy f o r a s u i t a b l e s t o p p i n g t e s t m i g h t be

u,

Again, i f A, i s n o t known, a d a p t i v e processes have been developed. Other methods s i m i l a r t o a Reduced System w i t h Chebyshev a c c e l e r a t i o n a r e a l s o a v a i l a b l e , and a r e o b t a i n a b l e as m o d i f i c a t i o n s o f t h e fundamental " r e d - b l a c k " o r d e r i n g concept. 8.

ALTERNATING-DIRECTION IMPLICIT (ADI) METHODS Peaceman-Rachford Method There a r e a number o f d i f f e r e n t AD1 methods. Only t h e Peaceman-Rachford method w i 11 b e c o n s i d e r e d h e r e (see B i r k h o f f , Varga and Young (1962) ; o r Varga (1962) ) . Consider t h e s i m p l e r model problem corresponding t o h = +, and t h e r e s u l t i n g l i n e a r system A! = b. Assume we a r e u s i n g t h e n a t u r a l o r d e r i n g .

Len Colgan

652

Then t h e m a t r i x A can b e w r i t t e n as A = Al 0 1 2 -1 0 1-1 2 - 1 0 0 -1 2 0 0 2 A l = i : 0 -1 0 0 0 0 0 0 0 0

;

0 0

0 0

0

0 1 2 1

0

0

0 0 0

0 0 0

-1 2

0

0

0

0

+ A Z , where 0 0 0

0 0 0

2 - 1

0

L! 0

0

01 OI 0

and 0 - 1 0 0 0 0 0’ 0 0 - 1 0 0 0 0 0 2 0 0 - 1 0 0 0 0 2 0 0 - 1 0 0 0 0 2 0 0 - 1 0 ‘ -1 0 0 2 0 0 - 1 0 - 1 0 0 2 0 0 0 0 - 1 0 0 2 0 0 0 0 - 1 0 0 2

l o

The e l l i p t i c p a r t i a l d i f f e r e n t i a l e q u a t i o n f r o m which t h i s arose was

when u(x,y) = 0 on t h e boundary o f t h e u n i t square 0

5

x

5

1, 0

5

y 5 1.

A f t e r u s i n g t h e f i n i t e d i f f e r e n c e d i s c r e t i z a t i o n and m u l t i p l y i n g by ( - h 2 ) , we can see t h a t A l i s d e r i v e d f r o m t h e a2u/ax2 term, and A; i s d e r i v e d f r o m a2u/ay2. I n a more general s i t u a t i o n , suppose we d i s c r e t i z e

on, say, a r e c t a n g u l a r r e g i o n . L e t t h e m a t r i c e s H, V, a

ax

+

c

N, ax

1respectively

c o n t a i n t h e terms corresponding t o

a2u + d au b 7 aY aY

and e u,

I

a f t e r m u l t i p l y i n g by ( - h 2 ) . Then A = H + V + 1. Assume i s a non-negative diagonal m a t r i x , and H, V a r e symmetric p o s i t i v e d e f i n i t e m a t r i c e s w i t h nonp o s i t i v e o f f - d i a g o n a l elements. Let

A1

= H +

%I, A Z = V + %,I, so t h a t A

I t i s w o r t h n o t i n g t h a t most o f t h e r i g o r o u s j u s t i f i c a t i o n o f t h e f o l l o w i n g method r e l i e s on H, V b e i n g s y m n e t r i c p o s i t i v e d e f i n i t e , being a m u l t i p l e o f t h e i d e n t i t y m a t r i x , and HV = VH. T h i s t h e n i m p l i e s t h a t A 1 and A P a r e symmetric p o s i t i v e d e f i n i t e , and AIAz = A2A1.

I

The e q u a t i o n Au_ =

i s c l e a r l y e q u i v a l e n t t o each o f t h e e q u a t i o n s

I t e r a t i v e Methods

(A1 +

un+lI)!

=

b-

(A2

-

o,+~

+

I)!

=

b -

(Ai

-

o,+~ I)!,

653

I)!

and (A2

f o r any

w"+~.

This suggests t h e Peaceman-Rachford AD1 method:

1

b-

( A , + W , + ~ I ) U _ ( " +=~ )b

(A2

( A 2 + O , + ~ I ) U _ ( " + ~ =)

(A,

-

W~+~I)U_(~) W,+~I)U("+~).

The % a r e a c c e l e r a t i o n parameters. Since A, + u,,+~I i s u s u a l l y a t r i d i a g o n a l matrix, and A 2 + w , , + ~ Ican be made t r i d i a g o n a l a f t e r a s u i t a b l e permutation, each of t h e i m p l i c i t s t e p s i n t h e i t e r a t i o n s can be c a r r i e d o u t u s i n g t h e wellknown Gaussian elimination a l g o r i t h m f o r tridiagonal systems. In e f f e c t , t h e f i r s t equation i s equivalent t o horizontal sweeps, whereas t h e second equation should be regarded a s v e r t i c a l sweeps; hence t h e category "A1 t e r n a t i n g Direction I m p l i c i t Methods". T h e a u x i l i a r y v e c t o r u _ ( " + ~ )i s not r e t a i n e d from one complete i t e r a t i o n t o t h e next. In p r a c t i c e , i t i s e f f i c i e n t t o use a c c e l e r a t i o n parameters in t h e manner y , .., 4 , wl, q , .., 3 ,o ~ w,z , .., i y , f o r a s u i t a b l e k . I f k=l, the process i s s t a t i o n a r y . rhe Peaceman-Rachford AD1 method can t h e o r e t i c a l l y be expressed i n t h e form

...

%,

-

U_("+l) -

B , + l ~( n ) +

kn+l.

where t h e i t e r a t i o n matrix B,+l B,+l

=

(A2

+

11 2

0,

i s given by

W , + ~ I ) - ' ( A I - U , + ~ I ) ( A+~o,+,I)-'(Az

-

W,+~I).

Fi

If p ( i i ) i s t h e s p e c t r a l r a d i u s of t h i s matrix product, then the average r a t e of convergence a f t e r n i t e r a t i o n s i s R, = - 1 ln(P(l!l B, 1). Although optimal a c c e l e r a t i o n parameters O, t o minimise t h i s average r a t e of convergence a r e not a v a i l a b l e , good ones, t h a t a r e simple t o c a l c u l a t e , have been found t o work q u i t e s a t i s f a c t o r i l y . Let a , 13 be found so t h a t 0 < s p 5 B and 0 < a I A s 6 f o r a l l eigenvalues of A l and eigenvalues A of A z . S u c h bounds on the eigenvalues can e a s i l y be determined using t h e r e l e v a n t t h e o r i e s . Then two commonly used parameters a r e t h e following: (1) Peaceman-Rachford Parameters Find the smalle:;ilr,t,egfr k such t h a t (&?-l)zk s a / B . Then, f o r i = l , Z , . . , k , l e t wi = a ( a / ~ ) These parameters a r e c l e a r l y n o t evenly spaced, but a r e geometrically spaced. ( I n t h e s p e c i a l case of t h e s t a t i o n a r y Peaceman-Rachford method, k = 1, and so o = / a @ . )

.

I t can then be shown t h a t

Len CoZgan

654

say, and so t h e average r a t e o f convergence

%

1

-

2

6.

The c o n d i t i o n (J2-1)2k s a/B i s used t o determine k because i t i s a consequence o f t h e r e s u l t t h a t Rk i s maximised when IS = J2-1. For example, c o n s i d e r t h e second model problem. Hence, = 4 s i n 2 n h / 2 and B = 4 cos1?h/2.

~1

E

a

i:

tan2

9,

and so we r e q u i r e (42-1)” Then we o b t a i n 0.01385,

WI z

(2)

~p

2

s tan

0.06986,

Wachspress Parameters Find t h e smallest integer k

6

a, = a(,)

(I-l)/(k-l)

2

F. For h 03

1

= 1/40,

0.35245,

10.006165,

, f o r i=1,2

w5 =

= 3.9938.

W?

1

04

this yields k

0.03110,

=

4.

1.7781.

2

2 such t h a t (J2-1)z‘k-”

s

.;

Then

,...,k .

For t h e second model problem w i t h h = 1/40, obtain WI = a

I t can b e shown t h a t

~g

t h i s y i e l d s k = 5.

0.1569,

~4

Then we

10.7916,

Numerical experiments suggest t h a t t h e Wachspress parameters a r e s u p e r i o r t o t h e Peaceman-Rachford parameters by a f a c t o r o f about 2, p r o v i d e d t h a t k i s determined as d e s c r i b e d . For b o t h s e t s o f parameters, t h e number o f i t e r a t i o n s r e q u i r e d t o o b t a i n t h e s p e c i f i e d accuracy v a r i e s a c c o r d i n g t o l n ( l / h ) as t h e mesh-size h tends t o zero. T h i s i s t o be compared w i t h l / h f o r b o t h (i) (ii)

t h e Successive O v e r r e l a x a t i o n method, and t h e Jacobi method w i t h Chebyshev a c c e l e r a t i o n

and w i t h 1/Jb f o r t h e a c c e l e r a t e d SSOR method. However, these o t h e r methods a r e perhaps a p p l i c a b l e o v e r a w i d e r range o f problems, and a l s o t h e e x t r a computations i n v o l v e d i n t h e i m p l i c i t s o l u t i o n s needed i n each i t e r a t i o n o f t h e Peaceman-Rachford method a r e a p p r e c i a b l e . Hence, t h i s AD1 method i s r e a l l y o n l y advantageous p r o v i d e d h i s v e r y small. I n f a c t , e x p e r i m e n t a t i o n suggests t h a t , p r o v i d i n g t h e A01 method i s a p p l i c a b l e , t h e r e e x i s t s a v a l u e o f h, say h*, such t h a t f o r h > h* t h e SOR-type methods a r e p r e f e r a b l e , b u t f o r h < h*, t h e Peaceman-Rachford method uses l e s s computer time. 9.

COMPARISON OF ITERATIVE METHODS Comments made a b o u t t h e convergence r a t e s o f t h e v a r i o u s i t e r a t i v e methods d e s c r i b e d i n t h e preceding s e c t i o n s g i v e a reasonable i n d i c a t i o n as t o t h e s u i t a b i l i t y o f each. O f course, t h e amount o f t i m e and machine s t o r a g e space necessary t o perform t h e i t e r a t i o n s themselves must be taken i n t o account. A number o f references i n t h e b i b l i o g r a p h y q u o t e numerical r e s u l t s f o r a range o f problems, and t h e r e i s a reasonable degree o f u n i f o r m i t y i n t h e i r conclusi o n s . For example, see B i r k h o f f , Varga and Young (1962); o r Westlake (1968).

I t e r a t i v e Methods

655

Most o f t h e i t e r a t i v e methods c o n s i d e r e d r e l y on t h e m a t r i x A b e i n g symmetric p o s i t i v e d e f i n i t e , o r else require ( t h e o r e t i c a l l y ) the existence o f a s y m m e t r i z a t i o n m a t r i x W. I f A i s n o t symmetric, newer methods t h a t a r e b a s i c a l l y e x t e n s i o n s o r v a r i a t i o n s o f t h e c l a s s i c methods have been f o r m u l a t e d . D e t a i l s can be found i n Hageman and Young (1981), o r M a n t e u f f e l (1975), o r Axelsson and Gustafsson (1977), o r Widlund ( 1 9 7 8 ) .

O f t h e methods o u t l i n e d i n t h e p r e c e d i n g s e c t i o n s , those t h a t can s e r i o u s l y be c o n s i d e r e d f o r implementation t o s o l v e l a r g e sparse l i n e a r systems a r e (1)

Successive O v e r r e l a x a t i o n method, which may need t o b e a d a p t i v e i f wont i s unknown.

(2)

Chebyshev a c c e l e r a t i o n o f t h e b a s i c Jacobi o r SSOR methods. Again, a d a p t i v e techniques w i l l be needed i f t h e eigenvalues o r r e l a x a t i o n parameters a r e n o t known beforehand.

(3)

Conjugate g r a d i e n t a c c e l e r a t i o n o f t h e Jacobi o r SSOR methods.

(4)

The reduced system concept, u s i n g a " r e d - b l a c k ' ' o r d e r i n g , i n c o n j u n c t i o n w i t h e i t h e r Chebyshev o r c o n j u g a t e g r a d i e n t a c c e l e r a t i o n .

(5)

Peaceman-Rachford A1 t e r n a t i n g D i r e c t i o n method.

I f t h e reduced system methods a r e a p p l i c a b l e , numerical r e s u l t s i n d i c a t e t h a t t h e y f r e q u e n t l y use fewer i t e r a t i o n s , l e s s t i m e and l e s s s t o r a g e than t h e o t h e r techniques. I n p a r t i c u l a r , t h e reduced system combined w i t h c o n j u g a t e g r a d i e n t a c c e l e r a t i o n has o f t e n proven t o be o p t i m a l i n each o f these c l a s s i f i c a t i o n s . The SSOR method, w i t h e i t h e r Chebyshev o r c o n j u g a t e g r a d i e n t a c c e l e r a t i o n a l s o produces f a v o u r a b l e r e s u l t s , and covers a w i d e r c l a s s o f problems. However, t h e a c t u a l programming i m p l e m e n t a t i o n i s more c o m p l i c a t e d . The Successive O v e r r e l a x a t i o n method i s p r e f e r r e d by many because o f i t s s i m p l i c i t y i n t h i s r e g a r d , b u t u s u a l l y c o n s i d e r a b l y more i t e r a t i o n s , and somewhat more e x e c u t i o n time, a r e r e q u i r e d . The Peaceman-Rachford method i s r e a l l y o n l y s u p e r i o r when t h e r e g i o n i s r e c t a n g u l a r ( o r s i m i l a r ) , t h e mesh-size i s s m a l l , and k > 1. I t i s a d v i s a b l e t o s e l e c t a v a l u e o f k, t h e number o f i t e r a t i o n parameters, which i s t o o l a r g e than t o use one which i s t o o s m a l l . Savings have been found by a c t u a l l y choosing a v a l u e o f k l a r g e r t h a n t h e ones produced by t h e f o r m u l a e d e s c r i b e d i n t h a t s e c t i o n .

Len Colgan

65 6

APPENDIX A LIST OF SYEIBOLS : c o e f f i c i e n t matrix of l i n e a r system : A = Al t A,

i n Peaceman-Rachford Plethod

: general i t e r a t i o n matrix : Jacobi i t e r a t i o n matrix

: Gauss-Seidel i t e r a t i o n matrix : Successive Overrelaxation i t e r a t i o n matrix

: Backward Successive Overrelaxation i t e r a t i o n matrix : Synmetric Successive Overrelaxation i t e r a t i o n matrix : Reduced system i t e r a t i o n matrix

: Peaceman-Rachford i t e r a t i o n matrices : diagonal p a r t o f A : submatrix of A i n "red-black" ordering : estimate f o r s p e c t r a l radius of BsoR

: submatrix of A i n "red-black" ordering : lower t r i a n g u l a r p a r t of A : l a r g e s t eigenvalue of B : estimate of El : order of matrix A

: transformed Chebyshev polynomials : general s p l i t t i n g matrix : v a r i a b l e used i n solving Chebyshev equation : variable used i n solving Chebyshev equation : asymptotic r a t e of convergence : average r a t e of convergence a f t e r n i t e r a t i o n s : r a t i o o f consecutive pseudo-residuals i n SOR method

: Chebyshev polynomials : upper t r i a n g u l a r p a r t of A : symnetrization matrix

: v a r i a b l e used i n solving Chebyshev equation : variable used i n solving Chebyshev equation : v a r i a b l e used i n s o l v i n g Chebyshev equation : right-hand s i d e of l i n e a r system

: mesh s i z e : vector i n general i t e r a t i o n : smallest eigenvalue of B : estimate of m : the number of the i t e r a t i o n

Iterative Methods

E(

")

: d i r e c t i o n v e c t o r s i n c o n j u g a t e g r a d i e n t method

r

: parameter i n Chebyshev a c c e l e r a t i o n

.U.

: s o l u t i o n o f l i n e a r system : estimate o f u a f t e r n i t e r a t i o n s

-r( n ) u( n)

I

l i

: residual vector a f t e r

n iterations

: lower e i g e n v a l u e bound i n Peaceman-Rachford method

(i

: c o n j u g a t e g r a d i e n t parameters

b

: upper e i g e n v a l u e bound i n Peaceman-Rachford method

B"

: c o n j u g a t e g r a d i e n t parameters

B Y

: Chebyshev a c c e l e r a t i o n c o n s t a n t

-

yn

- n)

: upper bound f o r (p(D"Ld'U)

: c o n j u g a t e g r a d i e n t a c c e l e r a t i o n parameters

(p(

: pseudo-residual v e c t o r

c x

: e r r o r vector

< n)

XJ

LJ

i n SSOR method

: general e i g e n v a l u e : l a r g e s t p o s i t i v e e i g e n v a l u e o f BJ : upper bound f o r XJ, i n SSOR method : eigenvalue o f BsoR

u

: s p e c t r a l r a d i u s o f BsoR

p(B)

: spectral radius o f B

0

: Chebyshev a c c e l e r a t i o n constant

w

: r e l a x a t i o n parameter

Ill

OPl

: o p t i m a l value f o r w i n SOR method

o*

: good v a l u e f o r w i n SSOR method

[I)

: Chebyshev a c c e l e r a t i o n parameters

n

: c o n j u g a t e g r a d i e n t a c c e l e r a t i o n parameters : Peaceman- Rac hf o r d acce 1e r a t i on parame t e r s

om

: l i m i t i n g value o f o

i n Chebyshev a c c e l e r a t i o n

657

Len Colgan

658

APPENDIX B SUBROUTINE JAC C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C

C C C C C

20

SUBROUTINE TO SOLVE LARGE SPARSE L I N E A R SYSTEMS USING JACOB1 I T E R A T I O N

THE COEFFICIENT MATRIX IS STORED IN 'UNSYMMETRIC SPARSE FORM SUBROUTINE JAC( A,JA,ISTART,N ,NP1 ,IADIM,B ,U ,ZETA, tSR, IADAPT, ITMAX,NUMITS ,V,D) GLOSSARY OF TERMS ~NON-ZERO ELEMENTS OF SPARSE COEFF. MATRIX (INPUT) A( ) :COLUMN POSITIONS OF ELEMENTS I N A ( I N P U T ) JA( I S T A R T ( I) :POSITION I N A & J A OF THE F I R S T ENTRY FROM :ROW I . LAST ENTRY IS I A D I M + l ( I N P U T ) :NO. OF EQUATIONS ( I N P U T ) N :N+1 ( I N P U T ) NP1 :DIMENSION OF A & J A ( I N P U T ) IADIM :RIGHT HAND S I D E OF EQUATIONS ( I N P U T ) B( 1 :SOLUTION. ON I N P U T CONTAINS I N I T I A L GUESS U( 1 :WHICH MAY B E ZERO VECTOR (INPUT,OUTPUT) :STOPPING TEST. RELATIVE ERRORcZETA ( I N P U T ) ZETA :SPECTRAL RADIUS OF I T E R A T I O N MATRIX SR :(MUST B E I N P U T I F ADAPT=O) IADAPT :=0 DOES NOT ADAPT TO ESTIMATE SR :=1 WILL ADAPT ( I N P U T ) :MAX NO. OF I T E R A T I O N S ( I N P U T ) I TMAX :NO. OF ITERATIONS PERFORMED (OUTPUT) NUMITS :TEMPORARY STORAGE WORK AREA OF S I Z E N V( 1 :DIAGONAL STORAGE WORK AREA D F S I Z E N D( :PREVIOUS PSEUDO-RESIDUAL NORM PSPREV PSNDRM :CURRENT PSEUDO-RESIDUAL NORM :CURRENT SOLUTION NORM VNORM :STOPPING TEST TEST REAL A( I A D I M ) ,B(N) ,U(N) ,ZETA,SR INTEGER J A ( I A D I M .) ..ISTART( .NP1). .N I A D I M , IADAPT, _ ,NP1, . +NUMITS ,ITMAX REAL V(N) ,PSPREV,PSNDRM,D(N) ,VNDRM,SUM,TEST

****

SET UP DIAGONAL ELEMENTS F I N D O R I G I N A L PSEUDO-RESIDUAL NORM CALCULATE F I R S T ITERATES

****

PSPREV=O. DO 10 I = l . N SUM=B( I ) ' DO 20 J=ISTART( I ) ,ISTART( I+l)-l IF(JA(J) .EQ.I)THEN D( I)=A(. J). ELSE SUM=SUM-A( J)*u(JA(J ) ) ENDIF CONTINUE

Iterative Methods

10

30 C C C C C 40

60

50

70 C C C C C C

C C C

659

V ( I)=SUM/D( I ) PSPREV=PSPREV+( V( I)-U( I ) )**2 CONTINUE PSPREV=SQRT( PSPREV) NUMITS=l DO 30 I = I , N U( I)=v( I) CONTINUE

****

PERFORM THE NEXT I T E R A T I O N CALCULATE THE NORMS OF THE ESTIMATE AND THE PSEUDO-RESIDUAL

****

NUMITS=NUMITS+l IF( NUMITS.GT.ITMAX)RETURN PSNORM=O. VNORM=O. DO 50 I = I , N SUM=B( I) 00 60 J = I S T A R T ( I),ISTAR I F ( J A ( J ) .NE.I)SUM=SUM=A CONTINUE V( I ) = S U M / D ( I ) PSNORM=PSNORM+( V( I ) - U ( I ) **2 VNORM=VNORM+V( I ) **2 CONTINUE PSNORM=SQRT( PSNORM) VNORM=SQRT( VNORM) DO 70 I = l . N U( I ) =v( I ) . CONTINUE

****

STOPPING TEST I F IADAPT=I,USE R A T I O OF CONSECUTIVE PSEUDO-RESIDUALS AS AN ESTIMATE FOR THE SPECTRAL RADIUS OF THE I T E R A T I O N MATRIX

****

I F ( I A D A P T . EQ. l)SR=PSNORM/PSPREV TEST=SR*PSNORM/( 1.-SR)/VNORM I F ( A B S ( T E S T ) .LT.ZETA)RETURN

****

PERFORM ANOTHER I T E R A T I O N

****

PSPREV=PSNORM GOT040 END

OUTPUT

The subroutine JAC, w i t h a p p r o p r i a t e P R I N T statements, y C o ’ = Q , and a c a l l i n g PROGRAM t o implement the s o l u t i o n of the f i r s t model problem ( N = 1 6 ) , y i e l d s t h e following o u t p u t : (showing t h e number o f t h e i t e r a t i o n , NUMITS; progressive estimate of the s p e c t r a l radius, SR; and t h e 16 v a r i a b l e values, U1, . . . , UI6. Only some o f the i t e r a t i o n s have been e x h i b i t e d . ) NUMITS = 1 .0100000 .0100000 .0100000 .0100000 .0100000 .0100000 SR = 0 (has n o t changed)

.0100000 .0100000

.0100000 .0100000

.0100000 .0100000

.0100000 .0100000

.0100000 .0100000

Len C o l g m

660

NUMITS = 2 .0150000 .0175000 .0175000 .0200000 SR = 770552

.0175000

.0200000 .0175000

.0150000

.0175000 .0150000

.0200000 .0175000

.0200000 .0175000

.0175000 .0150000

NUMITS = 3 .0187500 .0231250 .0231250 .0287500 SR = .805001

,0231250 .0287500

.0187500 .0231250

.0231250 .0187500

.0287500 .0231250

.0287500 .0231250

.0231250 .0187500

NUMITS = 4 .0215625 .0276563 .0276563 .0359375 SR = .a08623

.0276563 .0359375

.0215625 .0276563

.0276563 .0215625

.0359375 .0276563

.0359375 .0276563

.0276563 .0215625

NUMITS = 5 .0238281 .0312891 .0312891 .0417969 SR = .808978

.0312891 .0417969

.0238281 .0312891

.0312891 .0238281

.0417969 .0312891

.0417969 .0312891

.0312891 .0238281

NUMITS = 10 .0300404 .0413385 .0413385 .0580456 SR = .809017

.0413385 .0580456

.0300404 ,0413385

.0413385 .0300404

.0580456 .0413385

.0580456 .0413385

.0413385 .0300404

NUMITS = 15 ,0321921 .0448201 .0448201 .0636789 SR = .809017

.0448201 .0636789

.0321921 .0448201

.0448201 .0321921

.0636789 .0448201

.0636789 .0448201

.0448201 .0321921

NUMITS = 20 .0329378 .0460267 .0460267 .0656312 SR = .809017

.0460267 .0656312

.0329378 .0460267

.0460267 ,0329378

.0656312 .0460267

.0656312 .0460267

.0460267 .0329378

NUMITS = 60 ,0333333 .0466665 .0466665 .0666665 SR = .809017

.0466665 .0666665

.0333333 .0466665

.0466665 .0333333

.0666665 .0466665

.0666665 .0466665

.0466665 .0333333

NUMITS = 66 ,0333333 .0466666 .0466666 .0666666 SR = .809017

.0466666 .0666666

.0333333 .0466666

.(I466666 .0333333

.0666666 .0466666

.(I666666 .0466666

.0466666 .0333333

I t e r a t i v e Methods

APPENDIX C SUBROUTINE GS C C

SUBROUTINE TO SOLVE LARGE SPARSE L I N E A R SYSTEMS USING GAUSS-SEIDEL I T E R A T I O N

C C C C

THE COEFFICIENT SPARSE FORM

C

C C C C C C C C C C C C C C C C C C C C C C C C C C C C C

C C C C C

20

GLOSSARY

OF

MATRIX

IS STORED NI

'UNSYMMETRIC

TERMS

:NON-ZERO ELEMENTS OF SPARSE COEFF. MATRIX ( I N P U T ) A( ) :COLUMN POSITIONS OF ELEMENTS I N A ( I N P U T ) JA( ) I S T A R T ( ) : P O S I T I O N I N A & J A OF THE F I R S T ENTRY FROM :ROW I. L A S T ENTRY IS I A D I M + l ( I N P U T ) N :NO. OF EQUATIONS ( I N P U T ) NP1 :N+1 ( I N P U T ) IADIM :DIMENSION OF A & J A ( I N P U T ) :RIGHT HAND S I D E OF EQUATIONS ( I N P U T ) B( ) :SOLUTION. ON INPUT CONTAINS I N I T I A L GUESS u( 1 :WHICH MAY B E ZERO VECTOR (INPUT,OUTPUT) ZETA :STOPPING TEST. R E L A T I V E ERROR
****

SET UP DIAGONAL ELEMENTS F I N D O R I G I N A L PSEUDO-RESIDUAL NORM CALCULATE F I R S T ITERATES

****

PSPREV=O. DO 10 I = l , N SUM=B( I ) DO 20 J = I S T A R T ( I) ,ISTART( I+l)-I I F ( J A ( J) .EQ. 1)THEN D( I ) = A ( J ) ELSE SUM=SUM-A(J)*U(JA(J)) ENDIF CONTINUE

661

Len CoZgan

662

10 C C C C C 40

60

50 C C C C C C

C C C

TEMP=SUM/D( 1) PSPREV=PSPREV+(TEMP-U(I ) ) * * 2 U( I)=TEMP CONTINUE PSPREV=SQRT(PSPREV) NUMITS=l

****

PERFORM THE NEXT ITERATION CALCULATE THE NORMS OF THE ESTIMATE AND THE PSEUDO-RESIDUAL

****

NUMITS=NUMITS+l I F ( NUM I T S . GT I TMAX RETURN PSNORM=O. UNORM=O. DO 50 I=l,N SUM=B ( I) 00 60 J=ISTART( I ) ,ISTART( I + l ) - 1 I F ( JA ( 3 ) .NE .I ) SUM=SUM-A( J)*U( JA( 3 ) ) CONTINUE TEMP=SUM/D( I ) PSNORM=PSNORM+(TEMP-U( I ) ) * * 2 UNORM=UNORM+TEMP**2 U( I)=TEMP CONTINUE PSNORM=SQRT(PSNORM) UNORM=SQRT(UNORM)

.

****

STOPPING TEST I F IAOAPT=l,USE RATIO OF CONSECUTIVE PSEUDO-RESIDUALS AS AN ESTIMATE FOR THE SPECTRAL RAOIUS OF THE ITERATION MATRIX

****

I F ( 1ADAPT.EQ. l)SR=PSNORM/PSPREV TEST=SR*PSNORM/( 1.-SR)/UNORM IF(ABS(TEST) .LT .ZETA)RETURN

****

PERFORM ANOTHER ITERATION

****

PSPREV=PSNORM GOT040 EN0

OUTPUT

The subroutine GS yields the following output, under the same conditions as in Appendix 8. I t i s found t h a t a b o u t half as many i t e r a t i o n s are required f o r the f i r s t model problem being considered. NUM ITS = 1 .0100000 .0125000 .0131250 .0173438 SR=0.000000

.0131250 .0186719

.0132813 .0190820

.0125000 .0132813

,0162500 .0176563

.0173438 .0190820

.0176563 .0195410

NUMITS= 2 .0162500 .0214063 .0230078 .0321777 SR=.684742

.0230078 .0356299

.0201660 .0307568

.0214063 .0201660

.0293750 .0278564

.0321777 .0307568

.0278564 .0253784

I t e r a t i v e Methods

663

NUMITS=3 .@07031 .0282715 .0301538 .0434662 SR=.700422

.0301538 .0471115

.(I245026 .0367928

.0282715 .0245026

.(I402246 .0346814

.0434662 .0367928

,0346814 .0283964

NUMITS=4 .0241357 .(I336285 .0353993 .0514349 SR=.696357

.0353993 .0541139

.0275202 .0403618

.0336285 .0275202

.0485474 .0389370

.0514349 .0403618

.0389370 .0301809

NUMITS=5 .0268143 .0376902 .0391613 .0566937 SR=.686504

.0391613 .0585277

.0295246 .0425884

.0376902 .0295246

.0545626 .0416450

.0566937 .0425884

.0416450 .0312942

NUMITS=10 .0324713 .(I455369 .0457519 .0654688 SR=.655921

.0457519 .0656973

.(I328758 .0461820

.0455369 .0328758

.0651866 .0460676

.(I654688 .0461820

.0460676 .0330910

NUMI TS= 15 .0332296 .0465308 .0465568 .0665228 SR=. 654544

.0465568 .0665503

.0332784 .0466085

.0465308 .(I332784

.0664888 .0465947

.0665228 .0466085

.0465947 .0333042

NUMITS=20 .0333209 .0466503 .0466535 .(I666494 SR=.654509

.0466535 .OM6527

.0333267 .0466597

.0466503 .0333267

.0666453 .0466580

.0666494 .0466597

.0466580 .0333298

NUMI TS = 2 5 .0333318 ,0466647 .0466651 .0666646 SR=.654509

.0466651 .0666650

.0333325 .0466658

.0466647 .0333325

.0666641 .0466656

.(I666646 .0466658

.0466656 .0333329

NUMITS=30 .0333332 .0466664 .0466665 .Of566664 SR=.654508

.0466665 .0666665

.0333332 .0466666

.0466664 .0333332

.0666664 .0466665

.0666664 .0466666

.0466665 ,0333333

NUM ITS.34 .0333333 .0466666 .0466666 .0666666 SR=.654508

.0466666 .0666666

.0333333 .0466666

.0466666 .0333333

.OM6666 .0466666

,0666666 .0466666

.0466666 .0333333

Len CoZgan

684

APPENDIX 0 SUBROUTINE SOR C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C

C C C C C

SUBROUTINE TO SOLVE LARGE SPARSE LINEAR SYSTEMS USING SOR I T E R A T I O N WRITTEN BY L.COLGAN 1983 SOLVES SYSTEMS OF N SPARSE L I N E A R EQUATIONS WHERE THE COEFFICIENT MATRIX I S STORED I N UNSYMMETRIC SPARSE FORM SUBROUTINE SOR( A ,JA, I S T A R T ,N ,NP1 ,IAOIM,B ,U ,ZETA, +SR,OMEGA, IAOAPT, ITMAX ,NUMITS ,0) GLOSSARY OF TERMS

SR OMEGA IAOAPT ITMAX NUMITS

O( 1

PSPREV PSNORM UNORM TEST

H

TEMP SRJ T(I) IEST NEWITS CHANGE

:NON-ZERO ELEMENTS OF SPARSE COEFF. MATRIX ( I N P U T ) :COLUMN P O S I T I O N S OF ELEMENTS I N A ( I N P U T ) I ) : P O S I T I O N I N A & J A OF THE F I R S T ENTRY FROM :ROW I . LAST ENTRY IS I A O I M t l ( I N P U T ) :NO. OF EQUATIONS ( I N P U T ) : N t l (INPUT) :DIMENSIONS OF A & J A ( I N P U T ) :RIGHT HANO S I D E OF EQUATIONS ( I N P U T ) :SOLUTION. ON I N P U T CONTAINS I N I T I A L GUESS :WHICH MAY BE ZERO VECTOR (INPUT,OUTPUT) :STOPPING TEST. RELATIVE ERRORcZETA ( I N P U T ) :SPECTRAL RAOIUS OF I T E R A T I O N MATRIX :(MUST BE I N P U T I F IADAPT=O) :OVER-RELAXATION PARAMETER : (MUST BE INPUT I F IAOAPT=O) :=0 DOES NOT ADAPT TO ESTIMATE SR.OMEGA :=1 W I L L ADAPT ( I N P U T ) :MAX NO. OF I T E R A T I O N S ( I N P U T ) :NO. OF I T E R A T I O N S PERFORMED (OUTPUT) :DIAGONAL STORAGE WORK AREA O i S I Z E N :PREVIOUS PSEUDO-RESIDUAL NORM :CURRENT PSEUDO-RESIDUAL NORM :CURRENT SOLUTION NORM :STOPPING TEST :USED I N STOPPING TEST :TEMPORARY STORAGE OF NEW SOLUTION ELEMENT :SPECTRAL RADIUS OF JACOB1 MATRIX :UPPER BOUND FOR THE I ' T H ESTIMATE OF OMEGA :NUMBER OF DIFFERENT ESTIMATES OF OMEGA :NUMBER OF ITERATIONS USING THE NEW ESTIMATE :TEST I F SR IS REASONABLE ESTIMATE

REAL A ( I A D I M ) ,B(N) ,U(N) ,ZETA,SR,OMEGA INTEGER J A ( I A D I M ) , I S T A R T ( N P l ) ,N,NPl,IAOIM,IADAPT, +NUMITS.ITMAX REAL PSPREV ,PSNORM, O ( N) ,UNORM,SUM,TEST ,TEMP ,H ,SRJ, ~ ( 9 ,CHANGE ) INTEGER IEST,NEWITS DATA T/1.5,1.8,1.85,1.9,1.94,1.96,1.975,1.985,1.992/

****

SET UP DIAGONAL ELEMENTS F I N D O R I G I N A L PSEUDO-RESIDUAL NORM CALCULATE F I R S T ITERATES

****

I t e r a t i v e Methods

20

10

C C C C C 40

60

50 C C C C C C C C

I F ( IADAPT .EQ . 1)THEN IEST=l OMEGA=MIN ( OMEGA ,T ( 1) ) ENDIF PSPREV=O. DO 10 I = l , N SUM=B ( I ) DO 20 J=ISTART(I),ISTART(I+l)-1 I F ( J A ( J) .EQ. 1)THEN D( I ) = A ( J ) ELSE SUM=SUM-A( J)*U(JA( J ) ) ENDIF CONTINUE TEMP=OMEGA*SUM/D( I ) + ( 1.0-OMEGA)*U( PSPREV=PSPREV+(TEMP-U(1))**2 u( I )=TEMP CONTINUE PSPREV=SQRT( PSPREV ) NUMITS=l NEWITS-1

I)

****

PERFORM THE NEXT I T E R A T I O N CALCULATE THE NORMS OF THE ESTIMATE AND THE PSEUDO-RESIDUAL

****

NUMITS=NUMITStl IF(NUM1TS. GT. 1TMAX)RETURN NEWITS=NEWITS+l PSNORM=O. UNORM=D . DO 5 0 I = l , N SUM=B ( I ) DO 60 J = I S T A R T ( I ) ,I S T A R T ( I + l ) - l IF(JA( J ) .NE. I)SUM=SUM-A( J)*u(JA( J) ) CONTINUE TEMP=OMEGA*SUM/D( I ) + I 1.O-OMEGA)*U( 1) PSNORM=PSNORM+(TEMP-U ( I) ) **2 UNORM=UNORM+TEMP**Z U ( I )=TEMP CONTINUE PSNORM=SQRT( PSNORM) UNORM=SQRT( UNORM)

****

STOPPING TEST I F IADAPTZ1,USE R A T I O OF CONSECUTIVE PSEUDO-RESIDUALS AS AN ESTIMATE FOR THE SPECTRAL RADIUS OF THE I T E R A T I O N MATRIX MUST PERFORM A T LEAST TWO ITERATIONS WITH NEW VALUE OF OMEGA BEFORE STOPPING

****

I F ( IADAPT .EQ.O)THEN TEST=PSNORM/( 1.-SR)/UNORM I F ( A B S ( TEST) :LT .ZETAIRETURN GOT0 80 ENDIF I F ( NEW I T S . L T . 2)GOTO 80 IF(ABS(SR-PSNORM/PSPREV 1. GT. 0.01 ITHEN CHANGELO. ELSE

665

Len Colgan

666

C C C C

C C C 80

CHANGE.1. ENOIF SR=PSNORM/PSPREV I F ( SR .GE. l.O)GOT080 I F ( SR.LT.OMEGA-1.O)THEN H-OMEGA-1.0 ELSE H=SR ENDIF TEST=PSNORM/( 1.-H)/UNORM IF(ABS( TEST) LT .ZETA)RETURN

****

’

.

CHANGE OMEGA I F NECESSARY MUST BE AT LEAST FIVE ITERATIONS USING PREVIOUS OMEGA

****

I F ( NEW I T S .LT. 5.0R.SR.LT. (OMEGA-1. )**O. 75 .OR.CHANGE .EQ.O. )GOTO 80 IEST=I E S T + ~ NEWITS=O SRJ=(SR+OMEGA-l.O)/OMEGA/SQRT(SR) OMEGA=P.O/( l.O+SQRT( l.O-SRJ**Z)) I F ( IEST. LE .9 ITHEN OMEGA=MIN(OMEGA,T( IEST) ) ELSE OMEGA=MIN(OMEGA ,T( 9 ) ) ENDIF

****

PERFORM ANOTHER ITERATION

****

PSPREV=PSNORM GOT040 EN0

OUTPUT

The subroutine SOR y i e Is t h e fol.Jwing output f o r the f i r s t model problem. I n i t i a l l y t h e relaxation parameter w i s given t h e value 1, and the subroutine changes i t a t appropriate times only. T h i s occurs a f t e r t h e f i f t h i t e r a t i o n only. The value of H, used in t h e stopping t e s t , i s a l s o shown. NUM ITS=1 w=1.000000 .0100000 .0125000 .0131250 .0131250 .0173438 .0186719 H not y e t a p p l i c a b l e NUMITS=2 w=1.000000 .0162500 .0214063 .0230078 .0321777 H= A84742

.0132813 .0190820

.0125000 .0132813

.0162500 .0176563

.0173438 .0190820

.0176563 .0195410

.0230078 .0356299

.0201660 .0307568

.0214063 .0201660

.0293750 .0278564

.0321777 .0307568

.0278564 .0253784

,0301538 ,0471115

.0245026 ,0367928

.0282715 .0245026

.0402246 .0346814

.0434662 .0367928

.0346814 .0283964

NUMITS.3

w= 1.000000

.0207031 .0282715 .0301538 ,0434662 H=. 700422

I t e r a t i v e Methods

NUMITS=4 w=1.000000 .0241357 .0336285 .0353993 .0514349 H=.696357 NUMITS=5 w=1.000000 ,0268143 .0376902 .0391613 .0566937 H=. 686504

NUMITS=7 ~1.282128 .0312275 .0441691 .0449078 .0648243 H= .534703

NUMITS=9 u=1.282128 .0329194 .0462823 .0464766 .0665409 H=.434433

NUMITS=11 ~1.282128 .0333109 .0466583 .0466609 .0666648 H=. 282128

NUMITS=13 0=1.282128 0.333320 .0466667 .0466671 .0666673 H = . 282128

NUMITS= 15 u=1.282128 .0333335 .0466668 .0466668 .Of366667 H=. 367594

NUMI TS= 17 w=1.282128 .0333333 .0466667 ,0466667 .0666667 H=. 300329

667

.0353993 .0541139

.0275202 .04036 18

.03362a5 .0275202

.04a5474 .0389370

.0514349 .0403618

.0389370 ,0301809

.0391613 .0585277

.0295246 .0425884

.0376902 .0295246

.0545626 .0416450

.0566937 .0425884

.0416450 .0312942

.0449078 .0656558

.0326485 .0464129

.0441691 .0326485

.0638370 .0460275

.0648243 .0464129

.0460275 .0332224

,0464766 .0665959

,0332860 .0466478

.0462823 .0332860

.0663611 .0466257

.0665409 .0466478

.0466257 .0333289

.0466609 .0666674

.0333322 .0466675

.04665a3 .0333322

.0666587 .0466671

.0666648 .0466675

.0466671 .0333338

.0466671 .0666670

.0333335 .0466668

.0466667 .0333335

.0666673 .0466669

.0666673 .0466668

.0466669 .0333334

.0466668 .0666667

.0333334 .0466667

.0466668 .0333334

.0666668 .0466667

.0666667 .0466667

.0466667 .0333333

.0466667 .0666667

.0333333 .0466667

.0466667 .0333333

-0666667 .0466667

.0666667 .0466667

.0466667 .0333333

Len Colgan

668

APPENDIX E SUBROUTINE J S I C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C

C

C

SUBROUTINE TO SOLVE LARGE SPARSE LINEAR SYSTEMS USING JACOBI I T E R A T I O N WITH CHEBYSHEV ACCELERATION WRITTEN BY L.COLGAN 1983 SOLVES SYSTEMS OF N SPARSE L I N E A R EQUATIONS WHERE THE COEFFICIENT MATRIX I S STORED I N UNSYMMETRIC SPARSE FORM SUBROUTINE J S I(A, JA, ISTART ,N ,NP1, IADIM ,B ,U ,ZETA, +SRJ ,SRNJ, ICASE, IADAPT, ITMAX,NUMITS ,V,W ,0) GLOSSARY OF TERMS A( 1 JA( ) ISTART(I) N NP1 IADIM B( ) U( 1 ZETA SRJ SRNJ ICASE IADAPT ITMAX NUMITS V( ) W( 1 D( 1 PSORIG PSNORM UNORM TEST GAMMA SIGMA R OMEGA TEMP QA QT NEWITS

:IION-ZERO ELEMENTS OF SPARSE COEFF. MATRIX ( I N P U T ) :COLUMN POSITIONS OF ELEMENTS I N A ( I N P U T ) :POSITION IN A & JA OF THE FIRST ENTRY FROM :ROW I.LAST ENTRY I S I A D I M + l ( I N P U T ) :NO. OF EOUATIONS [ I N P U T ) :~ (INPUT) +i :DIMENSION OF A & J A ( I N P U T ) :RIGHT HAND S I D E OF EQUATIONS ( I N P U T ) :SOLUTION. ON INPUT CONTAINS I N I T I A L GUESS :WHICH MAY B E ZERO VECTOR (INPUT.OUTPUT1 :STOPPING TEST. RELATIVE ERRORCZETA (INPUT) :SPECTRAL RADIUS OF B A S I C JACOBI I T E R A T I O N MATRIX :(MUST BE INPUT) :LEAST EIGENVALUE ;MOST NEGATIVE ( INPUT) :=D DO NOT CHANGE SRNJ : = 1 L E T SRNJ=-SRJ ( I N P U T ) :=0 DOES NOT ADAPT TO ESTIMATE SRJ,SRNJ : = 1 WILL ADAPT ( I N P U T ) :MAX NO. OF ITERATIONS (INPUT) :NO. OF ITERATIONS PERFORMED (OUTPUT) :TEMPORARY STORAGE WORK AREA OF SIZE' N :STORAGE FOR PREVIOUS I T E R A T I O N OF S I Z E N :DIAGONAL STORAGE WORK AREA OF S I Z E N :ORIGINAL PSEUDO-RESIDUAL NORM WITH CURRENT SRJ :CURRENT PSEUDO-RESIDUAL NORM :CURRENT SOLUTION NORM :STOPPING TEST :PARAMETER :PARAMETER :PARAMETER :ACCELERATION PARAMETER :VARIOUS TEMPORARY STORAGES :ACTUAL RESIDUAL OUOTIENT :THEORETICAL RESIDUAL QUOTIENT :NUMBER OF ITERATIONS U S I N G THE NEW ESTIMATE

REAL GAMMA, SIGMA,R ,OMEGA,QA,QT INTEGER NEWITS

I t e r a t i v e Methods

C C C

20 10 C C C 25

**** SET UP DIAGONAL ELEMENTS

****

DO 10 I = I , N 00 20 J=ISTART(I),ISTART(I+l)-1

I F ( J A ( J) .EQ. 1)THEN O( I ) = A ( J ) GOT0 10 END1 F CONTINUE CONTINUE NUMITS=O

****

CALCULATE OTHER PARAMETERS

****

G A W A = 2 . 0 / ( 2 .O-SRJ-SRNJ) S I GMA=( S RJ-S RNJ ) / ( 2.0-SRJ-S RNJ ) R=SORT( 1.O-SIGMA**Z)

R=(i.O:R)/(l.O+R)

C

C 30

40

60

50 C C C

'

NEW ITS=O

****

CALCULATE THE ACCELERATION PARAMETER PERFORM THE NEXT I T E R A T I O N I F N E W I T S = l , F I N D O R I G I N A L PSEUDO-RESIDUAL NORM USING THE NEW VALUES OF SRJ,SRNJ

****

NUMITS=NUMITS+l I F ( NUMITS. GT. 1TMAX)RETURN NEWITS=NEWITS+l PSNORM=O. UNORM=O. I F ( NEWITS. EQ. 1)GOT040 I F ( NEWITS .EQ .PITHEN OMEGA=1.0/( 1.0-0.5*SIGMAX*2) ELSE OMEGAz1 . O / ( l.O-O.Z5*0MEGA*SIGMA**2) ENOIF 00 50 I = l . N SUM=B(I) ' DO 60 J=ISTART(I),ISTART(I+l)-l IF(JA(J) .NE. I)SUM=SUM-A(J)*U(JA( J)) CONTINUE TEMP=SUM/O( I) - u ( I I PSNORM=PSNORM+TEMP**Z IF( NEWITS. EQ. ITH HEN V( I)=GAMMA*TEMP+U( I ) ELSE V( I)=OMEGA*(GAMMA*TEMP+U(I))+(l.O-OMEGA)*W(I) E N 0 1F UNORM=UNORM+V( 1)**2 CONTINUE UNORM=SQRT( UNORM) PSNORM=SQRT( PSNORM)

****

STOPPING TEST

****

TEST=PSNORM/( 1.O-SRJ)/UNORM I F ( A B S ( TEST) . L T .ZETAIRETURN IF( NEWITS. EQ. 1 )PSORIG=PSNORM

669

Len Colgan

670

C C C

70 C

c

C

**** INTERCHANGE ARRAYS

****

00 70 I=l,N W( I)=U(I) U ( I )=Vi I ) CONTINUE

****

CHANGE PARAMETERS

****

IF NECESSARY

IF(IAOAPT.EQ.O.OR.NEWITS.EQ.1)GOT030

QA=PSNORM/PSORIG TEMP=R**( NEWITS-1) QT=2.0*SQRT(TEMP)/( l.O+TEMP) C C C C

IF(QA.LT.QT**0.75)GOTO30

****

SOLVE CHEBYSHEV EQUATION TO FINO A NEW ESTIMATE FOR SRJ

****

TEMP=0.5*( l.O+TEMP)*( QA+SQRT(QA*QA-QT*QT)) TEMPTEMP**( 1.O/(NEWITS-l .O). . TEMP=(TEMP+R/TEMP)/( 1.o+R)

SRJ=0.5*(SRJtSRNJtTEMP*(2.0-SRJ-SRNJ)) I F( ICASE .EQ .1) SRNJ=-SRJ GOT025 EN 0

OUTPUT The subroutine JSI y i e l d s t h e f o l l o w i n g o u t p u t f o r t h e f i r s t model problem. I t usesthe n o t a t i o n SRJ, SRNJ f o r M, m r e s p e c t i v e l y , s t a r t i n g w i t h SRJ=O, SRNJ=-1, adapting when convenient t o f i n d a b e t t e r SRJ, and t h e r e a f t e r l e t t i n g SRNJ=-SRJ. This occurred a f t e r the second and s i x t h i t e r a t i o n s o n l y . NUMITS=1 .0066667 .0066667 .0066667 .0066667 SRJ=0.000000

.0066667 .0066667

.0066667 .0066667

.0066667 .0066667

.0066667 .0066667

.0066667 .0066667

.0066667 .0066667

NUMITS=2 .0117467 .0129412 .0129412 .0141176 SRJ=O .OOOOOO

.0129412 .0141176

.0117647 .0129412

.0123412 .0117647

.0141176 .0123412

.0141176 .0129412

.0129412 .0117647

NUMITS=3 .0164706 ,0137059 .0197059 .0235294 SRJ=.762438

.0197053 .0235234

,0164706 .0187059

.0197059 .0164706

.0235294 .0197059

‘0235234 .0137059

.0197059 .0164706

NUMITS=4 .0231671 .0298375 .0298375 .0387883 SRJ=.762438

.0238375 .0387883

.0231671 .0298375

.0238375 .0231671

.0387883 .0238375

.0387883 .0298375

.0298375 .0231671

NUMITS.5 .0270955 .0363603 .0363603 .0496681 SRJ=.762438

.0363603 .0436681

.0270955 .0363603

.0363603 .0270355

.0436681 .0363603

.0436681 .0363603

.0363603 .0370955

Iberative Methods

671

NUMITS=6 .0293013 .0401694 .0401694 .0293013 .0401694 .0561959 .0561959 .0401694 .0401694 .0561959 .0561959 .0401694 .0293013 .0401694 .0401694 .0293013 SRJ=.762438

NUM ITS= 7

.0300847 .0414167 .0414167 .0300847 .0414167 .0581826 .0581826 .0414167 .0414167 .0581826 .0581826 .0414167 .0300847 .0414167 .0414167 .0300847 SRJ=.808666

NUMITS=8 .0313919 .0435149 .0435149 .0313919 .0435149 .0615504 .0615504 .0435149 .0435149 .0615504 .0615504 .0435149 .0313919 .0435149 .0435149 .0313919 SRJ=.808666 NUMITS=9 .0322941 .0449798 .0449798 .0322941 .0449798 .0639283 .0639283 .0449798 .0449798 .0639283 .0639283 .0449798 .0322941 .0449798 .0449798 .0322941 SRJ=.808666 NUMITS-10 .0327923 .0457924 .0457924 .0327923 .0457924 .0652538 .0652538 .0457924 .0457924 .0652538 .0652538 .0457924 .0327923 .0457924 .0457924 .0327923 SRJ=.808666 NUMITS=11 0.330548 .0462177 .0462177 .0330548 .0462177 .0659432 .0659432 .0462177 .0462177 .0659432 .0659432 .0462177 .0330548 .0462177 .0462177 .0330548 SRJ=.808666

N UMITS= 12

.0331913 .0464371 .0464371 .0331913 .0464371 .0662957 .0662957 .0464371 .0464371 .0662957 .0662957 .0464371 .0331913 .0464371 .0464371 .0331913 SRJ=.808666

NUMITS= 13 .0332610 .0465493 .0465493 .0332610 .0465493 .0664762 .0664762 .0465493 .0465493 .0664762 .0664762 .0465493 .0332610 .0465493 .0465493 .0332610 SRJ=.808666 NUMITS=14 .0332963 .0466065 .0466065 .0332963 .0466065 .0665691 .0665691 .0466065 .0466065 .0665691 .0665691 .0466065 .0332963 .0466065 .0466065 .0332963 SRJ=.808666 NUMITS=15 .0333142 -0466358 .0466358 .0333142 .0466358 .0666167 .0666167 .0466358 .0466358 -0666167 .0666167 .0466358 .0333142 .0466358 .0466358 .033142 SRJ=.808666 NUMITS=16 .0333235 .0466508 .0466508 .0333235 .0466508 .0666411 .0666411 .0466508 .0466508 .0666411 .0666411 .0466508 .0333235 .0466508 .0466508 .0333235 SRJ=.808666 NUMITS=17 .0333283 .0466585 .0466585 ,0333283 .0466585 .0666535 .0666535 .0466585 .0466585 .0666535 .0666535 ,0466585 .0333283 .0466585 .0466585 .0333283 SRJ=.808666

Len CoZgan

672

NUMITS=18 .0333308 .0466625 .0466625 .0666599 SRJ=.808666

.0466625 .0666599

.0333308 .0466625

.0466625 .0333308

.0666599 .0466625

.0666599 .0466255

,0466625 ,0333308

NUMITS=19 .0333320 ,0466645 .0466645 .Of566632 SRJ=. 808666

,0466645 .Of566632

.0333320 .0466645

.0466645 .0333320

.0666632 .0466645

.0666632 .0466645

.0466645 .03333?0

NUMITS=20 .0333326 .0466656 .0466656 .0666649 SRJ=.808666

.0466656 .0666649

.0333326 .0466656

.0466656 .0333326

,0666649 .0466656

.0666649 .0466656

.0466656 .0333326

NUM ITS=2 1 .0333330 ,0466661 .0466661 .0666657 SRJ=.808666

.0466661 .Of366657

.0333330 .0466661 .046666 1 .0333330

.0466661 .0466661 ,0333330

,0666657

.0666657

.0466661

NUMITS=22 .0333332 .0466664 .0466664 .Of366662 SRJ=.808666

.0466664 .Of366662

.0333332 .0466664

.0466664 .0333332

.0666662 .0466664

.0666662 .0466664

.0466664 .0333332

NUMITS=23 .0333332 .0466665 .0466665 .0666664 SRJ=.808666

.0466665 .0666664

.0333332 .0466665

.0466665 .0333332

.Of566664 .0466665

.Of366664 .0466665

.0466665 .0333332

NUM ITS=24 .0333333 .0466666 .0466666 .0666665 SRJ=.808666

.0466666 .Of566665

.0333333 .0466666

.0466666 .0333333

.0666665 .0466666

,0666665 .0466666

.0466666 .0333333

NUMITS=25 .0333333 .0466666 .0466666 .0666666 SRJ=.808666

.0466666 .0466666

.0333333 .0466666

.0466666 .0333333

.0666666 .0466666

.0666666 .0466666

.0466666 .0333333

I t e r a t i v e Methods

APPENDIX F SUBROUTINE JCG C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C

SUBROUTINE TO SOLVE LARGE SPARSE L I N E A R SYSTEMS U S I N G JACOBI I T E R A T I O N W I T H CONJUGATE GRADIENT ACCELERATION WRITTEN BY L.COLGAN 1983 SOLVES SYSTEMS OF N SPARSE L I N E A R EQUATIONS WHERE THE COEFFICIENT MATRIX I S STORED I N UNSYPMETRIC SPARSE FORM SUBROUTINE JCG(A,JA, I S T A R T ,N ,NP1 ,IADIM,B ,U,ZETA, +SRJ, IADAPT, ITMAX,NUMITS .W ,D ,PSU ,PSV ,PSW) GLOSSARY OF TERMS :NON-ZERO ELEMENTS OF SPARSE COEFF. MATRIX ( I N P U T ) A( ) :COLUMN POSITIONS O F ELEMENTS I N A ( I N P U T ) JA( ) I S T A R T ( I) : P O S I T I O N I N A & J A OF THE F I R S T ENTRY FROM :ROW I. LAST ENTRY I S I A D I M + l ( I N P U T ) N :NO. OF EQUATIONS ( I N P U T ) NP1 :N+1 ( I N P U T ) IADIM :DIMENSION OF A & JA (INPUT) :RIGHT HAND S I D E O F EQUATIONS ( I N P U T ) B( 1 :SOLUTION. ON INPUT CONTAINS I N I T I A L GUESS U( 1 :WHICH MAY BE ZERO VECTOR (INPUT,OUTPUT) ZETA :STOPPING TEST. R E L A T I V E ERRORcZETA ( I N P U T ) :SPECTRAL RADIUS O F B A S I C JACOBI I T E R A T I O N MATRIX SRJ :(MUST BE I N P U T ) IADAPT :=0 DOES NOT ADAPT TO ESTIMATE S R J :=1 W I L L ADAPT ( I N P U T ) ITMAX :MAX NO. OF ITERATIONS ( I N P U T ) NUMITS :NO. OF ITERATIONS PERFORMEO (OUTPUT) :STORAGE FOR PREVIOUS I T E R A T I O N OF S I Z E N W( 1 :DIAGONAL STORAGE WORK AREA OF S I Z E N D( PSNORM :CURRENT PSEUDO-RESIDUAL NORM UNORM :CURRENT SOLUTION NORM TEST :STOPPING TEST :STORAGE FOR CURRENT PSEUDO-RESIDUAL VECTOR PSU( ) :STORAGE FOR NEW PSEUDO-RESIDUAL VECTOR PSV( ) STORAGE FOR PREVIOUS PSEUDO-RESIDUAL VECTOR PSW( ) PSUPSV : I N N E R PRODUCT PSUPSU :INNER PRODUCT PSWPSW : I N N E R PRODUCT GAMMAU :PARAMETER GAMMAW :PARAMETER OMEGAU :CURRENT ACCELERATION PARAMETER OMEGAW :PREVIOUS ACCELERATION PARAMETER ALPHA( ) :DIAGONAL ELEMENTS OF TRIDIAGONAL MATRIX BETA( ) :SQUARE OF OFF-DIAGONAL ELEMENTS CHANGE :=0 I F SRJ IS ACCURATE ENOUGH TEMP :VARIOUS TEMPORARY STORAGES REAL A ( I A D I M ) , B ( N ) ,U(N) ,ZETA,SRJ INTEGER J A ( I A D I M ) , I S T A R T ( N P l ) ,N,NPl,IADIM,IADAPT, +NUMITS,ITMAX REAL PSNORM,W(N) ,D(N) ,PSU(N) ,PSV(N) ,PSW(N) ,UNORM,SUM, +TEST,TEMP ,PSUPSV ,PSUPSU ,PSWPSW ,GAMMAU ,GAMMAW, +OMEGAU ,OMEGAW,ALPHA( 200) ,BETA( 200)

67 3

Len CoZgan

674

INTEGER CHANGE

****

SET UP DIAGONAL ELEMENTS F I N O O R I G I N A L PSEUDO-RESIDUAL AND I T S INNER PRODUCT

****

20 10

C C C 30

PsuPsu=o.o DO 10 I = l , N SUM=B( I) DO 20 J = I S T A R T ( I ) , I S T A R T ( I + l ) - 1 I F ( J A ( J ) . E Q . 1)THEN D ( I ) = A ( J) ELSE SUM=SUM-A(J)*U(JA(J)) ENDIF CONTINUE PSU(I)=SUM/D(I)-U(I) PSUPSU=PSUPSU+D( I ) * P S U ( I)**2 CONTINUE CHANGE= 1 NUMITS=O GAMMAU=O. 0 GOT050

****

ESTIMATE SRJ

****

I F ( CHANGE .EQ .O.OR. TADAPT. EQ.O)GOT040

IF"JMITS.GT.2ITHEN CALL NEWTON( S R J ,CHANGE ,ZETA ,ALPHA,BETA,NUMITS) E L S E I F ( NUMITS .EQ. 1)THEN SRJ=ALPHA( 1) ELSE SRJ=ALPHA( l ) + A L P H A ( 2 ) TEMP=ALPHA( 1)-ALPHA( 2) SRJ=0.5*(SRJ+SQRT(TEMP**2+4.0*BETA( 2 ) ) ) ENDIF P R I N T 35,NUMITS,SRJ FORMAT( l X , ' N U M I T S = ' , I 4 / 1 X , , ' S R J =',E12.6/) .

35 C C C 40 C C C C 50

70 60

I

****

STOPPING TEST

****

TEST=PSNORM/( 1.O-SRJ)/UNORM I F ( ABS(TEST) . L T . ZETA) RETURN

****

CALCULATE NEW PARAMETERS NEEOS TO F I N D ( I T E R A T I O N M A T R I X ) * ( PSU)

****

NUMITS=NUMITS+l I F ( N U M I T S . G T . 1TMAX)RETURN PSUPSV=O .o DO 60 I = l , N SUM=O. 0 DO 70 J=ISTART( I 1, ISTART( I+ 1)- 1 IF( JA(J) .NE.I)SUM=SUM-A( J)*PSU( JA( J)) CONTINUE PSV( I ) = S U M / D ( I ) P s u P s v = P s u P s v + P s u ( I )*D( I )*PSV( I ) CONTINUE

I t e r a t i v e Methods

C C C

C C C C C

80

90

C C C L

GAMMAW=GAMMAU GAMMAUzl .O/( 1.0-PSUPSV/PSUPSU)

****

CALCULATE NEW ACCELERATION PARAMETERS

****

I F ( NUMITS. EQ. 1)THEN OMEGAU=1.0 ELSE TEMP=l .O/ ( 1.O-PSUPSU*GAMMAU/(OMEGAU*GAMMAW*PSWPSW)) OMEGAW=OMEGAU OMEGAU=TEMP END1 F

****

F I N D NEW SOLUTION ESTIMATE AND NORM F I N D NEW PSEUDO-RESIDUAL AND NORM INTERCHANGE ARRAYS FOR NEXT I T E R A T I O N

****

UNORM=O .O PSNORM=O.O PswPsw=PsuPsu PsuPsu=o.o I F ( N U M I T S . E Q . 1)THEN DO 80 I = l , N TEMP=GAMMAU*PSU( I )+u( I ) W(I)=U(I) U( I)=TEMP UNORM=UNORM+TEMP**Z TEMP=GAMMAU*PSV( I ) + ( l . O - G A M M A U ) * P S U ( I ) PSW( I ) = P S U ( I ) PSU( I )=TEMP TEMP=TEMP**Z PSNORM=PSNORM+TEMP PSUPSU=PSUPSU+D( I )*TEMP CONTINUE ALPHA( 1 ) z l . O - 1.O/GAMMAU ELSE DO 90 I = l , N TEMP=OMEGAU*( GAMMAU*PSU( I ) + U ( I ) ) + ( 1.O-OMEGAU)*W( I ) W( I ) = U ( I ) U( I ) = T E M P UNORM=UNORM+TEMP**Z TEMP=OMEGAU*(GAMMAU*PSV( I ) + (1.O-GAMMAU)*PSU( I ) ) t+( 1.0-OMEGAU)*PSW( I ) PSW(I)=PSU( I ) PSU( I )=TEMP . TEMP=TEMP**Z PSNORM=PSNORM+TEMP PSUPSU=PSUPSU+D( I )*TEMP CONTINUE ALPHA( N U M I T S ) = l . 0- 1.O/GAMMAU B E T A ( N U M I T S ) = ( OMEGAU-1 .O)/( GAMMAW*GAMMAU*OMEGAW*OMEGAU) ENDIF UNORM=SQRT( UNORM) PSNORM=SQRT( PSNORM) GOT030 END

****

SUBROUTINE TO F I N D NEW SRJ ESTIMATE LARGEST EIGENVALUE OF A TRIDIAGONAL MATRIX USES NEWTON'S METHOD

****

675

Len CoZgan

676

100

110

SUBROUTINE NEWTON( SRJ ,CHANGE ,ZETA,ALPHA,BETA,NUMITS) REAL SRJ ,ZETA,ALPHA( 200) ,BETA( 200) +,X,FX( 200) ,DFX(200) ,DELTAX INTEGER CHANGE,NUMITS X=SRJ FX( l)=ALPHA( 1)-X DFX( 1)=-1.0 FX(P)=FX( l)*(ALPHA(2)-X)-BETA(2) OFX(2)=DFX( l)*(ALPHA(2)-X)-FX( 1 ) 00 110 I=3,NUMITS FX( I ) = F X ( I-l)*(ALPHA( I)-X)-FX( I-2)*BETA( I ) DFX( I)=DFX( I-l)*(ALPHA( I ) - X ) - F X ( I-1)-DFX( 1-2)*BETA( I ) CONTINUE OELTAX=FX ( NUMITS) /DFX ( NUMITS) X=X-DELTAX IF( ABS(D E L T A X ~ G TZETA)GOTO~OO . IF(ABS (X-SRJ) .LT. ZETA) CHANGE=O SRJ=X RETURN END

OUTPUT The s u b r o u t i n e JCG y i e l d s t h e f o l l o w i n g o u t p u t f o r t h e f i r s t model problem. I t does n o t r e q u i r e any parameters d u r i n g t h e a c t u a l i t e r a t i o n s , b u t must e s t i m a t e M ( c a l l e d SRJ i n t h e s u b r o u t i n e ) f o r t h e s t o p p i n g t e s t . T h i s i s done i n s i d e t h e a d d i t i o n a l s u b r o u t i n e NEWTON. NUMITS=l .0400000 .0400000 .0400000 .0400000 SRJ=.750000

.0400000 .0400000

.0400000 .0400000

.0400000 .0400000

.0400000 .0400000

.0400000 .0400000

.0400000 .0400000

NUMITS=2 .0320000 .0480000 .0480000 .0640000 SRJ=.806186

.0480000 .0640000

.0320000 .0480000

.0480000 .0320000

.0640000 .0480000

.0640000 .0480000

.0480000 .0320000

NUMITS=3 .0333333 .0466667 .0466667 .0666667 SRJ=.809017

.0466667 .0666667

.0333333 .0466667

.0466667 .0333333

.0666667 .0466667

.0666667 .0466667

.0466667 .0333333

Iterative Methods

677

BIBLIOGRAPHY Axelsson, O . , (1977) " S o l u t i o n o f l i n e a r systems o f equations: i t e r a t i v e methods", i n L e c t u r e Notes i n Math.: Sparse m a t r i x techniques 572 (V.A. B a r k e r , ed.) Springer-Verlag, N.Y. Axelsson, 0. and Gustafsson, I . , (1977) "A m o d i f i e d upwind scheme f o r c o n v e c t i v e t r a n s p o r t e q u a t i o n and t h e use o f t h e c o n j u g a t e g r a d i e n t method f o r t h e s o l u t i o n o f nonsymmetric systems o f equations", Computer Sciences Dept. Report 77-12R, Chalmers Univ.Tech., Goteborg, Sweden. B a r t e l s , R . and D a n i e l , J.W., (1973) "A c o n j u g a t e g r a d i e n t approach t o n o n l i n e a r e l l i p t i c boundary v a l u e problems i n i r r e g u l a r r e g i o n s " , L e c t u r e Notes i n Math.363, S p r i n g e r - V e r l a g , N.Y. B i r k h o f f , G. and Varga, R.S., Trans.Amer.Math.Soc. 92.

(1959)

" I m p l i c i t a l t e r n a t i n g d i r e c t i o n methods",

B i r k h o f f , G., Varga, R.S. and Young, D., methods", Advances i n Computers 3 .

(1962)

"Alternating direction i m p l i c i t

Buzbee, B.L., Golub, G.H. and N i e l s o n , C.W., (1970) "The method o f odd/even r e d u c t i o n and f a c t o r i z a t i o n w i t h a p p l i c a t i o n t o P o i s s o n ' s e q u a t i o n " , Siam J.Numer. Anal. 7. Carre, B.A., (1961) "The d e t e r m i n a t i o n o f t h e o p t i m a l a c c e l e r a t i n g f a c t o r f o r successive o v e r r e l a x a t i o n " , t h e Computer J o u r n a l 4. Chandra, R., (1977) "Conjugate g r a d i e n t methods f o r p a r t i a l d i f f e r e n t i a l e q u a t i o n s " , d o c t o r a l t h e s i s , Yale U n i v e r s i t y , New Haven, C o n n e c t i c u t . Concus, P. and Golub, G.H., (1976) "A g e n e r a l i z e d c o n j u g a t e g r a d i e n t method f o r n o n s y m e t r i c systems o f l i n e a r equations", Stan-CS-76-536, Computer Science Dept., Stanford University, Palo Alto, C a l i f o r n i a . Concus, P . , Golub, G.H. and O'Leary, D.R., (1976) "A g e n e r a l i z e d c o n j u g a t e g r a d i e n t method f o r t h e numerical s o l u t i o n o f e l l i p t i c p a r t i a l d i f f e r e n t i a l e q u a t i o n s " , i n Sparse M a t r i x Computation (J.R. Bunch and D.J. Rose, eds.), Academic Press.

a2u a2u - au Douglas, J . Jr., (1955) "On t h e numerical i n t e g r a t i o n o f 7+ 7 - - by ay at a x i m p l i c i t methods , J.Soc.Ind.App1 .Math.3. Douglas, J . Jr., (1961) "A survey o f numerical methods f o r p a r a b o l i c d i f f e r e n t i a l e q u a t i o n s " Advances i n Computers, Vol.2, (F.L. A l t , ed.), Academic Press, N.Y. Douglas, J. J r . and Rachford, H.H. Jr., (1956) "On t h e numerical s o l u t i o n o f h e a t c o n d u c t i o n problems i n two o r t h r e e space v a r i a b l e s " , Trans.Amer.Math.Soc.82. E h r l i c h , L.W., (1963) "The B l o c k symmetric successive o v e r r e l a x a t i o n method", d o c t o r a l t h e s i s , Univ.Texas a t A u s t i n , A u s t i n , Texas. E n g e l i , M., Ginsburg, T., Rutishauser, H. and S t i e f e l , E., (1959) "Refined i t e r a t i v e methods f o r t h e computation o f t h e s o l u t i o n and t h e eigenvalues o f s e l f - a d j o i n t boundary v a l u e problems", M i t t . I n s t . f . a n g e w . Math.ETH Z u r i c h , No. 8. Evans, D.J., (1967) "The use o f p r e c o n d i t i o n i n g i n i t e r a t i v e methods f o r s o l v i n g l i n e a r equations w i t h symmetric p o s i t i v e d e f i n i t e m a t r i c e s " , J.Inst.Math.Appl.4. Forsythe, G.E. and Wasow, W.R., (1960) " F i n i t e D i f f e r e n c e Methods f o r P a r t i a l D i f f e r e n t i a l Equations", Wiley, N.Y.

Len CoZgan

678

Ginsburg, T.,

(1963)

"The c o n j u g a t e g r a d i e n t method", Numer.Math.5.

Golub, G.H., (1959) "The use o f Chebyshev m a t r i x p o l y n o m i a l s i n t h e i t e r a t i v e s o l u t i o n o f 1 i n e a r systems compared w i t h t h e methods o f successive o v e r r e l a x a t i o n " , d o c t o r a l t h e s i s , U n i v . I l l i n o i s , Urbana, I l l i n o i s . Golub, G.H. and Varga, R.S., (1961) "Chebyshev s e m i - i t e r a t i v e methods, successive o v e r r e l a x a t i o n methods, and second-order Richardson i t e r a t i v e methods, P a r t I and P a r t I I " , Numer.Math.3. Grimes, R.G., K i n c a i d , D.R., MacGregor, W . I . and'Young, D.M., (1978) "ITPACK Report: Adaptive I t e r a t i v e A l g o r i t h m s u s i n g Symmetric Sparse Storage", CNA-139, Center f o r Numerical A n a l y s i s , U n i v . o f Texas a t A u s t i n , A u s t i n , Texas. H a b e t l e r , G.J. and Wachspress, E.L., (1961) "Symmetric successive o v e r r e l a x a t i o n i n s o l v i n g d i f f u s i o n d i f f e r e n c e equations" , Math .Comp. 15. Hageman, L.A., (1972) "The e s t i m a t i o n o f a c c e l e r a t i o n parameters f o r t h e Chebyshev polynomials and t h e successive o v e r r e l a x a t i o n methods", WAPD-TM-1038, B e t t i s Atomic Power Laboratory, P i t t s b u r g h , Pennsylvania. Hageman, L.A. and Young, D.M., N.Y.

(1981)

" A p p l i e d i t e r a t i v e methods", Academic Press,

Hageman, L.A., Luk, F. and Young, D.M., (1977) "The a c c e l e r a t i o n o f i t e r a t i v e methods: p r e l i m i n a r y r e p o r t " , CNA-129, Center f o r Numerical A n a l y s i s , U n i v . o f Texas a t A u s t i n , A u s t i n , Texas. H e l l e r , J., (1960) "Simultaneous successive and a1 t e r n a t i n g d i r e c t i o n schemes", J. SOC. Indust.App1 .Math.8. H e n r i c i , P., (1962) " D i s c r e t e v a r i a b l e methods i n O r d i n a r y D i f f e r e n t i a l Equations", John W i l e y & Sons I n c . , N.Y. Hestenes, M.R. and S t i e f e l , E.L., (1952) s o l v i n g l i n e a r systems", NBS J.Res. 49.

"Methods o f c o n j u g a t e g r a d i e n t s f o r

Kahan, W., (1958) "Gauss-Seidel methods o f s o l v i n g l a r g e systems o f l i n e a r equations", d o c t o r a l t h e s i s , U n i v . o f Toronto, T o r o n t o . K i n c a i d , D.R., (1971) "An a n a l y s i s of a c l a s s o f norms o f i t e r a t i v e methods f o r systems of l i n e a r equations", d o c t o r a l t h e s i s , U n i v . o f Texas a t A u s t i n , A u s t i n , Texas. K i n c a i d , D.R. and Young, D.M. (1978) "Survey o f i t e r a t i v e methods", CNA-135, Center f o r Numerical A n a l y s i s , Univ. o f Texas a t A u s t i n , A u s t i n , Texas. M a n t e u f f e l , T.A., (1975) "An i t e r a t i v e method f o r s o l v i n g nonsymmetric l i n e a r systems w i t h dynamic e s t i m a t i o n o f parameters", d o c t o r a l t h e s i s , Univ. I l l i n o i s , Urbana , I1 1i n o i s . Peaceman, D.W. and Rachford, H.H. Jr., (1955) "The numerical s o l u t i o n of p a r a b o l i c and e l l i p t i c d i f f e r e n t i a l e q u a t i o n s " J.Soc.Indus.App1. Math. 3 . P r i c e , H. and Varga, R.S., (1962) "Recent numerical experiments comparing successive o v e r r e l a x a t i o n i t e r a t i v e methods w i t h i m p l i c i t a1 t e r n a t i n g d i r e c t i o n methods", Report No. 91, Gulf Research and Development Co., P i t t s b u r g h , Pennsylvania. Schwarz, H.R., Rutishauser, H. and S t i e f e l , E., Symmetric M a t r i c e s " , P r e n t i c e - H a l l Inc., N.J.

(1973) "Numerical a n a l y s i s o f

I t e r a t i v e Methods

679

Stone, H.L., (1968) " I t e r a t i v e S o l u t i o n s o f i m p l i c i t a p p r o x i m a t i o n s o f m u l t i dimensional p a r t i a l d i f f e r e n t i a l equations", Siam J.Num.Ana1. 5. Tee, G.J., (1963) " E i g e n v e c t o r s o f t h e successive o v e r r e l a x a t i o n process and i t s c o m b i n a t i o n w i t h Chebyshev s e m i - i t e r a t i o n " , C0mput.J. 6. Varga, R.S.,

(1962)

" M a t r i x I t e r a t i v e A n a l y s i s " , P r e n t i c e - H a l l , N.J.

Wachspress, E . L . , (1966) " I t e r a t i v e s o l u t i o n o f e l l i p t i c systems and a p p l i c a t i o n s t o t h e Neutron d i f f u s i o n equations o f r e a c t o r p h y s i c s " , P r e n t i c e - H a l l . N.J. Westlake, J.R., (1968) "A handbook o f numerical m a t r i x i n v e r s i o n s and s o l u t i o n o f l i n e a r e q u a t i o n s " , John W i l e y & Sons I n c . , N.Y. Widlund, O., (1978) "A Lanczos method f o r a c l a s s o f n o n s y m e t r i c systems o f l i n e a r equations", Siam J.Numer.Ana1. No. 4. Wrigley, H.E., (1963) "On a c c e l e r a t i n g t h e Jacobi method f o r s o l v i n g simultaneous equations by Chebyshev e x t r a p o l a t i o n when t h e eigenvalues o f t h e i t e r a t i o n m a t r i x a r e complex", Comput J. 6.

.

Young, D.M. (1950) " I t e r a t i v e methods f o r s o l v i n g p a r t i a l d i f f e r e n c e e q u a t i o n s o f e l l i p t i c type", d o c t o r a l t h e s i s , Harvard Univ., Cambridge, Massachusetts. Young, D.M. (1954) "On R i c h a r d s o n ' s method f o r s o l v i n g l i n e a r systems w i t h p o s i t i v e d e f i n i t e m a t r i c e s " , J .Ma t h .Phys X X X I I .

.

Young, D.M. (1970) "Convergence p r o p e r t i e s o f t h e s y m n e t r i c and unsymmetric successive o v e r r e l a x a t i o n methods and r e l a t e d methods", Math.Comp. 24. Young, D.M. N.Y.

(1971)

" I t e r a t i v e s o l u t i o n o f l a r g e l i n e a r systems" Academic Press,

Young, D.M. (1977) "On t h e a c c e l e r a t e d SSOR method f o r s o l v i n g l a r g e l i n e a r systems", i n Advances i n Math. 23.