Algebraic aspects of finite-element solutions

385

Computer Physics Reports 6 (1987) 385-414 North-Holland, Amsterdam

ALGEBRAIC ASPECTS OF FINITE-ELEMENT SOLUTIONS J.K. REID HanveIl

Laboratory,

Oxfordshire,

UK

In these six lectures, we consider the solution of sets of linear equations Ax = b and eigenvalue problems Ax = hMx that have arisen from finite-element models. The emphasis is on direct methods for symmetric and positive-definite problems, but iterative methods and more general problems are also discussed.

J. K. Reid / Algebraic aspects of finIte-element

386

solutrons

Contents

1. Introduction.

2.

3.

4.

5.

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

1.1. Scope ............................................ 1.2. Gaussian elimination for full matrices. ..................... 1.3. Numerical stability ................................... ................... 1.4. Equivalence with triangular factorization ................................ 1.5. The effects of roundoff 1.6. The symmetric case. .................................. 1.7. Interchanges for unsymmetric cases ....................... 1.8. The symmetric nondefinite case .......................... 1.9. Monitoring stability .................................. ..................................... 1 .lO. Ill-conditioning .................................. 1.11. Iterative refinement ................. 1.12. Forward substitution and back-substitution. Band, profile and frontal methods for symmetric positive-definite systems ....................................... 2.1. Introduction 2.2. Band matrices ...................................... 2.3. Variable-band matrices ................................ 2.4. Variable-band matrices in auxiliary storage .................. 2.5. Assembling variable-band matrices ........................ 2.6. Frontal method ..................................... .................................. 2.7. Static condensation. 2.8. Ordering .......................................... 2.9. Comparison between band and frontal solutions .............. ...................... Substructuring or the multifrontal method. 3.1. The multifrontal method ............................... ...................................... 3.2. Substructuring ..................... 3.3. Choosing assembly trees by dissection ...... 3.4. Choosing assembly trees by the minimum-degree algorithm ........ 3.5. Ordering the operations associated with an assembly tree 3.6. Substructuring large problems ........................... Iterative methods ........................................ ....................................... 4.1. Introduction 4.2. Classical stationary iterative methods ...................... .................................. 4.3. Conjugate gradients ....................... 4.4. Preconditioned conjugate gradients ................................... 4.5. Multigrid iteration Eigenvalue problems ...................................... ....................................... 5.1. Introduction 5.2. Sylvester’s inertia theorem .............................. 5.3. Inverse iteration ..................................... ................................ 5.4. Simultaneous iteration 5.5. The algorithm of Lanczos ..............................

. . . . . . . .

. .

. . .

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

. . . . .

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

. . . . . .

388 388 388 389 389 390 391 391 392 392 392 393 393 393 393 394 394 394 394 395 396 396 397 398 398 399 399 399 400 400 401 401 401 402 403 404 405 405 405 405 406 407

J.K. Reid / Algebraic aspects of finite-element

387

solutions

5.6. The block Lanczos algorithm ............................... 5.7. Application of Lanczos’ method ............................. 6. Nondefinite and unsymmetric matrices ............................ 6.1. Introduction .......................................... 6.2. Band methods ......................................... 6.3. The frontal method. ..................................... 6.4. The multifrontal method .................................. 6.5. Classical stationary iterative methods ......................... 6.6. Conjugate gradients ..................................... 6.7. Multigrid iteration ...................................... References ..................................................

4%

409 409 409 410 410 411 411 411 412 412

J. K. Reid / Algehrurc aspects of finite-element

.solutions

1. Introduction I. I. Scope

In these six lectures, equations

we will be concerned

principally

with the solution

of sets of n linear

Ax=h

(1.1)

that have arisen from finite-element models. The rz X n matrix A is usually symmetric (AT = A) and positive definite (zTAz > 0 unless z = 0) because an energy minimizing principle underlies the mathematical model. These properties are very helpful for both direct and iterative methods and for most of these lectures we assume that they hold, but do not neglect the symmetric nondefinite and the unsymmetric cases. The matrix A and the vector b always originate as sums A

=

CK”‘, e

b =

&&)

(1.2)

t’

of contributions from the finite elements. Each matrix K (e) has nonzeros only in a few rows and columns and may therefore be stored in the computer in packed form in a small dense array whose rows and columns are labelled by index lists that indicate their positions in A. If the matrix pattern is symmetric, the row list and the column list are identical, so only one need be stored. Similarly cCe) may be stored as a packed vector labelled by a row-index list. The matrix A is sparse and may be treated by any sparse-matrix technique, but often there are advantages to be gained from the particular form (1.2) and we will stress these. The eq. (1.1) is usually solved by a variant of Gaussian elimination. We spend most of this lecture introducing the technique and then in lectures 2 and 3 explain various ways of exploiting the sparsity. Lecture 4 is devoted to a discussion of iterative methods for solving eq. (1.1) and lecture 5 to a discussion of the eigenvalue problem Ax=hMx. Finally in lecture 6 we discuss the solution unsymmetric sets of linear equations.

(1.3) of nondefinite

symmetric

sets of linear equations

and

1.2. Gaussian elimination for full matrices

All direct methods for solving eq. (1.1) are based on variants of Gaussian therefore worthwhile to review the properties of the basic method as applied before considering how to accommodate sparsity. In particular, we will explain of Wilkinson [l], which tells us when we can rely on the computed results; it is though not widely understood by nonspecialists.

elimination. It is to full matrices, the error analysis essentially simple

389

J.K. Reid / Algebraic aspects of finite-element solutions

If no interchanges elimination are a(k+l)

are in use, the operations

= a(k) _ aj;’ ‘I

‘J

/,,(k+l)

b(k)

=

I

(ai;))-lajj;)

al’k”‘(

_

I

performed

at the

kth

major

step

of the

(1.4)

(i, j> k),

a,$‘)-1b:k) (i > k),

(1.5)

commencing with A (l) = A , b(l) = 6. After (n - 1) such steps, we have a system of linear equations whose matrix is upper triangular and which can be solved by the back-substitution n

c

xk=(a:“,‘)-’ bik’-

j=k+l

i

1.3. Numerical

all;‘xj ,

k=n,

n-l

,...,

1.

(1.6)

i

stability

The operations (1.4) are numerically stable in the presence of roundoff provided that the numbers do not grow drastically in size. That growth can cause loss of accuracy may be seen by considering the four-decimal example 1.234 + 402.1 = 403.3,

(1 J)

where it is immediately apparent that the last two decimals of 1.234 have been ignored. often supposed that numerical stability requires that no “pivot” d,,

It is

(14

= ai;’

be small and no “multiplier” ljk=af~)(a$))-‘,

i>k,

be large, but actually neither condition is necessary. Instead, a ,‘,“‘( a II;,‘) - ‘ai,) whole product become large; controlling multiplier is a possible means to this end.

(1.9) it is important to avoid letting the the size of the pivot or of the

with triangular factorization

I. 4. Equivalence

If we use the notation u

kJ

=

d-la(k.) kk

kJ

3

j
(1 .lO)

390

J. K. Reid / Aigebralc aspects of finite-element

together with (1.8) and (1.9) yield the relations

we may collect

all the operations

solutions

(1.4) that are applied

to u,, to

J-1 a!!) = 1.1

1 .d = ‘J JJ

a

‘J -

a(‘) II = d,, = arr -

c k=l

c

‘rkdkkUkJ,

i

(1

>j?

lrkdkkukr ,

.lla)

(lllb)

k=l

and 1-l

a!‘) = d;,u,, 'J

=

C

aij-

'rkdkkUk~,

i Cj.

(l.llc)

k=l

Since (1.9) and (1.10) imply that left-hand side may be incorporated

