Black box multigrid for systems

Black Box Multigrid for Systems J. E. Dendy,

Jr.

Theoretical Division, MS B284 Los Alarms National Laboratory Los Alamos, New Mexico

ABSTRACT This paper extends to systems the black box multigrid method, in which one provides a (logically rectangular) matrix problem and in which the coarser problems necessary to achieve a multigrid solulion are automatically derived. The extension is straightforward if the components of the vector solution of the system are defined at the same grid locations. Otherwise, as in the case of staggered grids, the situation is less satisfactory.

I.

INTRODUCTION

This paper box multigrid nonsymmetric method is that

is an extension to systems of [7] and [8], which consider black for a single scalar positive definite equation and a single equation, respectively. The idea behind the black box multigrid one has to provide a discrete approximation

LMU”

= FM

(1.1)

to a continuous equation LU = F; coarser grid problems necessary to achieve a multigrid solution are then automatically derived. One constraint in the black box multigrid method is that (1.1) must be logically rectangular. Algebraic multigrid [6], which is an extension of black box multigrid, does not have such a constraint and thus can handle a more general class of problems; the price it pays for this generality is a loss in efficiency. In [7] we claimed that the extension to systems would be straightforward. This claim is true if the components of the vector solution of the system live at the same grid locations, and this situation is discussed in Section II. But if this situation does not hold, as in the case of staggered grids, then the

J. E. DENDY, JR.

58

approach of Section II will not work, and an alternative approach is discussed in Section III.

II. NONSTAGGERED

GRIDS

In this section we consider the system LU= F, i.e.,

iclLijUj =F',

i=l

>..., Pa

and its discrete approximation on a grid G” with mesh spacing h,,,:

i.e.,

5 L;(w)*

= (F’ y,

i=l,...,p,

(2.1)

j=l where we assume each (Uj)“‘, j = 1,. . . , p, is defined on G,“. We also assume det L = det( Li j) f 0 and det L”’ = det( L;\:) # 0. The algorithms will be the same as the scalar case with a reinterpretation of the symbols. In its simplest form one constructs a sequence of grids G’ , . . . , G” with corresponding mesh sizes h,, . . . , h, where hi 1 = 2h,. One does a fixed number, IM, of relaxation sweeps on (2.1) and then drops down to grid G”-’ and the equation

where V,“- ’ is to be the (G MP‘)P approximation to V” = 7J” - uAf, where M = U&f is the last iterate on G.“, and where Ii;-’ :(G”)p -+ (G”‘-r)P; the V notation (G k)P means Gkx

... xGk. p times

Black Box Multigrid forSystems

59

To solve (2.2) approximately, one resorts to recursion, taking M-l>k>2, sweeps on G ’ before dropping down to Gk-‘, equation Lk-mlvk-l=

fk-l=

If-‘(fk

_

ID relaxation

and

the

Lkuk)

When grid G’ is reached, the equations L’V’ = f’ can be solved directly and v2 +- v2 + 1FV’ performed. Then one does IU relaxation sweeps on GkP ’ before forming vk + vk + I;_ rvk- ‘, 3 5 k _I M. Let us first consider the case p = 1. In that case [7], Lkpl is defined by Lk-1 = Ik-‘LkIk k_l, where Z~_l:Gk-l+Gk and Ii-‘:Gk-+GkP1. For k simplicity let us first consider the case of symmetric L” in one space dimension. If IF E Gk is the same point as IC E Gk-i, then I[_ 1 is just given by replacement:

(z;_Ivk-‘),F= vk-l IC

If

1) E Gk lies between 1 the template

(IF+

IF+

IC E Gk-’

[-w

’

and (ret

c

1) E Gk-’

and Lk has at

-E],

then

=

(I;_lvk_‘)

WV:;,_’ + Evk:l

IF+1

In other words, one uses the difference that L is the operator

-v.(

c

equation

.

for interpolation.

(2.3)

In the case

Do.)

and L”U

then ukk-

(if 10

k-l

z.z

1

-D

r-1,2L+P3-1,2+D~+1,2)~rD/+1,2~+1~

U,_ 1 = Vlk,’ ),,+

1 =

.z

so

and U,,, = V,“,,‘,), (2.3) is the same as defining that the flux D,_ 1,2( .z - U,_ l) equals the flux

