path: root/src/mg59p.F

#include "cctk.h"
c----------------------------------------------------------------------
      subroutine mgparm(m,ifd59,id5,id9,idi,idg,imx,jmx)
      implicit CCTK_REAL (a-h,o-z)
c----------------------------------------------------------------------
c  Given imx, jmx and ifd59 (See comments in mgsu2), mgparm calculates
c  the number of grids that will be needed, and the dimensions
*  needed for the coefficient, right hand side and solution arrays
*  to store values on all grid levels.
c .....................................................................
*
****  Parameters
*
*  INTEGERS:
*     m :
*         This is the number of grid levels that the multigrid
*         routine will use.
*
*     ifmg :
*          ifmg = 0 - The full multigrid algorithm is not used to
*                     obtain a good initial guess on the fine grid.
*                     (use this if you can provide a good initial guess)
*          ifmg = 1 - The full multigrid algorithm is used to obtain a
*                     good initial guess on the fine grid.
*
*     id5 :
*          Dimension of the arrays ac,aw,as,ae,an,q and f. id5 is the
*          total number of grid points on the finest grid and all
*          coarser grids.
*
*     id9 :
*          Dimension of the arrays asw,ase,ane,anw. If ifd59=5 then
*          id9=idi.  If ifd59=9 then id9=id5.
*          (NOTE: This routine specifically written for 9-point)
*
*     idi :
*          Dimension of the work arrays pu and pd. idi is the total
*          number of grid points on all of the coarser grids.
*
*     idg :
*          Dimension of the work array gam. It is set to the value im,
*          the number of grid points in the i-direction on the finest
*          grid.
*
*     imx,jmx :
*          These are the number of points in the i and j directions
*          including the two ficticious points.
*
      parameter(n5=1,n9=2,ni=3,jm=4,i9=5,j9=6,ifd=7,jred=8)
      dimension np2(20,8)
      iq5=1
      iq9=1
      iqi=1
      m=1
      np2(m,1)=jmx
      np2(m,2)=3
   10 if(np2(m,1).le.3) go to 20
      m=m+1
      np2(m,1)=np2(m-1,1)/2+1
      if(np2(m-1,2).eq.2.and.mod(np2(m-1,1),2).eq.1)
     + np2(m,1)=np2(m,1)+1
      np2(m,2)=2
      go to 10
   20 do 30 k=1,m
      np2(m-k+1,jm)=np2(k,1)
   30 np2(m-k+1,jred)=np2(k,2)
      do 40 k=m,1,-1
      ktot=imx*np2(k,jm)
      np2(k,n5)=iq5
      iq5=iq5+ktot
      np2(k,n9)=iq9
      if(k.lt.m.or.ifd59.eq.9) iq9=iq9+ktot
      np2(k,ni)=iqi
   40 if(k.lt.m) iqi=iqi+ktot
      do 50 k=1,m
      np2(k,i9)=imx
      np2(k,j9)=np2(k,jm)
   50 np2(k,ifd)=9
         if(ifd59.eq.5) then
      np2(m,i9)=1
      np2(m,j9)=1
      np2(m,ifd)=5
         endif
      id5=iq5-1
      id9=iq9-1
      idi=iqi-1
      idg=imx
      return
      end
*
*
*
************************************************************************
*
      subroutine mg9 (idim,ilower,iupper,jdim,jlower,jupper,
     &             cc,cn,cs,cw,ce,cnw,cne,csw,cse,u,rhs,
     &             id5,id9,idi,idg,ifmg,eps,rmax,ier)
      implicit CCTK_REAL (a-h,o-z)