the equalities I,, = ukk = 1 hold, in each of the cases the into the sum on the right to give the relation

min(r. J) 0

=

aij-

C

(1.12)

lrkdkkUkJ,

k=l

which in turn may be expressed

in the matrix

form

A = LDU,

(1.13)

where L is the lower-triangular matrix with lower-triangular entries I,J, D is the diagonal matrix with diagonal entries di,, and U is the upper-triangular matrix with upper-triangular entries u,,. 1.5. The effects of roundoff If all the symbols refer to computed quantities, an additional term
yn = min(i,

(1.14)

j) - 1

k=l

and we see that we have actually

factorized

exactly

a perturbed

matrix (1.15)

A+E=LDU,

where E has entries (1.14). Inspection of (1.4) shows at once that the associated error e$Jk)will be In fact, a detailed analysis (Reid comparable with a roundoff on the larger of a!J”’ and a!k+l). ‘J [2b] yields the bound e,, < (3.0l)rnm;x

( a:;’ 1,

(1.16)

J.K. Reid / Algebraic aspects of finite-element

where e (assumed less than 10e3) is the relative precision. of avoiding growth in the size of the matrix entries. 1.6. The symmetric

solutions

391

This result emphasises

the importance

case

If A is symmetric ( ajj = aji), it is readily verified from (1.4) that symmetry is preserved is , a(k) = a!$) i, j > k) and, h ence, that the final factorization may be written I’ ’ ‘J A = LDLT.

(1.17)

If A is also positive definite, Wilkinson [la] has shown that no overall matrix elements can occur, in the sense that the inequality

I aij (k) I G “”

‘3J

growth

in the size of

I aiJ I

(1.18)

holds. This highly satisfactory result means symmetric positive-definite case. I. 7. Interchanges

(that

for unsymmetric

that the algorithm

is unconditionally

stable in the

cases

In the unsymmetric case, it is usual to include row or column interchanges (or both) with the aim of limiting growth. For instance, if row interchanges are introduced to make the inequalities lark’1 >ulaj,k’l,

i>k

hold for a preset parameter

(1.19) u in the range

O
(1.20)

then an immediate max

deduction

1a!;+‘) 1 <

from (1.4) is that the inequality

(1 + u-l) max I a!;’ )

i

holds. A similar result holds if column laii)l

(1.21)

i

>,ulai;‘l,

j>k

hold. For full matrices, there is choosing the largest entry in the advantages in using another value grounds. Note also that for a sparse

interchanges

are used to ensure

that the inequalities (1.22)

little disadvantage in taking u = 1, which corresponds to row or column as pivot. For sparse matrices, there are of u to permit some freedom to choose pivots on sparsity matrix there are fewer occasions for growth to occur. If there

J.K. Rend/ Algehruic uperrs of finite-elemer~~

392

are p, nonzero overall bound

entries

u, ,, i cj,

then columri

j is modified

solurions

in only p, steps and therefore

the

(1.23) may be deduced

from the bound

I. 8. The symmetric

nondefinite

(1.21)

cuse

The symmetric nondefinite case may be treated by including the use of 2 x 2 pivots ‘I:%’ (Bunch [3]). The steps k and li + 1 are combined into one, u:t’ is a 1 x 2 matrix, ai:’ is a 2 x 2 matrix and ai:’ is a 2 X 1 matrix. By using symmetric interchanges, we may either obtain a 1 x 1 iivot that satisfies the threshold test (1.19) or a 2 X 2 pivot that satisfies the corresponding test (Duff et al. [4]): (1.24) and maintains the validity of the bound (1.23). Note that either a large diagonal entry is available to be used as a 1 X 1 pivot or a large off-diagonal entry may be permuted to become the off-diagonal entry of a 2 X 2 pivot that has both diagonal entries small. In the latter case, the pivot will not have large entries in its inverse and therefore is likely to satisfy (1.24). Duff et al. [4] showed that it is always possible to find a pivot providing u < 0.5. In practice, smaller values of u are used for good sparsity preservation. 1.9. Monitoring

stability

It is straightforward p = max 1aCk) f, I. I. 1. h

to monitor

the stability

by computing

the quantity (1.25)

Although this is not often done in codes, it is advisable (except in the symmetric and positive-definite case which is unconditionally stable) because if the bound (1.21) is actually attained, we are likely to have lost all accuracy. I. IO. Ill-conditioning Possible instability in Gaussian elimination should not be confused with possible ill-conditioning of the original problem (1.1). Informally, the condition number (ratio of largest eigenvalue to smallest eigenvalue in the symmetric case and ratio of singular values in the unsymmetric case) bounds the amount by which relative errors in the problem are magnified in the solution. If a solution has been obtained stably, we will have solved a nearby problem exactly and whether this has a solution near the wanted solution depends on the conditioning.

J.K. Reid / Algebruic aspects of finite-element

solutions

393

1.1 I. Iterative refinement Iterative refinement may be used to mitigate these effects and estimate the error in the solution. Given an approximate solution X, a correction to it may be found by solving the equation Ai$=b-AZ.

(1.26)

Note that if the factorization (1.13) is stored, this equation is easy to solve. Note also that if A is ill-conditioned, the residual b - AZ may be small although S is large. It is traditional (see Forsythe and Moler [5], for example) to use double-length arithmetic to compute the right-hand side of (1.26) with the aim of repeatedly correcting until the solution settles down, at which point it will have virtually full single-precision accuracy. However, such a procedure makes no sense unless b and A are represented exactly in single precision or are available to more than single-precision accuracy. We prefer Wilkinson’s suggestion (private communication) to use single-precision arithmetic and regard the size of S as an indicator of the likely error in X. Actually, one can go further (Marlow and Reid [6]) and introduce into the right-hand side of eq. (1.26) pseudo-random perturbations at the level of the uncertainties in A and b so that S indicates the effect of these uncertainties on the solution. 1.12. Forward substitution

and back-substitution

We have emphasized stability in connection with the operations on the matrix, that is the calculation of the triangular factorization (1.13). This is because this is the calculation that may be unstable. The operations performed on the right-hand side during forward substitution and back-substitution give rise to errors that lead to a solution that is exact for slightly perturbed factors L and U and no possibility of numerical instability is present.

2. Band, profile and frontal methods for symmetric positive-definite

systems

2. I. Introduction In this lecture,

we consider

the solution

of the equation

Ax=b.

P-1)

where A is an n x n matrix of the form A =

CK’“’

and K@) is a symmetric columns.

(2.4 and positive-definite

matrix whose entries are confined

to a few rows and

394

J. K. Reid / Algehruic aspects of finite-element

solutions

2.9. Band matrices If the variables can be so ordered that those associated with any one element matrix Kc”) have indices that differ by at most m, it is apparent from (2.2) that u,, = 0 if 1i -j 1 > m, that is that the matrix A is a band matrix of bandwidth 2m + 1 and semi-bandwidth m. It is also apparent by inspection of (1.4) that if Gaussian elimination without interchanges is applied, the band property is preserved in the reduced matrices Ack’. 2.3. Variable-band

matrices

For practical problems, the bandwidth is rarely fixed and Jennings [7] has demonstrated that it can be very worthwhile to exploit a variable bandwidth. It is clear from (1.4) that a zero coefficient ahead of the first nonzero in its row (or ahead of the first nonzero in its column) remains zero throughout the elimination. However, between the first nonzero of a row (or column) and the diagonal there is usually total fill-in. In fact, George and Liu [8a] have shown that if all rows (columns) after the first have a nonzero to the left of (above) the diagonal, they all fill in totally between this nonzero and the diagonal. Code to work with a varying bandwidth is scarcely any more complicated than that for a fixed bandwidth, and both can exploit vector hardware in a straightforward fashion since the inner loop can be expressed as the addition of a multiple of one full row vector to another. 2.4. Variable-band