J. E. DENDY, JR.

60

D )+ 1,,2(U,, 1 - z ), and in this case both equal

2Dl-1/2 D1+1,2 D I~-1/e + D, +l/2 thus the “average”

- 17,_1);

G’,+,

diffusion coefficient

is given by I

___

___

+D

1+1/2

The obvious analogue of (2.3) in the one space dimensional p > 1, is as follows: let Lfj have at IF + 1 the template

[

-W,j

Cij

case, but for

-E,j];

then C( I;_l~k-l)IF+l

if C is nonsingular,

(E,j).

let

where

= Wt$

’ + EV”,C+ ’1

C, W, and E are the matrices

(2.4) (Ci i),

( W;j )T and

In the case that

A4 be the matrix (DIj 1,2) and P be the matrix (D;! 1,2), and rewrite

(2.4)-(2.5)

(where

Cl_, = I$-’

and Lr,+ 1= t$>‘,)

as

M(z-U,_,)=P(U,+,-4, which expresses

equality of flues,

or

(2.6)

61

= [ M ( 1\/I+P)

I

‘AI-(M+P)(M+P)

‘A4 CT,-,

+M(n4+P)~‘Pq+, = -P(M+P)~~‘Mq,+M(M+P)

k,

1

=zM(M+P)-‘P~(C.5,,-~7,~,), since P( M + P)- ‘M = M( M + P)

‘P. [To see this, let S = AI + P. Then

=S-P-ll/I+MF’l’

=M(M+

P)

‘P.]

Thus, analogous to the scalar case, the average diffusion coefficient matrix is [i(,Mm’+P ‘)I-‘if MP1 and P ’ exist. The description of the symmetric case is completed by taking 1: ’ = (1: I )*. In case L” is nonsymmetric, we still use (2.4), but we replace C by symm’C = l[( C, j) + (C,;)], W hy symm’ IT’, and E by symm 1;:. Analogous to IS], we define

JL 1 by

if IF + G’ lies between IC E G’ ’ and IC: + 1 E GA ‘, where C* = (C,T ), etc., and we define If-’ = (jt__,)‘, where(A,j)t=(AT,). III two space dimensions we use the reduction to oue space dimension as iu 17, 81. For example, in the symmetric case, if (IF+ 1, JF) E G” lies between

J. E. DENDY, JR.

62

(IC, JC) E Gk -r and (IC + 1, JC) E Gk-‘, - NWjj

and Z,pj has at (IF + 1, JF) the template

- Ntj

- wij

qj

4Wij

- Sj

- NE,, -Eij

,

- SEij

I

then

= (NW+W+SW)U,~~T;,+(NE+E+SE)I$;&, if C - S - N is invertible, where C - S - N is the matrix (Cii - Sji,. N,j), etc. A similar formula holds for (IF, JF + 1) E G” between (IC, JC) E G ’ and (IC, JC + 1) E Gk- ‘. Enough information is now present to use the difference template at points (IF + 1, JF + 1) E Gk which are the centers of rectangles formed by (IC,JC), (IC+~,JC), (IC,JC+~), (ret l,~c+l)~G~~r. reduction to one space dimension is also used in the extension three

space

dimensions One

I”-ILkI” k

dimensions,

but

we

will content

ourselves

with

This idea of [5] of [7] to two

space

in this paper.

simple k_ 1.

detail

of

implementation concerns the that for matrices A, B, C,

computation

of

It is the observation

thus some operations can be saved by forming, for each k and j, I.IZ3klC,j = ( Z3C)ki and then computing A,k(BC)ki for each i. Although it is still a nontrivial calculation, the coding for each piece is already present in the codes in [7, 81 and needs only trivial modification for application in this situation. One other issue is how necessary this definition of the coarse grid operators is. For p = 1, we know that of the possibilities considered in [l], it is the only robust choice. For p > 1, let I;;“_ r be denoted by the matrix Z with block entries Zij; that is,

v; = clijv;-‘. How important

is I, j, i # j? Consider

the one dimensional

(2.7)

case, where L’ii is

Black Box Multigrid

given a)P]which if the should nearly

forSystems

63

by (2.5). Suppose M = aP. Then in (2.6), (M + P)-'M = [(I+ ‘aP = (~(1 + a)-‘Zd, where Id is the identity matrix. In this case, includes the constant coefficient case, Z,j = 0 for i f j. By continuity, jumps in Dk’, k = 1,. ..,p, I= 1,. ..,p,are not severe, then Zij, i # j, be small. (Even if the jumps are large, if they are proportional or so, i.e., M G cxP,then I,,, i # j, is small.) Hence, a reasonable

alternative fine

to Z,k_l:(Gk-l)p-+(Gk)p

and Z/-l:(Gk)p+(Gk-l)J’

1,”_ , based

on

is to de-

block diag( Lyj), (2.8a)

fi

1 based on block diag ( Lyj )t

and to define

Lk-I

=

ik&lLkik k

k-m,;

(2.8b)

this construction presupposes at least that the pointwise diagonal entries of blockdiag(Lfj) are nonzero. In most cases (2.8) gives as good convergence rates as (2.7), but in Section IV we give a pathological example in which this is not the case. Note, however, that both choices avoid special definitions of interpolation near the boundary; see Figure 3 in [7]. The basic relaxation scheme we consider is collective point Gauss-Seidel, with lexicographic ordering, i.e., G k is swept in lexicographic order, with U k being updated at q E Gk so that the residual is zero at q. This requires solution of a p x p system at q. As in the case of p = 1, collective point Gauss-Seidel gives an acceptable smoothing rate unless there is a severe anisotropy (as when AX K Ay), in which case collective line Gauss-Seidel must be used. Also, as in [7, 81 and for the same reasons, the use of the right hand side in interpolation is important. This is implemented as in [7, 81; at points where Zk is not given by just simple injection, C- ‘rk is added to If_ lvk- ‘, k-1 vk-l where C = diag(symm’( Lk)).

III.

STAGGERED

GRIDS

In this section we examine the extent to which the black box approach is feasible for systems defined on staggered grids. The approach used in Section II does not appear to be feasible in this case. We explain one possible

64 approach

J. E. DENDY, JR. by example,

beginning

with the Can&y-Riemann

r:, + VY = I;, UY -v,

in

equations

G?= (0, 1) X (0,l).

= F9‘,,

(3.1)

(r:(X),V(X)),,=G(X),

XE XI,

where (II, V I,, denotes the outward normal of (I:, C’ ), and where C: is subject to the compatibility condition

/0 The grid and positioning is approximated

IT1dxdy =

I

c:ds.


of discrete

variables is shown in Figlae

1, and (:3.1)

by 1,

+

fp\7’l

=

J’l[i’l

-

a’;v’,

= F”

at interior vertices 2,

= c;”

at midpoints of bolmdary links.

!/

(py’),,

I

q

cell centers

qy”’

at

2

The relaxation scheme suggested by (41 for (3.2) is Kaczmarz tive relaxation. Kaczmarz relaxation for a system

(:3.2a) (3.21))

or dist&n-

consists

of, successively

for each i, solving for 6 in

and then replacing xi by Xi + u,,6. If (3.3) is rewritten as Ax = J; then Kaczmarz relaxation is equivalent to Gauss-Seidel relaxation on the system ;\A 7‘y = f. In [4], one does Kaczmarz relaxation first on (:32a), then on (3.21,), and then drops down to the next coarser grid, where this same relaxation pattern is repeated. The process of coarsening is clone by cells, four fine cells constituting one coarse cell. In interpolation, bilinear interpolation is performed separately in the [J’s and V ‘s. It is noted in [4] that Kaczmarz relaxation of (3.2a) does not change the residuals of (3.21)) and vice versa. As an alternative to the above let u” and 0” he the current approximatiom of CT” and V” with /I = h,,. Then we want to solve for correctioils SC”’ and SV” so that

To solve (3.4) approximately,

let us form “AA”“’ globally: _ Al’F /J = ‘I

r

1,

1

that is, foresides (35)

with zero Nemnami boundary conditions; (X5) is an equation \ve can solve approximately by one or more cycles of black 1)0x multigrid. Add the resultant approximate solution, E”, to Su” and 6~” (which can he guessed zero initially) L la Kaczmarz; i.e.,

(3.6)

66

J. E. DENDY, JR.

and then form u’ + uh + au”,

The “AA

” equations

0’ + oh + 6~;“. Now consider

here are

_ @E/I = 2

(3.7)

r;

with zero Dirichlet conditions, which can also he solved approximately by one or more cycles of black box multigrid, and the approximate solution ei is added a la Kaczmarz to Su’* and Sol’ (again guessed initially to be zero). Since the Cauchy-Riemann equations are particularly simple, it worthwhile

considering

another example,

F,

v.U= -vU+vP=F where

the steady-state

in Q = (0,l)

X (0, 1))

(3.8)

in Sz,

U = (U,, U,), and

U(x) = G(x) = (G,b),G,(+, G is subject to the compatibility positioning of discrete variables mated

is

Stokes equations:

XE

ail.

condition j,Fdx = jcjaG. da: The grid and is shown in Figure 2, and (3.8) is approxi-

by

(3.9a)

at cell centers,

_ A”u!’ 1

+ a”p” 1

= F!’ I

at centers of j-faces,

j = 1,2,

(3.9b)

where A” is the usual five point Laplacian. For a point x near the boundary A”UF(x) may involve an exterior value Up. In contrast to [4], for convenience we define this value by linear instead of quadratic interpolation; i.e., &[Up(x’)+ U,l’(r)] = G:,(x”), where xb is the boundary point halfway between r and x’. This relation is then used to eliminate U/(x’) and thus modify the definition of A”Uj’,(x) and the right hand side. The relaxation scheme for (3.9) in [4] is as follows: with p” frozen, perform Gauss-Seidel sweeps on (3.9b) with right hand sides r/’ = FF - dyp” to

Black Box Multigrid for Systems

67

b

FIG. 2.

Discretization

of twdimensiod

52

Stokes eqwtions.

update U,‘, j = 1,2. Then perform Kaczmarz relaxation on (3.9a) in each cell, also adding changes to p” so that the residuals of (3.9b) are unchanged. If xt is the characteristic function at a cell center X, this can be written as

p” t p” - 6AII xx,h where 6 = (1/4h2)r0 = (1/4h’)(F: - E~_,J~U~). At cells next to the boundary A” with zero Neumann boundary conditions is used (see [4]); in these cells it is not true that the residuals rI and r, are unchanged. To use the black box approach here, freeze p” and set up the correction equation

solve this equation approximately by one or more cycles of black box multigrid, and add the approximate solution 6:’ to U/; i.e., U: + u:’ + c!.

J. E. DENDY, JR.

6X Now consider

as in the Cauchy

Kiemann

with zero Neumann box multigrid, Kaczmarz,

case, solve approximately

boundary

conditions,

and add the resultant

with one or more cycles of black

approximate

solution

as in (3.6). But in this case we also replace

the appropriate

modifications

6:; to SUB k la

p” by p” - A%“, with

near the boundary.

To what extent can the above be automated? First, there is the issue of storage of the coefficients of the operators, which must be available for the computation of residuals. The storage scheme we use is as follows: all the coefficients

are stored in one array; for each eqlration a number is assigned to

each mlknown which indicates the type of operator (e.g. A”‘, at, etc.) and the coefficients of the operator are stored in lexicographic order in the array; also required for each unknown is the beginning location. Such a storage scheme at least does not waste storage as would setting aside nine storage locations for each unknown

in each equation.

With the above storage scheme one can also, for the nth equation, indicate which unknowns are to be included in the operator A,,; it is then straightforward to write a routine

to compute

A,< A’,‘, and store its coefficients

in the

format expected by black box multigrid. Also needed are pointers to the operators that indicate how the corrections are to be added. These can coincide with the above storage pointers but may also include additional pointers and operators, as in the Stokes case above, where in addition to the Kaczmarz-like corrections, there is also the p” +- p” - A%” correction. There is a lot to find fault with here. First, in the simple cases considered, where the coefficients are readily computable, the storage penalty is stiff, since one would not ordinarily store the coefficients at all and since one must also store A ,,A7,i for each equation, an additional five to nine locations per grid point. The second complaint is also a complaint against the state of the art of multigrid for systems on staggered grids, and that is that the relaxation schemes appear to lieed hand tailoring. Indeed, the user needs so much more sophistication than in the single equation case that perhaps the name “grey box multigrid” should be used instead. More general systems, such as those considered in 121, would appear to need a combination of the techniques of this section and of Section II. In 121,

delrsity, pressure, aud intennal energy are all stored at cell centers a!lrl must I)e relaxed collectively: this would rule out treating each- equation separately as we were able to do for the Cauchy-Fiemann equations and for the Stokes equations. Some progress, however, has been made since [3] and [2]. In particular, the ideas set forth in Section 3.7 of [S] may lead to a machine tailored approa&.

IV.

NUMERICAL

EXAMPLES

The first example

is the biharmonic

- J[T’ = F

eqllation

in Q = (0,l)

written

as a system:

X (0,l).

311 ati, [:“=,3e’(>~

Yx(l -x)y(l

-y)

(4.1)

011m,

\c-here I; is chosen so that the solution is V’ = :3~‘r, !!u(1 - s)y( 1 - y) and I . ’ = II ’ “. S is discretized by A ,, , the five point Laplacian with mesh spacing h. The results are summarized in Table I. TOL indicates the tolerance to which the system is solved; i.e., the convergence criterion is that the discrete I,” norm of the residuals equals

The Cole \~a\ 11111with the option of I)eginning on the coarsest grid and “ I)ootstrapping” up to the finest grid before cy&ng [7]. The rrnml)ers p are conlprlted I)y the formula log[err(Zh )/err(1l)]/log 2, and L>J’-“w, j = (I/‘211 )( “: , ,, , - n; , ,). Table 1 reports the results when (2.7) is nsed: the c~otlvergence factors are the same to two decimal placles when (2.8) is used. TI~v second example is the hiharmonic equation with more realistic and harder I)ormdary conditions: - A( A[!‘) = F

where

I; is chosen

in 52,

so that I’“(r, y) = sin” ~.rsin’)~~.

This can be discretized

=

10

10

”

”

TOL=

10

”

39 x 39 ( M = 5),

TOL=

lexlR(~~4=4),

TOL

9X9(W==3),

Problem size and parameters

k ” j.2

max

’ f)

