A short introduction to the multi grid packages MG2 and MG3. 1. Introduction. The FORTRAN-90 multi grid packages mg2 and mg3 have been written to solve (large) systems of equations that are discretizations of Laplace-type equations. (strongly elliptic). They may or may not work for other types of equations. The package has two strong points: a. It will work for all type of boundary conditions. (Black box mg) b. It will work for any number of unknowns in all directions. So no need to work in powers of 2 The package has one weak point: The geometry of the region must be a block, if necessary by coordinate transformation. It can be (and has been) used in a multi-block environment. The problem has to be defined by giving the matrix as a 9-point stencil in 2D and a 27-point stencil in 3-D. (There is a way to handle a 19-point stencil and a 7-point stencil is under way). See section 2 Data. The boundary conditions must be given as part of the matrix and the right hand side. 2. Data. All real type data are DOUBLE PRECISION (real (kind(1D0)). There are two types of data involved in the package: matrix data and vector data. A. Vector data Both the right hand side and the solution are structured as 2 (3) dimensional arrays and all indices start at 1. So having nx*ny*nz unknowns all arrays have dimension (1:nx, 1:ny, 1:nz) *corresponding to the unknowns* of the problem. Slack space. Sometimes it is necessary to add some slack space to a vector for instance to accomodate Dirichlet boundary values or to let the array conform to some other data structure. For that reason there are 4 (6) slack variables sl1 .. sl6 and the actual dimensions of the arrays are (1-sl1:nx+sl2, 1-sl3:ny+sl4, 1-sl5:nz + sl6) No assumptions are made about the contents of the slack space and the slackspace is left invariant by the package. B. Matrix Data. The matrix has four indices, of which the first three have to conform to the vector data. (So if the vectors have slack space, so has the matrix) The fourth index contains the stencil map. There is a default mapping: 7 8 9 4 5 6 1 2 3 for 2D 7 8 9 16 17 18 25 26 27 4 5 6 z = k-1 13 14 15 z = k 22 23 24 z = k+1 1 2 3 10 11 12 19 20 21 for 3D It is possible, however for the user to define his own stencil mapping. See subroutine pmap and pmap3. For the 19 point stencil the default numbering is: 5 12 13 14 19 2 3 4 z=k-1 9 10 11 z = k 16 17 18 z = k+1 1 6 7 8 15 but when a 19 point stencil is used the user always must call the remap routine named pmap319. Symbolic names. It is quite difficult to remember the stencil map and for that reason the package offers predefined *symbolic* names, that always are the same: For 2D: p02 p12 p22 p01 p11 p21 p00 p10 p21 First index for x, second index for y, third index for z. So in 3D we have pijk. p020 p120 p220 p021 p121 p221 p022 p122 p222 p010 p110 p210 z=k-1 p011 p111 p211 z=k p012 p112 p212 z=k+1 p000 p100 p200 p001 p101 p201 p002 p102 p202 These quantities are always defined, also in the 19 point stencil, since on coarser grids this stencil also becomes a 27 point stencil. These symbolic names have been defined in the module mgpreset. To include them just use the clause use mgpreset in the calling program. According to the FORTRAN-90 syntax this clause must precede all specification statements. The coefficients of the matrix stencil of unknowns at the boundary *that work on unknowns outside the region* are not used in the package and could contain anything. On the coarser operators they contain zeros conventionally, but on the finest level you might want to put in real matrix data, for instance to calculate the contributions to the right hand side in a Dirichlet boundary condition. The contents of the stencil at points in the slack space are not used. 3. Usage. There are 4 subroutines that form the core of the package. For 2D these are pmap, autosize, prelude and v1step for 3D: pmap3, autosiz3, prelude3 and v1step3. For the 19 point molecule the first is replaced by pmap319 The package is for 3D much the same as for 2D. From now on we only describe the 2D package, mentioning 3D idiosyncrasies only when necessary. If we don't there should be a straightforward analogy. subroutine v1step ( mat, rhs, nsmooth, sol, rbuf, ibuf) mat : matrix of form as described in section 2 rhs : right hand side vector of form as described in section 2 nsmooth : number of smoothing steps on each level sol : solution vector rbuf : real buffer vector containing coarse grid info ibuf : integer buffer vector containing administration info Given the matrix mat and the right hand side rhs, v1step generates an approximate solution based on a multigrid V cycle, that could be used as preconditioner in a gradient method like BiCg-Stab or just in its own right as defect correction. In general it is somewht better to use it as a preconditioner, that saves something like 25% computation time. Nsmooth should be taken 2 in defect correction and 1 or 2 in BiCg-Stab. Higher values are useless. Alternating line Jacobi is used as smoother, 2 directions in 2D, 3 directions in 3D. For v1step to work properly, rbuf and ibuf must have been filled with the proper data. To do that, subroutine prelude must be called: subroutine prelude ( mat, nx, ny, acorsnin, ahilevel, rbuf, + DIMRBUF, sl1, sl2, sl3, sl4, ibuf) or 3D subroutine prelude3 ( mat, nx, ny, nz, acorsnin, ahilevel, rbuf, + DIMRBUF, sl1, sl2, sl3, sl4, sl5, sl6, ibuf) mat i matrix (must have been filled) nx, ny, nz i number of unknowns in various directions acorsnin i method of coarsening. Available are FULLCORS or SEMICORS, defined in the module mgpreset. FULLCORS recommended. ahilevel i highest level of coarsening. As a rule of thumb the next coarser dimension is determined by the formula nxc = nxf/2 + 1 with integer division. We usually coarsen until the smallest dimension becomes 2. After that further coarsening serves no purpose but will not hurt either. rbuf o flat array containing coarse grid operators and LU factorizations on exit DIMRBUF i declared dimension or allocated dimension of rbuf sl1..sl6 i slack parameters (see section 2) ibuf o integer array of at least size 300, containing various data for administrative purposes If the package is used in a multi block environment prelude must be called for every block with a different rbuf and ibuf for each block. There is a piece of data that prelude requires that may pose some difficulties and that is DIMRBUF. The minimum size of the buffer is some very complicated function of nx, ny, nz, ahilevel, acoarsening and the slacks. The easiest way out is to make rbuf an allocatable array and calculate the size dynamically using autosize (autosiz3). subroutine autosize ( nx, ny, acorsnin, ahilevel, + bufsize, sl1, sl2, sl3, sl4) or subroutine autosiz3 ( nx, ny, nz, acorsnin, ahilevel, + bufsize, sl1, sl2, sl3, sl4, sl5, sl6) On exit bufsize contains the minimum size rbuf must have to accomodate all data. Finally the stencil can be remapped using pmap (pmap3 or pmap319). This is a system wide operation, in other words, in a multi-block environment has to be done only once. (Does not hurt to do it more than once though, but the thing to avoid is to use various remaps for different blocks) It has to be done *before* any other call to the package. subroutine pmap (mapvec) subroutine pmap3 (mapvec) subroutine pmap319 (mapvec) Mapvec is a permutation of 1..9 (29, 19) such that the index in ith place will correspond to i in the default stencil. Hence data mapvec /2,3,4,5,1,6,7,8,9/ call pmap (mapvec) will result in the following stencil 7 8 9 5 1 6 2 3 4 and correspondingly p00 will have value 2, p11 will have value 1 and p01 value 6. data mapvec /2, 3, 4, 5, 6, 7, 8, 9,10, + 11,12,13,14, 1,15,16,17,18, + 19,20,21,22,23,24,25,26,27/ call pmap3 (mapvec) will result in the following stencil: 8 9 10 16 17 18 25 26 27 5 6 7 14 1 15 22 23 24 2 3 4 11 12 13 19 20 21 with corresponding values for the pijk. Finally, when a 19 point molecule has to be used, it is mandatory to call pmap19, even if the default numbering is used. Omitting a call to pmap in the 3D case assumes a 27 point stencil with default numbering. To get 1 in diagonal place use the following mapvector: data mapvec /2,3,4,5,6, + 7, 8, 9, 10, 1, 11, 12, 13, 14, + 15, 16, 17, 18, 19/ So the order will be call pmap SYSTEMWIDE, optional call autosize ) for each block, optional but recommended call prelude ) for each block, compulsory after which v1step may be called for this block as often as necessary. For 2D also f1step is available (that implements an F-cycle), but experience has shown that the V-cycle is quite competitive with that and is faster in most cases. (Wall clock time). 4. Libraries and example programs. In the Isnas package all subroutines and the module mgpreset are included. To use the MG package stand-alone there are 4 libraries: MG2 BICG2 MG3 BICG3 The subroutines needed to use the multi grid package are all in the MGi.f (i=2,3) libraries. The BICGi.f libraries contain some extra routines to implement a gradient method, like calculating an inner product and calculating a max norm. (The matrix vector multiplication and vector update are part of MGi.f) There are four example programs: dfc.f bicg.f dfc3.f and bicg3.f that implement a defect correction method and a bicg-Stab method in both 2D and 3D. To trim this to your needs you only have to supply your own fillmat(3) and/or fillvec(3) to feed data into the matrix and the right hand side respectively. Example fillmat and fillvec routines are included.