matrices in auxiliary storage

It frequently happens that the variable-band matrix is too large to fit into main storage and therefore has to be held in auxiliary storage (for example, on disk). If the greatest semi-bandwidth is m and we have main storage for (m + I)( m + 2)/2 variables and a reasonably large input-output buffer, it is straightforward to organize the elimination without writing any entries of the factors to auxiliary storage more than once. This is possible because stage k of the elimination involves only rows k to k + m. The earlier rows can have been already written to auxiliary storage and later rows are not yet required. For bigger problems, Jennings and Tuff [9] suggested grouping the rows into segments, each of which fits into half the available main storage. Each segment in turn is overwritten by the corresponding segment of the lower-triangular factor L and then written to auxiliary storage. To process a particular segment in this way requires reading in turn some of the previously written segments of L. The method works in principle for single-row segments, but will be efficient only if many rows are collected into each segment. 2.5. Assembling

variable-band

matrices

The assembly (2.2) of the element matrices K (e) into the overall variable-band matrix A is reasonably straightforward. If auxiliary storage is not in use, we may hold the rows in sequence, with storage allocated for all the positions between the first nonzero in the row and the diagonal and pointers to the positions of the diagonal entries. The matrices Kce) may be added easily into this structure in any order. If auxiliary storage is in use, it is usual to read a file containing the

solutions

395

accumulating as to compute the being assembled need not be read

many rows of A as will fit largest and smallest index does not overlap the range on this occasion.

J.K. Reid / Algebraic aspects of finite-element

packed matrices K”’ several times, on each occasion into the main store at once. It can be very helpful associated with each matrix K (e): if the range of rows and of indices for K”‘, K”’ mak es no contribution 2.6. Frontal method

Although used in computer programs in the early 196Os, the frontal method was first explained in the published literature by Irons [lo] and the method is hence firmly linked with his name. It depends on the observation that the assembly operations (see (2.2)) a IJ :=a,_+k!‘) ‘J

(2.3)

‘J

do not have to be all completed

before

elimination

operations

(see (1.4))

(2-4)

a 'J := a 'I -a,ka,-,lakj

are commenced. Because addition and subtraction are associative (apart from roundoff effects), the operation (2.4) may be performed before all assembly operations (2.3) for a,, have been completed, although it is necessary for a+, akk and akJ to be fully assembled or the wrong quantity will be subtracted. Irons eliminates each variable as soon as its row and column are fully summed, that is performs all operations (2.4) as soon as row and column k are fully summed. After this, row and column k play no further part in the elimination and may be written to auxiliary storage if this is in use. At a typical intermediate stage, we work only with those rows and columns that have been partially surnmed. These correspond to those variables that are associated with one or more assembled elements and with one or more unassembled elements. These variables lie on the border between the set of assembled elements and the rest. As the assembly and elimination progresses, they form a front that moves across the region. For example, suppose the elements labelled “a” in fig. 2.1 have been assembled. The front consists of variables 4, 5, 6 and 8. When element (4, 5, 8) is assembled, variable 4 becomes fully assembled and can be eliminated so the front reduces to 5, 6 and 8. When element (6, 8, 9) is assembled, variable 9 is added to the front but variable 6 can be eliminated. It is usual to hold the active rows and columns, corresponding

6

Fig. 2.1. A triangular

9

cigar-shaped

region.

to variables in the front, in a full matrix labelled with an index list of associated variables. In the fig. 2.1 example, the lists would successively be (4, 5, 6, 8), (5, 6. 8). (5. 6, 8, 9), (5. 8, 9). In fact. the frontal matrix is very similar to one of the original element matrices K”“, except that its size is much larger in practical problems. Notice that the variables usually leave the front in a different order from that in which they enter it. This is important to the efficiency of the method because it permits the front to remain small. but it does lead to slight programnling complications. An initial pass through the index lists associated with the matrices Kfel .rs needed to determine when each variable appears for the last time (that is, when it can be eliminated). There is also some complication in organizing the actual assemblies and eliminations within the frontal array. Normally, the freshly entering rows and columns are placed at the end of the array and the eliminations are performed from the middle with the revised matrix closed up to fill the gaps.

If a variable is associated with only a single element matrix K(“). it will enter the front with the assembly of Kc’) and will be immediately eliminated because its row and column are fully assembled. In fact, we can take the frontal principal one stage further and perform the elimination operations before assembling the element matrix into the frontal matrix. This is and means that the true cost of adding internal variables to known as “‘static condensation” finite elements is very slight indeed. For instance, the algebraic cost of working with 9-variable biyuadratic elements as shown in fig. 2.2 is virtually identical with that of working with the g-variable element shown in fig. 2.3.

So far, we have assumed that a good ordering of the variables is available for a variable-band solution or that a good ordering of the elements is available for a frontal solution. For an underlying geometry that is long and thin, as in fig. 2.1, it may be straightforward to set up such an order manually. In general. an automatic renumbering scheme is obviously desirable. Available algorithms are based on considering the graph that has a node for each variable and an “level sets”. The first level set is a single edge (i. j) for each nonzero entry a,, and constructing node and each successive set consists of all the neighbours of the nodes of the previous level set that are not already members of that or the previous set. A simple example is shown in fig. 2.4 where the level sets are (1). (2, 3, 4), (5, 6, 7), (8. 9. lo), (11. 12, 13). (14, 15, 16) and (17). If the matrix is permuted (by rows and columns) to correspond to the level sets being taken in order, it has a block tridiagonal form. as illustrated in fig. 2.5. Clearly, it is desirable to have a large

liIIll 0

Fig. 2.2. 9-node biquadratic

element.

u

Fig. 2.3. X-node element

(biquadratic

less term in .x’v’ ).

J.K. Reid / Algebraic aspects of finite-element

solutions

> > I :xxx :xX :xxx

:

good node order.

x x

xx

X

Fig, 2.4. A graph with an obviously

397

X XX X 1 X XXX X X XX X X xx x X xxx X X XX X X XX X X XXX X X XX X X XX X X XXXX X XXX XXXX

Fig. 2.5. The matrix corresponding

to the fig. 2.4 graph.

number of levels with few nodes in each level. This led Gibbs, Poole and Stockmeyer [II] to which is a pair of nodes such that if level sets are propose seeking a “pseudo-diameter”, generated from one, the final level set includes the other and restarting from any node in the final set does not yield a greater number of levels. A good implementation of this algorithm has been provided by Lewis [12]. There remains the problem of choosing the order within each level set. Lewis [12] implements two alternative and reasonably successful strategies, but recent work by Sloan [13] is very promising. He chooses each variable successively to minimize a weighted sum of the level set index and the resulting front size. He recommends that the weight for the size of the front should be double that for the level index. This heuristic is based on the notion of taking account both of the global information present in the level sets and the local objective of narrow front-width. For the frontal method, the element ordering may be deduced from the node ordering by associating each element with its earliest node in the order. 2.9. Comparison

between band and frontal solutions

For a given pivotal sequence, the variable-band and frontal methods do the same arithmetic and the storage demands are often similar. Each requires main storage for the symmetric active

b

Fig. 2.6. A small problem.

Fig. 2.7. The matrix from the fig. 2.6 problem.

398

J. K. Reid / Algehruic aspects of finite-element

solutions