” -C+(ih,jlr)/,

- z

a

8.90,

-

5;

7.61, - 4: p = 1.95

2.95, - 3; p = 1.92

-

Irk( ih, jh)

8.97, - 2; p = 1.61

2.74, - 1 1.12, - 2

I\!,)”

2.58, - 2; p = 1.80

#I(

p = 1.99

I; = 1.2

max ,,j

5.45, ~ 4; p = 1.99

2.17, - 3; p = 1.98 3.54, - 4; p = 1.99

8.54, - 3 1.41, ~ 3

I(

[J”

TABLE 1 FOR (4.1)

RESULTS

r,

,

.22, .12;

.21, .11:

.18, .09:

NCYC

NCYC

NCYC

=

=

=

10

9

9

Reduction in discrete L’ norm of the residual in first and last cycles and number of cycles

Block Box Multigrid

71

for Systems

as follows [S]: _ A’iu’ + M”U2 = F 0

,..., (N-l)h)X(h

on!J,,=(h

,...,

(N-l)h),

(4.2a) cl’ - A$12 = 0

A’; is the five point Laplacian

where h = l/N, tions, and i

(4.2b)

on Q2,,,

-2k1

if

(i,j)=(l,j),

j=2

(i,j)=(N-l,j),

,v(; =

*V- 2,

(i,j)=(i.l),

i=2

(i,j)=(i,N-l), - 4h _A

\O

,..., j=2

with zero boundary

condi-

N-2, ,..., N-2,

,...,

N-2,

i=2

,...,

N-2,

if

(i,j)=(l,1),(1,N-1),(N-l,1),(N-1,N-1), otherwise.

The results are summarized in Table 2, where D):,“U,, j = (q+ ,, j The degradation in convergence factors is disappointing. The c: ,, j)/2h. fourth order approximation to ( 6’/&r)U2 is an accident; in other examples we have tried, only second order approximation to (d/&r) U2 is achieved. The third example is an artificial example illustrating that (2.7) is better than (2.8). The problem is

where

-W(L),,W’)

-v.(D,,vU2)

= F’

-v.(D,,vU’)

-v.(D~~vU~)

= F”,

in

U,=U,=O

0n

&?= (0,48)X dQ2:

a, = [0,23] X [0,23] U 123,481 x [23,48],

q,(x.

Y)

1000

if

1

otherwise,

= 4,(X> Y) =

F(x, Y) =

i

(x,y)EQr

999

if

0

otherwise,

1

if

(r,y)EGr,

0

otherwise.

(x,y)EG,,

(0,48)

/-

” =-

If both (2.7) and (2.8) are used for this proMem, and both are i7m for ten cycles, the reduction in the discrete I,’ norm of the residual in the last cycles is 0.45 for (2.7) and 0.77 for (2.8); thus (2.7) is three times faster than (2.8) ill this case. We emphasize that (2.8) is generally nearly as good as (2.7) and that it was, in fact, difficult to concoct this example. The next example pertains to Section III. U7e will omit any example for the Cauchy-Kiemann equations for the following reason: if f;? = 0 in (:3.1) and if 11” and t.” are guessed zero initially, then the residuals of (3.21)) are zero initially; since u ” and u” are always changed in (3.2a) so that the residuals of (:3.21)) do not change, then there is no point in relaxing (3.21)) at all; hence, the Catchy-Riemann equations can be treated as a single equation. If & f 0 in (3.1), (t3.5) can be solved and the u)rrections (:3X) performed, yielding r, = 0, then (3.7) can be sclved and its corrections performed, yielding 1; = 0 and not changing r,. For (3.8) \ve take F so that the solution is I:, = i sins c’os y, Cr, = J C’OS s sin y, P = c’os x’cos y. For ear.11 equation, j = 1,2,3, to he solved by l)lack 1)0x multigrid we take IM = 1, IU = 0, In = 1, and we solve only until iir’il ih less than 0.2 times the initial 11 r ‘lj. it being pointless to solve each equation to great precision. The results are summarized in Table 3. The convergence factors ol)tained compare approximately with the ones in [4]. (Achi Bra& has commented that subsequent munerical experiments indicate that the Its:e of \2’-cycles greatly improves the results in [4]. \V-cycles are of no help in the present sitllatiotl, however.) We remark that the equation obtained from (X91)) is singular: hence, relaxation rather than direct solution was llsed to solve it on the coarsest grid. For a genera1 code oue should use some sort of generalized inverse solution on the coarsest grid, but since the staggered grid &orithnl lacks generality anyway, \ve have not invested this effort.