*
************************************************************************
*
*
*     This routine is a wrapper for the multigrid solver.  The input
*  arrays are the finite difference coefficients on the 2D grid with
*  the boundary conditions absorbed into them.
*
****  Parameters
*
*  INTEGERS:
*     idim,jdim :
*          These define sizes of the coefficient arrays passed to this
*          routine.
*
*     ilower,iupper,
*     jlower,jupper :
*          These are the indices of the computational grid that
*          correspond to upper and lower index limits for points on
*          which computations are actually done.
*
*     id5 :
*          Dimension of the arrays ac,aw,as,ae,an,q and f. id5 is the
*          total number of grid points on the finest grid and all
*          coarser grids.
*
*     id9 :
*          Dimension of the arrays asw,ase,ane,anw. If ifd59=5 then
*          id9=idi.  If ifd59=9 then id9=id5.
*          (NOTE: This routine specifically written for 9-point)
*
*     idi :
*          Dimension of the work arrays pu and pd. idi is the total
*          number of grid points on all of the coarser grids.
*
*     idg :
*          Dimension of the work array gam. It is set to the value im,
*          the number of grid points in the i-direction on the finest
*          grid.
*
*     ifmg :
*          ifmg = 0 - The full multigrid algorithm is not used to
*                     obtain a good initial guess on the fine grid.
*                     (use this if you can provide a good initial guess)
*          ifmg = 1 - The full multigrid algorithm is used to obtain a
*                     good initial guess on the fine grid.
*
*     ier :
*          This is an error flag with the following possible return
*          values:
*              ier =  0 => solver converged without error
*              ier = -1 => solver did not converge
*
*  REALS:
*     cc(idim,jdim),cn(idim,jdim),
*     cs(idim,jdim),cw(idim,jdim),ce(idim,jdim),
*     cnw(idim,jdim),cne(idim,jdim),
*     csw(idim,jdim),cse(idim,jdim) :
*          These are the finite difference coefficient arrays for a
*          nine-point stencil on the two dimensional grid as follows:
*
*
*               cnw    cn    cne
*            ^
*            |  cw     cc    ce
* increasing |
*  j values  |  csw    cs    cse
*   (theta)  |
*             --------> increasing i values (eta)
*
*  Ie. the coresspondence is : nw : i-1,j+1
*                              n  : i,j+1
*                              ne : i+1,j+1
*                              w  : i-1,j
*                              c  : i,j
*                              e  : i+1,j
*                              sw : i-1,j-1
*                              s  : i,j-1
*                              se : i+1,j-1
*
*     u :
*          Input: this contains the initial guess to the solution of the
*                 equation
*          Output: This contains the final approximation to the solution
*                 determined by the multigrid solver.
*     
*     rhs :
*          This array contains the values of the right hand side of the
*          equation at every point on the rwo dimensional grid.
*
*     eps :
*          eps > 0.  The maximum norm of the residual is calculated at
*                    the end of each multigrid cycle. The algorithm is
*                    terminated when this maximum becomes less than eps
*                    or when the maximum number of iterations (see ncyc)
*                    is exceeded.  It is up to the user to provide a
*                    meaningfull tolerance criteria for the particular
*                    problem being solved.
*
*     rmax:
*          This is the final value of the residual calcu;ated.
*
************************************************************************
*
*

      integer idim,ilower,iupper,jdim,jlower,jupper,ifmg
      integer id5,id9,idi,idg,ier
      CCTK_REAL cc(idim,jdim),cn(idim,jdim),cs(idim,jdim),cw(idim,jdim),
     &     ce(idim,jdim),cnw(idim,jdim),cne(idim,jdim),csw(idim,jdim),
     &     cse(idim,jdim),u(idim,jdim),rhs(idim,jdim)
      CCTK_REAL eps,rmax