matrix. The maximum front size cannot exceed the maximum bandwidth, but is likely to be similar. Furthermore, the total number of entries stored by the frontal method cannot exceed the total number stored by the band method, but is likely to be similar. A simple illustration of how they can differ is shown in figs. 2.6 and 2.7. Variables 1, 2, and 5 can be eliminated by static condensation within their elements so prior to each main set of eliminations, the fronts are (3, 4, 8, 9, lo), (3, 4, 6, 7, 8, 9, 10) and (6, 7, 8, 9, lo), so the maximum front size is 7. The maximum semi-bandwidth is 10, in row 10. The zeros in positions (8,2) and (9,2) are stored explicitly in the band method, but not in the frontal method, thanks to static condensation. The frontal method has the advantage of working with a full matrix during factorization. The band method has the advantage of working with contiguous row vectors during forward substitution and back-substitution.

3. Substructuring

or the multifrontal method

3.1. The multifrontal

method

It is desirable in the frontal method to keep the size of the front small and that is the principal aim of the automatic ordering scheme that we discussed in section 2.8. However, other orderings may be able to reduce the amount of computation. For example, the algorithm of minimum degree (Scheme 2 of Tinney and Walker [14]) has a good reputation. This chooses each pivot ai:’ successively to minimize the number of off-diagonal entries in the pivot row (the degree in the associated graph). Here, many of the early pivotal steps will often involve disjoint sets of variables and the front size may be large. However, if the front is ordered so that variables involved in a pivotal step are grouped, the frontal matrix will be block diagonal. Great savings are available by storing the blocks separately, each with its own index list. This is the essence of the multifrontal method. We have many fronts, rather than just one. Whereas the basic frontal method adds a single element matrix into the frontal matrix and then eliminates any rows and columns that become fully summed, the multifrontal method adds two or more of the frontal and element matrices together to make a fresh frontal matrix and then eliminates rows and columns that become fully summed. We have already noted that a frontal matrix is very like an element matrix. Both can be stored as full matrices labelled with an index list and we may perform eliminations within an element matrix (static condensation). In fact, the frontal matrices are often called generated elements or generalized elements (Speelpenning [15]). The sequence of additions can conveniently be represented as a tree in which the sons of a node correspond to the elements and generated elements whose matrices are added at a single step. As an example, consider the problem shown in fig. 3.1. The frontal method is represented by the tree shown in fig. 3.2 and a multifrontal ordering is shown in fig. 3.3.

Fig. 3.1. A finite element

problem

J.K. Reid / Algebraic aspects of finite-element

solutions

399

5

A+c?l 3

9

1 2

1 Fig. 3.2. Assembly tree for frontal elimination on the fig. 3.1 example, using the frontal method.

4

10

3

4

Fig. 3.3. Assembly

11

12

5

L 8

6

I

tree for another

ordering.

3.2. Substructuring Any assembly tree can be regarded as a representation of the division of the problem into successively finer sets of substructures. For example, fig. 3.3 actually represents nested dissection (George [16]) applied to the problem shown in fig. 3.1. The first dissection divides the problem into its left and right-hand halves. The next divides into quarters and the final level corresponds to the elements themselves. Any node of the tree corresponds to a substructure consisting of the elements associated with the leaves of the subtree of which the node is the root; for instance, node 14 in fig. 3.3 corresponds to the substructure consisting of elements 5, 6, 7 and 8. The variables of the associated frontal matrix (or generalized element) are those on the boundary of the substructure; algebraically they are the variables in common between the structure and the rest of the problem. 3.3. Choosing assembly trees by dissection George [16] considered the regular nested dissection of a square grid of q2 square elements, with each level being a division into four quarters having the same structure. He showed that 0( q3) operations and 0( q* log q) storage locations are needed in place of the 0( q4) operations and O( q3) storage locations needed by bandwidth methods. Actual gains for practical values of q, though worthwhile, are not enormous because the multiplying factors favour the band technique. For general problems, George and Liu [8b] suggest using the algorithm of Gibbs, Poole and Stockmeyer [ll] to find a pseudo-diameter (see section 2.8) and then take the middle level set as a separator. They then repeat on the halves, and so on. On the regular square problem, they found results that were very nearly as good as those of the dissection of George [16] that he analysed theoretically. 3.4. Choosing assembly trees by the minimum-degree

algorithm

The whole elimination process may be modelled symbolically by working only with index lists. Corresponding to adding element and generated element matrices, we merge the index lists. Corresponding to eliminations, we remove indices from the list. These symbolic operations are

400

J. K. Reid / Aigehraic uspects of .finite-elemmt

solution.,

an order of magnitude faster than the corresponding actual operations because (with careful coding) the work at a tree node is linear rather than quadratic in the lengths of the index lists. Furthermore, the length of the resulting list cannot exceed the sum of the lengths of the lists being merged, so storage will always be available. This allows assembly trees to be built efficienctly from the bottom up, using the algorithm of minimum degree. We start with the index lists of all the elements and calculate the degree of each variable, which is the number of other variables with which it appears in one or more lists. A variable of minimum degree is chosen and its elimination is simulated symbolically. If it appears in only one list, it is removed from that list (static condensation). If it appears in more than one list, the lists are merged, the variable is eliminated from the merged list and father-son tree links are recorded between a new node corresponding to the merged list and the nodes corresponding to the lists that were merged. The degrees of the variables are revised and the process is repeated. In this way, the assembly tree is built from the bottom up. The Harwell codes MA27 (Duff and Reid [17] and TREESOLVE (Reid [2d]) use this approach. 3.5. Ordering the operutions associuted with un ussemh!v tree Given an assembly tree, there is considerable scope for reordering the associated operations. A static condensation must be performed ahead of the assembly of the element in which it appears, but apart form this the static condensations may be performed in any order. It is sensible to perform them when the element is assembled. Similarly, it is sensible to perform all the eliminations associated with a node that is not a leaf immediately its generalized element has been assembled. The assemblies may be reordered, provided no node is assembled before all its sons have been assembled into it. A good ordering is obtained by “postordering” during a depth-first search: each node is ordered when a single backtrack for it is performed. For example, a depth-first search of the tree in fig. 3.3 with priority to the left leads to the order (9, 10, 13, 11, 12, 14, 15). The advantage of such a ordering is that the generated elements required at each stage are the most recently generated ones of those so far unused, which means that a stack may be used for temporary storage. For instance, the order we have just obtained for the tree of fig. 3.3 would lead to successive stack contents (9). (9, lo), (13). (13, 1 l), (13, 11, 12). (13. 14). 3.6. Substructuring Iurge problems Several advantages accrue from substructuring large problems. Two subtrees without any in common may be processed in parallel. If any change is made in one of the substructures then the eliminations associated with it and with the overall structure must be recomputed. but not those associated with the other substructures. Really enormous problems may need several disk packs if treated as a whole, whereas data associated with each substructure and with the overall interconnection may each fit on a pack. If several substructures are identical then the elimination operations may be performed just once on a “master” copy. though the forward and backward substitution operations must be performed for each since they are likely to be different on each substructure.

nodes

J. K. Reid / Algebraic aspects of finite-element

solutions

401

4. Iterative methods 4.1. Introduction We have chosen to emphasize direct methods in these lectures because they have traditionally been used to solve finite-element equations. However, for very large problems they can be very expensive, and an iterative method is an alternative. Indeed, iterative methods have traditionally been used to solve finite-difference equations. 4.2. Classical stationary Classical N#k’

#)

stationary =

b

iterative methods iterative

methods

require

a starting

vetor X(O) and take the form

A#-‘),

_

= &+l)

(4.6)

+ #U,

(4.7)

for k = 1, 2,. . . , where N is chosen so that eq. (4.6) is particularly easy to solve. Note that iterative refinement (section 1.11) fits into this pattern if N is taken as the computed triangular factorization of A. By subtracting (4.6) from the corresponding equation with k replaced by k + 1 we find the equation = (N - A)Sck’,

N,j’!=+”