TAHLE 3 RESULTS

FOR (3.8)

I,’ error, kllltl pXdIllekrS 10 x

over

10

1’

2.94, - 2 2.94. ~ 2

6.30. -- 2

1.17, ~ 1

’

4.42, ~ 3: ,I = 1.48 4.42, 3; p = 1.48

7.95. - 3; p = 1.88 7.9.5, ~ 3: p = 1.88

2.72, ~ 2; p = 1.21

3.63,

1.33, ~ 3; p = 1.73 1.33, - 3; p = 1.73

2.03, - 3: p = 1.96 2.0.3, - 3; p = 1.96

1.05, ~ 2; p = 1.37

1.00, - 2; 1, = 1.8.3

10

10

antI

1.23, - 2

40 x 40 TOL=

r;,

1.23, ~ 2

20 x 20, TOL=

r.;,

’

10,

TOL=

meshpoint

’

last

cycle

ant1 nmnl,er

of cycles

,665:

NCYC

=

25

,718:

NCYC

=

“9

,757:

NCYC

=

30

- 2: p = 1.68

74 V.

J. E. DENDY, JR. CONCLUSION

In this paper we have investigated the extension to systems of the black box multigrid method. For nonstaggered grids, the extension is straightforward, but for staggered grids it is not very satisfactory. In both cases we do not yet know of any physically