*
************************************************************************
*
*  Variable definitions:
*
*  Integers:
*     ifd59 :
*          ifd59 = 5 - means a 5-point finite difference stencil
*                      (ac,an,as,aw,ae) is defined on the finest grid.
*          ifd59 = 9 - means a 9-point finite difference stencil
*                      (ac,an,as,aw,ae,anw,ane,asw,ase) is defined on
*                      the finest grid by the user.
*                      (NOTE: This routine specifically written
*                             for 9-point)
*
*     ncyc :
*          The maximum number of multigrid "v"-cycles to be used. If
*          the maximum norm of the residual is not less than tol at the
*          end of ncyc cycles, the algorithm is terminated.
*          (NOTE: ncyc <= 40 )
*          
*     np(20,8) :
*          Input: When the iskip=1,-1 or -2 option is used, np2 is
*                 assumed to contain the grid information for umgs2.
*          Output: When the iskip=0 option is used, the grid
*                 information for umgs2 is returned in np2.
*          (NOTE: This is only useful for multiple instance problems)
*
*     iskip :
*          iskip = 0 - The coarse grid information, coarse grid
*                      operators and interpolation coefficients are
*                      calculated by umgs2.  This information is stored
*                      in the arrays ac, aw, as, asw, ase, pu, pd, np2
*                      and the variable.
*          iskip = 1 - The calculation of the coarse grid information,
*                      coarse grid operators and interpolation
*                      coefficients is skipped.  This option would be
*                      used when umgs2 has been called with iskip=0 and
*                      is being called again to solve a system of
*                      equations with the same matrix. This would be
*                      the case in, say, parabolic problems with time
*                      independent coefficients.
*          iskip =-1 - The set up of pointers (ugrdfn) is skipped.
*                      Coarse grid operators and interpolation
*                      coefficients are calculated and the given matrix
*                      equation is solved.  This option would be used
*                      when umgs2 has been called with iskip=0 and is
*                      being called again to solve a system of equations
*                      with a different matrix of the same dimensions.
*                      This would be the case for, say, parabolic
*                      problems with time dependent coefficients.
*          iskip =-2 - The set up of pointers (ugrdfn) is skipped.
*                      Coarse grid operators and interpolation
*                      coefficients are calculated and returned.
*                      No matrix solve.
*
*     ipc :
*          ipc = 0 or 1.
*          ipc is a multigrid parameter which determines the type of
*          interpolation to be used.  Usually ipc=1 is best.  However,
*          if the boundary condition equations have been absorbed into
*          the interior equations then ipc=0 can be used which results
*          in a slightly more efficient algorithm.
*
*     nman :
*          nman = 0 usually.
*          nman = 1 signals that the fine grid equations are singulari
*          for the case when homogeneous Neumann boundary conditions are
*          applied along the entire boundary.  In this case, the
*          difference equations are singular and the condition that the
*          integral of q over the domain be zero is added to the set of
*          difference equations.  This condition is satisfied by adding
*          the appropriate constant vector to q on the fine grid.  It is
*          assumed, in this case, that a well-defined problem has been
*          given to mgss2, i.e. the integral of f over the domain is
*          zero.
*
*     im :
*          The number of grid points in the x-direction (including two
*          ficticious points)
*     jm :
*          The number of grid points in the y-direction (including two
*          ficticious points)
*
*     linp :
*          This is a dummy argument left over from the authors
*          development of the code
*             Use:  common /io/ linp,lout
*
*     lout :
*          lout = unit number of output file into which the maximum norm
*          of the residual after each multigrid v-cycle" is printed.
*             Use:  common /io/ linp,lout
*
*     iscale :
*          Flag to indicate whether problem can be diagonally scaled to
*          speed convergence of the multigrid solver.
*
*     REALS:
*
*     ac(id5),an(id5),as(id5),aw(id5),ae(id5),
*     anw(id9),ane(id9),asw(id9),ase(id9) :
*          Input: ac, an, as, aw, ae, anw, ane, asw and ase contain the
*                 stencil coefficients for the difference operator on
*                 the finest grid. When the iskip=1 option is used,
*                 these arrays also are assumed to contain the coarse
*                 grid difference stencil coeficients.
*          Output: when the iskip=0 option is used, the coarse grid
*                 stencil coeficients are returned in ac, an, as, aw,
*                 ae, anw, ane, asw and ase.
*
*     ru(idi),rd(idi),rc(idi) :
*          Real work arrays.
*
*     pu(idi),pd(idi),pc(idi) :
*          Real work arrays.
*          Input: when the iskip=1 option is used, these arrays are
*                 assumed to contain the interpolation coefficients used
*                 in the semi-coarsening multigrid algorithm.
*          Output: when the iskip=0 option is used, the interpolation
*                 coeficients are returned in pu and pd.
*
*     f(id5) :
*          f contains the right hand side vector of the matrix
*          equation to be solved by umgs2.
*
*     q(id5) :
*          If ifmg=0, q contains the initial guess on the fine grid.
*          If ifmg=1, the initial guess on the fine grid is determined
*                     by the full multigrid process and the value of
*                     q on input to umgs2 not used.
*
*     tol :
*          tol > 0.  The maximum norm of the residual is calculated at
*                    the end of each multigrid cycle. The algorithm is
*                    terminated when this maximum becomes less than tol
*                    or when the maximum number of iterations (see ncyc)
*                    is exceeded.  It is up to the user to provide a
*                    meaningfull tolerance criteria for the particular
*                    problem being solved.
*          tol = 0.  Perform ncyc multigrid cycles.  Calculate and print
*                    the maximum norm of the residual after each cycle.
*          tol =-1.  Perform ncyc multigrid cycles.  The maximum norm of
*                    the final residual is calculated and returned in
*                    the variable rmax in the calling list of umgs2.
*          tol =-2.  Perform ncyc multigrid cycles.  The maximum norm of
*                    the residual is not calculated.
*
*     rmax :
*          If tol.ge.-1., the final residual norm is returned in rmax.
*
************************************************************************
*
*

      integer ncyc,ifd59
      parameter (ncyc=40,ifd59=9)

      integer np2(20,8)
      integer iskip,ipc,nman
      parameter (iskip=0,ipc=1,nman=0)
      integer irc,irurd,im,jm
      integer linp,lout
      common /io/ linp,lout

      CCTK_REAL ac(id5),an(id5),as(id5),aw(id5),ae(id5),
     &     anw(id9),ane(id9),asw(id9),ase(id9),
     &     q(id5),f(id5),gam(idg)

      CCTK_REAL ru(idi),rd(idi),rc(idi),pu(idi),pd(idi),pc(idi)

      CCTK_REAL tol

      integer iscale,itry


      ier=0