which can be rewritten 8(k+l)

= TS’k’

in the form (4.9)

>

where T is the iteration T = N-‘(N

matrix

-A).

(4.10)

It is equally easy to show that similar equations x

_

X(k+l)

(4.8)

= T(

x

_

#))

apply to the errors (4.11)

and residuals b -

AXck+l) = T(b - A#)).

(4.12)

The rate of convergence depends on the size of the largest eigenvalue of T. Intuitively, the requirement is for N to be a good approximation to A. It is the presence of T in eqs. (4.9) (4.11) and (4.12) as a fixed multiplier that leads to the methods being called “stationary”. The simplest iteration of this kind is Jacobi’s, in which N is the diagonal of A (point Jacobi) or a block diagonal submatrix (block Jacobi).

402

qf finite-element

J. K. Reid / Algebraic aspects

solutions

The next simplest is Gauss-Seidel iteration, where N is the lower-triangular part of A (point Gauss-Seidel) or a block lower-tri~gular submatrix (block Gauss-Seidel). If the diagonal (or block-diagonal) part of N is divided by the scalar factor o, this has the effect of scaling up each component (or block component) of the change as it is calculated, and the resulting method is called successive over-relaxation (SOR). If A is symmetric and positive definite, SOR can be shown to converge for any w in the range 0 < ti -C 2 (see Young [lg]). 4.3. Conjugate gradients Like the stationary methods, requires a starting approximation

the method of conjugate gradients (Hestenes and Stiefel X(O) to the solution. If the corresponding residual is

r(@ = b _ AX’@, the kth iterate

(4.13)

xCk) differs

from x (O) by .~ a vector

S, = span( p(O), Ad*‘,

y(k):,&AX(k),

from the Krylov

A*r”“, . . _ , Ak- ‘do) ).

The choice is such that the corresponding

are orthogonal,

[19])

k=O * I , 2 ,...

subspace (4.14)

residuals (4.15)

that is satisfy the relations

( r(r))T+)

= 0

for

i -f,j_

(4.16)

Again, we assume that A is symmetric and positive definite. Because of this property, the iterates xfk) can be calculated in a very straightforward fashion with the aid of auxilary vectors p” ). with (‘I = r(O), that are “conjugate”, that is satisfy the relations P (p”‘)TAp(J’ At iteration

= 0,

i #j.

k, we calculate (

\

yiii)TrW)

ak= ( p'k')TAp(k)

’

_&k+l)

= X(k)

+ akp(k’,

,.W+V

= p(k)

_ akAp’k’,

(ktl)

T

(kill

= r(k+l)

)

k+l)

) y(

(rW)TrW P

(4.17)

’

+ ~~~‘k’_

(4.18)

J.K. Reid / Algebraic aspects

offinite-element

solutions

403

That the formulae (4.18) yield vectors rCk), pCk) and x(k) - x(O) that lie in the space S, is readily verified. For proofs of the relations (4.16) and (4.17), we refer the reader to Reid [2c]. An important property (see Reid [2c] for a proof) is that of all possible changes to x(‘) that lie in S,, this choice is such the corresponding residual rCkf mini~~es

Q(r)= rTA-‘r. The residual

(4.19)

may be rewritten

in the form

dkf = P,(A)r’O’, where Pk is a polyno~al

(4.20) of degree k satisfying

the equation

P,(O) = 1.

(4.21)

If A has eigenvalues A,, k = 1, 2,. . . , n, and a corresponding set of orthogonal k = 1, 2,. . . , n, then the original residual r(O) may be expressed in the form

4’) = C ykek

eigenvectors

ek,

(4.22)

k and the i th residual

is

r(I)= CykPj(Xk)ek.

(4.23)

k

If there is a polynomial Qi of degree i with Qj(0) = 1 and Qj(hk), k = 1, 2,. . . , ~1, all small, the method of conjugate gradients will be able to choose a polyno~al Pj that gives a small residual r(‘) since it minimizes the quadratic form (4.19). There will be such a polynomial if the range of eigenvalues is small, or if all but a few eigenvalues lie in a small range. Notice that in such cases the method automatically works well without any need to estimate the eigenvalues or any other parameters. 4.4. Preconditioned conjugate gradients Preconditio~ng seeks to transform the problem to one with eigenvalues in a small range so that rapid convergence is thereby obtained. This is usually done by calculating an approximate factorization A = LLT,

(4.24)

and then working (L-lAL-T)y whose solution LTx=y.

with the problem = L-lb,

is related

(4.25) to x by the equation (4.26)

J. K. Reid / Algebraic aspects of finite-element

404

solutions

Meijerink and Van der Vorst 1201 and others have proposed various approximate factorizations. all of which can be viewed as o~tting some of the operations (1.4) of Gaussian eli~nation. They range from making no change to any off-diagonal entry through making no fill-ins and making prescribed fill-ins to making fill-ins only if their numerical value is large. It is hard to draw firm conclusions about which is preferable since numerical experience is of faster convergence being associated with more work per iteration. On the basis of “Occam’s razor”, I prefer the simpler procedures that do not change any off-diagonal entries or make no fill-in. The code of Munksgaard 1211, which is MA31 in the Hat-well library, contains a parameter that can be used to adjust the amount of fill-in from permitting none to permitting all so that the factorization is complete. 4.5. Multigrid

iteration

Multigrid iteration is based on the observation that many traditional iterative methods (see section 4.2) are good at reducing the high-frequency components of the error but poor at reducing low-frequency components. Typically, after a few iterations the residual F(k) =: b _ A#“

(4.27)

is smooth and so can be represented on a coarser grid, say as rj”‘. If A, is a matrix that represents the same differential operator on the coarse grid as A represents on the given grid, then we may solve the equation A #k’ = +Q 1 1

(4.28)

on the coarse grid and then interpolate 6;“’ onto the original grid to provide a correction for xCk) there. This provides a means of speeding the decay of the smooth components of the error. The idea can be nested, that is eq. (4.28) can itself be solved iteratively with the help of an even coarser grid, etc., and very economical overall solution is obtained. The total amount of work performed on all the coarse grids is often less than for one fine grid iteration and the overall convergence rate can be very fast and independent of the fineness of the original grid. Notice that the basic iteration does not necessarily have to be convergent. It has to be able to reduce the high-frequency errors, but does not necessarily have to reduce the low-frequency ones, since those are handled on the coarser grid. It is therefore often called a “smoothing” iteration. The Gauss-Seidel iteration (see section 4.2) in either its point or block form has been widely used as a smoother with success, and there is usually no advantage in replacing it with SOR. The Jacobi iteration can also be used, but Brandt [22] reports that it often has to be under-relaxed (that is the change is scaled down) and it is not as successful as Gauss-Seidel. Another possibility (Wesseling [23]) is to use an incomplete factorization of A (as in section 4.4) for the matrix N in eq. (4.6) and thereby construct a smoothing iteration. Other choices that have to be made at each level of the multigrid method are the number of smoothing iterations to be performed before going to the next grid, the number of multigrid iterations to be performed there, and the number of smoothing iterations to be performed afterwards. The method seems to be very tolerant of the exact choice. Indeed Wesseling [23] concludes in favour of about the simplest possible, namely 0, 1, 1 for the three numbers.

J.K. Reid / Algebraic aspects of finite-element

solutions

405

For fuller discussions of the technique, we refer the reader to Brandt [22], Stiiben and Trottenberg [24], and Hackbusch [25a]. Nicolaides [26] considers grids obtained by repeated dissection of the region and shows the method in a very favourable light. It has been used with success by Bank and Sherman [27] in connection with selective mesh refinement near singularities.

5. Eigenvalue problems 5.1. Introduction In this lecture we consider

the eigenvalue

problem