meaningful

systems for which the black box

method is desperately needed, that is, problems like the second example in Section IV with highly discontinuous coefficients and strong enough coupling between vector components that they must be treated as a system and not component

by component.

Perhaps we are just ignorant of such problems,

or

perhaps

when they have arisen in the past they have had to be abandoned

because

of the formidable

numerical

difficulties

Thanks and a tip of the hat to Achi Bran&

that they present. Blair Swartz, John Ruge,

and

Burt Wendroff REFERENCES 1

2

3

4

5 6

7 8 9

R. Alcouffe, A. Bran&, J. Dendy, Jr., and J. Painter. The multi-grid method for the diffusion equation with strongly discontinuous coefficients, SIAM J. Sci. Stutist. Comput. 2: 430-454 (1981). A. Brandt, Multigrid solutions to steady-state compressible Navier-Stokes equations. I, in Proceedings of the 5th International Symposium on Computing lMethot1.s in Applied Science and Engineering, pp. I-16. A. Brandt, Multigrid Techniques: 1984 Guide with Applications to Fluid Dynamics, lecture notes for the Computational Fluid Dynamics Lecture Series, von K&m&n Institute for Fluid Dynamics, Rhode-Saint-Genkse, Belgium, 26-30 Mar. 1984. A. Brandt and N. Dinar, Multigrid solutions to elliptic flow problems, in Numerical Methods for Partial Differential Equations (S. Parter, Ed.), Academic, 1979, pp. 53- 147. A. Behie and P. Forsyth, Jr., Multigrid solution of three-dimensional problems with discontinuous coefficients, Appl. Math. Comput. 13: 229-240 (1983). A. Brandt, S. McCormick, and J. Ruge, Algebraic multigrid (AMG) for automatic multigrid solution with application to geodetic computations, SIAM J. Sci. Statist. Comput., submitted for publication. J. Dendy, Jr., Black box multigrid, J. Conput. Phys. 48: 366-386 (1982). J. Dendy, Jr., Black box multigrid for nonsymmetric problems, Appt. Math. Comput. 13: 261-283 (1983). L. W. Ehrlich, Solving the biharmonic equation as coupled finite difference equations, SIAM J. Numer. Anal. 8: 278-303 (1971).