- Home
- Documents
*An iterative solution method for solving sparse ... An iterative solution method for solving sparse*

prev

next

out of 14

View

0Download

0

Embed Size (px)

Journal of Computational and Applied Mathematics 15 (1986) 339-352

North-Holland 339

An iterative solution method for solving sparse nonsymmetric linear systems

G.J. MAKINSON and A.A. SHAH Mathematical Institute, University of Kent, Canterbu?. Kent CT2 71VF, United Kingdom

Received 20 January 1985 Revised 7 Augustus 1985

Abstract: An algorithm is presented for the general solution of a set of linear equations Ax = b. The method Lvorks exceptionally well for the solution of large sparse systems of linear equations, the co-efficient matrix A of which need

not be symmetric but should have workable splits. The method can be applied to problems which arise in convection-diffusion, flow of fluids and oil reservoir modeling. The difference of the upper secondary diagonals (super

diagonals) and the lower secondary diagonals (sub diagonals) of the matrix A leads to a decomposition of A into a difference of a symmetric matrix, having the same lower structure as that of A. and a strictly upper triangular matrix. The symmetric matrix is made diagonally dominant and the system is solved iteratively.

1 Introduction

In order to solve the following sparse system

Ax=b (1.1)

where A is n x n matrix and b is a n X 1 column matrix, five well known basic iterative methods can be considered (i) The RF Method (ii) The Jacobi Method (iii) The Gauss-Seidel Method (iv) The Successive Over-Relaxation (SOR) method (v) The Symmetric Successive Over-Relaxation

(SSOR) Method. The developments presented by Lanczos [12], Golub and Varga [7], Hestenes and Stiefel [lo] are also of particular interest.

Though in some cases the Jacobi, Gauss-Seidel, SOR, SSOR methods are defined relative to the fixed partitioning imposed on the co-efficient matrix A, methods are known to converge if the following conditions are satisfied

(1) The co-efficient matrix A is symmetric and positive definite.

(2) The diagonal elements of A are non-zero. An iterative solution method for the linear system of which the co-efficient matrix is a

symmetric M-matrix is proposed by Meijerink and Van der Vorst [14]. As pointed out by Kershaw [ll] on a typical hard problem Meijerink and Van der Vorst method [14] is about 8000 times faster than the point Gauss-Seidel method, 200 times faster than the alternating direction implicit method and 30 times faster than the block successive overrelaxation method with optimum relaxation factor. However, the algorithm is not so effective if the co-efficient matrix is not symmetric and positive definite. Some authors tackle the unsymmetric case by forming AATx = ATb but this approach leads to ill conditioning and increases the condition number.

0377-0427/86/$3.50 Q 1986, Elsevier Science Publishers B.V. (North-Holland)

We encounter problems in which A is not necessarily positive definite nor even symmetric. For the solution of this class of problems we have Fletcher’s [5] ‘ bi-CG algorithm’. Fletcher’s [S] method is based on the Lanczos method which resembles the two-term form of the conjugate gradient method. Wong [17] has also done some work, which is very effective for the problems he considers, using his ‘row-sums agreement factorization’. In fact very few methods for tackling

such class of problems are presently available in the literature. The Manteuffel [13] adaptive procedure is also worth mentioning. Duff [4] has pointed out that for the more general case of unsymmetric linear systems there are no highly refined and efficient techniques. The method which is presented in this paper is very similar to Widlund’s [16] method which has been given recently. Widlund’s [16] method is closely related to the work by Concus and Golub [3]. They derive the algorithm using the Krylov sequence

Y(O) KV’O’, . . . ( Ck- ‘+P, . . . . It is a Lanczos, conjugate gradient-like method for the solution of the equation Ax = b. The co-efficient matrix is broken in to a difference

A=P-Q

where P is the symmetric part of A and Q is the skew-symmetric part of A. Thus

P=&4 +/IT) and Q = $(A -AT)

where A* is the transpose of A. The Cholesky factors of the symmetric part are computed to save storage and arithmetic operations. Young and Jea [18,19] have studied the acceleration of Widlund’s [16] method. At present this method is the best available method for the nonsym- metrizable cases. The method works only if the coefficient matrix A is positive real.

Our algorithm, the generalized Widlund’s method or the GW method, splits the unsymmetric matrix A in such a way that the symmetric part becomes diagonally dominant which can be decomposed into the Incomplete Cholesky factors [ll] or the Cholesky factors. The other matrix in the splitting becomes a triangular matrix. The method converges in the case when the coefficient matrix is nearly symmetric, for simplicity we call it e-semisymmetric (see Definition 3.3). The method is guaranteed to converge if the coefficient matrix contains workable splits (Definition 3.2). The method presented in this paper has a universal scope of application. The only restriction is that A should have workable splits.