Ax = XMx,

(5.1)

where A is symmetric and M is symmetric and positive definite. Usually only a few of the smallest eigenvalues A and the corresponding eigenvectors x (the fundamental modes) are wanted, so we concentrate on these. The problem is nonlinear, so some form of iterative method is certainly needed. 5.2. Sylvester’s inertia theorem Several of the methods (A -pM)

that we will discuss involve

the factorization

= LDLT

(5.2)

for a suitable scalar shift p, where L is unit lower triangular and D is block diagonal with blocks of size 1 and 2. Sylvester’s inertia theorem shows that the number of eigenvalues of the problem (5.1) that are greater than, equal to, and less than p are respectively equal to the number of eigenvalues of D that are positive, zero, and negative. The eigenvalues of D are the diagonal blocks of size 1 and the eigenvalues of the diagonal blocks of size 2. These are always easy to calculate, so whenever we form a factorization (5.2) it is trivial to discover how many eigenvalues are less than p. This property is the basis of the method of bisecton, where an interval known to contain the kth eigenvalue is halved by performing a factorization with p at the midpoint, which tells us whether the k th eigenvalue is to the left or right of it. However, in this simple form the algorithm would be too expensive for our application. 5.3. Inverse iteration The eigenvalue iteration

nearest

(A -pM)~(~+‘)q~+~

to p and corresponding

= Mx’“‘,

eigenvector

may be calculated

by the inverse

(5.3)

J. K. Reid / Algebraic aspects of finite-element

406

where qk+l (1

is a normalizing

(k+i))TMX(k+i)

factor chosen,

for instance,

solutions

so that _t~(~+i) satisfies

the condition

= I.

(5.4)

This works well if the wanted eigenvalue is simple and p is much nearer to it than any other eigenvalue, since it can be shown that, if the eigenvalues in order of nearness to p are A,, h 2,...r each iteration reduces the error asymptotically by the factor

A, -P

I I h-P

(5.5)

.

Note that the factorization (5.2) may be formed once and then used repeatedly. Indeed it can be used to find several eigensolutions. Once a approximation y, has been found for the that includes eigenvector corresponding to eigenvalue X1, a new iteration may be commenced periodic M-orthogonalization with respect to yi, thus T(k) = n(k) - ( y:Mx(k))Y,f The resulting

iteration

will converge

(5.6)

to the second

eigensolution

asymptotically

at the rate

A, -P

I I

(5.7)

X3-P.

In general, we M-orthogonalize with respect to all known eigenvectors. Jensen [28] uses this technique, restarting whenever the rate of convergence is unsatisfactory and using the progress of the iterate to judge a more appropriate position for p. He also makes use of Sylvester’s law of inertia (section 5.2) to ensure that all wanted eigensolutions are found. 5.4. Simultaneous

iteration

In order to cope with multiple eigenvalues and clusters of nearly equal eigenvalues, Jennings [29] and Rutishauser [30] generalized inverse iteration by working with simultaneously. Iteration (5.3) is replaced by (A - PWX

(k+l)@+l)

where now Qckil)

= Mx’k’,

is chosen

(xW+l))TMx(k+l)

=

I,

so that Xck+i) satisfies

Clint and r vectors

(5.8) the condition (5.9)

so that the columns of Xck+‘) are M-orthonormal. Note that condition (5.9) is a generalization of condition (5.4). To extract individual eigenvalues and eigenvectors requires the application of a Rayleigh-Ritz procedure using the columns of Xck) as a test space. This involves solving the eigenvalue problem (X’k’)TAX(k)u = yu,

(5.10)

401

J.K. Reid / Algebraic aspects of finite-element solutions

whose order is r, the number of simultaneous iterates. Each y approximates an eigenvalue of the original problem and the corresponding approximate eigenvector is XCk)u. For A,, each iteration now reduces the error asymptotically by the factor A,-P

I

x r+l -P

I

(5 .ll)

.

It is clear that a good rate of convergence is available by choosing expression (5.11) to be small for all eigenvalues of interest. 5.5. The algorithm

Y sufficiently

large

for

of Lanczos

The method of Lanczos [31] is closely related to the method of conjugate gradients that we described in section 4.3. In this section we describe its application to the simple problem (5.12)

Ax = xx. Starting

with an initial vector

V(O),an orthonormal

Sk = span( U(O),Ad”, A2do), . . . , A%(‘)), is generated P ktl

set of vectors

from the Krylov

k = 0, 1, . . .

spaces (5.13)

by the recurrence

Z)(k+l)

=

Au(k)

_

ak27 (k)

_

bkD(k-l),

kc

0

)

1

(5.14)

)...

where ak

PO=0

=

k = 0, 1,. . . ,

(u~~))~Au~~),

(5.15)

1, . ..) is chosen so that the relation

and Pk+i(k=O, ( VW+iJ)TU(k+i) = 1

holds. The recurrences

(5.14) can be collected

AV = VT + ,Ok+lEk,

(5.16) into the matrix

equation (5.17)

where V is the matrix having columns o(O), &), . . . , dk) E, is the matrix that is zero except for the last column which is the vector yCk+i) and T is the kdiagonal matrix

p2. a2

(5.18)

J. K. Reid / Algebraic osprcts

408

offinrte-element solutions

If T has an eigensolution (5.19)

Tu = pu,

and we multiply AVu

eq. (5.17) by u, we find the equation

= VTu +

&+,~p(~+~

= pVu i &+,uxdk +‘I.

(S.20)

where uA- is the k th component of u. If pk+ ,uA is small, it follows that p is an approximate eigenvalue of A and Vu is the corresponding approximate eigenvector. It is relatively cheap to compute the required eigensolutions of T since it is tridiagonal. To obtain corresponding eigenvectors of A we must either store all the columns of V, perhaps on auxiliary storage, or recalculate them by running through the recurrence (5.14) again with the known values for crk and Pk. If we multiply (5.17) by VT and use the orthogonality of the columns of V (so that VTV = I) we find the relation (5.21)

VTAV = T.

This implies that the eigensolutions of T give Rayleigh-Ritz approximations to the eigensolutions of A, just as in simultaneous iteration (section 5.4). If it were not for roundoff errors, the orthonormal set { Y(‘). i = 0, 1, 2..} must terminate with of T correspond to eigensolutions of A. P k+l = 0 and at this point all the eigensolutions Unfortunately roundoff errors mean that the set of vectors zi(‘) IS not accurately orthogonal. However we did not use the orthogonality of V to obtain our result (5.20). This remains accurately valid in the presence of roundoff (see Paige [32]), so whenever / ,Ox+ luk 1 is small we can obtain an accurate eigensolution of A. What we cannot ensure is that we obtain all the eigensolutions of A and that we do not obtain the same eigensolution more than once. Practical experience (see Paige [32], for example) with the method is that the end eigenvalues are found first and are often found more than once before all the middle eigenvalues are found. That the end eigenvalues are likely to be found first is a natural consequence of the Rayleigl-Ritz relation (5.21). For instance, we will have found the vector in S, that maximises the Rayleigh quotient (5.22)

vTAu/vTz;.

That spurious additional copies of end eigenvalues appear is a consequence of orthogonality in the sequence (u(O), &), . . .) not being maintained. The fact that the end eigenvalues are found first is advantageous in that it is usually these that are wanted. 5.6. The block Lanczos algorithm Cullum and Donath algorithm that replaces

[33] and Underwood [34] have suggested a generalization of Lanczos’ columns. The the vectors uck) by n x Y matrices Vck’ with orthogonal

409

J. K. Reid / Algebraic aspects of finite-element solutions

scalars (Ye and Pk become YX r matrices and T becomes Rayleigl-Ritz principle still applies, now with test space spar-0

a block

tridiagonal

matrix.

(0), AV’O’ ,...,AkVco)).