*
* Set some parameters for multigrid solver
*
      lout=6
c      rewind(unit=lout)
      irc=0
      irurd=0
      im=iupper-ilower+3
      jm=jupper-jlower+3
      tol=eps

*
*  Set up coefficients into vectors with correct indexing
*
      do 110 j=jlower,jupper
      do 100 i=ilower,iupper
      n=(j+1-jlower)*im + i+1-(ilower-1)
      ac(n)=cc(i,j)
      an(n)=cn(i,j)
      as(n)=cs(i,j)
      aw(n)=cw(i,j)
      ae(n)=ce(i,j)
      anw(n)=cnw(i,j)
      ane(n)=cne(i,j)
      asw(n)=csw(i,j)
      ase(n)=cse(i,j)
      q(n)=u(i,j)
      f(n)=rhs(i,j)
  100 continue
  110 continue


*
*  Determine whether we can diagonal scale the problem to speed
*  convergence. Can only be done if there are no zeros on the main
*  diagonal (ie. central difference coefficient).
*
      iscale=1
      do 200 j=jlower,jupper
      do 205 i=ilower,iupper
      n=(j+1-jlower)*im + i+1-(ilower-1)
      if (ac(n) .eq. 0.) then
        iscale=0
      endif
  205 continue
  200 continue

*
*  Do the diagonal scaling if we can.
*
      if (iscale.eq.1) then

      do 210 j=jlower,jupper
      do 215 i=ilower,iupper
        n=(j+1-jlower)*im + i+1-(ilower-1)
        f(n)=f(n)/ac(n)
        ase(n)=ase(n)/ac(n)
        asw(n)=asw(n)/ac(n)
        ane(n)=ane(n)/ac(n)
        anw(n)=anw(n)/ac(n)
        ae(n)=ae(n)/ac(n)
        aw(n)=aw(n)/ac(n)
        as(n)=as(n)/ac(n)
        an(n)=an(n)/ac(n)
        ac(n)=1.
  215 continue
  210 continue

      endif

c 
c   Now call the multigrid routine

      itry=0
      
 1122 call umgs2(
     . ac,aw,as,ae,an,asw,ase,ane,anw,q,f,pu,pd,pc,ru,rd,rc,gam,np2,
     . ifd59,ifmg,ncyc,tol,nman,im,jm,id5,id9,idi,m,iskip,rmax,
     . ipc,irc,irurd)

      if ((rmax.gt.tol).and.(itry.le.5)) then
        itry=itry+1
        print*,"Retry #",itry
        goto 1122
      endif

      if (rmax.gt.tol) then
        ier = -1
        print*,"Did not converge"
        print*," maximum residual    = ",rmax
        print*," tolerance           = ",tol
      endif

*
*  Convert the solution back to the 2D array form
*
      do 510 j=jlower,jupper
      do 500 i=ilower,iupper
      n=(j+1-jlower)*im + i+1-(ilower-1)
      u(i,j)=q(n)
  500 continue
  510 continue

      return
      end



*
*
*
************************************************************************
*
      subroutine mg5 (idim,ilower,iupper,jdim,jlower,jupper,
     &             cc,cn,cs,cw,ce,u,rhs,
     &             id5,id9,idi,idg,ifmg,eps,rmax,ier)
      implicit CCTK_REAL (a-h,o-z)
