C For compiling, type f77 mg2d.f to get an executable file a.out C For running the executable file, type a.out C C MULTIGRID SOLVER FOR 2D POISSON EQUATION C----------------------------------------------------------------------- C Multigrid Solver for solving the 2D Poisson equation C Included here are the following files: C mgtest.f Test driver (main program) C mgsolv.f Multigrid solver (subroutines) C mgcv.f V-cycle control routine (subroutine) C C Send questions/comments to: C Jun Zhang C Department of Computer Science C University of Kentucky C Lexington, KY 40506-0046 C phone: (606) 268-2379 C e-mail: jzhang@cs.uky.edu C C Disclaimer: C This program is provided as it is. You use it at your C risk. The author is not responsible for the correctness C of the program and is not responsible for any loss or C damage as a result of using this software. C C last modified on January, 16, 1999 C----------------------------------------------------------------------- c program mgtest C C Program to test multigrid algorithm to solve 2D Poisson equation C implicit none integer ierr, io, j, k, lq, lqq, nx, ny, ipar(1:10) integer n1, n2, n3 parameter ( n1 = 1, n2 = 1, n3 = 1 ) parameter ( nx = 64, ny = 64, lq = 3*(nx+1)*(ny+1) ) real*8 a, b, c, d, ffun, h, rn, tol, ufun parameter ( a = 0.0, b = 1.0, c = 0.0, tol = 1.0e-10 ) real*8 f(0:nx,0:ny), u(0:nx,0:ny), r(0:nx,0:ny), q(lq) real*8 error, errmx real t1, t2, t, etime integer jmx, kmx C C The following statements are for testing the small storage option: C real*8 qq(0:nx,0:ny,2) equivalence ( q, qq ), ( qq(0,0,1), u ), ( qq(0,0,2), f ) lqq = -lq C C Set up right hand side and boundary values on the grid C h = (b - a)/nx d = c + ny*h do 10 k = 1, ny-1 do 10 j = 1, nx-1 u(j,k) = 0.0 f(j,k) = ffun(a+j*h, c+k*h) 10 continue do 20 k = 0, ny u( 0,k) = ufun( a, c+k*h) u(nx,k) = ufun( b, c+k*h) 20 continue do 30 j = 1, nx-1 u(j, 0) = ufun( a+j*h, c) u(j,ny) = ufun( a+j*h, d) 30 continue C C Solve the linear system using the multigrid solver C io = 2 C C Set up parameters C ipar(1) = nx ipar(2) = ny ipar(3) = n1 ipar(4) = n2 ipar(5) = n3 ipar(6) = io ipar(7) = lqq ipar(10) = ierr c t1 = etime(t) c call mgsolv( h, f, u, r, tol, q, rn, ipar ) c t2 = etime(t) c c compute the maximum absolute error c errmx = 0.0 do 50 k = 1, ny-1 do 50 j = 1, nx-1 error = u(k,j) - ufun(a+j*h,c+k*h) if ( abs(errmx) .lt. abs(error) ) then errmx = error jmx = j kmx = k end if 50 continue C print * print * print *, ' Total elapsed time =', t2 - t1 write (6,90) error, jmx, kmx 90 format(//,'maximum error =', e15.8,' at i=', i5,' , j=', i5,//) if ( ierr .eq. 0 ) write (*,*) ' ieer =', ipar(10) C stop C----------------- end of driver -------------------------------- end C double precision function ffun( x, y ) implicit none real*8 x, y, dcos, dexp C C The right hand side function f(x,y) C ffun = 5.2d1*dcos( 4.0d0*x + 6.0d0*y ) C ffun = -( x*x + y*y )*dexp( x*y ) C ffun = -( x*x*(1.0d0 - x*x)*(2.0d0 - 1.2d1*y*y) C & + y*y*(1.0d0 - y*y)*(2.0d0 - 1.2d1*x*x) ) return C-------------- end of right hand side function ----------------- end C double precision function ufun( x, y ) implicit none real*8 x, y, dcos, dexp C C Exact solution u(x,y), also prescribes boundary values C ufun = dcos( 4.0d0*x + 6.0d0*y ) C ufun = dexp( x*y ) C ufun = x*x*y*y*( 1.0d0 - x*x )*( 1.0d0 - y*y ) return C------------- end of exact solution ---------------------------- end C subroutine mgsolv( h, f, u, r, tol, q, rn, ipar ) implicit none integer ipar(1:10) real*8 h, f(0:ipar(1),0:ipar(2)), u(0:ipar(1),0:ipar(2)) real*8 r(0:ipar(1),0:ipar(2)), tol, q(*), rn C********************************************************************** C* subroutine mgsolv * C* MultiGrid Solver for the Poisson Equation * C* (Two dimensions, Dirichlet boundary conditions) * C********************************************************************** C C Description: C C This routine solves the two-dimensional Poisson equation C C -u (x,y) - u (x,y) = f(x,y) (on a rectangle) C xx yy C C u(x,y) = g(x,y) (on the boundary) C C using second order contral finite differences on the grid C x(j) = x0 + j*h, j=0, 1, ..., nx, C y(k) = k0 + k*h, k=0, 1, ..., ny, C where _h_ is the mesh size and (x0,y0) is the grid origin. C The user supplies values of _f_ and _g_ on this grid, and this C routine returns the (approximate) corresponding values of _u_. C C Arguments C C on entry: C C h: real*8, grid mesh spacing in both _x_ and _y_ directions C C f: real*8, array containing the values of the right hand side C _f_ of the Poisson equation on interior of the grid: C f(j,k) = f( x0+j*h, y0+k*h ) (00: Iterate until the residual norm is <= tol C tol=0: Do one multigrid cycle C tol<0: Do abs( tol ) multigrid cycles C C Work space: C C q: real*8, vector or array for work space. Approximately C 8*(nx+1)*(ny+1)/3 storage locations are needed. C If _q_ is not large enough (as indicated by _lq_), C the solution will not be computed, and the exact C length needed for _q_ will be returned in _ierr_. C ***** See also the Small Storage Option below ***** C C ipar: integer, array of length 10, used for passing parameters C C On exit: C C u: real*8, contains the computed solution C C rn: real*8, contains the final residual norm C C********************************************************************* C C ipar(1) = nx, number of grid intervals in _x_ (must be divisible C by 4 C C ipar(2) = ny, number of grid intervals in _y_ (must be divisible C by 4 C C ipar(3) = n1, number of presmoothing sweeps C C ipar(4) = n2, number of postsmoothing sweeps C C ipar(5) = n3, number of coarsest grid smoothing sweeps C C ipar(6) = io, number of postsmoothing sweeps C C io: integer, specifies whether or not output is desired: C io= 0: no output C io=23: multigrid convergence trace C (see comments in mgcv for other options) C C ipar(7) = lq: integer, length of the work space _q_ C C ipar(10) = ierr, error flag, ipar(10) = 0 on normal return C C Error conditions: C C ipar(10) = 1: Called with _nx_ not divisible by 4 C ipar(10) = 2: Called with _ny_ not divisible by 4 C ipar(10) = 3: Called with _h_ zero or negative C ipar(10) = 4: Specified tolerance tol>0.0 not achieved C ipar(10) >= 5: The supplied work space _q_ is too small: the C returned value of IERR is the length needed. C C Small storage option: C C To use the minimum possible storage, force _u_ and _f_ to occupy C the first and second sets of (nx+1)*(ny+1) locations of _q_, C respectively (e.g., using equivalence or common ), and flag this C use by setting _lq_ as minus the length of _q_. In this case the C actual arguments referenced by _u_ and _f_ are not accessed. C CAUTION: With this small storage option, the values of _f_ C stored in _q_ are altered on return (multiplied by h**2). C C Method: C C Solves the problem iteratively by repeated multigrid V-cycles C using Gauss-Seidel relaxation with natural ordering. C C Input/Output: C C Produces printed output as described above for argument _io_ C C Language: Fortran 77 C C Component Routines C C vcopy, vscale (vector copy and scaling routines) C mgrlax, mgftoc, mgctof, mgsolc (interface to mgcv) C relax, ftoc ctof (relaxation and grid transfers) C mgcv (V cycle conrtol routine) C C Author: C C Jun Zhang C Department of Computer Science C University of Kentucky C 773 Anderson Hall C Lexingtin, KY 40506-0046 C E-mail: jzhang@cs.uky.edu C C Last modified on January 16, 1998 C C This code was modified from the one provided to the author by C C Scott R. Fulton C Department of Mathematics and Computer Science C Clarkson University, Potsdam, NY 13676 C C********************************************************************** C Disclaimer: C This software is provided as it is. There is no guarantee associated C with this software. You use it at your own risk. The author(s) is C (are) not responsible for any loss or damage as a result of using C this software. C********************************************************************** C integer j, k, lastq, msize, nc, n1, n2, n3, io integer lq, ierr, nx, ny real*8 d, pi, rhojac, rn3, rntol, w, hh integer mm, m, mx, my, ipu, ipf parameter ( mm = 10 ) common /mgcp2d/ hh(mm), m, mx(mm), my(mm), ipu(mm), ipf(mm) save /mgcp2d/ C nx = ipar(1) ny = ipar(2) n1 = ipar(3) n2 = ipar(4) n3 = ipar(5) io = ipar(6) lq = ipar(7) C C Check consistency of the grid size requested. C ierr = 0 if ( mod( nx, 4 ) .ne. 0 .or. nx .lt. 4 ) ierr = 1 if ( mod( ny, 4 ) .ne. 0 .or. ny .lt. 4 ) ierr = 2 if ( h .le. 0.0 ) ierr = 3 if ( ierr .ne. 0 ) then ipar(10) = ierr return end if C C Set up the grid structure: grid sizes, meshes, and locations in q C lastq = 0 j = nx k = ny d = h m = 0 10 m = m + 1 mx(m) = j my(m) = k hh(m) = d msize = ( j+1 )*( k+1 ) ipu(m) = lastq + 1 lastq = lastq + msize ipf(m) = lastq + 1 lastq = lastq + msize j = j/2 k = k/2 d = 2*d if ( mod(j,2) .eq. 0 .and. mod(k,2) .eq. 0 .and. m .lt. mm ) & go to 10 C C Check to see if the work space q is large enough C if ( lastq .gt. abs( lq ) ) then ierr = lastq return end if C C Compute the number of sweeps needed on the coarsest grid C if ( mx(m) .le. 2 .and. my(m) .le. 2 ) then n3 = 1 else pi = acos( -1.0 ) rhojac = 0.5*(cos( pi/mx(m) ) + cos( pi/my(m) )) rn3 = (n1 + n2)*log( 0.25 )/(2.0*log( rhojac )) n3 = rn3 + 0.99 end if C C Copy the input data into q (if needed) and scale f by h*h C if ( lq .gt. 0 ) then call vcopy( nx, ny, f, q(ipf(1)), 1, nx-1, 1, ny-1 ) call vcopy( ny, ny, u, q(ipu(1)), 0, nx, 0, ny ) end if call vscale( h*h, nx, ny, q(ipf(1)), 1, nx-1, 1, ny-1 ) C C Solve the linear system by multigrid cycling C if ( tol .gt. 0.0 ) then rntol = tol nc = 100 else rntol = 0.0 nc = nint( abs( tol ) ) if ( nc .le. 0 ) nc = 1 end if w = 0.0 call mgcv( m, 1, q, r, n1, n2, n3, rntol, nc, io, rn, w ) if ( tol .gt. 0.0 .and. rn .gt. tol ) ierr = 4 C C Copy computed solution into u (unless using small storage option) C if ( lq .gt. 0 ) call vcopy( nx, ny, q(ipu(1)), u, 0, nx, 0, ny ) ipar(10) = ierr return C--------------------- end of mgsolve ---------------------------- end C subroutine vcopy( nx, ny, v1, v2, j1, j2, k1, k2 ) implicit none integer nx, ny, j1, j2, k1, k2 real*8 v1(0:nx,0:ny), v2(0:nx,0:ny) C C copy values between two 2D arrays C integer j, k C do 10 k = k1, k2 do 10 j = j1, j2 v2(j,k) = v1(j,k) 10 continue return C------------------ end of vcopy ------------------------- end C subroutine vscale( value, nx, ny, v, j1, j2, k1, k2 ) implicit none integer nx, ny, j1, j2, k1, k2 real*8 value, v(0:nx,0:ny) C C Scales a 2D array _v_ by the constant _value_ C integer j, k C do 10 k = k1, k2 do 10 j = j1, j2 v(j,k) = value*v(j,k) 10 continue return C---------------------- end of scale2 ------------------------ end C subroutine mgsolc( l, q, rn, w ) implicit none integer l real*8 q(*), rn, w C C Dummy routine for the coarest grid solution C return C---------------------- end of mgsolc ------------------------- end C subroutine mgrlax( l, q, rn, w ) implicit none integer l real*8 q(*), rn, w C C Front end for relaxation on a single grid C integer mm, m, mx, my, ipu, ipf real*8 hh parameter ( mm = 10 ) common /MGCP2D/ hh(mm), m, mx(mm), my(mm), ipu(mm), ipf(mm) save /MGCP2D/ C call relax( mx(l), my(l), hh(l), q(ipf(l)), q(ipu(l)), rn ) w = w + real( mx(l)*my(l) )/real( mx(1)*my(1) ) return C-------------------- end of mgrlax --------------------------- end C subroutine relax( nx, ny, h, f, u, rn ) implicit none integer nx, ny real*8 h, f(0:nx,0:ny), u(0:nx,0:ny), rn C C Relaxation for 2D Poisson equation: Gauss-Seidel C integer is, j, js, k real*8 r, unew C C Regular Gauss-Seidel sweeps C rn = 0.0 do 10 k = 1, ny-1 do 10 j = 1, nx-1 unew = 0.25*( (u(j+1,k) + u(j-1,k) + u(j,k+1) & + u(j,k-1)) + f(j,k) ) C r = unew - u(j,k) rn = rn + r*r u(j,k) = unew 10 continue rn = sqrt( 16.0*rn/(h*h) ) return C---------------------- end of mgrelax --------------------------- end C subroutine mgftoc( lf, lc, q, r, w ) implicit none integer lf, lc real*8 q(*), r(*), w C C Front end for Fine-TO-Coarse grid transfer C integer mm, m, mx, my, ipu, ipf real*8 hh parameter ( mm = 10 ) common /MGCP2D/ hh(mm), m, mx(mm), my(mm), ipu(mm), ipf(mm) save /MGCP2D/ C call ftoc( mx(lf), my(lf), q(ipu(lf)), q(ipf(lf)), & mx(lc), my(lc), q(ipu(lc)), q(ipf(lc)), r ) return C--------------------- end of mgftoc ----------------------------- end C subroutine ftoc( nxf, nyf, uf, ff, nxc, nyc, uc, fc, r ) implicit none integer nxf, nyf, nxc, nyc real*8 uf(0:nxf,0:nyf), ff(0:nxf,0:nyf), r(0:nxf,0:nyf) real*8 uc(0:nxc,0:nyc), fc(0:nxc,0:nyc) C C Fine-TO-Coarse transfer of 2D Poisson equation C real*8 res integer jc, j, kc, k C C Compute residual on the fine grid C do 10 k = 1, nyf-1 do 10 j = 1, nxf-1 r(j,k) = ff(j,k) - ( 4.0*uf(j,k) - & (uf(j+1,k) + uf(j-1,k) + uf(j,k+1) + uf(j,k-1)) ) 10 continue C C Transfer the interior residual by full weighting scheme C do 20 kc = 1, nyc-1 do 20 jc = 1, nxc-1 k = 2*kc j = 2*jc res = 4.0*r(j,k) + 2.0*( r(j,k+1) + r(j,k-1) + r(j+1,k) & + r(j-1,k) ) + r(j+1,k+1) + r(j-1,k+1) & + r(j+1,k-1) + r(j-1,k-1) fc(jc,kc) = 0.25*res 20 continue C C Set the initial approximate solution on the coarse grid to zero C do 30 kc = 0, nyc do 30 jc = 0, nxc 30 uc(jc,kc) = 0.0 return C------------------------ end of ftoc -------------------------------- end C subroutine mgctof( lc, lf, q, w ) implicit none integer lc, lf real*8 q(*), w C C Front end for Coarse-TO-Fine grid transfer C integer mm, m, mx, my, ipu, ipf real*8 hh parameter ( mm = 10 ) common /MGCP2D/ hh(mm), m, mx(mm), my(mm), ipu(mm), ipf(mm) save /MGCP2D/ C call ctof( mx(lc), my(lc), q(ipu(lc)), & mx(lf), my(lf), q(ipu(lf)) ) return C----------------------- end of mgctof ----------------------------- end C subroutine ctof( nxc, nyc, gc, nxf, nyf, gf ) implicit none integer nxc, nyc, nxf, nyf real*8 gc(0:nxc,0:nyc), gf(0:nxf,0:nyf) C C Coarse-TO-Fine transfer of 2D Poisson equation C integer jc, jf, kc, kf real*8 avg C if ( nxf .ne. 2*nxc .or. nyf .ne. 2*nyc ) go to 900 C C Transfer the correction by bilinear interpolation and add C gf(0,0) = gf(0,0) + gc(0,0) do 10 jc = 1, nxc jf = 2*jc gf(jf-1,0) = gf(jf-1,0) + 0.5*(gc(jc-1,0) + gc(jc,0)) gf(jf ,0) = gf(jf ,0) + gc(jc,0) 10 continue do 30 kc = 1, nyc kf = 2*kc gf(0,kf-1) = gf(0,kf-1) + 0.5*(gc(0,kc-1) + gc(0,kc)) gf(0,kf ) = gf(0,kf ) + gc(0,kc) do 20 jc = 1, nxc jf = 2*jc avg = 0.5*(gc(jc,kc-1) + gc(jc,kc)) gf(jf-1,kf-1) = gf(jf-1,kf-1) + & 0.5*(0.5*(gc(jc-1,kc-1)+gc(jc-1,kc)) + avg) gf(jf ,kf-1) = gf(jf ,kf-1) + avg gf(jf-1,kf) = gf(jf-1,kf) + 0.5*(gc(jc-1,kc) + gc(jc,kc)) gf(jf ,kf) = gf(jf ,kf) + gc(jc,kc) 20 continue 30 continue C----------------------- end of ctof -------------------------------- return C C Error exit C 900 write (*,9000) nxf, nyf, nxc, nyc 9000 format(' ctof error: grid sizes are not compatable:', & /' nxf =',i5,' nyf =',i5,' nxc =',i5,' nyc =',I5) stop end C----------------------------------------------------------------------- subroutine mgcv( lc, lf, q, r, n1, n2, n3, tol, nc, io, rn, w) implicit none integer lc, lf, n1, n2, n3, nc, io real*8 q(*), r(*), tol, rn, w C C MultiGrid V-Cycle control routine C From MGCPAC version 1.07 (12/02/88) C C Description C C This routine controls the execution of one or more V-cycles C of a multigrid method. The user must define the grid structure C and supply subroutines for the relaxation and grid transfers. C C Arguments C C Input C C ************ Grid Structure Definition ************* C C LC, LF Specify the indices of the coarsest and finest C levels, respectively. The level indices may be C defined as desired, as long as: (1) "adjacent" C levels have indices which differ by one, and (2) C there are at least 2 levels defined. Level indices C are passed to user-supplied routines (see below). C C Q Work array to be passed to other routines C (may be used to store data and variables as desired) C C *********** V-Cycle Structure Definition *********** C C N1 Number of relaxation sweeps to make before C going to the next coarser level (nonnegative) C C N2 Number of relaxation sweeps to make before C going to the next finer level (nonnegative) C C N3 Specifies what to do on the coarsest level: C N3 >= 1: make N3 relaxation sweeps C N3 <= 0: call routine MGSOLC to solve directly C C **************** Control Parameters **************** C C TOL Specifies convergence tolerance/stopping criterion. C The problem is said to have converged when the C residual norm on the finest level is less than: C TOL (if TOL>0.0) C abs(TOL)*TERR (if TOL<0.0) C (TERR is supplied through the common block /MGTERR/ C and can be set dynamically--see FAS comments below). C Multigrid cycles are repeated until the problem C converges (precisely NC cycles are made if TOL=0.0). C C NC Specifies the maximum number of cycles allowed C C IO Specifies printed output desired (IO=0 for none). C For printed output, use IO as the sum of values C specifying the desired output options as follows: C 1: level, residual norm, work (each sweep) C 2: convergence factors per cycle C 4: overall convergence factors C 8: convergence tolerance when reset (TOL<0.0) C 16: page header C For example, with IO=18 the output consists of a C page header and convergence factors for each cycle. C C RN Variable used to return the final residual norm C C W Variable used to return the amount of work done C (must be initialized by the user--usually to zero) C C Output C C Q May have been modified by the user-supplied routines C MGRLAX, MGSOLC, MGFTOC, and MGCTOF (see below) C C RN Contains the last residual norm computed on the C finest level (see Required Routines below) C C W Has been incremented by the amount of work done C (in work units--see Required Routines below) C C Error Condition C C If a positive convergence tolerance specified by TOL is not C met within NC V-cycles, an error message is printed. C C Required Routines (to be supplied by the user) C C SUBROUTINE MGRLAX( L, Q, RN, W ) C C ReLAXation: performs one relaxation sweep on level L. C Returns a residual norm RN associated with the sweep C and increments W in work units (if desired). C C SUBROUTINE MGSOLC( LC, Q, RN, W ) C C SOLution on Coarsest level: solves on coarsest level LC. C Returns a residual norm RN associated with solving C and increments W in work units (if desired). C NOTE: This routine is not called unless N3.le.0. C C SUBROUTINE MGFTOC( LFINE, LCOAR, Q, W ) C C Fine-TO-Coarse transfer: performs the transfer C from the fine level LFINE to the coarse level LCOAR C and increments W in work units (if desired). C C SUBROUTINE MGCTOF( LCOAR, LFINE, Q, W ) C C Coarse-TO-Fine transfer: performs the transfer C from the coarse level LCOAR to the fine level LFINE C and increments W in work units (if desired). C C Usage Notes C C The array Q is intended to be used to store the various C variables required for the problem. It is not accessed, but C simply passed to the routines described above. The user can C define the storage within Q as desired. Any pointers needed C to access data in Q must be passed in COMMON (or in Q). C C Use with the FAS method C C In FAS (Full Approximation Scheme) multigrid methods, a good C estimate of the actual local truncation error can be generated C during the solution process by computing the relative local C truncation error between two levels. To use this to test for C convergence, set TOL < 0.0 and include the COMMON block C REAL*8 TERR COMMON /MGTERR/ TERR SAVE /MGTERR/ C C in your routine MGFTOC. Within this routine, set TERR to be C an approximation to the fine level truncation error (norm) as C described above. Then each time the problem is transferred from C the finest level using MGFTOC, the convergence tolerance on the C finest level will be reset to abs(TOL)*TERR. C C Input/Output C C Produces printed output as described above for argument IO C C Language C C FORTRAN 77 C C Author C C Scott R. Fulton C Department of Mathematics and Computer Science C Clarkson University, Potsdam, NY 13676 C integer icycle, inc, isweep, l, ns0, ns1 real*8 cfact, eps, rnfine, rn0, rn1 real*8 sfact, wfact, w0, w1, wu logical first, lsweep, lcycle, lovall, lreset, lphead C C Check the input arguments C rn = 0.0 if ( abs( lf - lc ) .le. 0 ) then write (*,*) ' mgcv error: too few levels:' write (*,*) ' lc = ', lc, ' lf = ', lf return end if C C Initializations C if ( lc .lt. lf ) then inc = +1 else inc = -1 end if eps = tol rn = 0.0 first = .true. lsweep = mod( io , 2 ) .ne. 0 lcycle = mod( io/2 , 2 ) .ne. 0 lovall = mod( io/4 , 2 ) .ne. 0 lreset = mod( io/8 , 2 ) .ne. 0 lphead = mod( io/16, 2 ) .ne. 0 if ( lphead ) write (*,2000) lf, lc, n1, n2, n3, tol, nc C C ************************ V-Cycle Algorithm ************************** C do 100 icycle = 1, nc if ( lsweep ) write (*,*) C C Downward branch: fine to coarse C do 20 l = lf, lc+inc, -inc do 10 isweep = 1, n1 call mgrlax( l, q, rn, w ) if ( l .eq. lf ) then if ( first ) then rn0 = rn ns0 = -1 w0 = w first = .false. rn1 = rn0 ns1 = ns0 w1 = w0 end if rnfine = rn ns0 = ns0 + 1 ns1 = ns1 + 1 end if if ( lsweep ) write (*,2100) l, rn, w 10 continue call mgftoc( l, l-inc, q, r, w ) if ( l .eq. lf .and. tol .lt. 0.0 ) then eps = abs( tol )*terr if ( lreset ) write (*,2400) eps end if 20 continue C C Solve on the coarsest level C l = lc if ( n3 .gt. 0 ) then do 30 isweep = 1, n3 call mgrlax( l, q, rn, w ) if ( lsweep ) write (*,2100) l, rn, w 30 continue else call mgsolc( l, q, rn, w ) if ( lsweep ) write (*,2100) l, rn, w end if C C Upward branch: coarse to fine C do 50 l = lc+inc, lf, inc call mgctof( l-inc, l, q, w ) do 40 isweep = 1, n2 call mgrlax( l, q, rn, w ) if ( l .eq. lf ) then if ( first ) then rn0 = rn ns0 = -1 w0 = w first = .false. rn1 = rn0 ns1 = ns0 w1 = w0 end if rnfine = rn ns0 = ns0 + 1 ns1 = ns1 + 1 end if if ( lsweep ) write (*,2100) l, rn, w 40 continue 50 continue C C Print convergence factors for this cycle C if ( lcycle ) then if ( rn1 .gt. 0.0 ) then cfact = rnfine/rn1 sfact = 0.0 if ( ns1 .gt. 0 ) sfact = cfact**(1.0/ns1) wu = w - w1 wfact = 0.0 if ( wu .gt. 0.0 ) wfact = cfact**(1.0/wu) write (*,2200) icycle, cfact, sfact, wfact end if rn1 = rnfine ns1 = 0 w1 = w end if C C Test for convergence C rn = rnfine if ( tol .ne. 0.0 .and. rn .le. eps ) go to 200 100 continue C C Finished: print overall convergence factors and quit C 200 if ( lovall .and. rn0 .gt. 0.0 ) then cfact = rnfine/rn0 sfact = 0.0 if ( ns0 .gt. 0 ) sfact = cfact**(1.0/ns0) wu = w - w0 wfact = 0.0 if ( wu .gt. 0.0 ) wfact = cfact**(1.0/wu) write (*,2300) cfact, sfact, wfact end if if ( tol .eq. 0.0 ) return if ( rn .le. eps ) return write (*,*) ' mgcv error: convergence tolerance not achieved' write (*,*) ' Tolerance: tol = ', tol write (*,*) ' Res. norm: rn = ', rn write (*,*) ' Cycles : nc = ', nc return C C _format_ statements C 2000 format(///' MultiGrid V-Cycle trace', & //' Levels:',i3,' is finest and',i3,' is coarsest', & /' Sweeps: n1 =',i3,' n2 =',i3,' n3 =',i3, & /' Stop : tol =',1pg11.3,' nc=',i4) 2100 format(' Level',i3,': Residual Norm =',1pe10.3,' Work =',0pf7.3) 2200 format(' Cycle',i3,': Convergence =',1pe10.3,'/V-cycle', & 0pf7.3,'/sweep',f7.3,'/work unit') 2300 format(/' overall : convergence =',1pe10.3,' ', + 0pf7.3,'/sweep',f7.3,'/work unit') 2400 format(11x,' new tolerance =',1pe10.3) C------------------ end of mgcv ------------------------------------- end