The

(5.23)

It is usual to stop at quite a small value of k and restart using the latest estimates for the wanted eigenvectors as columns of the new V (‘) . There is no difficulty over calculating a cluster of eigenvalues provided they number less than r in total, including multiplicities, whereas the simple Lanczos method is slow to resolve a cluster and (if roundoff is ignored) finds only one representative of a multiple eigenvalue. 5.7. Application

of Lanczos’

Lanczos’ method factorized as

method

may be applied

to the generalized

M = UTU by working

(5.1) provided

M is

problem

= yy.

In this case, the eigenvalues

(5.25) are related

(A-p)-‘=y

and the eigenvectors

problem

(5.24)

with the transformed

U(A -PM)-‘US

eigenvalue

by the equation (5.26)

by the equation

ux=y.

(5.27)

The eigenvalues of the original problem near to p now corresponding to the end eigenvalues of the transformed problem, which the Lanczos algorithm finds the most rapidly. Ericsson and Ruhe [35] use several shifts p, choosing a new shift only when the Lanczos method is making poor progress and using the eigenvalue estimates provided by the Lanczos method to choose the next value of p. They use Sylvester’s inertia theorem (section 5.2) to check exactly where each shift p lies within the whole set of eigenvectors.

6. Nondefinite and unsymmetric

matrices

6.1. Introduction In lectures 2, 3, and 4, we assumed that we were solving a system of linear equations whose matrix was symmetric and positive definite. Most finite-element problems are of this nature, but not all. In this lecture, we explain extensions to the cases where the matrices do not have these properties.

410

J. K. Reid / Algebrmc aspects of finite-element

solutiom

6.2. Bund methods Band and variable-band methods (sections 2.2 and 2.3) may be used at the risk of instability for symmetric nondefinite matrices. It is perfectly practical to monitor the stability by calculating max ] a:,“’ ] (see section 1.9). If the nondefiniteness is because we are handling a matrix A - pM in the course of an eigenvalue calculation (lecture 5) it is perfectly feasible to use a slightly different shift p. Similarly. band and variable-band methods may be used for unsymmetric matrices at the risk of instability, but max 1a:,“) 1 may be monitored. Of course, we no longer have the advantage of being able to store and calculate only half the matrix. If interchanges are included for numerical stability, the form of the matrix is no longer preserved. If only row interchanges are used, the variable (or fixed) band form is preserved in the lower-triangular part, but it widens in the upper-triangular part. For a fixed bandwidth matrix, the upper-triangular semibandwidth at most doubles, and many codes have been written which allocate storage for exactly double and hold explicit zeros at the heads of columns whose lengths do not increase as much. It would be possible to write code that dynamically allocated more space for the columns of a variable-band matrix whenever necessary, but we do not know of one. 6.3. The frontal method Similarly, the frontal method can be used for a symmetric nondefinite or unsymmetric matrix if no interchanges are included. We would again recommend monitoring the stability. However, it turns out to be relatively easy to include interchanges, as was first proposed by Hood [36]. Hood treats unsymmetric matrices, but assumes that each element matrix Kce) has a symmetric structure, that is its row index list is identical to its column index list. He suggests continuing assembly until the front exceeds a critical size, and then performing eliminations until the front no longer exceeds the critical size. Each pivot is chosen as the largest entry in the submatrix of fully assembled rows and columns. Because pivots are not necessarily on the diagonal, unsymmetric permutations are introduced and the row index list for the front in general differs from the column index list. The idea has been developed by Cliffe et al. [37] and Duff [38a,b]. Since separate index lists are needed for the front, there is no difficulty over allowing separate lists for the element matrices K”‘, thereby allowing their structure to be unsymmetric. Indeed the Harwell code MA32 (Duff [38a]) even allows finite-difference equations to be solved since it permits each row of A to be treated as an element matrix K”‘. Furthermore, Duff obtains much better sparsity preservation by using threshold pivoting (section 1.7), that is requiring the pivot uj&’ to satisfy one of the inequalities (1.19) and (1.22) for a given threshold U. Duff [38c] reports that, in practice, it is usually possible to perform eliminations in most of the fully-summed rows and columns so that the front is little bigger than it would have been in the symmetric positive-definite case. Some rows and columns may be fully summed but not taken as pivotal because they contain large entries away from the submatrix of fully summed rows and columns. At a later stage in the elimination, there will be more choice within such a row or column, so that it may be possible to choose a pivot satisfying the threshold test. It is certainly possible eventually because in the end the whole frontal matrix is fully summed.

J.K. Reid / Algebraic aspects of finite-element

411

solutions

In the symmetric nondefinite case, we may preserve symmetry while obtaining a stable factorization by using a mixture of 1 x 1 and 2 X 2 pivots (see section 1.8). Again, we must restrict the pivot to the submatrix of fully-summed rows and columns so the front size may be greater than it would have been in the definite case. 4.4. The multifrontal

method

The multifrontal method (lecture 3) may be generalized in just the same way as the frontal method. Now the generated matrices may be a little bigger than they would have been in the symmetric positive-definite case, but the organizational aspects are otherwise identical. The Harwell code MA27 (Duff and Reid [17a]) is able to treat nondefinite symmetric matrices in this way and a variant MA37 treats unsymmetric matrices with a symmetric pattern. Both codes choose their assembly trees by the algorithm of minimum degree (see section 3.4). 6.5. Classical stationary

iterative methods

The classical iterative methods (section 4.2) can be applied when A is unsymmetric. Gauss-Seidel iterations are convergent (Varga ref. [39], page 73) for diagonally matrices, that is when the condition

holds. The condition

that the iteration

is convergent

if all the eigenvalues

Jacobi and dominant

of the iteration

T=N-‘(N-A)

matrix (6.2)

are less than unity in modulus useful in practice.

holds in the general

situation,

but this result

may not be very

6.6. Conjugate gradients If A is unsymmetric, breaks down. Instead,

the theoretical foundation of the conjugate gradient Vinsome [40] suggested minimizing the quantity

method

of section 4.3

rTr

for residuals ,.(k) =

b

_

AX(k)

that still correspond

64 to iterates

x(“) differing

from

x(O) by a vector

from the Krylov

subspace

412

J. K. Reid / Algebraic aspects of finite-element

Vinsome orthogonal equation

P

uses a set of vectors p’“’ which set. To maintain this orthogonality

(k +I) = r(k+l) -

.volutions

are such that {Ap”‘, i = 0, 1, . . } forms an he replaces the last equation of (4.18) by the

(6.6)

,~"&,P(ll.

where

(6.7)

The second and third equations

of (4.18) remain

