- Email: [email protected]
Projector preconditioning and domain decomposition methods
Projector Preconditioning and Domain Decomposition Methods ZdenEk Dostal Hurnicktj Btav &AV A. l&ma 1768 708 00 Ostrava, Czechoslovakia
ABSTRACT
We
show
exploited
that the method
in the context
of preconditioning
of domain decomposition
by conjugate methods.
method is assessed by model analysis of finite difference equation on a rectangular
1.
projector
The performance
discretization
may be of the
of the Poisson
region.
INTRODUCTION
The emergence of parallel processors has stimulated the development of algorithms which can be decomposed into large subtasks with as little communication between them as possible. From this point of view, among the most attractive algorithms for numerical solution of systems of linear equations arising from discretization of elliptic partial differential equations are those based on domain decomposition. As the name of these algorithms suggests, the domain of the problem is decomposed into subdomains, which imposes a certain block structure on the matrix A of the discretized system. Then either this structure is used for an efficient preconditioning of A by its approximate factorization, or A is reduced to a small, well-conditioned system with unknowns corresponding to the meshpoints on boundaries between the subregions. For some recent algorithms and their analysis, we refer to [l, 21, and for a survey of domain decomposition papers, we refer to [3]. In this paper, we show that a variant of the domain decomposition technique may be obtained as a special case of a method which was introduced independently by B. A. Nicolaides 141 and the author [5], and which may be considered a particular implementation of the conjugate gradient method in a subspace [S]. Namely, we use the preconditioning by an
75
APPLZED MATHEMATICS AND COMPUTATZON 37:75-81(1990) 0 Elsevier Science Publishing Co., Inc., 1990 655 Avenue of the Americas, New York, NY 10010
0096-3003/90/$03.50
76
ZDENkK DOSTiL
A-conjugate projector on the subspace which is formed by the nodal variables corresponding to nodes inside the subregions. In this case, the matrix which must be inverted (or rather factorized) in the algorithm is block diagonal, so that we can use the parallel processor efficiently for this most laborious part of the solution process. We illustrate the performance of the method by analysis of a model problem. We believe that our approach provides a simple geometric insight into the domain decomposition methods. 2.
ALGORITHM We shall be concerned
with the solution of a system of linear equations Ax=b,
where
A is an n X n positive definite
(1)
matrix arising from the discretization
of a system of elliptic partial differential equations in the bounded region 1R in R” with the boundary JR. Suppose that CI is decomposed into the subregions R i, . . . , R, and the dividing boundary r, so that R is the union of the disjoint sets R,, . . , iIlk, I’. We shall also suppose that if the variables xi and xj correspond to the meshpoints of R, u ~~,\ r and CI, u dR,\ r, respectively, and if p < q, then i < j and the entry aij of A is equal to zero. For convenience, we suppose that the largest indices belong to m variables which correspond to nodes on r. All these requirements may be easily met if A is assembled by finite element discretization of Ri, i = 1,. . . , k, and then modified to satisfy the boundary conditions. The decomposition induces a block structure on A, A = (Aij). Notice that only the blocks Aii and Ai,k+ 1 are different from zero, Aii being positive definite for i = 1,. . . , k + 1. The algorithm for domain decomposition by conjugate projector exploits an auxiliary subspace which is generated by variables corresponding to nodes outside r. Though the specification of the subspace is the only difference from the general algorithm for preconditioning by conjugate projector, for convenience let us describe the whole algorithm. A.
Definition
of Decomposition
and Initial Co-rrection
(i) U = A , where 1 is the identity matrix of order n - m, and 0 is the ( 1 m X n zero matrix. (ii) Find the factorization
77
Domain Decomposition Methods
Notice that L = diag(L,,, . . . , L,,), where Lii are lower triangular matrices of the same order as Aii, i = 1,. . , k. Notice that P = ULwTL-‘UTA is a conjugate (iii) x0 = UL-TL-‘UTb. projector and that x0 = PA-lb. For convenience B.
put Q = I-
P and denote by W the range of Q.
Iterations (iv) i = 0, p,, = r0 = b - Ax,. (v) If ri = 0, then put r = xi and stop.
(vi) qi = Qpi. (vii) oi = rT9bi/9:A~q~. (viii) xi+r = xi + oi9i. (ix) ri+l=ri-~,A9i(=b-Axi+l). (x) Pi = r:+ rA9, /9:-9i. (xi) pi+i=ri+i-pipi’ (xii) Put i = i + 1 and return to (v). Notice
that the iterations
generate
the same residuals
as the conjugate
gradient method for solution of the system AQy = ra with y0 = 0, where only the positive definite restriction AQlAW of the symmetric matrix AQ (considered here as a linear mapping) to AW takes part in the process of solution. Hence the error of iterations may be estimated by the spectral condition number k(AQ]AW) of AQ]AW. We shall do this in detail for a model problem in the next section. Some other results concerning this algorithm may be found in [4, 51. We shall see that the matrix of AQ]AW in the auxiliary subset of the standard basis is the Schur complement corresponding to nodes on I. 3.
ANALYSIS
OF MODEL
PROBLEM
Consider
-Au=b u=o
in R, on alR,
where LR = (0,2) X (0,l) is an open rectangle. First we define a rectangular grid for discretization of (2). Put h a natural number, xi = ih, yj = jh, i = 0, 1,. . .,2p, j = 0, 1,. . , p. use a vertical meshline ordering (Figure 1) and take points in the order: first the points of Kl, = (0,l) X (0, l), then the points of 0, (0, l), and finally the points on F = (1) X (0,l).
(2)
= l/p, p We shall following = (1,2) X
78
ZDENkK
FIG. 1.
DOSTiiL
Model problem.
The central difference approximation with enhancing of the boundary conditions yields the system of linear equations Ax = b, where
C
D
DT=(O,...,
B = tridiag( - 1,4, -l),
C=tridiag(-l,B,-I), The matrices
0
1 and B are of order p - 1; the matrices
-I).
C and D have p -
1
block rows, each with p - 1 rows. Now put
where the order of I is that of diag(C, C) and 0 computations show that
has p - 1 rows. Simple
AQ=A[Z-U(UTAU)-‘UrA] =hP2diag(0,0,B-2DTC-‘D). Since AW is the orthogonal
complement
k(AQIAW)
of the range of U, we have
= k(B -2DTC-‘D).
Domuin Decomposition Methods
79
The matrix M = B - 2 DTC-‘D
is the Schur complement
of A corresponding
to the nodal variables of r, whose condition number k(M) is known to be O(h-‘), h + 0 [2]. Since k(A) is O(h-“), this gives us an idea of the preconditioning effect of the method. Moreover, we can use the approximate formula
k(M) = ;(l+ of [2] to assess more accurately
e-sr)(l-
the performance
e-*“)p
of the presented
method.
In
particular, for p = 100 we get k(M) = 90, which, substituted into the simplified Chebyshev error bound for the Euclidean norm of the residua
implies that for reduction of lra] by a factor of lop4 one needs about 47 steps of the above algorithm. In this case, the execution of the step A requires about i(2 X 100*)X lOO* = 10s flops for the factorization of the diagonal blocks and not more than 2 X lOO* X 120 = 2.4 X lo6 flops for the initial approximation, while 47 steps of iterations of the part B do not require more than 47X lOO* X 120 2: 6 X lo7 flops for iterations. If we realize that the factorization of A requires about 2 X 10s flops (apart from the diagonal blocks there is also the nonzero last block column), we see that our algorithm compares favorably with the standard Gauss elimination even without the parallel processor. As the time for both the factorization and the iterations may be halved on parallel processors, we can conclude that there are problems for which the method presented may be successful. We have not compared the performance of our algorithm with that of fast solvers [6], as our aim is not to solve the model problem but to get an idea of the performance of the method in a more general case when no fast solver can be applied. 4.
CONCLUSIONS
AND COMMENTS
In our presentation, the preconditioning effect of our algorithm may be explained also by the fact that the eigenvectors which correspond to small eigenvalues of A may be approximated by the vectors of the range of U (see [5] for such an estimate). This gives us some intuitive idea of the performance of the algorithm in the general case.
80
ZDENkK
FIG. 2.
DOSTAL
Numbering option.
Algorithms of this type may also be considered a tool for the solution of numbering problems, as their performance depends on numbering of nodes in the interior of the subregions only, which may lead reduction of the effective bandwidth (see Figure 2).
to considerable
Let us mention that the method may be described also as a special case of the well-known multigrid algorithms [9]. In this interpretation, the “coarse grid” is the subspace generated by the range of U, the operator of restriction is the conjugate projector, the operator of prolongation is the identity, and the smoothing is performed on the conjugate complement of the coarse grid. Our choice of smoothing enables us to enjoy the nice self-accelerating properties of the standard conjugate gradient method [lo]. Finally, we recall that the Schur complement may be again efficiently preconditioned [2]. REFERENCES 0. Axelsson and B. Polman, Methods
Block Preconditioning
I, Report 8735 (1987),
0. Axelsson and B. Polman, Method
Block Preconditioning
II, Report 8807 (1988),
P. E. Bjorstad tioned
and 0.
Problem
Solvers
Decomposition
Univ. of Nijmegen.
Solving elliptic problems
in Elliptic
Decomposition
Univ. of Nijmegen.
and Domain
Dept. of Mathematics,
B. Widlund,
into substructures,
and Domain
Dept. of Mathematics,
II (G.
on regions
parti-
Birkhoff and A.
Schoenstadt, Eds.), Academic, 1984, pp. 245-256. B. A. Nicolaides, value problems,
Deflation of conjugate SZAM J. Numer.
Z. Dost&l, Conjugate
gradient
gradients
with applications
Anal. 24(2):355-365
method with preconditioning
nut. J. Comput. Math. 23:315-324
(1988).
to boundary
(1987). by projector,
Inter-
Domain Decomposition Methods Methods of Numerical Mathematics, Nauka, Moscow, 1980. and C. Reinsch, Handbook for Automatic Computation-Linear
6
G. Marchuk,
7
J. H. Wilkinson Algebra,
8
Springer,
G. I. Marchuk a subspace
9 10
New York, 1970.
and Yu. A. Kuznetsov,
and its applications,
Mathematical Moscow,
81
Mod&Zing
(G.
A generalized
conjugate
gradient
in Problems of Computational I.
Marchuk
and
V.
P.
Dymnikov,
Eds.),
Mir,
1985, pp. 10-32.
W. Hackbusch,
Multi-grid Methods and Applications,
Springer,
Berlin,
H. A. van der Vorst and A. van der Sluis, The rate of convergence gradients,
method in
Mathematics and
Preprint
354, Dept.
of Mathematics,
Univ. of Utrecht,
1985.
of conjugate 1984.
ABSTRACT
We
show
exploited
that the method
in the context
of preconditioning
of domain decomposition
by conjugate methods.
method is assessed by model analysis of finite difference equation on a rectangular
1.
projector
The performance
discretization
may be of the
of the Poisson
region.
INTRODUCTION
The emergence of parallel processors has stimulated the development of algorithms which can be decomposed into large subtasks with as little communication between them as possible. From this point of view, among the most attractive algorithms for numerical solution of systems of linear equations arising from discretization of elliptic partial differential equations are those based on domain decomposition. As the name of these algorithms suggests, the domain of the problem is decomposed into subdomains, which imposes a certain block structure on the matrix A of the discretized system. Then either this structure is used for an efficient preconditioning of A by its approximate factorization, or A is reduced to a small, well-conditioned system with unknowns corresponding to the meshpoints on boundaries between the subregions. For some recent algorithms and their analysis, we refer to [l, 21, and for a survey of domain decomposition papers, we refer to [3]. In this paper, we show that a variant of the domain decomposition technique may be obtained as a special case of a method which was introduced independently by B. A. Nicolaides 141 and the author [5], and which may be considered a particular implementation of the conjugate gradient method in a subspace [S]. Namely, we use the preconditioning by an
75
APPLZED MATHEMATICS AND COMPUTATZON 37:75-81(1990) 0 Elsevier Science Publishing Co., Inc., 1990 655 Avenue of the Americas, New York, NY 10010
0096-3003/90/$03.50
76
ZDENkK DOSTiL
A-conjugate projector on the subspace which is formed by the nodal variables corresponding to nodes inside the subregions. In this case, the matrix which must be inverted (or rather factorized) in the algorithm is block diagonal, so that we can use the parallel processor efficiently for this most laborious part of the solution process. We illustrate the performance of the method by analysis of a model problem. We believe that our approach provides a simple geometric insight into the domain decomposition methods. 2.
ALGORITHM We shall be concerned
with the solution of a system of linear equations Ax=b,
where
A is an n X n positive definite
(1)
matrix arising from the discretization
of a system of elliptic partial differential equations in the bounded region 1R in R” with the boundary JR. Suppose that CI is decomposed into the subregions R i, . . . , R, and the dividing boundary r, so that R is the union of the disjoint sets R,, . . , iIlk, I’. We shall also suppose that if the variables xi and xj correspond to the meshpoints of R, u ~~,\ r and CI, u dR,\ r, respectively, and if p < q, then i < j and the entry aij of A is equal to zero. For convenience, we suppose that the largest indices belong to m variables which correspond to nodes on r. All these requirements may be easily met if A is assembled by finite element discretization of Ri, i = 1,. . . , k, and then modified to satisfy the boundary conditions. The decomposition induces a block structure on A, A = (Aij). Notice that only the blocks Aii and Ai,k+ 1 are different from zero, Aii being positive definite for i = 1,. . . , k + 1. The algorithm for domain decomposition by conjugate projector exploits an auxiliary subspace which is generated by variables corresponding to nodes outside r. Though the specification of the subspace is the only difference from the general algorithm for preconditioning by conjugate projector, for convenience let us describe the whole algorithm. A.
Definition
of Decomposition
and Initial Co-rrection
(i) U = A , where 1 is the identity matrix of order n - m, and 0 is the ( 1 m X n zero matrix. (ii) Find the factorization
77
Domain Decomposition Methods
Notice that L = diag(L,,, . . . , L,,), where Lii are lower triangular matrices of the same order as Aii, i = 1,. . , k. Notice that P = ULwTL-‘UTA is a conjugate (iii) x0 = UL-TL-‘UTb. projector and that x0 = PA-lb. For convenience B.
put Q = I-
P and denote by W the range of Q.
Iterations (iv) i = 0, p,, = r0 = b - Ax,. (v) If ri = 0, then put r = xi and stop.
(vi) qi = Qpi. (vii) oi = rT9bi/9:A~q~. (viii) xi+r = xi + oi9i. (ix) ri+l=ri-~,A9i(=b-Axi+l). (x) Pi = r:+ rA9, /9:-9i. (xi) pi+i=ri+i-pipi’ (xii) Put i = i + 1 and return to (v). Notice
that the iterations
generate
the same residuals
as the conjugate
gradient method for solution of the system AQy = ra with y0 = 0, where only the positive definite restriction AQlAW of the symmetric matrix AQ (considered here as a linear mapping) to AW takes part in the process of solution. Hence the error of iterations may be estimated by the spectral condition number k(AQ]AW) of AQ]AW. We shall do this in detail for a model problem in the next section. Some other results concerning this algorithm may be found in [4, 51. We shall see that the matrix of AQ]AW in the auxiliary subset of the standard basis is the Schur complement corresponding to nodes on I. 3.
ANALYSIS
OF MODEL
PROBLEM
Consider
-Au=b u=o
in R, on alR,
where LR = (0,2) X (0,l) is an open rectangle. First we define a rectangular grid for discretization of (2). Put h a natural number, xi = ih, yj = jh, i = 0, 1,. . .,2p, j = 0, 1,. . , p. use a vertical meshline ordering (Figure 1) and take points in the order: first the points of Kl, = (0,l) X (0, l), then the points of 0, (0, l), and finally the points on F = (1) X (0,l).
(2)
= l/p, p We shall following = (1,2) X
78
ZDENkK
FIG. 1.
DOSTiiL
Model problem.
The central difference approximation with enhancing of the boundary conditions yields the system of linear equations Ax = b, where
C
D
DT=(O,...,
B = tridiag( - 1,4, -l),
C=tridiag(-l,B,-I), The matrices
0
1 and B are of order p - 1; the matrices
-I).
C and D have p -
1
block rows, each with p - 1 rows. Now put
where the order of I is that of diag(C, C) and 0 computations show that
has p - 1 rows. Simple
AQ=A[Z-U(UTAU)-‘UrA] =hP2diag(0,0,B-2DTC-‘D). Since AW is the orthogonal
complement
k(AQIAW)
of the range of U, we have
= k(B -2DTC-‘D).
Domuin Decomposition Methods
79
The matrix M = B - 2 DTC-‘D
is the Schur complement
of A corresponding
to the nodal variables of r, whose condition number k(M) is known to be O(h-‘), h + 0 [2]. Since k(A) is O(h-“), this gives us an idea of the preconditioning effect of the method. Moreover, we can use the approximate formula
k(M) = ;(l+ of [2] to assess more accurately
e-sr)(l-
the performance
e-*“)p
of the presented
method.
In
particular, for p = 100 we get k(M) = 90, which, substituted into the simplified Chebyshev error bound for the Euclidean norm of the residua
implies that for reduction of lra] by a factor of lop4 one needs about 47 steps of the above algorithm. In this case, the execution of the step A requires about i(2 X 100*)X lOO* = 10s flops for the factorization of the diagonal blocks and not more than 2 X lOO* X 120 = 2.4 X lo6 flops for the initial approximation, while 47 steps of iterations of the part B do not require more than 47X lOO* X 120 2: 6 X lo7 flops for iterations. If we realize that the factorization of A requires about 2 X 10s flops (apart from the diagonal blocks there is also the nonzero last block column), we see that our algorithm compares favorably with the standard Gauss elimination even without the parallel processor. As the time for both the factorization and the iterations may be halved on parallel processors, we can conclude that there are problems for which the method presented may be successful. We have not compared the performance of our algorithm with that of fast solvers [6], as our aim is not to solve the model problem but to get an idea of the performance of the method in a more general case when no fast solver can be applied. 4.
CONCLUSIONS
AND COMMENTS
In our presentation, the preconditioning effect of our algorithm may be explained also by the fact that the eigenvectors which correspond to small eigenvalues of A may be approximated by the vectors of the range of U (see [5] for such an estimate). This gives us some intuitive idea of the performance of the algorithm in the general case.
80
ZDENkK
FIG. 2.
DOSTAL
Numbering option.
Algorithms of this type may also be considered a tool for the solution of numbering problems, as their performance depends on numbering of nodes in the interior of the subregions only, which may lead reduction of the effective bandwidth (see Figure 2).
to considerable
Let us mention that the method may be described also as a special case of the well-known multigrid algorithms [9]. In this interpretation, the “coarse grid” is the subspace generated by the range of U, the operator of restriction is the conjugate projector, the operator of prolongation is the identity, and the smoothing is performed on the conjugate complement of the coarse grid. Our choice of smoothing enables us to enjoy the nice self-accelerating properties of the standard conjugate gradient method [lo]. Finally, we recall that the Schur complement may be again efficiently preconditioned [2]. REFERENCES 0. Axelsson and B. Polman, Methods
Block Preconditioning
I, Report 8735 (1987),
0. Axelsson and B. Polman, Method
Block Preconditioning
II, Report 8807 (1988),
P. E. Bjorstad tioned
and 0.
Problem
Solvers
Decomposition
Univ. of Nijmegen.
Solving elliptic problems
in Elliptic
Decomposition
Univ. of Nijmegen.
and Domain
Dept. of Mathematics,
B. Widlund,
into substructures,
and Domain
Dept. of Mathematics,
II (G.
on regions
parti-
Birkhoff and A.
Schoenstadt, Eds.), Academic, 1984, pp. 245-256. B. A. Nicolaides, value problems,
Deflation of conjugate SZAM J. Numer.
Z. Dost&l, Conjugate
gradient
gradients
with applications
Anal. 24(2):355-365
method with preconditioning
nut. J. Comput. Math. 23:315-324
(1988).
to boundary
(1987). by projector,
Inter-
Domain Decomposition Methods Methods of Numerical Mathematics, Nauka, Moscow, 1980. and C. Reinsch, Handbook for Automatic Computation-Linear
6
G. Marchuk,
7
J. H. Wilkinson Algebra,
8
Springer,
G. I. Marchuk a subspace
9 10
New York, 1970.
and Yu. A. Kuznetsov,
and its applications,
Mathematical Moscow,
81
Mod&Zing
(G.
A generalized
conjugate
gradient
in Problems of Computational I.
Marchuk
and
V.
P.
Dymnikov,
Eds.),
Mir,
1985, pp. 10-32.
W. Hackbusch,
Multi-grid Methods and Applications,
Springer,
Berlin,
H. A. van der Vorst and A. van der Sluis, The rate of convergence gradients,
method in
Mathematics and
Preprint
354, Dept.
of Mathematics,
Univ. of Utrecht,
1985.
of conjugate 1984.