*
************************************************************************
*
*
*     This routine is a wrapper for the multigrid solver.  The input
*  arrays are the finite difference coefficients on the 2D grid with
*  the boundary conditions absorbed into them.
*
****  Parameters
*
*  INTEGERS:
*     idim,jdim :
*          These define sizes of the coefficient arrays passed to this
*          routine.
*
*     ilower,iupper,
*     jlower,jupper :
*          These are the indices of the computational grid that
*          correspond to upper and lower index limits for points on
*          which computations are actually done.
*
*     id5 :
*          Dimension of the arrays ac,aw,as,ae,an,q and f. id5 is the
*          total number of grid points on the finest grid and all
*          coarser grids.
*
*     id9 :
*          Dimension of the arrays asw,ase,ane,anw. If ifd59=5 then
*          id9=idi.  If ifd59=9 then id9=id5.
*          (NOTE: This routine specifically written for 5-point)
*
*     idi :
*          Dimension of the work arrays pu and pd. idi is the total
*          number of grid points on all of the coarser grids.
*
*     idg :
*          Dimension of the work array gam. It is set to the value im,
*          the number of grid points in the i-direction on the finest
*          grid.
*
*     ifmg :
*          ifmg = 0 - The full multigrid algorithm is not used to
*                     obtain a good initial guess on the fine grid.
*                     (use this if you can provide a good initial guess)
*          ifmg = 1 - The full multigrid algorithm is used to obtain a
*                     good initial guess on the fine grid.
*
*     ier :
*          This is an error flag with the following possible return
*          values:
*              ier =  0 => solver converged without error
*              ier = -1 => solver did not converge
*
*  REALS:
*     cc(idim,jdim),cn(idim,jdim),
*     cs(idim,jdim),cw(idim,jdim),ce(idim,jdim):
*          These are the finite difference coefficient arrays for a
*          five-point stencil on the two dimensional grid as follows:
*
*
*                      cn       
*            ^
*            |  cw     cc    ce
* increasing |
*  j values  |         cs       
*   (theta)  |
*             --------> increasing i values (eta)
*
*  Ie. the coresspondence is : n  : i,j+1
*                              ne : i+1,j+1
*                              c  : i,j
*                              e  : i+1,j
*                              s  : i,j-1
*
*     u :
*          Input: this contains the initial guess to the solution of the
*                 equation
*          Output: This contains the final approximation to the solution
*                 determined by the multigrid solver.
*     
*     rhs :
*          This array contains the values of the right hand side of the
*          equation at every point on the rwo dimensional grid.
*
*     eps :
*          eps > 0.  The maximum norm of the residual is calculated at
*                    the end of each multigrid cycle. The algorithm is
*                    terminated when this maximum becomes less than tol
*                    or when the maximum number of iterations (see ncyc)
*                    is exceeded.  It is up to the user to provide a
*                    meaningfull tolerance criteria for the particular
*                    problem being solved.
*     rmax:
*          This is the final value of the residual calcu;ated.
*
************************************************************************
*
*

      integer idim,ilower,iupper,jdim,jlower,jupper,ifmg
      integer id5,id9,idi,idg,ier
      CCTK_REAL cc(idim,jdim),cn(idim,jdim),cs(idim,jdim),cw(idim,jdim),
     &     ce(idim,jdim),u(idim,jdim),rhs(idim,jdim)
      CCTK_REAL eps,rmax