with (Yechosen

as

(,.(b)TAp”“

(6.X) The principle difficulty is that all the previous vectors p”’ are needed in eq. (6.6) so must be stored. However, Vinsome proposed its use in conjunction with preconditioning by incomplete factorization, that is working with the problem (LU)-‘Ax

= (LU) - ‘6,

which often leads to a small number discarding some of the vectors p(l).

(6.9)

of iterations.

The algorithm

may also be modified

by

6.7. Multigrid iterution Multigrid iteration is also applicable to unsymmetric problems. A possible source of difficulty lies in the choice of a suitable smoothing iteration. Most of the available theoretical results assume that A is diagonally dominant so that Gauss-Seidel iteration is convergent and therefore suitable as a smoother.

References [I] a. b. [2] a. b. c.

J.H. Wilkinson, Error analysis of direct methods of matrix inversion. J. ACM 8 (1961) 281-330. J.H. Wilkinson, The Algebraic Eigenvalue Problem (Oxford University Press, London, 1965). J.K. Reid, ed., Large Sparse Sets of Linear Equations (Academic Press. New York. London. 1971). J.K. Reid, A note on the stability of Gaussian elimination. J. Inst. Maths. Appl. 8 (1971) 374-375. J.K. Reid, On the method of conjugate gradients for the solution of large sparse systems of linear equations. In ref. [2a], 231-254. d. J.K. Reid, TREESOLVE. A Fortran package for solving large sets of linear finite-element equations. Report CSS 155. Computer Science and Systems Division, AERE Harwell (1984).

J.K. Reid / Algebraic aspects of finite-element

solutions

413

[3] J.R. Bunch, Partial pivoting strategies for symmetric matrices, SIAM J. Numer. Anal. 11 (1974) 521-528. [4] I.S. Duff, J.K. Reid, N.Munksgaard and H.B. Neilsen, Direct solution of sets of linear equations whose matrix is sparse, symmetric and indefinite, J. Inst. Maths. Applies. 23 (1979) 235-250. f5] G. Forsythe and C.B. Moler, Computer Solution of Linear Algebraic Equations (Prentice-Hall, New Jersey, 1967). [6] S. Marlow and J.K. Reid, Fortran subroutines for the solution of linear equations, inversion of matrices and evaluation of determinants. AERE R6899, HMSO, London (1971). [7] A. Jennings, A compact storage scheme for the solution of symmetric linear simultaneous equations. Computer J. 9 (1966) 351-361. [8] a. A. George and J.W.H. Liu, Some results on fill for sparse matrices, SIAM J. Numer. Anal. 12 (1975) 452-455. b. A. George and J.W.H. Liu, An automatic nested dissection algorithm for irregular finite-element problems, SIAM J. Numer. Anal. 15 (1978) 1053-1069. [9] A. Jennings and A.D. Tuff, A direct method for the solution of large sparse symmetric simultaneous equations, in: ref. [2a], 97-104. [lo] B.M. Irons, A frontal solution program for finite-element analysis. Intern. J. Numer. Meth. Engng. 2 (1970) 5-32. [II] N.E. Gibbs, W.G. Poole Jr. and P.K. Stockmeyer, An algorithm for reducing the bandwidth and profile of a sparse matrix, SIAM J. Numer. Anal. 13 (1976) 236-250. [12] J.G. Lewis, Implementation of the Gibbs-Poole-Stockmeyer and Gibbs-King algorithms, ACM Trans. Math. Softw. 8 (1982) 180-189 and 190-194. (131 S.W. Sloan, An algorithm for profile and wavefront reduction of sparse matrices, Intern. J. Numer. Meth. Engng. (1986) 239-251. 1141 W.F. Tinney and J.W. Walker, Direct solutions of sparse network equations by optimally ordered triangular factorization. Proc. IEEE 55 (1967) 1X01-1809. [15] B. Speelpenning, The generalized element method. Report UIUCDCS-R-78-946, Department of Computer Science. University of Illinois at Urbana-Champaign, Illinois (1978). [16] A. George, Nested dissection of a regular finite-element mesh, SIAM J. Numer. Anal. 10 (1973) 345-363. [17] a. I.S. Duff and J.K. Reid, MA27 - A set of Fortran subroutines for solving sparse symmetric sets of linear equations, AERE R10533, HMSO, London (1982). b. IS. Duff and J.K. Reid, The multifrontai solution of indefinite sparse symmetric linear systems, ACM Trans. Math. Softw. 9 (1983) 302-325. [18] D.M. Young, Iterative Solution of Large Linear Systems (Academic Press, New York, London, 1971). [19] M.R. Hestenes, and E. Stiefel, Methods of conjugate gradients for solving linear systems, National Bureau of Standards J. Research 49 (1952) 409-436. [20] J.A. Meijerink and H.A. van der Vorst, An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix, Math. Comput. 31 (1977) 1488162. [Zl] N. Munksgaard, Solving sparse sysmmetric sets of linear equations by preconditioned conjugate gradients, ACM Trans. Math. Softw. 6 (1980) 206-219. 122) A. Brand& Guide to multigrid development, in: ref. [25b] 220-312. [23] P. Wesseling, Theoretical and practical aspects of a multigrid method, SIAM J. Sci. Stat. Comput. 3 (1982) 387-407. [24] K. Stiiben and U. Trottenberg, Multigrid methods: fundamental algorithms, model problem analysis and applications, in: ref. [25b] l-176. [25] a. W. Hackbusch, Multigrid convergence theory, in: ref. [25b] 177-219. b. W. Hackbusch and U. Trottenberg, eds., Multigrid methods. Lecture notes in mathematics, vol. 960 (Springer-Verlag, Berlin, Heidelberg, New York, Tokyo, 1981). [26] R.A. Nicolaides, On some theoretical and practical aspects of multigrid methods, Math. Comp. 33 (1979) 933-952. [27] R.E. Bank and A.H. Sherman, An adaptive multi-level method for elliptic boundary-value problems, Computing 26 (1981) 91-105. [28] P.S. Jensen, The solution of large symmetric eigenproblems by sectioning, SIAM J. Numer. Anal. 9 (1972) 534-545. [29] M. Clint and A. Jennings, The evaluation of eigenvalues and eigenvectors of real symmetric matrices by simultaneous iteration, Computer J. 13 (1970) 76-80.

414

J. K. Reid / Algebraic aspects of finite-element

solutions

[30] a. H. Rutishauser. Computational aspects of F.L. Bauer’s simultaneous iteration method, Numerische Math. 13 (1969) 4-13. b. H. Rutishauser, Simultaneous iteration methods for symmetric matrices, Numerische Math. 16 (1970) 205-223. [31] C. Lanczos, An iteration method for the solution of the eigenvalue problem of linear differential and integral operators, National Bureau of Standards J. Research 45 (1950) 255-281. [32] CC. Paige, Computational variants of the Lanczos method for the eigenproblem. J. Inst. Maths. Appl. 10 (1972) 373-381. [33] J.K. Cullum and W.E. Donath, A block Lanczos algorithm for computing the q algebraically largest eigenvalues and a corresponding eigenspace of large, sparse, real symmetric matrices. Proc. IEEE Conf. on Decision and Control, Phoenix. Arizona (1974). [34] R. Underwood, An iterative block Lanczos algorithm for the solution of large sparse symmetric eigenproblems. Report STAN-CS-75-496, Department of Computer Science. Stanford University, Stanford, California (1975). [35] T. Ericsson and A. Ruhe, The spectral transformation Lanczos method for numerical solution of large sparse generalized symmetric eigenvalue problems, Math. Comput. 35 (1980) 1251-l 268. [36] P. Hood, Frontal solution program for unsymmetric matrices. Intern. J. Numer. Meth. Engng. 10 (1976) 3799400. [37] K.A. Cliffe, C.P. Jackson, J. Rae and K.H. Winters, Finite element flow modelling using velocity and pressure variables, AERE-R.9202. HMSO, London (1978). [38] a. I.S. Duff, MA32 ~ A package for solving sparse unsymmetric systems using the frontal method. AERE R10079, HMSO, London (1981). b. I.S. Duff, Enhancements to the MA32 package for solving sparse unsymmetric equations. AERE R11009, HMSO, London (1983). c. I.S. Duff, Design features of a frontal code for solving sparse unsymmetric linear systems out-of-core. SIAM J. Sci. Stat. Comput. 5 (1984) 270-280. [39] R.S. Varga, Matrix Iterative Analysis (Prentice-Hall, New Jersey, 1962). [40] P.K.W. Vinsome, Orthomin, an iterative method for solving sparse sets of simultaneous linear equations. Sot. Pet. Eng. of A.I.M.E.. SPE 5729 (1976).