In Section 2 we present notation and terminology. In Section 3 we describe our algorithm. The Section 4 contains the proof of convergence of the algorithm and in Section 5 the performance of the algorithm is explained with the help of examples and comparisons with standard methods.

2. Notation and terminology

Let (n) denote the set of positive integers. Let R”,” denote the set of all n x n real matrices A = [u~,~], for all i, j E (n). Let R”,*” c R”*” denote the set of real symmetric matrices. Let RI*” c R”*” denote the set of nonsymmetric matrices. Let R, n*n c R”-” denote the set of real lower

trtangular matrices. Let R”,*” c R”.” denote the set of real upper triangular matrices, Let R” denote the real n-dimensional vector space of all real vectors r = [ri, r2, r,, . . . , r,]‘, where ri E R for all i E (~2). Similarly C”*“, C$n, Cl*‘, CrVn, C,“.“, C”, c are defined for complex matrices,

vector spaces and vectors.

G. Makinson, A.A. Shah / Linear systems 341

3. Description of algorithm

Let the coefficient matrix of system (1.1) be such that A = [a,.,] E R:*", for all i. j E (n) is a nonsymmetric sparse matrix. Let A = L, + DA + 0" be the sphtting of A into strictly lower, diagonal, and strictly upper triangular matrices. Define S = L, + DA + Ll and H = Ll - U,. For certain diagonal matrices A E R",." consider the splitting

A=(S+A)-(H+A) (3.1)

where for all i, j E (n), the matrices

S = [sij] E R",*", H= [hi,,] E R",," are such that

s,,i = ai,i 9 s;,~ = ai,j if i > j, si,j = sj,i if i ,j, hi,j = ai,j - ajVi if i

342 G. Makmson, A.A. Shah / Linear systems

Theorem 3.3 (Young [20, p. 801). S(G) < 1 if and onfy if there exists a positire definite matrix P such that the matrix M gioen by

M=P-GPGH

is positive definite. If G is real, then P can be taken to be real.

Definition 3.2. A nonsymmetric matrix A E R:.” of system (l.l), having the GW splitting A = (S + A) - (H + A), is said to contain a pair of workable splits ((S + A). (H + A)) if, for an arbitrary nonnegative diagonal matrix A, A( H -t A)T is positive real.

Such matrices arise in the solution of the linear systems resulting from the discretization of elliptic boundary value problems where the first order partial derivative terms have different signs. In general it is sufficient for matrices to have a pair of workable splits if for an arbitrary nonnegative diagonal matrix A the modulus of the inner product of (H + A) is less than the inner product of (S + A). In some cases the splitting

A = (DA + u, + UAr) - (c/AT - LA), (s+A)=(D,+A+u,+u~)

and

(H+A)=(A+U,T-LJ

makes A( H + A)T positive real.

Definition 3.3. A matrix A = [a,,j] E R”*” is called an e-semisymmetric if for a small positive number E, 1 aj,j -a,,,1

G. Makinson. A.A. Shah / Linear sperm 343

Proof. A direct result of Theorem 3.1. 0

Theorem 4.1. If the coefficient matrix A = [a,.j] E R:.“, of the system (1 .l) is r-semisymmetric and can be represented by A = S - CP where c -C A,/( n - 1) and A, is the smallest eigenvalue of S. a matrix defined in Section 3, for P = [pi.,] E R”,.“, 1 p,., 1 < 1, the method (4.1) with the iteration matrix G given by (4.2) converges.

Proof. A = [a,.j] E Rz”’ is

344 G. Makinson, A.A. Shah / Linear sysrems

Table 4.1 Approximate number of multiplication operations needed for one iteration in the solution of Ax = 15. (A is an N x N matrix)

Method No. of Operations

SOR Gauss-Seidel Jacobi OAHM

6N 5N

SN IN

But (S+A)=A +(H+A),

K= (S+A)-‘[(A + (H+A))(A~+ (H+A)~) - (H+A)(H+A)‘]((S+A)-l)T

=(s+~)-'[AA~+A(H+~)~+(H+~)A~+(H+A)(H+A)~

= (S-A)-‘[AAT+A(H+A)~+ (H+A)A~]((S+A)-~)~.

As AA is SPD, (AAT)l/* exists. Therefore,

K=(S+A)-1(AAT)*‘2[I+(AAT)-1’2(/d(H+A)T

+ (HfA)AT)((AAT)-1’2)T](~~T)*‘2((S+A)-1)T

=[(S+A)-1(AAT)1’2][I+(AAT)-1’2(A(H+A)T

+ (H+A)AT)((hIT)-“2)T][(S+A)-1(AAT)1’2]T.

Using Lemma 4.1 and Theorem 3.2 it can be proved that K is positive definite and Theorem 3.3 proves the statement. q

At the end of this section (see Table 4.1) we give the approximate number of multiplication operations needed for the solution of (1.1) by the d