*
************************************************************************
*
*  Variable definitions:
*
*  Integers:
*     ifd59 :
*          ifd59 = 5 - means a 5-point finite difference stencil
*                      (ac,an,as,aw,ae) is defined on the finest grid.
*          ifd59 = 9 - means a 9-point finite difference stencil
*                      (ac,an,as,aw,ae,anw,ane,asw,ase) is defined on
*                      the finest grid by the user.
*                      (NOTE: This routine specifically written
*                             for 5-point)
*
*     ncyc :
*          The maximum number of multigrid "v"-cycles to be used. If
*          the maximum norm of the residual is not less than tol at the
*          end of ncyc cycles, the algorithm is terminated.
*          (NOTE: ncyc <= 40 )
*          
*     np(20,8) :
*          Input: When the iskip=1,-1 or -2 option is used, np2 is
*                 assumed to contain the grid information for umgs2.
*          Output: When the iskip=0 option is used, the grid
*                 information for umgs2 is returned in np2.
*          (NOTE: This is only useful for multiple instance problems)
*
*     iskip :
*          iskip = 0 - The coarse grid information, coarse grid
*                      operators and interpolation coefficients are
*                      calculated by umgs2.  This information is stored
*                      in the arrays ac, aw, as, asw, ase, pu, pd, np2
*                      and the variable.
*          iskip = 1 - The calculation of the coarse grid information,
*                      coarse grid operators and interpolation
*                      coefficients is skipped.  This option would be
*                      used when umgs2 has been called with iskip=0 and
*                      is being called again to solve a system of
*                      equations with the same matrix. This would be
*                      the case in, say, parabolic problems with time
*                      independent coefficients.
*          iskip =-1 - The set up of pointers (ugrdfn) is skipped.
*                      Coarse grid operators and interpolation
*                      coefficients are calculated and the given matrix
*                      equation is solved.  This option would be used
*                      when umgs2 has been called with iskip=0 and is
*                      being called again to solve a system of equations
*                      with a different matrix of the same dimensions.
*                      This would be the case for, say, parabolic
*                      problems with time dependent coefficients.
*          iskip =-2 - The set up of pointers (ugrdfn) is skipped.
*                      Coarse grid operators and interpolation
*                      coefficients are calculated and returned.
*                      No matrix solve.
*
*     ipc :
*          ipc = 0 or 1.
*          ipc is a multigrid parameter which determines the type of
*          interpolation to be used.  Usually ipc=1 is best.  However,
*          if the boundary condition equations have been absorbed into
*          the interior equations then ipc=0 can be used which results
*          in a slightly more efficient algorithm.
*
*     nman :
*          nman = 0 usually.
*          nman =1 signals that the fine grid equations are singular for
*          the case when homogeneous Neumann boundary conditions are
*          applied along the entire boundary.  In this case, the
*          difference equations are singular and the condition that the
*          integral of q over the domain be zero is added to the set of
*          difference equations.  This condition is satisfied by adding
*          the appropriate constant vector to q on the fine grid.  It is
*          assumed, in this case, that a well-defined problem has been
*          given to mgss2, i.e. the integral of f over the domain is
*          zero.
*
*     im :
*          The number of grid points in the x-direction (including two
*          ficticious points)
*     jm :
*          The number of grid points in the y-direction (including two
*          ficticious points)
*
*     linp :
*          This is a dummy argument left over from the authors
*          development of the code
*             Use:  common /io/ linp,lout
*
*     lout :
*          lout = unit number of output file into which the maximum norm
*          of the residual after each multigrid v-cycle" is printed.
*             Use:  common /io/ linp,lout
*
*     iscale :
*          Flag to indicate whether problem can be diagonally scaled to
*          speed convergence of the multigrid solver.
*
*     REALS:
*
*     ac(id5),an(id5),as(id5),aw(id5),ae(id5),
*     anw(id9),ane(id9),asw(id9),ase(id9) :
*          Input: ac, an, as, aw, ae, anw, ane, asw and ase contain the
*                 stencil coefficients for the difference operator on
*                 the finest grid. When the iskip=1 option is used,
*                 these arrays also are assumed to contain the coarse
*                 grid difference stencil coeficients.
*          Output: when the iskip=0 option is used, the coarse grid
*                 stencil coeficients are returned in ac, an, as, aw,
*                 ae, anw, ane, asw and ase.
*
*     ru(idi),rd(idi),rc(idi) :
*          Real work arrays.
*
*     pu(idi),pd(idi),pc(idi) :
*          Real work arrays.
*          Input: when the iskip=1 option is used, these arrays are
*                 assumed to contain the interpolation coefficients used
*                 in the semi-coarsening multigrid algorithm.
*          Output: when the iskip=0 option is used, the interpolation
*                 coeficients are returned in pu and pd.
*
*     f(id5) :
*          f contains the right hand side vector of the matrix
*          equation to be solved by umgs2.
*
*     q(id5) :
*          If ifmg=0, q contains the initial guess on the fine grid.
*          If ifmg=1, the initial guess on the fine grid is determined
*                     by the full multigrid process and the value of
*                     q on input to umgs2 not used.
*
*     tol :
*          tol > 0.  The maximum norm of the residual is calculated at
*                    the end of each multigrid cycle. The algorithm is
*                    terminated when this maximum becomes less than tol
*                    or when the maximum number of iterations (see ncyc)
*                    is exceeded.  It is up to the user to provide a
*                    meaningfull tolerance criteria for the particular
*                    problem being solved.
*          tol = 0.  Perform ncyc multigrid cycles.  Calculate and print
*                    the maximum norm of the residual after each cycle.
*          tol =-1.  Perform ncyc multigrid cycles.  The maximum norm of
*                    the final residual is calculated and returned in
*                    the variable rmax in the calling list of umgs2.
*          tol =-2.  Perform ncyc multigrid cycles.  The maximum norm of
*                    the residual is not calculated.
*
************************************************************************
*
*

      integer ncyc,ifd59
      parameter (ncyc=40,ifd59=5)

*     integer id5,id9,idi,idg
*  This is for a 103x28 grid
*     parameter(id5=6695,id9=3811,idi=3811,idg=103)
*  This is for a 203x56 grid
*     parameter(id5=24969,id9=13601,idi=13601,idg=203)
*  This is for a 403x118 grid
*     parameter(id5=100750,id9=52793,idi=52793,idg=403)

      integer np2(20,8)
      integer iskip,ipc,nman
      parameter (iskip=0,ipc=1,nman=0)
      integer irc,irurd,im,jm
      integer linp,lout
      common /io/ linp,lout

      CCTK_REAL ac(id5),an(id5),as(id5),aw(id5),ae(id5),
     &     anw(id9),ane(id9),asw(id9),ase(id9),
     &     q(id5),f(id5),gam(idg)

      CCTK_REAL ru(idi),rd(idi),rc(idi),pu(idi),pd(idi),pc(idi)

      CCTK_REAL tol

      integer iscale,itry



      ier=0
*
* Set some parameters for multigrid solver
*
      lout=6
c      rewind(unit=lout)
      irc=0
      irurd=0
      im=iupper-ilower+3
      jm=jupper-jlower+3
      tol=eps


*
*  Set up coefficients into vectors with correct indexing
*
      do 110 j=jlower,jupper
      do 100 i=ilower,iupper
      n=(j+1-jlower)*im + i+1-(ilower-1)
      ac(n)=cc(i,j)
      an(n)=cn(i,j)
      as(n)=cs(i,j)
      aw(n)=cw(i,j)
      ae(n)=ce(i,j)
      anw(n)=0.
      ane(n)=0.
      asw(n)=0.
      ase(n)=0.
      q(n)=u(i,j)
      f(n)=rhs(i,j)
  100 continue
  110 continue


*
*  Determine whether we can diagonal scale the problem to speed
*  convergence. Can only be done if there are no zeros on the main
*  diagonal (ie. central difference coefficient).
*
      iscale=1
      do 200 j=jlower,jupper
      do 205 i=ilower,iupper
      n=(j+1-jlower)*im + i+1-(ilower-1)
      if (ac(n) .eq. 0.) then
        iscale=0
      endif
  205 continue
  200 continue

*
*  Do the diagonal scaling if we can.
*
      if (iscale.eq.1) then

      do 210 j=jlower,jupper
      do 215 i=ilower,iupper
        n=(j+1-jlower)*im + i+1-(ilower-1)
        f(n)=f(n)/ac(n)
        ae(n)=ae(n)/ac(n)
        aw(n)=aw(n)/ac(n)
        as(n)=as(n)/ac(n)
        an(n)=an(n)/ac(n)
        ac(n)=1.
  215 continue
  210 continue

      endif

c 
c   Now call the multigrid routine

      itry=0
      
 1122 call umgs2(
     + ac,aw,as,ae,an,asw,ase,ane,anw,q,f,pu,pd,pc,ru,rd,rc,gam,np2,
     + ifd59,ifmg,ncyc,tol,nman,im,jm,id5,id9,idi,m,iskip,rmax,
     + ipc,irc,irurd)


      if ((rmax.gt.tol).and.(itry.le.5)) then
        itry=itry+1
        print*,"Retry #",itry
        goto 1122
      endif

      if (rmax.gt.tol) then
        ier = -1
        print*,"Did not converge"
        print*," maximum residual    = ",rmax
        print*," tolerance           = ",tol
      endif

*
*  Convert the solution back to the 2D array form
*
      do 510 j=jlower,jupper
      do 500 i=ilower,iupper
      n=(j+1-jlower)*im + i+1-(ilower-1)
      u(i,j)=q(n)
  500 continue
  510 continue

      return
      end