      program mglab
c ********************** MULTIGRID LABORATORY *****************************
c *                                                                       *
c *            Tutorial software to accompany the book:                   *
c *                                                                       *
c *              AN INTRODUCTION TO MULTIGRID METHODS                     *
c *                      by P.Wesseling                                   *
c *                  Wiley, Chichester, 1992                              *
c *************************************************************************
c *                       Main program                                    *
c *************************************************************************
c
c  Programmer    P.Wesseling
c  Version 1.2   date 18-04-92
c  Version 1.1   date 16-12-91  (F.N. van de Vosse: first revision)
c  Version 1.0   date 10-12-91  (developed at HP Vectra RS/20C
c
c  Copyright (c) 1992  P.Wesseling
c  Address:  Delft University of Technology
c            Department of Technical Mathematics and Informatics
c            Mekelweg 4,  2628 CD Delft,  The Netherlands
c            Telefax +3115 787209
c            email witawes@dutinfh.tudelft.nl
c  Permission to copy or distribute this software in hard copy or soft copy
c  is hereby granted by the programmer.
c
c **************************************************************************
c
c  Multigrid method for the solution of the following boundary value problem:
c    -(ay')' + (by)' + cy = s,   0 < x < 1,  a > 0
c  Nonlinear multigrid algorithm as described in section 8.7 of the book.
c  The problem solved is only one-dimensional, and the purpose of this
c  software is purely tutorial.
c
c **************************************************************************
c
c                  INPUT / OUTPUT PARAMETERS
c
      implicit none
      integer  histry(2,1000), ib(13), icase,  igal, ihis,
     $         iout(4), ism, it, levels, nc, nmg, np(11), ns(2)
      double precision aa,bb, cc, delta, f0, f1, om, spar, ss, tol
      integer idim
      parameter (idim = 1000)
      double precision amat(idim,3),dmat(idim,3),r(idim),u(idim),
     $                 ue(idim), wu(10), rhs(idim), s(10), ut(idim),
     $                 v(idim), uce(idim)
      character cycle
c **************************************************************************
c
c  aa        i  governs the diffusion coefficient in the diff. equation
c               see function a
c  amat      o  matrix (discretization of differential equation)
c                 igal=0: matrices on all grids created by matrhs
c                 igal=1: matrix on finest grid created by subroutine matrhs
c                         and initialized to 0.d0 on the other grids, where the
c                         matrices will be made by subroutine galerk
c  bb        i  governs the velocity coefficient in the diff. equation
c               see function b
c  cc        i  governs the third coefficient in the diff. equation
c               see function c
c  cycle     i   governs type of multigrid schedule
c                 cycle = 'A' : adaptive schedule
c                 cycle = 'F' : F-cycle
c                 cycle = 'V' : V-cycle
c                 cycle = 'W' : W-cycle
c  delta     i  used in setting tolerance eps in adaptive cycle according
c                to equation (8.3.1) in the book
c  dmat      o  matrix used when defect correction is applied: more accurate
c                discretization than amat
c  f0        i  value for dirichlet or neumann boundary condition at x = 0
c  f1        i  value for dirichlet or neumann boundary condition at x = 1
c  histry    o  gives sequence of grids visited and number of smoothings
c                 carried out during each visit of a grid
c                 histry(1,i) : number of grid visited
c                 histry(2,i) : number of smoothings carried out at
c                   visit to grid(histry(1,i))
c  ib        i  ib(1) = 1: Dirichlet at x = 0
c               ib(1) = 2: Neumann   at x = 0
c               ib(2) = 1: Dirichlet at x = 1
c               ib(2) = 2: Neumann   at x = 1
c               ib(3) = 1: central
c               ib(3) = 2: upwind
c               ib(3) = 3: hybrid with diffusion term deleted when switched to
c                          upwind
c               ib(3) = 4: hybrid with diffusion term maintained when switched
c                          to upwind
c               ib(4) = 1: vertex-centered
c               ib(4) = 2: cell-centered
c               ib(5) governs updating of ut in nonlinear multigrid method
c                 implemented in subroutine mg
c               ib(6) chooses type of prolongation - see subroutine prol:
c                 ic = ib(6)
c               ib(7) chooses type of restriction for residual - see
c                 subroutine restri: ic = ib(7)
c               ib(8) chooses type of restriction for ut - see
c                 subroutine restri: ic = ib(8)
c               ib(9) = 1: direct solver is used instead of multigrid
c               ib(10) = n>0: defect correction is applied n times
c               ib(11) = n>0: nested iteration is applied with gamma~=n
c                 multigrid iterations (see program 1 in section 8.4)
c               ib(12) chooses type of prolongation in nested iteration
c               ib(13) chooses type of restriction for coarse right-hand-
c                 sides in nested iteration
c  icase     i  identification number given by user to his problem
c                icase = 1: interface problem: eq. (3.3.1) with
c                           discontinuous diffusion coefficient a(x)
c                           defined in function a
c                icase = 2: convection diffusion eq. (6.6.1) (1-dimensional)
c  idim      *  length of arrays reserved for matrix and grid vectors on
c                all grids together
c  igal      i  igal = 0: coarse grid discretization; matrices are made by
c                       discretization on all grids
c               igal = 1: Galerkin coarse grid discretization; only finest
c                       grid matrix is made and coarse grid matrices are
c                       initialized to 0.d0 before call of Galerkin
c  ihis      o   total number of entries in histry
c  iout      i   governs amount of output; iout larger -> more output
c                 iout(1) governs output in main program  mglab
c                 iout(2) and iout(3) govern output in subroutine mg
c                 iout(4) governs output in subroutines matrhs and galerk
c  ism       i   chooses smoothing method; see subroutine smoothing
c  it        o   gives current multigrid iteration number
c  levels    i   number of grids in multigrid method; < 10
c  nc        i   number of grid points on coarsest grid
c  nmg       i   maximum number of multigrid iterations
c  np        o   array of pointers to start of grids and to 1+end of finest grid
c                 np(1)=1     : start of coarsest grid
c                 np(levels)  : start of finest grid
c                 np(levels+1): end of finest grid + 1
c  ns        i   ns(1):  number of pre-smoothing iteraties
c                ns(2):  number of post-smoothing iterations
c  om        i  underrelaxation parameter in smoothing method
c  rhs       o  right-hand-side of discretization on finest grid made by
c                subroutine matrhs
c  s        i   underrelaxation parameters in the nonlinear multigrid
c                 algorithm. See algorithm TG in section 8.2 of the book.
c  spar      i  used in generating underrelaxation parameters s(k) in
c               the nonlinear multigrid algorithm. See algorithm TG in
c               section 8.2 of the book.
c  ss        i  governs the right-hand-side in the diff. equation
c               see function s
c  tol       i  accuracy requirement:  |Au-rhs)| < tol*|rhs| on finest grid.
c                Governs also termination ofcoarse grid correction in
c                adaptive cycle
c  u         o  new approximation of solution
c  uce       o  exact solution of differential equation
c  ue        o  exact discrete solution
c  ut        i  u~ in nonlinear multigrid algorithm.
c                 On coarse grids input is optional.
c                 It may or may not be updated during multigrid iteration
c                 depending on ib(5) = 1 : ut is updated else not.
c                 For linear problems the treatment of ut makes no difference.
c                 If ut is input on coarse grids it may have been obtained
c                 from nested iteration
c  v         o  initial approximation of solution
c  wu        o  wu(k) = work units expended on grid k
c
c ***************************************************************************
c
c                 LOCAL QUANTITIES
c
      integer i1, i2, ibb
      double precision anorm,bdef(idim),errmax,errl2,rdef(idim),ub(idim)
c
c  anorm     name of function that computes norms
c  bdef      scratch array for defect correction
c  errmax    scratch variable
c  err12     scratch variable
c  i1        input device number
c  i2        output device number
c  ibb       scratch variable
c  rdef      right-hand-side for second order discretization (def. corr.)
c  ub        scratch array for defect correction
c
c ***************************************************************************
c
c                   I/O
c
c  input file:  mglab.dat  device number = i1   free format
c  output file: mglab.res  device number = i2   free format
c
c **************************************************************************
c
c            SUBROUTINES CALLED BY MAIN PROGRAM
c
c  axplsy      u = a*ut + u;  a:real,  u, ut: gridfunctions
c  copy        copies gridfunctions
c  defres      computes residual for defect correction
c  dirsol      computes solution by direct method
c  err         puts out error messages
c  galerk      galerkin coarse grid approximation
c  info        user-provided subroutine that prints information about
c                his case; body may be empty
c  inout       machine dependent subroutine that opens input and output
c                files
c  input       reads and writes input
c  matrhs      computes matrix and right-hand-side on finest grid and matrices
c              on coarse grid in case of discretization coarse grid
c              approximation
c  mg          multigrid algorithm
c  nested      nested iteration: program 1 of section 8.4
c  pointr      computes pointers to starting addresses of grids
c  smooth      smoothing
c  uexact      computes exact solution of differential equation
c  vecout      prints gridfunctions
c
c ***************************************************************************
c
c                   FUNCTIONS CALLED BY MAIN PROGRAM
c
c                    none
c
c ***************************************************************************
c
c                 COMPLETE LIST OF SUBROUTINES
c
c       axplsy    chrfun    copy      defres    defrp1    diagn     dirsol
c       err       errmg     galerk    gener     ibchck    info      inout
c       input     matout    matrhs    mg        nested    pointr    prol
c       resid     restri    scal      smooth    vecout
c
c ***************************************************************************
c
c                 COMPLETE LIST OF FUNCTIONS
c
c       a    anorm   b     c    omega    s      u
c
c ***************************************************************************
c
c                METHOD
c
c  Multigrid method for the solution of the following boundary value
c    problem:
c    -(ay')' + (by)' + cy = s,   0<x<1, a>0
c
c  Nonlinear multigrid algorithm as described in section 8.7 of the book.
c
c  The diffusion coefficient a(x) is supposed to be possibly discontinuous
c    at grid points. The functions b(x) and c(x) are supposed to be possibly
c    discontinuous in the middle between at grid points. The user may put
c    discontinuities elsewhere (and will often have to on some coarse grids
c    if igal = 0) at his own risk.
c
c  Finite volume discretization of 1-dimensional second order boundary
c    value problem.
c  Upwind, central or hybrid discretization of convection (first order)
c    term, depending on ib(3).
c  Cell-centered or vertex-centered discretization, depending on ib(4).
c
c*************************************************************************
c
c                      STRUCTURE DIAGRAM
c
c  Statements that are not represented in the structure diagram are bracketed
c    by c1**** and c2*****
c
c #===========================================================================#
c #  matrhs                                                                   #
c #___________________________________________________________________________#
c # \                                                                       / #
c # F \                           ib (9) eq 1                             / T #
c #_____\_______________________________________________________________/_____#
c # \                                   /|                                    #
c # F \         igal eq 1             / T| exact solution of differential eq. #
c #_____\___________________________/____|                                    #
c #                   |  galerk          | direct solution of discrete system #
c #___________________|__________________|                                    #
c #   ut = v = 0                         | smoothing on finest grid only      #                                   |
c #______________________________________|                                    #
c # \                                   /|                                    #
c # F \         ib (11) > 0           / T|                                    #
c #_____\___________________________/____|                                    #
c #                   | nested iteration |                                    #
c #___________________|__________________|                                    #
c # \                                   /|                                    #
c # F \         ib (10) eq 0          / T|                                    #
c #_____\___________________________/____|                                    #
c # central discret.  | direct solution  |                                    #
c # defect correction | multigrid        |                                    #
c #===========================================================================#
c
c******************************************************************************
      parameter( i1=11, i2=12 )
c1*********************************
      write(*,*)
     $ '                ************ MULTIGRID LABORATORY ************'
      write(*,*)
     $ '                *   Tutorial software to accompany the book  *'
      write(*,*)
     $ '                *    AN INTRODUCTION TO MULTIGRID METHODS    *'
      write(*,*)
     $ '                *              by P.Wesseling                *'
      write(*,*)
     $ '                *          Wiley, Chichester, 1991           *'
      write(*,*)
     $ '                **********************************************'
      write(*,*) ' '
      write(*,*) ' '
c
c open the in- and output files
c
      call inout( i1, i2 )
c
c read the input
c
      call input( i1, i2, ib, icase, igal, iout, ism, levels, nc, nmg,
     $            ns, aa, bb, cc, delta, f0, f1, om, spar, ss, tol,
     $            cycle)
c
c close the input file
c
      close( i1 )
c
c call user-provided subroutine that prints information about his case
c
      call info( icase, ib )
c
c set pointers to starting adresses of grids
c
      call pointr( ib(4), iout(1), levels, nc, np, i2 )
c
c check of workspace reserved for amat
c
      if ( np(levels+1) .gt. idim) then
         call err( i1, i2, 4 )
      endif
c
c compute matrix and right-hand-side on finest grid and matrices
c on coarse grid in case of discretization coarse grid approximation
c
c2*********************************
      call matrhs( ib, icase, idim, igal, iout(4), i2, levels, nc,
     $             np, aa, amat, bb, cc, f0, f1, rhs, ss)
c
c if ib(9)=1  a single grid solution method is applied else multigrid solution
c
      if ( ib(9) .ne. 1 ) then
c
c galerkin coarse grid approximation if igal=1
c
         if ( igal .eq. 1 ) then
            call galerk( ib, idim, iout(4), i2, levels, np, amat, 1 )
         endif
c
c initialization of solution and exact solution
c
         do 5 ibb = 1, np(levels+1)
           ut(ibb) = 0d0
           ue(ibb) = 0d0
    5    continue
         do 6 ibb = np(levels), np(levels+1) - 1
           v(ibb) = ut(ibb)
    6    continue
c1*********************************
         ihis = 0
         do 7 ibb = 1, 10
           wu(ibb) = 0.d0
    7    continue
c2*********************************
c
c nested iteration if ib(11)>0
c
         if ( ib(11) .gt. 0 ) then
            call nested ( amat, rhs, cycle, delta, histry, ib, idim,
     $                    ihis, iout, it, ism, i2, levels, nmg, np, ns,
     $                    om, r, s, tol, u, ue, ut, v, wu )
c1***************************************************************
c
c Output of results
c
            if (iout(1) .gt. 1) then
        write(i2,*)'Multigrid solution on finest grid after nested it.'
              call vecout( i2, idim, levels, np, u )
            endif
c2***************************************************************
         endif
         if ( ib(10) .eq. 0 ) then
c
c  calculate the exact discrete solution on the finest grid
c
            call dirsol( idim, i2, levels, np, amat, ue, rhs )
c1***************************************************************
            if (iout(1) .gt. 0) then
              write(i2,*) 'Exact discrete solution on finest grid'
              call vecout( i2, idim, levels, np, ue )
            endif
c2***************************************************************
            call mg( amat, rhs, cycle, delta, histry, ib, idim, ihis,
     $               iout, it, ism, i2, levels, nmg, np, ns, om, r, s,
     $               tol, u, ue, ut, v, wu )

         else
c
c defect correction on central discretization
c
            ibb = ib(3)
            ib(3) = 1
            call matrhs( ib, icase, idim, 1, iout(4), i2, levels, nc,
     $                   np, aa, dmat, bb, cc, f0, f1, rdef, ss )
            ib(3) = ibb
c
c  calculate the exact discrete solution on the finest grid
c
            call dirsol( idim, i2, levels, np, dmat, ue, rdef )
c1***************************************************************
            if (iout(1) .gt. 0) then
              write(i2,*) 'Exact discrete solution on finest grid'
              call vecout( i2, idim, levels, np, ue )
            endif
c2***************************************************************
c
c defect correction according to section 4.6 of the book
c one multigrid iteration is applied
c
            call mg( amat, rhs, cycle, delta, histry, ib, idim, ihis,
     $               iout, it, ism, i2, levels, 1, np, ns, om, r, s,
     $               tol, ub, ue, ut, v, wu )
            do 10 ibb = 1, ib(10)
               call defres( idim, iout(4), i2, levels, np, amat, dmat,
     $          rhs, bdef, ub)
               call copy( idim, levels, np, v, ub )
               call mg( amat, bdef, cycle, delta, histry, ib, idim,
     $                  ihis, iout, it, ism, i2, levels, 1, np, ns, om,
     $                  r, s, tol, ub, ue, ut, v, wu )
   10       continue
            call copy( idim, levels, np, u, ub )
         endif
c1***************************************************************
c
c Output of results
c
            if (iout(1) .gt. 1) then
              write(i2,*) 'Multigrid solution on finest grid'
              call vecout( i2, idim, levels, np, u )
            endif
            call axplsy( idim, levels, np, -1d0, ue, u )
            if (iout(1) .gt. 0) then
              write(i2,*) 'Local error on finest grid'
              call vecout( i2, idim, levels, np, u )
            endif
            errmax = anorm( 1, idim, u, levels, np )
            errl2 = anorm( 2, idim, u, levels, np )
            write(i2,'(a,2e13.5)')
     $       'Maximum and L2 norm of the error = ', errmax, errl2
            write(*,'(a,2e13.5)')
     $       'Maximum and L2 norm of the error = ', errmax, errl2
            if (iout(1) .gt. 1) then
              write(i2,*) 'Work unit distribution over the grids'
              do 11 ibb = 1,levels
                 write(i2,'(i3,e10.2)') ibb, wu(ibb)
  11          continue
            endif
c2***************************************************************
      else
c
c computation of the exact solution of the differential equation
c
         call uexact( i2, icase, idim, levels, ib,np,uce,f0,f1,aa,bb)
c1***************************************************************
         if (iout(1) .gt. 0) then
           write(i2,*) 'Exact continuous solution on finest grid'
           call vecout( i2, idim, levels, np, uce )
         endif
c2***************************************************************
c
c computation of exact discrete solution by a direct method
c
         call dirsol( idim, i2, levels, np, amat, ue, rhs )
c1***************************************************************
         if (iout(1) .gt. 0) then
           write(i2,*) 'Direct solver solution on finest grid'
           call vecout( i2, idim, levels, np, ue )
         endif
         call axplsy( idim, levels, np, -1d0, ue, uce )
         if (iout(1) .gt. 0) then
           write(i2,*) 'Local direct solver error on finest grid'
           call vecout( i2, idim, levels, np, uce )
         endif
         errmax = anorm( 1, idim, uce, levels, np )
         errl2 = anorm( 2, idim, uce, levels, np )
         write(i2,'(a,2e13.5)')
     $         'Maximum and L2 norm of the error = ',errmax, errl2
         write(*,'(a,2e13.5)')
     $         'Maximum and L2 norm of the error = ',errmax, errl2
c2***************************************************************
c
c test the smoother and the discretization error
c
         do 15 ibb = np(levels), np(levels+1) -1
           v(ibb) = 0d0
   15    continue
         call smooth( i2, idim, iout, ism, levels, np, ns(1), amat,
     $                levels, om, rhs, v, u )
c1***************************************************************
         write(i2,*) '  Solution by smoother on finest grid'
         call vecout( i2, idim, levels, np, u )
         call axplsy( idim, levels, np, -1d0, ue, u )
        write(i2,*)' Error on finest grid after solution with smoother'
         call vecout( i2, idim, levels, np, u )
         errmax = anorm( 1, idim, u, levels, np )
         errl2 = anorm( 2, idim, u, levels, np )
         write(i2,'(a,2e13.5)')
     $         'Maximum and L2 norm of the error = ',errmax, errl2
         write(*,'(a,2e13.5)')
     $         'Maximum and L2 norm of the error = ',errmax, errl2
      endif
      write(*,*)'                         **************************'
      write(*,*)'                          results are in mglab.res'
      write(*,*)'                         **************************'
      close(i2)
c2***************************************************************
      end
      double precision function a(icase,aa,x)
c***************************************************************
c
c User defined function to set the diffusion coefficient
c
c***************************************************************
      integer icase
      double precision aa,x,y
      if ( icase .eq. 1 ) then
c
c interface problem according to eqs. (3.3.1) and (3.3.6); xstar=aa,
c epsilon = 1d-3
c
         if ( x. le. aa) then
            y = 1d-3
         else
            y= 1d0
         endif
         a = y
      else if ( icase .eq. 2 ) then
c
c convection-diffusion problem
c
         a = aa
      else
         a = aa
      endif
      end
      double precision function anorm( ic, idim, b, k, np )
c **************************************************************************
c
c     computes norms of vector b on grid k
c         ic = 1:    maximum norm
c         ic = 2:    l2-norm scaled by dimension of b
c ****************************************************************************
c
c                  INPUT  PARAMETERS
c
      implicit none
      integer ic, idim, k, np(10)
      double precision b(idim)
c ***************************************************************************
c
c                 LOCAL QUANTITIES
c
      integer i
      double precision a
c ***************************************************************************
c
c                 I/O
c                none
c
c ***************************************************************************
c
c            FUNCTIONS
c
c              none
c
c **************************************************************************
c
c            SUBROUTINES CALLED
c
c               none
c
c ***************************************************************************
c
c            METHOD
c
c   Can be done with BLAS subroutines IDAMAX, DNRM2 and SCAL; has not been
c   implemented
c
c ***************************************************************************
      a = 0.d0
      if ( ic .eq. 1 ) then
         do 10 i = np(k), np(k+1) - 1
            a = max(a, abs(b(i)))
  10     continue
         anorm = a
      else if ( ic .eq. 2 ) then
         do 20 i = np(k), np(k+1)-1
            a = a + b(i)*b(i)
  20     continue
         anorm = sqrt(a/(np(k+1)-np(k)))
      endif
c     end of anorm
      end
      subroutine axplsy( idim, k, np, a, ut, u )
c***********************************************************************
c
c     u = a*ut + u    on grid k
c                     structure analogous to BLAS routine DAXPY
c
c ****************************************************************************
c
c                  INPUT/OUTPUT PARAMETERS
c
      implicit none
      integer idim, k, np(11)
      double precision a, u(idim), ut(idim)
c ***************************************************************************
c
c                 LOCAL QUANTITIES
c
      integer i
c **************************************************************************
c
c            SUBROUTINES CALLED
c
c   BLAS routine DAXPY may be used on vector computers - not implemented
c
c ***************************************************************************
      do 10 i = np(k), np(k+1)-1
         u(i) = u(i) + a*ut(i)
  10  continue
c     end of axplusy
      end
      double precision function b(icase,bb,x)
c ***************************************************************************
c
c User defined function to set the velocity coefficient
c
c ***************************************************************************
      integer icase
      double precision bb,x
      if (icase .eq. 1) then
         b = bb
      else if ( icase .eq. 2 ) then
         b = bb
      else
         b = bb
      endif
      end
      double precision function c(icase,cc,x)
c ***************************************************************************
c
c User defined subroutine to set the third coefficient
c
c ***************************************************************************
      implicit none
      integer icase
      double precision cc,x
      if (icase .eq. 1) then
         c = cc
      else if ( icase .eq. 2 ) then
         c = cc
      else
         c = cc
      endif
      end
      subroutine chrfun( i2, ib, ca, cab, cp, cr, scpmtx )
c******************************************************************************
c
c  This subroutine defines the characteristic functions of A, A(bar), P and R,
c  and computes the lower and upper bound of the structure of A on coaser grids
c  according to the choice of P and R
c  Used in Galerkin coarse grid approximation
c
c******************************************************************************
c
c  Input/Output parameters
c
      implicit none
      integer i2, ib(10), ca(-1:1), cab(-1:1), cp(-2:1),
     $        cr(-2:1), scpmtx
c******************************************************************************
c
c  i2      i   logical output device number
c  ib      i   for further information, see subroutine galerk
c  ca      o   characteristic function of fine grid operator A _
c  cab     o   characteristic function of coarse grid operator A
c  cp      o   characteristic function of prolongation operator P
c  cr      o   characteristic function of restriction operator R
c  scpmtx  o/i -scpmtx and scpmtx are the assumed lower and upper bound of the
c              structure of the matrix A. In our cases here, for every possible
c              combination of P and R, scpmtx has only two possible values:
c              either 1 or 2.
c
c******************************************************************************
c
c  Local quantities
      integer i
c
c  i           do-loop variable
c
c******************************************************************************
c
c Initialization
c
      do 5 i=-2,1
         cr(i)=0
         cp(i)=0
 5    continue
      do 6 i=-1,1
         ca(i)=0
         cab(i)=0
 6    continue
c
c Define cr, cp, ca and cab according to ib(6) and ib(7)
c
      if ( ib(6) .eq. 1 ) then
         cp(-1)=1
         cp(0)=1
         cp(1)=1
      else if ( ib(6) .eq. 2 ) then
         cp(-1)=1
         cp(0)=1
      else if ( ib(6) .eq. 3 .or. ib(6) .eq. 4 ) then
         cp(-2)=1
         cp(-1)=1
         cp(0)=1
         cp(1)=1
      else
         write(i2,*) 'chrfun : no such restriction available'
         stop
      end if
c
      if ( ib(7) .eq. 1 .or. ib(7) .eq. 2 ) then
         cr(-1)=1
         cr(0)=1
         cr(1)=1
      else if ( ib(7) .eq. 3 ) then
         cr(-1)=1
         cr(0)=1
      else if ( ib(7) .eq. 4 .or. ib(7) .eq. 5 ) then
         cr(-2)=1
         cr(-1)=1
         cr(0)=1
         cr(1)=1
      else
         write(i2,*) 'chrfun: no such prolongation available'
         stop
      end if
c
      do 10 i=-scpmtx,scpmtx
         ca(i)=1
 10   continue
c
c If enlargement of stencil occurs, then scpmtx=2
c
      if ( ib(6) .eq. 3 .or. ib(6) .eq. 4 ) then
         if ( ib(7) .ne. 3 ) scpmtx=2
      end if
      do 15 i=-scpmtx,scpmtx
         cab(i)=1
 15   continue
c
c end of subroutine chrfun
c
      return
      end
      subroutine copy( idim, k, np, u, v )
c***********************************************************************
c
c     u = v    on grid k
c
c ****************************************************************************
c
c                  INPUT  PARAMETERS
c
      implicit none
      integer idim, k, np(11)
      double precision u(idim), v(idim)
c ***************************************************************************
c
c                 LOCAL QUANTITIES
c
      integer i
c **************************************************************************
c
c            SUBROUTINES CALLED
c
c   BLAS routine DCOPY may be used on vector computers - has been
c     commented out with cv
c
c ***************************************************************************
      do 10 i = np(k), np(k+1)-1
         u(i) = v(i)
  10  continue
c
c  vector version
c
cv      call DCOPY (np(k+1) - np(k),u(np(k)),1,v(np(k)),1)
c     end of copy
      end
      subroutine defres( idim, io, i2, k, np, amat, dmat, b, bdef, ub )
c **************************************************************************
c
c     residual computation for defect correction according to section 4.6
c     of the book
c
c     bdef = b - dmat*ub + amat*ub            on grid k
c
c     no attempt has been made to emulate BLAS routines
c
c ****************************************************************************
c
c                  INPUT  PARAMETERS
c
      implicit none
      integer idim, io, i2, k, np(11)
      double precision amat(idim,3), dmat(idim,3), b(idim), bdef(idim),
     $                 ub(idim)
c
c  amat    less accurate discretization matrix
c  dmat    more accurate discretization matrix
c  b       right-hand-side of more accurate discretization
c  bdef    defect correction residual
c  io      governs output
c  i2      output device number
c  ub      solution to which defect correction is applied
c
c ***************************************************************************
c
c                 LOCAL QUANTITIES
c
      integer i
c
c ***************************************************************************
      do 5 i = np(k), np(k+1)-1
         bdef(i) = b(i) - (dmat(i,2)-amat(i,2))*ub(i)
    5 continue
      do 6 i = np(k), np(k+1)-2
         bdef(i) = bdef(i) - (dmat(i,3)-amat(i,3))*ub(i+1)
    6 continue
      do 7 i = np(k)+1, np(k+1)-1
         bdef(i) = bdef(i) - (dmat(i,1)-amat(i,1))*ub(i-1)
    7 continue
      if ( io .gt. 2 ) then
         write(i2,*)'Right hand side for defect correction'
         call vecout( i2, idim, k, np, bdef )
      endif
c     end of defres
      end
      subroutine defrp1( i2, ib, p, r )
c
c******************************************************************************
c
c  Subroutine defrp1 defines operator-independent prolongation and restriction
c  Used in Galerkin coarse grid approximation
c
c******************************************************************************
c
c  Input/Output
c
      implicit none
      integer i2, ib(10)
      double precision p(-2:1,3), r(-2:1,3)
c******************************************************************************
c
c  i2     i   logical output device number
c  ib     i   for further information, see subroutine galerk
c  p      o   array (-2:1,1:3), stores the elements of prolongation operator;
c             the operator for the interior points are stored in p(*,2); that
c             for the left boundary in p(*,1) and that for the right boundary
c             in  p(*,3)
c  r      o   array (-2:1,1:3), stores the elements of restriction operator;
c             the storage arrangement is similar to p
c
c ***************************************************************************
c
c  LOCAL QUANTITIES
c
      integer idx1, idx2
c
c  idx1,idx2,  scratch variables
c
c ***************************************************************************
c
c Initialization
c
      do 6 idx1 = 1, 3
         do 7 idx2 = -2, 1
            p(idx2,idx1) = 0d0
            r(idx2,idx1) = 0d0
    7    continue
    6 continue
c
c Define prolongation operator
c
      if ( ib(6) .eq. 1 ) then
         p(-1,1)=0d0
         p(0,1)=1d0
         p(1,1)=.5d0
         p(-1,2)=.5d0
         p(0,2)=1.d0
         p(1,2)=.5d0
         p(-1,3)=.5d0
         p(0,3)=1.d0
         p(1,3)=0.d0
      else if ( ib(6) .eq. 2 ) then
         p(-1,1)=1.d0
         p(0,1)=1.d0
         p(-1,2)=1.d0
         p(0,2)=1.d0
         p(-1,3)=1.d0
         p(0,3)=1.d0
      else if ( ib(6) .eq. 3 ) then
         p(-2,1)=0.d0
         p(-1,1)=.5d0
         p(0,1)=.75d0
         p(1,1)=.25d0
         p(-2,2)=.25d0
         p(-1,2)=.75d0
         p(0,2)=.75d0
         p(1,2)=.25d0
         p(-2,3)=.25d0
         p(-1,3)=.75d0
         p(0,3)=.5d0
         p(1,3)=0.d0
      else if ( ib(6) .eq. 4 ) then
         p(-2,1)=0.d0
         p(-1,1)=1.d0
         p(0,1)=.75d0
         p(1,1)=.25d0
         p(-2,2)=.25d0
         p(-1,2)=.75d0
         p(0,2)=.75d0
         p(1,2)=.25d0
         p(-2,3)=.25d0
         p(-1,3)=.75d0
         p(0,3)=1.d0
         p(1,3)=0.d0
      else
        write(i2,*) 'defrp1: ib(6) wrong, this option not implemented'
        stop
      end if
c
c  Define restriction operator
c
      if ( ib(7) .eq. 1 ) then
         r(-1,1)=0.d0
         r(0,1)=.5d0
         r(1,1)=.25d0
         r(-1,2)=.5d0
         r(0,2)=1.d0
         r(1,2)=.5d0
         r(-1,3)=.25d0
         r(0,3)=.5d0
         r(1,3)=0.d0
      else if ( ib(7) .eq. 2 ) then
         r(-1,1)=0.d0
         r(0,1)=1.5d0
         r(1,1)=.5d0
         r(-1,2)=.5d0
         r(0,2)=1.d0
         r(1,2)=.5d0
         r(-1,3)=.5d0
         r(0,3)=1.5d0
         r(1,3)=0.d0
      else if ( ib(7) .eq. 3 ) then
         r(-1,1)=1.d0
         r(0,1)=1.d0
         r(-1,2)=1.d0
         r(0,2)=1.d0
         r(-1,3)=1.d0
         r(0,3)=1.d0
      else if ( ib(7) .eq. 4 ) then
         r(-2,1)=0.d0
         r(-1,1)=.75d0
         r(0,1)=.75d0
         r(1,1)=.25d0
         r(-2,2)=.25d0
         r(-1,2)=.75d0
         r(0,2)=.75d0
         r(1,2)=.25d0
         r(-2,3)=.25d0
         r(-1,3)=.75d0
         r(0,3)=.75d0
         r(1,3)=0.d0
      else if ( ib(7) .eq. 5 ) then
         r(-2,1)=0.d0
         r(-1,1)=1.d0
         r(0,1)=0.75d0
         r(1,1)=.25d0
         r(-2,2)=.25d0
         r(-1,2)=.75d0
         r(0,2)=.75d0
         r(1,2)=.25d0
         r(-2,3)=.25d0
         r(-1,3)=.75d0
         r(0,3)=1.d0
         r(1,3)=0.d0
      else
         write(i2,*) 'defrp1: ib(7) wrong, this option not implemented'
         stop
      end if
c
c end of subroutine defrp1
      return
      end
      subroutine diagn( i, idim, i2, j, k, np, v )
c***********************************************************************
c
c     performs diagnostic output
c
c ****************************************************************************
c
c                  INPUT  PARAMETERS
c
      implicit none
      integer i, idim, i2, j, k, np(10)
      double precision v(idim)
c
c   i     i    chooses text for output
c   j     i    determines amount of output
c   k     i    grid level
c
c ***************************************************************************
c
c                 LOCAL QUANTITIES
c
      double precision a1, a2, anorm
c ***************************************************************************
c
c                 I/O
c     output to logical unit number i2
c
c ***************************************************************************
c
c            FUNCTIONS
c
c              none
c
c **************************************************************************
c
c            SUBROUTINES CALLED
c
c               vecout
c
c ***************************************************************************
      write(i2,*)' diagnostic #',i
      if ( i .eq. 1 ) then
         if ( j .gt. 1 ) then
            write(i2,*) ' maxnorm, l2-norm of right-hand-side on',
     $                  ' coarsest grid (level',k,')'
            a1 = anorm( 1, idim, v, k, np )
            a2 = anorm( 2, idim, v, k, np )
            write(i2,'(2e11.3)') a1, a2
         endif
         if ( j .gt. 2 ) then
            write(i2,*)'  right-hand-side on coarsest grid (level',k,')'
            call vecout( i2, idim, k, np, v )
         endif
      else if ( i .eq. 2 ) then
         write(i2,*)' maxnorm, l2-norm of residual after direct ',
     $              'solution'
         write(i2,*)' on coarsest grid (level=',k,')'
         a1 = anorm( 1, idim, v, k, np )
         a2 = anorm( 2, idim, v, k, np )
         write(i2,'(2e11.3)') a1, a2
      else if ( i .eq. 3 ) then
         write(i2,*)' current solution',
     $              ' just before move to finer grid (level',k,')'
         call vecout( i2, idim, k, np, v )
      else if ( i .eq. 4 ) then
         write(i2,*)' correction',
     $              ' just before move to finer grid (level',k,')'
         call vecout( i2, idim, k, np, v )
      else if ( i .eq. 5 ) then
         write(i2,*)' prolongated correction',
     $              ' just after move to finer grid (level',k,')'
         call vecout( i2, idim, k, np, v )
      else if ( i .eq. 6 ) then
         write(i2,*)' maxnorm, l2-norm of residual',
     $              ' after coarse grid (level',k,')'
         write(i2,*)' correction'
         a1 = anorm( 1, idim, v, k, np )
         a2 = anorm( 2, idim, v, k, np )
         write(i2,'(2e11.3)') a1, a2
      else if ( i .eq. 7 ) then
         write(i2,*)' current solution',
     $              ' after coarse grid correction (level',k,')'
         call vecout( i2, idim, k, np, v )
      else if ( i .eq. 8 ) then
         write(i2,*)' maxnorm, l2-norm of residual',
     $              ' after post-smoothing (level',k,')'
         a1 = anorm( 1, idim, v, k, np )
         a2 = anorm( 2, idim, v, k, np )
         write(i2,'(2e11.3)') a1, a2
      else if ( i .eq. 9 ) then
         write(i2,*)' maxnorm, l2-norm of residual',
     $              ' before pre-smoothing (level',k,')'
         a1 = anorm( 1, idim, v, k, np )
         a2 = anorm( 2, idim, v, k, np )
         write(i2,'(2e11.3)') a1, a2
      else if ( i .eq. 10 ) then
         write(i2,*)' maxnorm, l2-norm of residual',
     $              ' after pre-smoothing (level',k,')'
         a1 = anorm( 1, idim, v, k, np )
         a2 = anorm( 2, idim, v, k, np )
         write(i2,'(2e11.3)') a1, a2
      else if ( i .eq. 11 ) then
         write(i2,*)' maxnorm, l2-norm of error',
     $              ' at moment of leaving finest grid (level',k,')'
         a1 = anorm( 1, idim, v, k, np )
         a2 = anorm( 2, idim, v, k, np )
         write(i2,'(2e11.3)') a1, a2
      else if ( i .eq. 12 ) then
         write(i2,*)' maxnorm, l2-norm of final error'
         a1 = anorm( 1, idim, v, k, np )
         a2 = anorm( 2, idim, v, k, np )
         write( i2,'(2e11.3)') a1, a2
      else if ( i .eq. 13 ) then
         write(i2,*)' maxnorm, l2-norm (level',k,')'
         a1 = anorm( 1, idim, v, k, np )
         a2 = anorm( 2, idim, v, k, np )
         write( i2,'(2e11.3)') a1, a2
      else
      endif
c     end of diagn
      end
      subroutine dirsol( idim, i2, k, np, amat, u, b )
c **************************************************************************
c
c     computes solution of amat*u = b   on grid k by direct solution method
c
c ****************************************************************************
c
c                  INPUT  PARAMETERS
c
      implicit none
      integer idim, i2, k, np(11)
      double precision amat(idim,3), u(idim), b(idim)
c ***************************************************************************
c
c                 LOCAL QUANTITIES
c
      integer j0, j, n, nmax
      parameter (nmax = 1000)
      double precision bet, gam(nmax)
      if ( idim .gt. nmax ) then
         write(i2,*)' dirsol: nmax < idim '
         stop
      endif
      j0 = np(k)-1
      n = np(k+1) - np(k)
      bet = amat(j0+1,2)
      if ( bet .eq. 0d0 ) then
         write(i2,*) 'dirsol: zero pivot found in row ',1
         write(i2,*) '        take a look at your matrix'
         stop
      endif
      u(j0+1) = b(j0+1)/bet
      do 10 j = j0+2, j0+n
         gam(j) = amat(j-1,3)/bet
         bet = amat(j,2) - amat(j,1)*gam(j)
         if ( bet .eq. 0d0 ) then
            write(i2,*) 'dirsol: zero pivot found in row ',j-j0
            write(i2,*) '        take a look at your matrix'
            stop
         endif
         u(j) = ( b(j) - amat(j,1) * u(j-1) ) / bet
   10 continue
      do 20 j = j0+n-1, j0+1, -1
         u(j) = u(j) - gam(j+1)*u(j+1)
   20 continue
c     end of dirsol
      end
      subroutine err(i1,i2,i)
c***********************************************************************
c
c   Write error message for mglab program and stop the program
c
c***********************************************************************
      integer i

      write( i2, '(a)' ) '********** ERROR ***********'
      if (i .eq. 1) then
         write( i2, '(a)' )
     $          ' levels too large: increase dimension of np to '
         write( i2, '(a)' )' levels+1'
         goto 10
      else if (i .eq. 3) then
         write( i2, '(a)' )' ib(4) out of range '
         goto 10
      else if (i .eq. 4) then
         write( i2, '(a)' )' dimension of amat too small'
         goto 10
      else if (i .eq. 5) then
         write( i2, '(a)' )' igal out of range'
         goto 10
      endif
  10  write(*,*)'******************************************'
      write(*,*)' !!! error  !!!! results are in discr.res'
      write(*,*)'******************************************'
      write( i2, '(a)' ) '***** PROGRAM STOPPED ******'
      close(i2)
      stop
c     end of err
      end
      subroutine errmg (i2,i)
c***********************************************************************
c
c   Write error message for mglab program and stop the program
c
c***********************************************************************
      implicit none
      integer i2,i
      if (i.eq.1) then
          write(i2,*) ' igamma out of range in subroutine smooth'
          write(*,*)' aborted - error in smoothing'
          stop
       else
          write(i2,*) ' i out of range in subroutine errmg'
          write(*,*)' aborted - error in errmg'
       endif
c      end of errmg
       end
      subroutine galerk( ib, idim, iout, i2, levels, np, amat, scpmtx )
c **************************************************************************
c
c  computation of coarse grid matrices by Galerkin approximation
c
c **************************************************************************
c
c  INPUT / OUTPUT PARAMETERS:
c
      implicit none
      integer ib(10), idim, iout, i2, levels, np(11), scpmtx

      double precision amat(idim,-1:1)
c
c   ib      i  ib(6) chooses type of prolongation - see subroutine prol:
c               ic = ib(6)
c   ib(6) = 1:  vertex-centered linear interpolation: equation (5.3.1) of
c               the book: (1/2)*[1 2 1]
c   ib(6) = 2:  cell-centered piecewise constant interpolation:
c               equation (5.3.5) of the book: [1 1]
c   ib(6) = 3:  cell-centered linear interpolation with Dirichlet boundary
c               modification:  equation (5.3.6) of the book with w=0
c               in equation (5.3.24) at boundaries: (1/4)*[w 2+w 2+e e]
c   ib(6) = 4:  cell-centered linear interpolation with Neumann boundary
c               modification:  equation (5.3.6) of the book with stencil
c               [P*] = (1/4)*[ 0  4  3  1 ] at boundary x=0 and similar at x=1
c   ib(6) = 5:  vertex-centered operator dependent prolongation
c               equations (5.4.2), (5.4.3) of the book: NOT IMPLEMENTED
c
c               ib(7) chooses type of restriction for residual - see
c               subroutine restri: ic = ib(7)
c   ib(7) = 1:  vertex-centered restriction according to the one-
c               dimensional version of equation (5.3.19) of the book
c               stencil = 0.5 * [1 2 1]
c               boundary modification type Dirichlet:
c               at x=0: stencil = 0.5 * [0 2 1];
c               at x=1: stencil = 0.5 * [1 2 0]
c   ib(7) = 2:  vertex-centered restriction according to the one-
c               dimensional version of equation (5.3.19) of the book
c               stencil = 0.5 * [1 2 1]
c               boundary modification type Neumann:
c               at x=0: stencil = 0.5 * [0 3 1];
c               at x=1: stencil = 0.5 * [1 3 0]
c   ib(7) = 3:  cell-centered restriction according to equation (5.3.23)
c               of the book
c               stencil = [1  1]
c   ib(7) = 4:  cell-centered restriction according to equation (5.3.24)
c               of the book
c               stencil = (1/4)*[1 3 3 1]
c               boundary modification type Dirichlet:
c               at x=0: stencil = (1/4)*[0 3 3 1];
c               at x=1: stencil = (1/4)*[1 3 3 0]
c   ib(7) = 5:  cell-centered restriction according to equation (5.3.24)
c               of the book
c               stencil = (1/4)*[1 3 3 1]
c               boundary modification type Neumann:
c               at x=0: stencil = (1/4)*[0 4 3 1];
c               at x=1: stencil = (1/4)*[1 3 4 0]
c   ib(7) = 6:  vertex-centered operator-dependent restriction according
c               to equation (5.4.22) of the book
c               stencil = [a1 a2 a3]
c               boundary modification type Neumann:
c               at x=0: a1 = 0, a2:= 2*a2-a3
c               at x=1: a3=0, a2:=2*a2-a1: NOT IMPLEMENTED
c   ib(7) = 7:  vertex-centered operator-dependent restriction according
c               to equation (5.4.22) of the book
c               stencil = [a1 a2 a3]
c               boundary modification type Dirichlet:
c               at x=0: a1 = 0; at x=1 : a3 = 0: NOT IMPLEMENTED
c  idim   *   parameter that passes the first dimension of amat
c  iout   i   governs amount of output; iout larger -> more output
c  i2     i   output device number
c  levels i   number of grids in multigrid method; < 16
c  np     i   array of pointers to start of grids and to 1+end of finest grid
c             np(1)=1     : start of coarsest grid
c             np(levels)  : start of finest grid
c             np(levels+1): end of finest grid + 1
c  amat   o   matrix (discretization of differential equation) on coarse grids
c  scpmtx i/o -scpmtx and scpmtx are the assumed lower and upper bound of the
c             structure of the matrix A. The matrix A on the finest grid is
c             assumed to have three-point stencil, which implies scpmtx=1. It
c             may be changed to 2 if enlargement of stencils of A occurs.
c
c ***************************************************************************
c
c  LOCAL QUANTITIES:
c
      integer ca(-1:1), cab(-1:1), cp(-2:1), cr(-2:1),
     $     i, idx1, idx2, idx3, idx4, ip, k, level, m, n, q
      double precision p(-2:1,1:3), r(-2:1,1:3)
c ***************************************************************************
c
c  ca         characteristic function of fine grid operator A _
c  cab        characteristic function of coarse grid operator A
c  cp         characteristic function of prolongation operator P
c  cr         characteristic function of restriction R
c  i          do-loop variable over the coarse grid points
c  idx1,idx2, scratch variables
c  idx3,idx4
c  ip         ip=1, if vertex-centered coarsening is used;
c             ip=0, if cell-centered coarsening is used.
c             For vertex-centered grid, the starting number of grid points
c             is zero; the final number is the number of grid points
c             minus 1. But for cell-centered grid, the starting number
c             of grid points is one; the final number is the number of
c             grid points.
c  k          do-loop variable over the structure of matrix A on the
c             fine grid
c  level      do-loop variable over grid levels
c  m          do-loop variable over the structure of restriction operator
c  n          do-loop variable over the structure of matrix A on the
c             coarse grid
c  q          variable whose range is within the scope of the structure of
c             the transpose of prolongation operator
c  p          array (-2:1,1:3), stores the elements of prolongation operator;
c             the operator for the interior points are stored in p(*,2); that
c             for the left boundary in p(*,1) and that for the right boundary
c             in  p(*,3)
c  r          array (-2:1,1:3), stores the elements of restriction operator;
c             the storage arrangement is similar to p
c
c****************************************************************************
c
c  FUNCTIONS CALLED: none
c
c****************************************************************************
c
c  SUBROUTINES CALLED:
c
c  chrfun: defines the characteristic functions of A, A(bar), P and R, resp.,
c          and computes the lower and upper bound of the structure of A on
c          coarser grids according to the choice of P and R
c  defrp1: defines operator-independent prolongation and restriction operator
c  defrp2: defines operator-dependent prolongation and restriction operator
c  ibchck: checks whether the choices of ib(6) and ib(7) are consistent, and
c          computes the flag ip
c
c****************************************************************************
c
c  METHOD:
c
c  Galerkin coarse grid approximation according to algorithm CALRAP
c    in section 6.2 of the book
c
c****************************************************************************
c
c Check whether combination of ib(6) and ib(7) is correct and compute ip
c
      call ibchck( i2, ib, ip )
c
c Define operator-independent prolongation and restriction operator
c
      call defrp1( i2, ib, p, r )
c
c Define the characteristic functions of A, A(bar), P and R, and compute the
c possible lower and upper bound of the structure of A on the coaser grid.
c Here, cp and cr are independent of grid level, but cab is not (stencils of A
c could be enlarged). In our case, the bound does not exceed 2.
c
      call chrfun( i2, ib, ca, cab, cp, cr, scpmtx )
      do 100 level = levels-1,1,-1
c
c Initialize the matrix on the coarse grid to zero
c
        do 5 idx1 = np(level),np(level+1)-1
          do 10 idx2 = -1,1
            amat(idx1,idx2) = 0.d0
   10     continue
    5   continue
c
c CALRAP procedure:
c compute the matrix on the coarse grid
c loop over the structure of A on the coarse grid
c
        do 15 n = -1,1
          if ( cab(n) .eq. 1 ) then
c
c loop over the structure of restriction operator
c
            do 20 m = -2,1
              if ( cr(m) .eq. 1 ) then
c
c loop over the structure of A on the fine grid
c
                do 25 k = -1,1
                  if ( ca(k) .eq. 1 ) then
                    q = m+k-2*n
                    if ( q .ge. -2 .and. q .le. 1 ) then
c
c q belongs to the structure of the transpose of prolongation operator
c
                      if ( cp(q) .eq. 1 ) then
c
c Now idx1 and idx2 are the number of coarse gird and fine grid points,
c respectively
c
                        idx1 = np(level+1)-np(level)
                        idx2 = np(level+2)-np(level+1)
c
c treat vertex-centered case (grid points denoted by *)
c             numbering on fine grid        *  *  *  *  *  *  *
c                                           0  1  2  3  4  5  6
c             numbering on coarse grid      *     *     *     *
c                                           0     1     2     3
c
                        if ( ip .eq. 1 ) then
c
c loop over the interior coarse grid points
c i is local coarse grid coordinate; see picture above
c
                          do 30 i=1,np(level+1)-np(level)-2
c
c  i+n is a coarse grid point and 2i+m is a fine grid point
c
                            if ( (i+n.ge.0 .and. i+n.le.idx1-1) .and.
     $                          (2*i+m.ge.0 .and. 2*i+m.le.idx2-1)) then
                              idx3 = np(level)+i
                              idx4 = np(level+1)+2*i+m
c
c if p or r is operator-dependent, call defrp2 (not implemented)
c
c                             call defrp2(i, ib, im, idim, level, ip,
c                                         np, scpmtx, amat, solut, p, r)
c
                              amat(idx3,n) = amat(idx3,n) +
     $                                       r(m,2)*amat(idx4,k)*p(q,2)
                            end if
 30                       continue
c
c left boundary point
c
                          i = 0
                          if ( (i+n.ge.0.and.i+n.le.idx1-1) .and.
     $                       (2*i+m.ge.0.and.2*i+m.le.idx2-1) ) then
                            idx3 = np(level)
                            idx4 = np(level+1)+2*i+m
c
c if p or r is operator-dependent, call defrp2 (not implemented)
c
c                           call defrp2( i, ib, im, idim, level, ip, np,
c                                        scpmtx, amat, solut, p, r )
c
c                           operator p for the boundary point is used
c
                            if ( i+n .eq. 0 ) then
                              amat(idx3,n) = amat(idx3,n)+
     $                                       r(m,1)*amat(idx4,k)*p(q,1)
c
c operator p for the interior point is used
c
                            else
                              amat(idx3,n) = amat(idx3,n) +
     $                                       r(m,1)*amat(idx4,k)*p(q,2)
                            end if
                          end if
c
c right boundary point
c
                          i = idx1-1
                          if ( (i+n.ge.0.and.i+n.le.idx1-1) .and.
     $                       (2*i+m.ge.0.and.2*i+m.le.idx2-1)) then
                            idx3 = np(level+1)-1
                            idx4 = np(level+1)+2*i+m
c
c if p or r is operator-dependent, call defrp2 (not implemented)
c
c                           call defrp2( i, ib, im, idim, level, ip,
c     $                                 np, scpmtx, amat, solut, p, r )
c
c operator p for the boundary point is used
c
                            if ( i+n .eq. idx1-1 ) then
                              amat(idx3,n) = amat(idx3,n) +
     $                                       r(m,3)*amat(idx4,k)*p(q,3)
c
c operator p for the interior point is used
c
                            else
                              amat(idx3,n) = amat(idx3,n) +
     $                                       r(m,3)*amat(idx4,k)*p(q,2)
                            end if
                          end if
c
c treat cell-centered case (grid points denoted by *)
c               numbering on fine grid     | * | * | * | * | * | * |
c                                            1   2   3   4   5   6
c               numbering on coarse grid   |   *   |   *   |   *   |
c                                              1       2       3
c
                        else
c
c loop over the interior coarse grid points. i is local coarse grid
c coordinate; see picture above
c
                          do 35 i = 2,np(level+1)-np(level)-1
c
c i+n is a coarse grid point and 2i+m is a fine grid point
c
                            if ( (i+n.ge.1.and.i+n.le.idx1) .and.
     $                         (2*i+m.ge.1.and.2*i+m.le.idx2) ) then
                              idx3 = np(level)+i-1
                              idx4 = np(level+1)+2*i+m-1
c
c if p or r is operator-dependent, call defrp2 (not implemented)
c
c                             call defrp2( i, ib, im, idim, level, ip,
c    $                                  np, scpmtx, amat, solut, p, r )
c
c operator p for left boundary point is used
c
                              if ( i+n .eq. 1 ) then
                                amat(idx3,n) = amat(idx3,n) +
     $                                     r(m,2)*amat(idx4,k)*p(q,1)
c
c operator p for right boundary point is used
c
                              else if (i+n .eq. idx1 ) then
                                amat(idx3,n) = amat(idx3,n) +
     $                                      r(m,2)*amat(idx4,k)*p(q,3)
c
c operator p for interior points
c
                              else
                                amat(idx3,n) = amat(idx3,n) +
     $                                      r(m,2)*amat(idx4,k)*p(q,2)
                              end if
                            end if
 35                       continue
c
c left boundary point
c
                          i = 1
                          if ( (i+n.ge.1.and.i+n.le.idx1) .and.
     $                       (2*i+m.ge.1.and.2*i+m.le.idx2) ) then
                            idx3 = np(level)
                            idx4 = np(level+1)+2*i+m-1
c
c if p or r is operator-dependent, call defrp2 (not implemented)
c
c                           call defrp2( i, ib, im, idim, level, ip,
c     $                                 np, scpmtx, amat, solut, p, r )
c
c operator p for the boundary point is used
c
                            if ( i+n .eq. 1 ) then
                              amat(idx3,n) = amat(idx3,n) +
     $                                   r(m,1)*amat(idx4,k)*p(q,1)
c
c operator p for the interior point is used
c
                            else
                              amat(idx3,n) = amat(idx3,n) +
     $                                    r(m,1)*amat(idx4,k)*p(q,2)
                            end if
                          end if
c
c right boundary point
c
                          i = idx1
                          if ( (i+n.ge.1.and.i+n.le.idx1) .and.
     $                       (2*i+m.ge.1.and.2*i+m.le.idx2) ) then
                            idx3 = np(level+1)-1
                            idx4 = np(level+1)+2*i+m-1
c
c if p or r is operator-dependent, call defrp2 (not implemented)
c
c                           call defrp2( i, ib, im, idim, level, ip, np,
c     $                                  scpmtx, amat, solut, p, r )
c
c operator p for the boundary point is used
c
                            if ( i+n.eq.idx1 ) then
                              amat(idx3,n) = amat(idx3,n) +
     $                                    r(m,3)*amat(idx4,k)*p(q,3)
c
c operator p for the interior point is used
c
                            else
                              amat(idx3,n) = amat(idx3,n) +
     $                                   r(m,3)*amat(idx4,k)*p(q,2)
                            end if
                          end if
                        end if
                      end if
                    end if
                  end if
   25           continue
              end if
   20       continue
          end if
   15   continue
c
c Output:
c write the matrix for grid level=level
c
        write(i2,*) 'galerk: matrix at grid=',level
        do 40 i = np(level),np(level+1)-1
          write(i2,'(i3,3e11.3)')
     $          i-np(level)+1,(amat(i,idx1),idx1 = -scpmtx,scpmtx)
   40 continue
c
  100 continue
c
c end of subroutine galerkin
c
      return
      end
      subroutine gener( i, ib, idim, icase, np, amat, aa, bb, cc,
     $                  f0, f1, scal, rhs, ss )
c***********************************************************************
c
c  Generate matrix amat (tri-diagonal) on a given grid depending on
c  ib(4)
c
c***********************************************************************
c
c               INPUT / OUTPUT PARAMETERS
c
      implicit none
      integer i, ib(13), idim, icase, np(11)
      double precision aa, amat(1:idim,-1:1), bb, cc, f0, f1,
     $ rhs(1:idim), scal, ss
c
c  i       i  grid number
c  ib      i  see subroutine matrhs
c  idim    i  see subroutine matrhs
c  icase   i  see subroutine matrhs
c  np      i  see subroutine matrhs
c  a          diffusion coefficient (function)
c  aa      i  parameter for diffusion coefficient
c  amat    o  see subroutine matrhs
c  b          velocity coefficient (function)
c  bb      i  parameter for velocity coefficient
c  c          coefficient (function)
c  cc      i  parameter for c coefficient
c  f0      i  see subroutine matrhs
c  f1      i  see subroutine matrhs
c  rhs     o  see subroutine matrhs
c  s          right-hand-side of differential equation (function)
c  scal    i  scaling factor for dirichlet condition in vertex-centered case
c  ss      i  parameter for function s
c
c****************************************************************************
c
c            LOCAL QUANTITIES
c
      integer j
      double precision  h, a, al, ar, b, bl, br, c, fm, fp, om, omdif,
     $ omega, pe, s, w, x
c
c  j          grid point index
c  h          mesh size
c  a          diffusion coefficient (function)
c  al         diffusion coefficient
c  ar         diffusion coefficient
c  b          velocity coefficient (function)
c  bl         velocity coeficient
c  br         velocity coeficient
c  c          coefficient (function)
c  fm         coefficient in flux function: flux = fp*u(j+1) + fm*u(j)
c  fp         see fm
c  om         switching variable between upwind, central, hybrid
c  omdif      switching variable between upwind, central, hybrid
c  omega      switching variable between upwind, central, hybrid (function)
c  pe         mesh peclet number
c  s          right-hand-side of differential equation (function)
c  w          harmonic average of al, ar
c  x          x-coordinate
c
c ***************************************************************************
c
c            FUNCTIONS
c
c  a(icase,aa,x)  see subroutine matrhs
c  b(icase,bb,x)  see subroutine matrhs
c  c(icase,cc,x)  see subroutine matrhs
c  s(icase,ss,x)  see subroutine matrhs
c  omega(pe,type) switching function for hybrid scheme
c
c **************************************************************************
c
c                METHOD
c
c  Finite volume discretization of 1-dimensional second order boundary
c    value problem.
c    See chapter 3 and section 9.7 of the book
c  Upwind, central or hybrid discretization of convection (first order)
c    term, depending on ib(3).
c  Cell-centered or vertex-centered discretization, depending on ib(4)
c
c  Equation: -(ay')'+(by)'+cy = s
c
c  Boundary conditions:
c     Dirichlet: y(0) = f0  or  Neumann: -ay'(0) = f0
c     Dirichlet: y(1) = f1  or  Neumann: +ay'(1) = f1
c     The gridfunction values at the boundary points are not eliminated,
c     as would be possible in the vertex-centered case. Instead,
c     the corresponding equations are multiplied with scal. If scal >> 1
c     then the matrix is made symmetric by adding an artificial element.
c
c  The diffusion coefficient a(x) is allowed to be discontinuous at x = j*h
c    with h the stepsize. The user may put discontinuities elsewhere (and
c    will often have to on some coarse grids if igal = 0) at his own risk.
c
c ***********************************************************************
c
      if ( ib(4) .eq. 1 ) then
c
c=======================================================================
c  Vertex-centered
c=======================================================================
         h = 1d0/(dble( np(i+1) - np(i) ) - 1d0 )
c
         do 10 j = np(i), np(i+1)-1
               x = ( j - np(i) )*h
c
c  visit vertices
c
            amat(j,0) =  c( icase, cc, x ) * h
            rhs(j) =  s( icase, cc, x ) * h
  10     continue
c
c  the boundary finite volumes have half size
c
         amat(np(i),0) = 0.5d0*amat(np(i),0)
         amat(np(i+1)-1,0) = 0.5d0*amat(np(i+1)-1,0)
         rhs(np(i)) = 0.5d0*rhs(np(i))
         rhs(np(i+1)-1) = 0.5d0*rhs(np(i+1)-1)
c
c  visit interior cell interfaces
c
         do 11 j = np(i), np(i+1)-2
            x = (j - np(i) + 0.5)*h
            w = a(icase,aa,x)
            pe = h*b(icase,bb,x)/w
            om = omega (pe, ib(3))
            if ( ib(3).eq.2 .or. ib(3).eq.4 ) then
c
c  diffusion term not switched off
c
              omdif = 0d0
            else
              omdif = om
            endif
            bl = b(icase,bb,x-0.5d0*h)
            br = b(icase,bb,x+0.5d0*h)
c
c  flux = fp*u(j+1) + fm*u(j)
c
            fp =  -w*(1d0-omdif)/h + 0.5d0*((1d0-om)*br+om*(br-abs(br)))
            fm =   w*(1d0-omdif)/h + 0.5d0*((1d0-om)*bl+om*(bl+abs(bl)))
            amat(j,1) = fp
            amat(j+1,-1) =- fm
            amat(j,0) = amat(j,0)+ fm
            amat(j+1,0) = amat(j+1,0) - fp
  11     continue
c
c x = 0: boundary condition
c
         j = np(i)
         amat(j,-1) = 0d0
         if ( ib(1) .eq. 1 ) then
c
c Scale matrix with scal at Dirichlet boundary
c
            amat(j, 0) = scal*a(icase,aa,0d0)/h
            amat(j, 1) = 0d0
            rhs(j) = scal*a(icase,aa,0d0)*f0/h
c
c make matrix symmetric if scal large
c
            if ( scal .ge. 1d6) then
              amat(j,1) = amat(j+1,-1)
            endif
         else
c
c Neumann condition
c
            amat(j,0) = amat(j,0) - b(icase,bb,0d0)
            rhs(j) = rhs(j) + f0
         endif
c
c x = 1: boundary condition
c
         j = np(i+1) - 1
         amat(j,1) = 0d0
         if ( ib(2) .eq. 1 ) then
c
c Scale matrix with scal at Dirichlet boundary
c
            amat(j, 0) = scal*a(icase,aa,1d0)/h
            amat(j, -1) = 0d0
            rhs(j) = scal*a(icase,aa,1d0)* f1/h
c
c make matrix symmetric if scal large
c
            if ( scal .ge. 1d6) then
              amat(j, -1) = amat(j-1,1)
            endif
         else
c
c Neumann condition
c
            amat(j,0) = amat(j,0) + b(icase,bb,1d0)
            rhs(j) = rhs(j) + f1
         endif
      else if ( ib(4) .eq. 2 ) then
c
c=======================================================================
c  Cell-centered
c=======================================================================
c
         h = 1d0/dble( np(i+1) - np(i) )
c
c visit cell centers
c
         do 20 j = np(i), np(i+1)-1
            x = (j - np(i) + 0.5)*h
            amat(j,0) =  c( icase, cc, x ) * h
            rhs(j) =  s( icase, cc, x ) * h
  20     continue
c
c visit interior cell interfaces
c
         do 21 j = np(i), np(i+1)-2
c
c  consider right face of cell j
c
            x = (j - np(i) + 0.5)*h
            al = a(icase,aa,x)
            ar = a(icase,aa,x+h)
            w = 2*al*ar/(al+ar)
            pe = h*b(icase,bb,x+h/2d0)/w
            om = omega (pe, ib(3))
            if ( ib(3).eq.2 .or. ib(3).eq.4 ) then
c
c  diffusion term not switched off
c
              omdif = 0d0
            else
              omdif = om
            endif
            bl = b(icase,bb,x)
            br = b(icase,bb,x+h)
c
c  flux = fp*u(j+1) + fm*u(j)
c
            fp =  -w*(1d0-omdif)/h+ 0.5d0*((1d0-om)*br+om*(br-abs(br)))
            fm =  w*(1d0-omdif)/h+ 0.5d0*((1d0-om)*bl+om*(bl+abs(bl)))
            amat(j,1) = fp
            amat(j+1,-1) =- fm
            amat(j,0) = amat(j,0)+ fm
            amat(j+1,0) = amat(j+1,0) - fp
  21     continue
c
c x = 0: boundary condition
c
         j = np(i)
         amat(j,-1) = 0d0
         w = a(icase,aa,0.5d0*h)
         bl = b(icase,bb,0d0)
         br = b(icase,bb,0.5d0*h)
         pe = h*bl/w
         om = omega (pe, ib(3))
         if ( ib(3).eq.2 .or. ib(3).eq.4) then
c
c  diffusion term not switched off
c
              omdif = 0d0
            else
              omdif = om
            endif
         if ( ib(1) .eq. 1 ) then
c
c Dirichlet condition
c
c  flux = fm*f0 + fp*u(j)
c
            fm = 2d0*w*(1d0-omdif)/h+ (1d0-om)*bl+0.5d0*om*(bl+abs(bl))
            fp = -2d0*w*(1d0-omdif)/h+ 0.5d0*om*(br-abs(br))
         else
c
c Neumann condition
c
c  flux = fm*f0 + fp*u(j)
c
            fm =1d0-omdif+(1d0-om)*0.5d0*pe+om*0.25d0*h*(bl+abs(bl))/w
            fp = (1d0-om)*bl + om*0.5d0*(bl+abs(bl)+br-abs(br))
         endif
         amat(j,0) = amat(j,0) - fp
         rhs(j) = rhs(j) + fm*f0
c
c x = 1: boundary condition
c
         j = np(i+1)-1
         amat(j,1) = 0d0
         w = a(icase,aa,1d0-0.5d0*h)
         bl = b(icase,bb,1d0-0.5d0*h)
         br = b(icase,bb,1d0)
         pe = h*br/w
         om = omega (pe, ib(3))
         if ( ib(3).eq.2 .or. ib(3).eq.4 ) then
c
c  diffusion term not switched off
c
           omdif = 0d0
         else
           omdif = om
         endif
         if ( ib(2) .eq. 1 ) then
c
c Dirichlet condition
c
c  flux = fp*f1 + fm*u(j)
c
            fm = 2d0*w*(1d0-omdif)/h+0.5d0*om*(bl+abs(bl))
            fp = -2d0*w*(1d0-omdif)/h+ 0.5d0*om*(br-abs(br))+(1d0-om)*br
         else
c
c Neumann condition
c
            fp =-1d0+omdif+(1d0-om)*0.5d0*pe+om*0.25d0*h*(br+abs(br))/w
            fm = (1d0-om)*br + om*0.5d0*(bl+abs(bl)+br-abs(br))
         endif
            amat(j,0) = amat(j,0) + fm
            rhs(j) = rhs(j) - fp*f1
      endif
c     end of gener
      end
      subroutine ibchck( i2, ib, ip )
c*****************************************************************************
c
c  This subroutine checks whether combination of ib(6) and ib(7) is correct
c  and compute ip
c
c*****************************************************************************
c
c  Input/Output parameters
c
      implicit none
      integer i2,ib(10),ip
c*****************************************************************************
c
c  i2     i   logical output device number
c  ib     i   for more information, see subroutine galerk
c  ip     o   flag for type of grid used: ip=1, vertex-centered grid is used;
c             ip=0, cell-centered grid is used
c
c*****************************************************************************
c
c  Check the choice of prolongation
c
      if ( ib(6) .eq. 1 .or. ib(6) .eq. 5 ) then
         ip = 1
         if ( ib(7) .gt. 2 ) then
            write( i2, * )
     $      'ibchck: vertex-centered P matches cell-centered R'
            stop
         end if
      else if ( ib(6) .eq. 2 .or. ib(6) .eq. 3 .or. ib(6) .eq. 4 ) then
         ip = 0
         if ( ib(7) .eq. 1 .or. ib(7) .eq. 2
     $          .or. ib(7) .eq.6 .or. ib(7) .eq. 7 ) then
            write(i2,*)
     $      'ibchck: cell-centered P matches vertex-centered R'
            stop
         end if
      end if
c
c Usually, enlargement of stencils is not expected
c
      if ( ib(6) .eq. 3 .or. ib(6) .eq. 4 ) then
         if ( ib(7) .ne. 3 ) write( i2, * )
     $      'ibchck: stencil enlarged under RAP, you really want so?'
      end if
c
c end of subroutine ibchck
c
      return
      end
      subroutine info(icase,ib)
c***********************************************************************
c
c   User provided subroutine to print information concerning mglab
c
c***********************************************************************
      integer icase, ib(*)
      if (icase.eq.1) then
      else
      endif
c     end of info
      end
      subroutine inout( i1, i2 )
c***********************************************************************
c
c   Open the mglab input and output files
c
c***********************************************************************
      implicit none
      integer i1, i2

      open( i1, file='mglab.dat' )
      open( i2, file='mglab.res' )

      end
      subroutine input( i1, i2, ib, icase, igal, iout, ism, levels, nc,
     $                  nmg, ns, aa, bb, cc, delta, f0, f1, om, spar,
     $                  ss, tol, cycle )
c
c ***************************************************************************
c
c   read and write the user provided input for mglab
c
c ***************************************************************************

      implicit none
      integer ib(13), icase, igal, iout(4), ism, i1, i2, levels, nc,
     $        nmg, ns(2)
      double precision aa, bb, cc, delta, f0, f1, om, spar, ss, tol
      character cycle
c ***************************************************************************
c
c       Description of parameters in main program
c       This subroutine reads and prints input
c
c ***************************************************************************
      integer i
      character*32 string

      read( i1, '(a)' ) string
      read( i1, * ) ( ib(i), i=1,13 )
      read( i1, '(a)' ) string
      read( i1, * ) icase, igal, ( iout(i), i=1,4 ), ism, levels
      read( i1, '(a)' ) string
      read( i1, * ) nc, nmg, ns(1), ns(2)
      read( i1, '(a)' ) string
      read( i1, * ) aa, bb, cc, ss, f0, f1
      read( i1, '(a)' ) string
      read( i1, * ) om, spar, delta, tol
      read( i1, '(a)' ) cycle
      write( i2, '(a)' )
     $ '********** TUTORIAL MULTIGRID LABORATORY ********************'
      write( i2, '(a)' ) ' '
      write( 12, '(a)' )
     $ '   ib1  ib2  ib3  ib4  ib5  ib6  ib7  ib8  ib9 ib10 ib11 ib12  i
     $b13'
      write( i2, '(13i5)' ) (ib(i),i=1,13)
      write( i2, '(a)' ) 'icase igal iout              ism levels   nc'
      write( i2, '(9i5)' ) icase, igal, ( iout(i), i=1,4),ism,levels,nc
      write( i2, '(a)' ) 'nmg  ns(1) ns(2)'
      write( i2, '(3i5)' ) nmg, ns(1), ns(2)
      write( i2, '(a)' )
     $       '   aa        bb        cc        ss        f0       f1'
      write( i2, '(6E10.3)' ) aa, bb, cc, ss, f0, f1
      write( i2, '(a)' ) '   om        spar      delta     tol   '
      write( i2, '(4E10.3)' ) om, spar, delta, tol
      write( i2, '(a,a)' ) '  cycle = ', cycle
      if( levels .gt. 10 ) then
         call err( i1, i2, 1 )
      endif
c     end of input
      end
      subroutine matout( i2, idim, k, np, a )
c **************************************************************************
c
c     prints matrix a  on grid k
c
c ****************************************************************************
c
c                  INPUT  PARAMETERS
c
      implicit none
      integer idim, i2, k, np(11)
      double precision a(idim,3)
c ***************************************************************************
c
c                 LOCAL QUANTITIES
c
      integer i, n, m
c **************************************************************************
c
c            SUBROUTINES CALLED
c
c                none
c ***************************************************************************
c
c                 I/O
c    output goes to logical unit number i2
c
c ***************************************************************************
      write( i2, * ) ' matrix at grid', k
      n = np(k+1)-np(k)
      m = n/4
      do 10 i= np(k),np(k+1)-1
         write( i2, '(i3,3e11.3)' ) i-np(k)+1, a(i,1), a(i,2), a(i,3)
  10  continue
c     end of matout
      end
      subroutine matrhs(ib, icase, idim, igal, iout, i2, levels, nc,
     $                  np, aa, amat, bb, cc, f0, f1, rhs, ss )
c **************************************************************************
c
c  Discretization of 1-dimensional boundary value problem on finest or on
c   all grids for application of multigrid methods
c   domain: 0 < x < 1
c   -d/dx(adu/dx) + d/dx(bu) + cu = s
c
c **************************************************************************
c
c                  INPUT / OUTPUT PARAMETERS
c
      implicit none
      integer ib(10), icase, idim, igal, iout, i1, i2, levels, nc,
     $        np(11)
      double precision aa, amat(idim,3), bb, cc, f0, f1, rhs(idim), ss
c
c  aa      i  governs the diffusion coefficient in the diff. equation
c  amat    o  discretization matrix
c  bb      i  governs the velocity coefficient in the diff. equation
c  cc      i  governs the third coefficient in the diff. equation
c  f0      i  value for dirichlet or neumann boundary condition at x = 0
c  f1      i  value for dirichlet or neumann boundary condition at x = 1
c  s       o  right-hand-side
c  a       i  diffusion constant (function)
c  ss      i  governs the right-hand-side in the diff. equation
c  ib      i  governs type of boundary conditions and discretzation of
c             convection term and choice between vertex-centered and
c             cell-centered discretization
c             ib(1) = 1: Dirichlet at x = 0
c             ib(1) = 2: Neumann   at x = 0
c             ib(2) = 1: Dirichlet at x = 1
c             ib(2) = 2: Neumann   at x = 1
c             ib(3) = 1: central
c             ib(3) = 2: upwind
c             ib(3) = 3: hybrid with diffusion term deleted when switched to
c                        upwind
c             ib(3) = 4: hybrid with diffusion term maintained when switched
c                        to upwind
c             ib(4) = 1: vertex-centered
c             ib(4) = 2: cell-centered
c  icase   i  identification number given by user to his problem
c  idim    *  parameter that passes the first dimension of amat to gener
c  igal    i  igal = 0: matrices are made on all grids
c             igal = 1: matrix is made on finest grid only, and coarse
c                       grid matrices are initialized to 0d0
c  iout   i   governs amount of output; iout larger -> more output
c  i2     i   output device number
c  levels i   number of gids in multigrid method; < 16
c  nc     i   number of grid points on coarsest grid
c  np     i   array of pointers to start of grids and to 1+end of finest grid
c             np(1)=1     : start of coarsest grid
c             np(levels)  : start of finest grid
c             np(levels+1): end of finest grid + 1
c  amat   o   matrix (discretization of differential equation) on finest
c             grid  and initialized to 0.d0 on the other grids, or matrices
c             on all grids, depending on igal
c  f0     i   boundary value at x = 0
c  f1     i   boundary value at x = 1
c  rhs    o   right-hand-side on finest grid
c
c ***************************************************************************
c
c                 LOCAL QUANTITIES
c
      integer i, j
      double precision scal
c
c  i,j     general loop variables or scratch variables
c  scal    see subroutine gener
c
c ***************************************************************************
c
c                 I/O
c
c  output file: discr.res  device number = i2   free format
c
c **************************************************************************
c
c            SUBROUTINES CALLED
c
c  gener(i,ib,idim,icase,np,amat,aa,bb,cc,f0,f1)  generates matrix on grid i
c
c  info(icase,ib)  user-provided subroutine that prints information about
c                   his case; body may be empty
c
c  matout(i2,idim,levels,np,amat)  prints amat for inspection
c
c  err(i)         gives error messages
c
c ***************************************************************************
c
c                METHOD
c
c  See subroutine gener
c
c ************************************************************************
      scal = 2d0
      if ( igal .eq. 1 ) then
c
c Galerkin: generate finest grid matrix
c
         call gener( levels, ib, idim, icase, np, amat, aa, bb, cc,
     $               f0, f1, scal, rhs, ss )
c
c Initialize amat for coarse grids
c
         do 11 j = 1, 3
            do 10 i = 1, np(levels) - 1
               amat(i,j) = 0d0
   10       continue
   11    continue
      else if ( igal .eq. 0 ) then
c
c Generate grids for all levels
c
         do 12 i = 1, levels
            call gener( i, ib, idim, icase, np, amat, aa, bb, cc,
     $                   f0, f1, scal, rhs, ss )
   12    continue
      else
c
c Wrong value of igal
c
         call err( i1, i2, 5 )
      endif
      if ( iout .gt. 2 ) then
         do 13 i = 1, levels
            call matout( i2, idim, i, np, amat )
  13     continue
      endif
      if ( iout .gt. 2 ) then
         write(i2,*)'Right hand side'
         call vecout( i2, idim, levels, np, rhs )
      endif
c     end of matrhs
      end
      subroutine mg( amat, b, cycle, delta, histry, ib, idim, ihis,
     $               iout, it, ism, i2, levels, nmg, np, ns, om, r, s,
     $               tol, u, ue, ut, v, wu )
c
c **************************************************************************
c
c  Nonlinear multigrid algorithm as described in section 8.7 of the book.
c  The problem to be solved is assumed to be linear.
c
c **************************************************************************
c
c                  INPUT / OUTPUT PARAMETERS
c
      implicit none
      integer histry(2,1000), ib(10), idim, iout(4), ism, it, i2,
     $        levels, np(11), nmg, ns(2)
      double precision amat(idim,3), b(idim), delta, om, r(idim), s(10),
     $                 tol, u(idim), ue(idim), ut(idim), v(idim), wu(10)
      character cycle
c
c  amat     i   matrices on all grids
c  b        i   right-hand-sides on grid(levels)
c               Remark. The right-hand-sides on the coarse grids, called
c                     f in subroutine MG in section 8.7 in the book, are
c                     also stored in b
c  cycle    i   governs type of multigrid schedule
c                 cycle = 'A' : adaptive schedule
c                 cycle = 'F' : F-cycle
c                 cycle = 'V' : V-cycle
c                 cycle = 'W' : W-cycle
c  delta    i   used in setting tolerance eps in adaptive cycle according
c                 to equation (8.3.1) in the book
c  histry  i,o gives sequence of grids visited and number of smoothings
c                carried out during each visit
c                histry(1,i) : grid number
c                histry(2,i) : number of smoothings carried out at
c                  visit to grid(histry(1,i))
c  ib       i   ib(5) governs updating of ut -- see information on ut
c               ib(6) chooses type of prolongation
c               ib(7) chooses type of restriction for residual
c               ib(8) chooses type of restriction for ut
c  idim     i   size of matrix and grid vectors
c  iout     i   governs output
c                 iout(2) governs diagnostic output and output in
c                    subroutine smooth
c                 iout(3) governs performance output
c  ism      i   chooses smoothing method
c  it       o   gives current multigrid iteration number
c  i2       i   logical unit number
c  levels   i   total number of grids
c                 1 : coarsest grid;  levels: finest grid
c  nmg      i   maximum number of iterations
c  np       i   array of pointers to start of grids and to 1+end of finest
c                 grid
c                 np(1)=1     : start of coarsest grid
c                 np(levels)  : start of finest grid
c                 np(levels+1): end of finest grid + 1
c  ns       i   ns(1):  aantal pre-smoothing iteraties
c               ns(2):  number of post-smoothing iterations
c  om       i   underrelaxation parameter in smoothing method
c  r        *   residual b-Au (local)
c  s        i   underrelaxation parameters in the nonlinear multigrid
c                 algorithm. See algorithm TG in section 8.2 of the book.
c  tol      i   accuracy requirement:     |Au-b)| < tol*|b|
c  u        o   new approximation of solution
c  ue       i   exact discrete solution
c  ut       i   u~ in nonlinear multigrid algorithm.
c                 On coarse grids input is optional.
c                 It may or may not be updated during multigrid iteration
c                 depending on ib(5) = 1 : ut is updated else not.
c                 For linear problems the treatment of ut makes no difference.
c                 If ut is input on coarse grids it may
c                 have been obtained from nested iteration
c  v        i   initial approximation of solution.
c  wu       i,o wu(k) = work units expended on grid k
c
c ***************************************************************************
c
c                 LOCAL QUANTITIES
c
      logical goon, finer, test
      integer i, igamma, ihis, n(10), k
      double precision aa, eps(10), t(10), res0, res1, res2, rfac, anorm
c
c  ****       Remark. The right-hand-sides on the coarse grids, called
c                     f in subroutine MG in section 8.7 in the book, are
c                     stored in b
c  finer      true:   go to finer grid
c             false:  go to coarser grid
c  goon       if (goon) then continue multigrid iteration else terminate
c  test       governs output
c  aa         scratch space
c  i          general loop variable or scratch variable
c  igamma     maximum number of coarse grid corrections
c  ihis       counter for histry array
c  k          general loop variable or scratch variable
c  n          counter of coarse grid iterations
c  eps        bound for residual norm to be reached on finest grid
c               and bound for residual norm to be reached on coarse grids
c               with adaptive schedule
c  t          t(k) < 0 implies coarse grid convergence within tolerance
c                in adaptive schedule
c  res0      l2-norm of initial residual on finest grid
c  res1      l2-norm of old residual on finest grid
c  res2      l2-norm of new residual on finest grid
c  rfac      reduction factor of l2-norm of residual
c
c ***************************************************************************
c
c                 I/O
c
c  output goes to logical unit number i2
c
c ***************************************************************************
c
c            FUNCTIONS
c
c  anorm
c
c **************************************************************************
c
c            SUBROUTINES CALLED
c
c  axplsy copy diagn dirsol errmg prol resid restri scal smooth
c
c ***************************************************************************
c
c                METHOD
c
c  Nonlinear multigrid algorithm based on structure diagrams given in
c  section 8.7 of the book.
c  The structure diagram has been modified slightly in order to obtain a
c  correct F-cycle.
c  An additional termination criterion has been built in, namely that the
c    residual on the finest grid satisfies  |Au-b)| < tol*|b|.
c  Additional statements for output, testing etc. that are not included in the
c    structure diagram are bracketed by c1******* and c2*******.
c  The problem to be solved is assumed to be linear.
c  The rudiments of CPU time measurements are included and are commented out
c  by ctime
c ***************************************************************************
c
c         STRUCTURE DIAGRAM OF NONLINEAR MULTIGRID ALGORITHM
c         Implemented in subroutine mg
c         The while statement has been implemented in the standard way with
c           a goto
c
c #==========================================================================#
c #  Choose ut(k), k = 1,2,...K            /* K is index of finest grid */   #
c #  Initialize v(K) = ut(K)               /* K is called levels        */   #
c #__________________________________________________________________________#
c # \                                                                       /#
c # F \                          cycle eq A                               / T#
c #_____\_______________________________________________________________/____#
c # \                                /|   gamma = 7                          #
c # F \       cycle eq V           / T|   eps(K) = tol*|b|                   #
c #_____\________________________/____|   r(K) = A*v -b                      #
c #  gamma = 2      |  gamma = 1      |   t(K) = |r| - eps                   #
c #_________________|_________________|______________________________________#
c #  k = K                                                                   #
c #  n(K) = nmg                                                              #
c #__________________________________________________________________________#
c # \                                                                       /#
c # F \                      cycle eq A                                   / T#
c #_____\_______________________________________________________________/____#
c #  goon = n(K) ge 0                 |    goon = n(K) ge 0 and t(K) gt 0    #
c #___________________________________|______________________________________#
c #  while (goon) do                                                         #
c #  ________________________________________________________________________#
c #  |  finer = n(k) eq 0  or  k eq 1                                        #
c #  |-----------------------------------------------------------------------#
c #  | \                                                                   / #
c #  | F \                      cycle eq A                               / T #
c #  |_____\___________________________________________________________/_____#
c #  |                                |     finer = finer  or  t(K) lt 0     #
c #  |-----------------------------------------------------------------------#
c #  | \                                                                   / #
c #  | F \                        finer                                  / T #
c #  |_____\___________________________________________________________/_____#
c #  |  S ( v ,u, b )                   |\                                  /#
c #  | /*  S: smoothing          */     |F \            k eq 1            / T#
c #  | /*  v:old; u:new; b:rhs   */     |____\__________________________/____#
c #  |  r = b - A*u                     |    |         direct solution       #
c #  |  Choose ut                       |    |-------------------------------#
c #  |  r(k-1) = Rr                     |    | \                           / #
c #  |  b(k-1) = A*ut + s*r(k-1)        |    | F \      cycle eq F       / T #
c #  |  v(k-1) = ut                     |    |_____\___________________/_____#
c #  |----------------------------------|    |      |    gamma = 1           #
c #  | \                              / |__ _|______|________________________#
c #  | F \       cycle eq A         / T | \                                / #
c #  |___|________________________/_____| F \           k lt K           / T #
c #  |   | t(k-1)=|r(k)|-eps(k)         |_____\________________________/_____#
c #  |   | eps(k-1)=delta*s(k-1)*|r(k)| |     |                              #
c #  |___|______________________________|     |         k =k + 1             #
c #  |         k = k - 1                |     |         u=u+P(u-ut)/s        #
c #  |         n(k) = gamma             |     |         S ( u, v, b )        #
c #  |                                  |_____|______________________________#
c #  |                                  |   n(k) = n(k) - 1                  #
c #  |                                  |   u = v                            #
c #  |                                  |------------------------------------#
c #  |                                  | \                                / #
c #  |                                  | F \       cycle eq A           / T #
c #  |                                  |_____\________________________/_____#
c #  |                                  |      |   r = b - Au                #
c #  |                                  |      |   t = |r| - eps             #
c #  |                                  |------------------------------------#
c #  |                                  | \                                / #
c #  |                                  | F \    cycle eq F and k eq K   / T #
c #  |                                  |_____\________________________/_____#
c #  |                                  |      |        gamma = 2            #
c #  |__________________________________|______|_____________________________#
c #  | \                                                                   / #
c #  | F \                       cycle eq A                              / T #
c #  |_____\___________________________________________________________/_____#
c #  |  goon = n(K) ge 0                |  goon = t(K) gt 0  and  n(K) ge 0  #
c #==========================================================================#
c
c ***************************************************************************
c
c Reset the time
c
ctime      t0 = time
c
      test = .false.
      do 2 i = 1, 10
         s(i) = 1d0
    2 continue
      if ( cycle .eq. 'A' ) then
c
c Adaptive schedule : max number of coarse grid corrections allowed is 7
c
         igamma = 7
         eps(levels) = tol * anorm( 2, idim, b, levels, np )
         call resid( idim, levels, np, b, amat, v, r )
         t(levels) = anorm( 2, idim, r, levels, np ) - eps(levels)
      else if ( cycle .eq. 'V' ) then
         igamma = 1
      else if ( cycle .eq. 'W' .or. cycle .eq. 'F' ) then
         igamma = 2
      else
         call errmg ( i2, 1 )
      endif
      k = levels
      n(levels) = nmg

c1**************************
      call resid( idim, levels, np, b, amat, v, r )
      res0 = anorm( 2, idim, r, levels, np )
      res1 = res0
      res2 = res1
      write( i2, '(a,e11.3)' ) ' l2-norm of initial residual',res1
c2**************************

      if ( cycle .eq. 'A' ) then
         goon = ( n(levels) .ge. 0 .and.  t(levels) .gt. 0d0 )
      else
         goon = ( n(levels) .ge. 0 )
      endif

   10 if ( goon ) then
         finer = ( n(k) .eq. 0 .or. k .eq. 1 )
         if ( cycle .eq. 'A' ) then
            finer = ( finer .or. t(k) .le. 0d0 )
         endif
         if ( finer ) then

c1************
            if (test) then
             write(i2,*) ' GOING TO A FINER GRID (from k = ',k,' )'
            endif
c2***********

            if ( k .eq. 1 ) then

c1**************************
               if ( iout(2) .gt. 1 ) then
c
c norm of right-hand-side on coarsest grid
c
                  call diagn( 1, idim, i2, iout(2), k, np, b )
               endif
c2**************************

               call dirsol( idim, i2, k, np, amat, u, b )

c1*****************
cc possibility of non-acurate solving on coarse grid instead of solving
cc with dirsol
cc
cc             call smooth( i2, idim, iout(2), ism, k, np, ns(2),
cc   $                         amat, levels, om, b, v, u )
cc
               if (test) then
                write(i2,*)' k=',k,' solution  A'
                call vecout(i2,idim,k,np,u)
               endif
c2*****************
               if ( cycle .eq. 'F' ) then
                  igamma = 1
               endif
c1**************************
               if ( ihis .lt. 1000 ) then
                  ihis = ihis + 1
                  histry(1,ihis) = k
                  histry(2,ihis) = 1
c                  write( i2, * ) ' history: level = ',histry(1,ihis),
c     $                           ' smoothings = ',histry(2,ihis)
c               else
c                  write( i2, * ) ' array histry too small'
               endif
               if ( iout(2) .gt. 1 ) then
c
c norm of residual after solution on coarsest grid
c
                  call resid( idim, k, np, b, amat, u, r )
                  call diagn( 2, idim, i2, iout(2), k, np, r )
               endif
               wu(k) = wu(k) + 2d0**(k-levels)
ct1
               if (test) then
                call resid( idim, k, np, b, amat, u, r )
                write(i2,*)' k=',k,' residual  B'
                call vecout(i2,idim,k,np,r)
               endif
ct2
c2**************************
            endif
            if ( k .lt. levels ) then
c
c go to finer grid
c
c1**************************
               if ( iout(2) .gt. 3 ) then
c
c current solution just before move to finer grid
c
                  call diagn( 3, idim, i2, iout(2), k, np, u )
               endif
               if (test) then
                write(i2,*)' k=',k,' solution before move to finer  C'
                call vecout(i2,idim,k,np,u)
               endif
c2****************************
               k = k + 1
c
c start of part B of multigrid algorithm
c u = u -ut
c
               call axplsy( idim, k-1, np, -1d0, ut, u )
c1**************************
               if ( iout(2) .gt. 3 ) then
c
c correction just before move to finer grid
c
                  call diagn( 4, idim, i2, iout(2), k-1, np, u )
               endif
               if (test) then
                write(i2,*)' k=',k,' correction before move to finer  D'
                call vecout(i2,idim,k-1,np,u)
               endif
c2**************************
c
c  v = Pu
c
               call prol( amat, ib(6), idim, i2, k, np, u, v )
c1**************************
               if ( iout(2) .gt. 3 ) then
c
c prolongated correction just after move to finer grid
c
                  call diagn( 5, idim, i2, iout(2), k, np, v )
               endif
               if (test) then
                write(i2,*)' k=',k,' prolongated correction  E'
                call vecout(i2,idim,k,np,v)
               endif
c2**************************
c
c  u = v/s + u
c
               call axplsy( idim, k, np, 1d0/s(k-1), v, u )
c1**************************
               if ( iout(2) .gt. 1 ) then
c
c residual after coarse grid correction
c
                  call resid( idim, k, np, b, amat, u, r )
                  call diagn( 6, idim, i2, iout(2), k, np, r )
               endif
               if ( iout(2) .gt. 3 ) then
c
c current solution after coarse grid correction
c
                  call diagn( 7, idim, i2, iout(2), k, np, u )
               endif
               if (test) then
                write(i2,*)' k=',k,' current sol. after cg corr.  F'
                call vecout(i2,idim,k,np,u)
               endif
c2**************************
c
c now follows post-smoothing
c
ctime          t1 = time
c
c  u is old, v is new approximation
c
               call smooth( i2, idim, iout(2), ism, k, np, ns(2),
     $                         amat, levels, om, b, u, v)
c1******************************
               if (test) then
                write(i2,*)' k=',k,' sol. after post smoothing  G'
                call vecout(i2,idim,k,np,v)
               endif
ctime          t2 = time
ctime          t2 = (t2-t1)/ns(2)
ctime          write( i2, e11.3 ) ' CPU time for 1 WU =',t2
               if ( ihis .lt. 1000 ) then
                  ihis = ihis + 1
                  histry(1,ihis) = k
                  histry(2,ihis) = ns(2)
c                  write( i2, * ) ' history: level = ',histry(1,ihis),
c     $                           ' smoothings = ',histry(2,ihis)
c               else
c                  write( i2, * ) ' array histry too small'
               endif
               if ( iout(2) .gt. 1 ) then
c
c residual after post-smoothing
c
                  call resid( idim, k, np, b, amat, v, r )
                  call diagn( 8, idim, i2, iout(2), k, np, r )
               endif
               if (test) then
                write(i2,*)' k=',k,' residual after post smoothing  H'
                call resid( idim, k, np, b, amat, v, r )
                call vecout(i2,idim,k,np,r)
               endif
               wu(k) = wu(k) + ns(2)*(2d0**(k-levels))
c2**************************
c
c end of part B of multigrid algorithm
c
            endif
            n(k) = n(k) - 1
c
c  u = v
c
            call copy( idim, k, np, u, v )
c1**************************
c            if ( ib(5) .eq. 1 ) then
c
c ut = u
c
c               call copy( idim, k, np, ut, u )
c            endif
c2******************************
            if ( cycle .eq. 'A' ) then
               call resid( idim, k, np, b, amat, u, r )
               t(k) = anorm( 2, idim, r, k, np ) - eps(k)
            endif
            if (k.eq.levels .and. cycle.eq.'F') then
               igamma = 2
            endif
         else
c
c go to coarser grid
c
c1*************************
            if (test) then
              write(i2,*) ' GOING TO A COARSER GRID (from k = ',k,' )'
            endif
            if ( iout(2) .gt. 1 ) then
c
c              residual before pre-smoothing
c
               call resid( idim, k, np, b, amat, v, r )
               call diagn( 9, idim, i2, iout(2), k, np, r )
            endif
c2**************************
c
c now follows pre-smoothing
c
c1****************************
               if (test) then
              write(i2,*) ' grid number',k,' v before presmoothing  I'
                 call vecout( i2, idim, k, np, v )
                 call resid( idim, k, np, b, amat, v, r )
                 write(i2,*)' k=',k,' residual before presmoothing  J'
                 call vecout(i2,idim,k,np,r)
               endif
c2******************************
            call smooth( i2, idim, iout(2), ism, k, np, ns(1),
     $                      amat, levels, om, b, v, u )
c1******************************
            if (test) then
             write(i2,*) ' solution after presmoothing            K'
             call vecout( i2, idim, k, np, u )
            endif
c2*******************************
            call resid( idim, k, np, b, amat, u, r )
c1********************************
            if (test) then
             write(i2,*) ' residual after presmoothing            L'
             call vecout( i2, idim, k, np, r )
            endif
            if ( ihis .lt. 1000 ) then
               ihis = ihis + 1
               histry(1,ihis) = k
               histry(2,ihis) = ns(1)
c               write( i2, * ) ' history: level = ',histry(1,ihis),
c     $                        ' smoothings = ',histry(2,ihis)
c            else
c               write( i2, * ) ' array histry too small'
            endif
            if ( iout(2) .gt. 1 ) then
c
c print residual after pre-smoothing
c
               call diagn( 10, idim, i2, iout(2), k, np, r )
            endif
            wu(k) = wu(k) + ns(1)*(2d0**(k-levels))
c2**************************

            if ( ib(5) .eq. 1 ) then
c
c ut = restriction(u)
c
               call restri( ib(8), idim, i2, k, np, ut, u )
            endif
            call restri( ib(7), idim, i2, k, np, r, r )

c1**************************
            if (test) then
            write(i2,*) ' restriction of residual after presmoothing M'
             call vecout( i2, idim, k-1, np, r )
            endif
c2**************************

c
c  r = -s*r
c
            call scal( idim, k-1, np, -s(k-1), r )
c
c  b = r -amat*ut
c
            call resid( idim, k-1, np, r, amat, ut, b )
c
c  b = -b
c
            call scal( idim, k-1, np, -1d0, b )
            call copy( idim, k-1, np, v, ut)

c1****************************
            if (test) then
             write(i2,*) ' righthandside for coarser grid          N'
             call vecout( i2, idim, k-1, np, b )
             write(i2,*) ' end of part A of multigrid algorithm'
            endif
c2*****************************

c
c end of part A of multigrid algorithm
c
            if ( cycle .eq. 'A' ) then
               aa = anorm( 2, idim, r, k, np )
               t(k-1) = aa - eps(k)
               eps(k-1) = delta*s(k-1)*aa
            endif
            k = k - 1
            n(k) = igamma
         endif
         if ( cycle .eq. 'A' ) then
            goon = ( t(levels) .gt. 0d0 .and. n(levels). ge. 0 )
         else
            goon = (n(levels) .ge. 0 )
         endif

c1**************************
         if (k.eq.levels) then
               call copy(idim,k,np,ut,u)
               if (test) then
                  write(i2,*) ' new ut'
                  call vecout(i2,idim,k,np,ut)
               endif
               if ( iout(3) .gt. -1 ) then
                  call resid( idim, levels, np, b, amat, u, r )
                  res1 = res2
                  res2 = anorm( 2, idim, r, levels, np )
                  rfac = max(1d-20,res2)/max(1d-20,res1)
                  it = nmg - n(levels)
                  if (it .le. nmg) then
                    write( i2, '(a,i2,a,e11.3)' ) ' iteration number ',
     $                 it, ' reduction factor ', rfac
                    write( i2,'(a,e11.3)')' l2-norm of residual ',res2
                  endif
               endif
               if ( iout(2) .gt. 1 ) then
c
c error at moment of leaving finest grid
c
                  call axplsy( idim, k, np, -1d0 , ue, v )
                  call diagn( 11, idim, i2, iout(2), k, np, v )
                  call copy( idim, k, np, v, u )
               endif
         endif
c2**************************

         goto 10
      endif

c1**************************
      call resid( idim, levels, np, b, amat, u, r )
      res2 = anorm( 2, idim, r, levels, np )
      write(i2,'(a,e11.3)')' l2-norm of res. when leaving mg ',res2
      it = nmg - n(levels)-1
      rfac = (max(1d-12,res2)/max(1d-12,res0))**(1d0/it)
      write(i2,'(a,i2,a,e11.3)')' total iterations ',it,
     $                     ' average reduction factor ', rfac
      if (iout(2) .gt. 0) then
        write(i2,*)' Schedule history'
        do 20 i =1, ihis
          write(i2,'(2i3)') histry(1,i),histry(2,i)
  20    continue
      endif
      if ( iout(2) .gt. 1 ) then
c
c final error
c
         call copy( idim, k, np, v, u )
         call axplsy( idim, k, np, -1d0, ue, v )
         call diagn( 12, idim, i2, iout(2), k, np, v )
      endif
ctime t2 = time
ctime t2 = t2 - t0
ctime write(i2,'(a,e11.3)')' total CPU time = ',t2
c2**************************
c     end of mg
      end

      subroutine nested( amat, b, cycle, delta, histry, ib, idim, ihis,
     $                   iout,  it, ism, i2, levels, nmg, np, ns, om, r,
     $                   s, tol, u, ue, ut, v, wu )
c **************************************************************************
c
c  Nested iteration algorithm as described in section 8.4 of the book.
c
c **************************************************************************
c
c                  INPUT / OUTPUT PARAMETERS
c
      implicit none
      integer histry(2,1000), ib(13), idim, ihis, iout(4), ism, it, i2,
     $        levels, np(11), nmg, ns(2)
      double precision amat(idim,3), b(idim), delta, om, r(idim), s(10),
     $                 tol, u(idim), ue(idim), ut(idim), v(idim), wu(10)
      character cycle
c
c  amat     i   matrices on all grids
c  b        i   right-hand-side on grid(levels)
c           *   the right-hand-sides on coarse grids are computed by
c               restriction
c  cycle    i   passed to subroutine mg
c  delta    i   passed to subroutine mg
c  histry  o   output of subroutine mg
c  ib       i   passed to subroutine mg
c               ib(11) number of multigrid iterations used in nested
c                 iteration (=gamma~ in program 1 of section 8.4)
c               ib(12) chooses type of prolongation in nested iteration
c               ib(13) chooses type of restriction for coarse right-hand-
c                 sides in nested iteration
c  idim     i   size of matrix and grid vectors
c  ihis     i,o counter for histry array
c  iout     i   governs output
c                 iout(2) governs diagnostic output and output in
c                    subroutine smooth
c                 iout(3) governs performance output
c  ism      i   passed to subroutine mg
c  it       o   passed to subroutine mg
c  i2       i   logical unit number
c  levels   i   total number of grids
c                 1 : coarsest grid;  levels: finest grid
c  np       i   array of pointers to start of grids and to 1+end of finest
c                 grid
c                 np(1)=1     : start of coarsest grid
c                 np(levels)  : start of finest grid
c                 np(levels+1): end of finest grid + 1
c  ns       i   passed to subroutine mg
c  om       i   passed to subroutine mg
c  r        *   passed to subroutine mg
c  s        i   passed to subroutine mg
c  tol      i   passed to subroutine mg
c  u        o   new approximation of solution
c  ue       i   exact discrete solution
c  ut       i,o initial and final approximations on all grids
c  v        *   scratch space (local)
c  wu       i,o   wu(k) = work units expended on grid k
c
c ***************************************************************************
c
c                 LOCAL QUANTITIES
c
      integer  k
c
c ***************************************************************************
c
c                 I/O
c
c  output goes to logical unit number i2
c
c ***************************************************************************
c
c            FUNCTIONS none
c
c **************************************************************************
c
c            SUBROUTINES CALLED
c
c   copy  mg  prol  restri  smooth
c
c ***************************************************************************
c
c                METHOD
c
c  Nested iteration according to program 1 of section 8.4
c  Additional statements for output, testing etc. that are not included in the
c    description given in the book are bracketed by c1******* and c2*******.
c  The problem to be solved is assumed to be linear.
c  The rudiments of CPU time measurements are included and are commented out
c  by ctime
c
c ****************************************************************************
c
c Reset the time
c
ctime      t0 = time
c
c compute coarse right-hand-sides
c
      do 10 k = levels, 2, -1
         call restri( ib(13), idim, i2, k, np, b, b )
  10  continue
c
c1**************************
      if ( iout(2) .gt. 1 ) then
c
c norm of right-hand-side on coarsest grid
c
         write(i2,*)' nested iteration, grid 1'
         call diagn( 1, idim, i2, iout(2), 1, np, b )
      endif
c2**************************
c
      call smooth( i2, idim, iout, ism, 1, np, ns, amat, levels,
     $                   om, b, ut, v )
      call copy( idim, 1, np, ut, v )
c
c1**************************
      if ( iout(2) .gt. 1 ) then
c
c norm of residual after solution on coarsest grid
c
         call resid( idim, 1, np, b, amat, ut, r )
         write(i2,*)' nested iteration, grid 1'
         call diagn( 8, idim, i2, iout(2), 1, np, r )
      endif
      k = 1
      wu(k) = wu(k) + 2d0**(k-levels)
      if ( ihis .lt. 1000) then
         ihis=ihis+1
         histry(1,ihis)=k
         histry(2,ihis)=1
      else
         write(i2,*)' array histry too small'
      endif
c2**************************
c
      do 20 k = 2, levels
         call prol( amat, ib(12), idim, i2, k, np, ut, ut )
c
c1**************************
         if ( iout(2) .gt. 1 ) then
c
c norm of residual before multigrid improvement
c
            call resid( idim, k, np, b, amat, ut, r )
            write(i2,*)
     $          ' nested iteration, grid ',k, 'res. before mg corr.'
            call diagn( 13, idim, i2, iout(2), k, np, r )
         endif
c2**************************
c
         call copy( idim, k, np, v, ut )
         call mg( amat, b, cycle, delta, histry, ib, idim, ihis, iout,
     $               it, ism, i2, k, ib(11), np, ns, om, r, s, tol,
     $               u, ue, ut, v, wu )
         call copy( idim, k, np, ut, u )
c
c1**************************
         if ( iout(2) .gt. 1 ) then
c
c norm of residual after multigrid improvement
c
            call resid( idim, k, np, b, amat, ut, r )
            write(i2,*)
     $         ' nested iteration, grid ',k, 'res. after mg corr.'
            call diagn( 13, idim, i2, iout(2), k, np, r )
         endif
c2**************************
c
   20 continue
c
c end of nested
c
      end
      double precision function omega( p, t )
c***********************************************************************
c
c   Switching function between central, hybrid and upwind scheme
c
c***********************************************************************
c
c       INPUT/OUTPUT PARAMETERS
c
      implicit none
      integer  t
      double precision p
c
c    t      i  t = 1:  central scheme
c              t = 2:  upwind scheme
c              t = 3:  hybrid scheme
c              t = 4:  hybrid scheme
c              t = something else:  central scheme
c    p      i  mesh peclet number
c
c***********************************************************************
c
c       LOCAL QUANTITIES
c
      double precision om
c
c***********************************************************************
      if( t.eq.1 ) then
         om = 0d0
      else if( t.eq.2 ) then
         om = 1d0
      else if( t.eq.3 .or. t.eq.4 ) then
         om = 5d0*p - 9.5d0
         om = min(om,1d0)
         om = max(om,0d0)
      else
         om = 0d0
      endif
      omega = om
c  end of omega
      end

      subroutine pointr( idis, iout, levels, nc, np, i2 )
c***********************************************************************
c
c        computation of pointer array np
c
c***********************************************************************
c
c       INPUT/OUTPUT PARAMETERS
c
      implicit none
      integer idis, iout, levels, nc, np(11), i2
c
c    idis   i  idis = 1:  vertex-centered multigrid
c              idis = 2:  cell-centered multigrid
c    iout   i  governs output
c    levels i  description in main program
c    nc     i  description in main program
c    np     o  description in main program
c
c***********************************************************************
c
c       LOCAL QUANTITIES
c
      integer i, j
c
c***********************************************************************

      np(1) = 1
      j=nc
      if ( idis .eq. 1 ) then
c
c Vertex centered
c
         do 10 i = 2,levels
            np(i) = np(i-1) + j
            j = 2*j - 1
  10     continue
         np(levels+1) = np(levels) + j
      else if ( idis .eq. 2 ) then
c
c Cell centered
c
         do 11 i = 2,levels
            np(i)= np(i-1) + j
            j = 2*j
  11     continue
         np(levels+1) = np(levels) + j
      else
c
c Wrong value of idis = ib(4)
c
         write( i2, * ) ' idis wrong in subroutine pointer'
      endif
c
c Check for output
c
      if ( iout .gt. 2 ) then
c
c Output required
c
         write( i2, * ) ' grid pointer array'
         write( i2, '(2i5)' ) ( i,np(i), i=1,levels)
      endif
c
c end of pointer
c
      end
      subroutine prol( amat, ic, idim, i2, k, np, u, v )
c **************************************************************************
c
c     computes prolongation of u from grid k-1 to grid k
c       v = Pu
c
c ****************************************************************************
c
c                  INPUT/OUTPUT PARAMETERS
c
      implicit none
      integer ic, idim, i2, k, np(11)
      double precision amat(idim,3), u(idim), v(idim)

c
c     ic   i  governs choice of type of prolongation; see below
c     u    i  function on grid(k-1) to be prolongated
c     v    o  prolongation of u
c
c ***************************************************************************
c
c                 LOCAL QUANTITIES
c
      integer i, ipu, ipv, nu
c ***************************************************************************
c
c                 I/O
c    output goes to logical unit number i2
c
c ***************************************************************************
c
c            FUNCTIONS
c
c             none
c
c **************************************************************************
c
c            SUBROUTINES CALLED
c
c     none;  possible to build this with BLAS routines
c
c ***************************************************************************
c
c            METHOD
c
c   ic = 1:  vertex-centered linear interpolation: equation (5.3.1) of the book
c   ic = 2:  cell-centered piecewise constant interpolation:
c               equation (5.3.5) of the book
c   ic = 3:  cell-centered linear interpolation with Dirichlet boundary
c               modification:  equation (5.3.6) of the book with w=0
c               in equation (5.3.24) at boundaries
c   ic = 4:  cell-centered linear interpolation with Neumann boundary
c               modification:  equation (5.3.6) of the book with stencil
c               [P*] = (1/4)[ 0  4  3  1 ] at boundary x=0 and similar at x=1
c   ic = 5:  vertex-centered operator dependent prolongation
c                  equations (5.4.2), (5.4.3) of the book
c
c ***************************************************************************
      ipv = np(k) - 1
      ipu = np(k-1) - 1
      nu = ipv - ipu
c
c many of the following loops can be done with the BLAS routines DCOPY,
c DAXPY and DSCALE. So the loops are splitted into simple parts.
c
      if ( ic .eq. 1 ) then
c
c vertex centered linear interpolation:
c
c        v: 1   2   3   4   5   6   7   8   9
c           | /   \ | /   \ | /   \ | /   \ |  [0 2 1] [1 2 1] [1 2 0]  *.5
c        u: 1       2       3       4       5
c
         do 10 i = 1, nu
            v(ipv+2*i-1) = u(ipu+i)
  10     continue
         do 12 i = 1, nu
            v(ipv+2*i) = 0.5d0*u(ipu+i)
  12     continue
         do 14 i = 1, nu
            v(ipv+2*i) = v(ipv+2*i) + 0.5d0*u(ipu+i+1)
  14     continue
      else if ( ic .eq. 2 ) then
c
c cell centered piecewise constant interpolation:
c
c        v: 1   2   3   4   5   6   7   8   9
c            \ /     \ /     \ /     \ /
c        u:   1       2       3       4
c
         do 20 i = 1, nu
            v(ipv+2*i-1) = u(ipu+i)
   20    continue
         do 22 i = 1, nu
            v(ipv+2*i) = u(ipu+i)
   22    continue
      else if ( ic .eq. 3 .or. ic .eq. 4 ) then
c
c cell centered linear interpolation Dirichlet/Neumann boundary
c
c        v: 1   2   3   4   5   6
c              ____/ \ / \____        [ 0 2 3 1 ] [ 1 3 3 1 ] [ 1 3 2 0 ] *.25
c        u:   1       2       3
c
         do 30 i = 1, 2*nu
            v(ipv+i) = 0d0
   30    continue
         do 32 i = 1, nu
            v(ipv+2*i-1) = v(ipv+2*i-1) + .75d0*u(ipu+i)
   32    continue
         do 34 i = 1, nu
            v(ipv+2*i) = v(ipv+2*i) + .75d0*u(ipu+i)
   34    continue
         do 36 i = 2, nu
            v(ipv+2*i-2) = v(ipv+2*i-2) + .25d0*u(ipu+i)
   36    continue
         do 38 i = 1, nu - 1
            v(ipv+2*i+1) = v(ipv+2*i+1) + .25d0*u(ipu+i)
   38    continue
         if ( ic .eq. 3 ) then
c
c Dirichlet boundary correction
c
            v(ipv+1) = 0.5*u(ipv+1)
            v(ipv+2*nu) = 0.5*u(ipv+nu)
         else
c
c Neumann boundary correction
c
            v(ipv+1) = u(ipv+1)
            v(ipv+2*nu) = u(ipv+nu)
         endif
      else
         write(i2,'(a,i2,a)') 'prol: ic = ',ic,' out of range '
         write(*,*) ' aborted'
         stop
      endif
c
c     end of prol
c
      end
      subroutine resid( idim, k, np, b, amat, u, r )
c **************************************************************************
c
c     computes residual r = b - amat*u  on grid k
c
c ****************************************************************************
c
c                  INPUT  PARAMETERS
c
      implicit none
      integer idim, k, np(11)
      double precision b(idim), amat(idim,3), u(idim), r(idim)
c ***************************************************************************
c
c                 LOCAL QUANTITIES
c
      integer i
c ***************************************************************************
c
c                 I/O
c                none
c
c ***************************************************************************
c
c            FUNCTIONS
c
c             none
c
c **************************************************************************
c
c            SUBROUTINES CALLED
c
c   BLAS level 2 routine _GBMV for vector computers; has not yet been
c    implemented
c
c ***************************************************************************
      i=np(k)
      r(i) = b(i) - amat(i,2)*u(i) - amat(i,3)*u(i+1)
      do 10 i = np(k)+1, np(k+1)-2
         r(i) = b(i) - amat(i,1)*u(i-1) - amat(i,2)*u(i) -
     $                 amat(i,3)*u(i+1)
  10  continue
      i = np(k+1)-1
      r(i) = b(i) - amat(i,1)*u(i-1) - amat(i,2)*u(i)
c     end of resid
      end
      subroutine restri( ic, idim, i2, k, np, v, u )
c***********************************************************************
c
c     computes restriction of u from grid k to grid k-1
c       v = Ru
c
c ****************************************************************************
c
c                  INPUT/OUTPUT  PARAMETERS
c
      implicit none
      integer ic, idim, i2, k, np(11)
      double precision u(idim), v(idim)
c
c  ic   i    chooses type of restriction
c            = 1:  vertex-centered restriction according to the one-
c                    dimensional version of equation (5.3.19) of the book
c                    stencil = 0.5 * [1 2 1]
c                    boundary modification type Dirichlet: at x=0:
c                    stencil = 0.25 * [0 2 1];
c                    at x=1 :stencil = 0.25 * [1 2 0]
c            = 2:  vertex-centered restriction according to the one-
c                    dimensional version of equation (5.3.19) of the book
c                    stencil = 0.5 * [1 2 1]
c                    boundary modification type Neumann: at x=0:
c                    stencil = 0.5 * [0 3 1];
c                    at x=1 :stencil = 0.5 * [1 3 0]
c            = 3:  cell-centered restriction according to equation (5.3.23)
c                    of the book
c                    stencil =  [1  1]
c            = 4:  cell-centered restriction according to equation (5.3.24)
c                    of the book
c                    stencil = (1/4)*[1 3 3 1]
c                    boundary modification type Dirichlet: at x=0:
c                    stencil = (1/4)*[0 3 3 1];
c                    at x=1 :stencil = 0.5*[1 3 3 0]
c            = 5:  cell-centered restriction according to equation (5.3.24)
c                    of the book
c                    stencil = (1/4)*[1 3 3 1]
c                    boundary modification type Neumann: at x=0:
c                    stencil = (1/4) * [0 4 3 1];
c                    at x=1 :stencil = 0.5 * [1 3 4 0]
c            = 6:  vertex-centered operator-dependent restriction according
c                    to equation (5.4.22) of the book
c                    stencil = [a1 a2 a3]
c                    boundary modification type Neumann: at x=0:
c                    a1 = 0, a2:= 2*a2-a3
c                    at x=1 : a3=0, a2:=2*a2-a1
c            = 7:  vertex-centered operator-dependent restriction according
c                    to equation (5.4.22) of the book
c                    stencil = [a1 a2 a3]
c                    boundary modification type Dirichlet: at x=0: a1 = 0
c                    at x=1 : a3 = 0
c
c ***************************************************************************
c
c                 LOCAL QUANTITIES
c
      integer i, ipu, ipv, nv
c ***************************************************************************
c
c                 I/O
c    output goes to logical unit number i2
c
c ***************************************************************************
c
c            FUNCTIONS
c
c             none
c
c **************************************************************************
c
c            SUBROUTINES CALLED
c
c     version for vector computers with BLAS routines
c     has been commented out with cv
c
c ***************************************************************************
c
c            METHOD
c
c   See above description of input parameter ic
c   The algorithm has been designed to be easily vectorizable and can be
c     replaced immediately by calls of BLAS routines.
c   Calls to BLAS have been commented out with cv
c
c ***************************************************************************
      ipu = np(k) - 1
      ipv = np(k-1) - 1
      nv = ipu - ipv
c
c much of the following loops can be done with the BLAS routines DCOPY,
c DAXPY and DSCALE. So the loops are splitted into simple parts.
c
      if ( ic .eq. 1 .or. ic.eq.2 ) then
c
c vertex centered restriction Dirichlet/Neumann boundary:
c
c        v: 1     2     3     4     5
c           | \ / | \ / | \ / | \ / |     [0 2 1] [1 2 1] [1 2 0] *.5
c        u: 1  2  3  4  5  6  7  8  9
c
         v(ipv+1) = 0d0
         do 10 i = 2, nv
            v(ipv+i) = 0.5d0*u(ipu+2*i-2)
  10     continue
         do 12 i = 1, nv
            v(ipv+i) = v(ipv+i) + u(ipu+2*i-1)
  12     continue
         do 14 i = 1, nv - 1
            v(ipv+i) = v(ipv+i) + 0.5d0*u(ipu+2*i)
  14     continue
         if ( ic .eq. 1 ) then
c
c vertex centered restriction Dirichlet boundary:
c
            v(ipv+1) = 0.5d0*v(ipv+1)
            v(ipv+nv) = 0.5d0*v(ipv+nv)
         endif
         if ( ic .eq. 2 ) then
c
c vertex centered restriction Neumann boundary:
c
c           v: 1     2     3     4     5
c              | \ / | \ / | \ / | \ / |     [0 3 1] [1 2 1] [1 3 0] *.5
c           u: 1  2  3  4  5  6  7  8  9
c
            v(ipv+1) = v(ipv+1)+0.5*u(ipu+1)
            v(ipv+nv) = v(ipv+nv)+0.5*u(ipu+2*nv-1)
         endif
      else if ( ic .eq. 3 ) then
c
c cell centered restriction piecewise constant
c
c        v:   1       2       3       4
c            / \     / \     / \     / \      [ 1 1 ] [ 1 1 ] [ 1 1 ]
c        u: 1   2   3   4   5   6   7   8   9
c
         do 30  i = 1, nv
            v(ipv+i) = u(ipu+2*i-1)
   30    continue
         do 32 i = 1, nv
            v(ipv+i) = v(ipv+i) + u(ipu+2*i)
   32    continue
      else if ( ic .eq. 4 .or. ic .eq. 5 ) then
c
c cell centered restriction Dirichlet/Neumann boundary
c
c        v:   1   ___ 2 ___   3
c                /   / \   \          [ 0 3 3 1 ] [ 1 3 3 1 ] [ 1 3 3 0 ] *.25
c        u: 1   2   3   4   5   6
c
         v(ipv+1) = 0d0
         do 40 i = 2, nv
            v(ipv+i) = 0.25d0*u(ipu+2*i-2)
  40     continue
         do 42 i = 1, nv
            v(ipv+i) = v(ipv+i) + 0.75d0*u(ipu+2*i-1)
  42     continue
         do 44 i = 1, nv
            v(ipv+i) = v(ipv+i) + 0.75d0*u(ipu+2*i)
  44     continue
         do 46 i = 1, nv - 1
            v(ipv+i) = v(ipv+i) + 0.25d0*u(ipu+2*i+1)
  46     continue
         if ( ic .eq. 5 ) then
c
c cell centered restriction Neumann boundary
c
c           v:   1   ___ 2 ___   3
c                   /   / \   \       [ 0 4 3 1 ] [ 1 3 3 1 ] [ 1 3 4 0 ] *.25
c           u: 1   2   3   4   5   6
c
            v(ipv+1) = v(ipv+1) + 0.25d0*u(ipu+1)
            v(ipv+nv) = v(ipv+nv) + 0.25d0*u(ipu+2*nv)
         endif
      else
         write(i2,'(a,i2,a)')'restri: ic = ',ic,' out of range '
         write(*,*) ' aborted'
         stop
      endif
c
c     end of restri
c
      end
      double precision function s(icase,ss,x)
c***********************************************************************
c
c     user provided function to fill the right hand side as a function
c     of the x-coordinate
c
c***********************************************************************
      integer icase
      double precision ss,x
      if ( icase .eq. 1 ) then
         s = ss
      else if ( icase .eq. 2) then
         s = ss
      else if ( icase .eq. 3 ) then
         s = ss
      else
         s = ss
      endif
      end
      subroutine scal( idim, k, np, a, u )
c***********************************************************************
c
c     u = a*u    on grid k
c
c ****************************************************************************
c
c                  INPUT  PARAMETERS
c
      implicit none
      integer idim, k, np(11)
      double precision a, u(idim)
c ***************************************************************************
c
c                 LOCAL QUANTITIES
c
      integer i
c **************************************************************************
c
c            SUBROUTINES CALLED
c
c   BLAS routine DSCAL may be used on vector computers - has been
c    commented out with cv
c
c ***************************************************************************
c
      do 10 i = np(k), np(k+1)-1
         u(i) = a*u(i)
  10  continue
c
c   vector version
c
cv      call DSCAL(np(k+1)-np(k),a,u(np(k)),1)
c     end of scal
      end
      subroutine smooth( i2, idim, iout, ism, k, np, ns, a, levels,
     $                   om, b, u, v )
c **************************************************************************
c
c   smoothing of equation Au = b  on grid k
c   old approximate solution is u
c   new approximate solution is v
c
c **************************************************************************
c
c                  INPUT / OUTPUT PARAMETERS
c
      implicit none
      integer i2, idim, iout, ism, k, np(11), ns, levels
      double precision a(idim,-1:1), b(idim), om, u(idim), v(idim)
c
c  a        i   matrix
c  b        i   right-hand-side
c  om       i   underrelaxation parameter in smoothing method
c  u        i   old approximation of solution
c  v        o   new approximation of solution
c  i2       i   logical device number
c  idim     i   size of matrix and grid vectors
c  iout     i   governs output
c  ism      i   chooses smoothing method
c               ism = 1  Jacobi with damping parameter om; section 4.3
c                  of the book and equation (4.1.7); damping is not applied
c                  at boundaries
c               ism = 2  forward Gauss-Seidel; section 4.3 of the book
c               ism = 3  backward Gauss-Seidel; section 4.3 of the book
c               ism = 4  symmetric Gauss-Seidel; section 4.3 of the book
c               ism = 5  white-black; section 4.3 of the book
c               ism = 6  Jacobi with damping parameter om; section 4.3
c                  of the book and equation (4.1.7); damping is applied
c                  uniformly in the interior and at the boundaries.
c  k        i   grid level
c  np       i   array of pointers to start of grids and to 1+end of finest
c                 grid
c                 np(1)=1     : start of coarsest grid
c                 np(levels)  : start of finest grid
c                 np(levels+1): end of finest grid + 1
c  ns       i   number of smoothing iterations
c
c ***************************************************************************
c
c                 LOCAL QUANTITIES
c
      integer i, n, ifirst, ilast, itotal
      double precision  om1
      logical odd
c ***************************************************************************
c
c                 I/O
c
c  output goes to logical unit number i2
c
c ***************************************************************************
c
c            FUNCTIONS     none
c
c **************************************************************************
c
c            SUBROUTINES CALLED
c
c                 copy
c
c ***************************************************************************
c
c                METHOD
c
c   see information on ism above
c
c ****************************************************************************
      ifirst = np(k)
      ilast = np(k+1) - 1
      itotal = ilast - ifirst
      odd = .false.
      if ( mod(itotal,2) .eq. 1 ) odd = .true.
      if ( a(ifirst,-1) .ne. 0d0 .or. a(ilast,1) .ne. 0d0 ) then
         write(i2,*)' smooth: matrix a wrong at boudaries'
         write(*,*)' aborted'
         stop
      endif
c
c  first copy the initial solution u to v
c
      call copy( idim, k, np, v, u )
      if ( ism .eq. 1 ) then
c
c  damped Jacobi; no damping at boundaries
c
         do 10 n = 1, ns
             i = np(k)
             v(i) = (b(i)-a(i,1)*u(i+1))/a(i,0)
             i = np(k+1)-1
             v(i) = (b(i)-a(i,-1)*u(i-1))/a(i,0)
             do 5 i = np(k)+1, np(k+1) -2
               v(i) = (b(i)-a(i,-1)*u(i-1)-a(i,1)*u(i+1))/a(i,0)
               v(i) = om*v(i)+(1d0-om)*u(i)
c       write(i2,'(a,e11.3)')' v after',v(i)
    5       continue
c            write(i2,*)' v after smoothing',n
c            call vecout( i2, idim, k, np, v )
            call copy( idim, k, np, u, v )
   10    continue
c
      else if ( ism .eq. 2 ) then
c
c  forward Gauss-Seidel
c
         do 20 n = 1, ns
            i = ifirst
            v(i) = ( b(i) - a(i,1)*v(i+1) ) / a(i,0)
            do 15 i = ifirst + 1, ilast - 1
               v(i) = ( b(i) - a(i,1)*v(i+1) - a(i,-1)*v(i-1) ) / a(i,0)
   15       continue
            i = ilast
            v(i) = ( b(i) - a(i,-1)*v(i-1) ) / a(i,0)
   20    continue
      else if ( ism .eq. 3 ) then
c
c  backward Gauss-Seidel
c
         do 30 n = 1, ns
            i = ilast
            v(i) = ( b(i) - a(i,-1)*v(i-1) ) / a(i,0)
            do 25 i = ilast - 1, ifirst + 1, -1
               v(i) = ( b(i) - a(i,1)*v(i+1) - a(i,-1)*v(i-1) ) / a(i,0)
   25       continue
            i = ifirst
            v(i) = ( b(i) - a(i,1)*v(i+1) ) / a(i,0)
   30    continue
      else if ( ism .eq. 4 ) then
c
c  symmetric Gauss-Seidel ( first forward and then backward )
c
         i = ifirst
         v(i) = ( b(i) - a(i,1)*v(i+1) ) / a(i,0)
         do 40 n = 1, ns
            do 34 i = ifirst + 1, ilast - 1
               v(i) = ( b(i) - a(i,1)*v(i+1) - a(i,-1)*v(i-1) ) / a(i,0)
   34       continue
            i = ilast
            v(i) = ( b(i) - a(i,-1)*v(i-1) ) / a(i,0)
            do 36 i = ilast - 1, ifirst + 1, -1
               v(i) = ( b(i) - a(i,1)*v(i+1) - a(i,-1)*v(i-1) ) / a(i,0)
   36       continue
            i = ifirst
            v(i) = ( b(i) - a(i,1)*v(i+1) ) / a(i,0)
   40    continue
      else if ( ism .eq. 5 ) then
c
c  damped white-black
c
         if ( om .le. 1d-12 ) then
            write(*,*) 'smooth: om = ',om, 'has wrong value'
            stop
         endif
         om1 = (1d0 - om)/om
         do 50 n = 1, ns
c
c odd numbers first
c
            i = ifirst
            v(i) = ( om1*a(i,0)*v(i) - a(i,1)*v(i+1) + b(i) ) /a(i,0)
            do 42 i = ifirst + 2, ilast - 1, 2
               v(i) =  om*( - a(i,-1)*v(i-1) + om1*a(i,0)*v(i)
     $                               - a(i,1)*v(i+1) + b(i) ) /a(i,0)
   42       continue
            if ( odd ) then
               i = ilast
               v(i) =  om*( - a(i,-1)*v(i-1)
     $                      + om1*a(i,0)*v(i) + b(i) ) /a(i,0)
            endif
c
c even numbers next
c
            do 44 i = ifirst + 1, ilast - 1, 2
               v(i) =  om*( - a(i,-1)*v(i-1) + om1*a(i,0)*v(i)
     $                           - a(i,1)*v(i+1) + b(i) ) /a(i,0)
   44       continue
            if ( .not. odd ) then
               i = ilast
               v(i) =  om*( - a(i,-1)*v(i-1)
     $                      + om1*a(i,0)*v(i) + b(i) ) /a(i,0)
            endif
   50    continue
c
      else if ( ism .eq. 6 ) then
c
c  damped Jacobi
c
         do 60 n = 1, ns
c            write(i2,*)' b'
c            call vecout( i2, idim, k, np, b )
c            write(i2,*)' u before smoothing'
c            call vecout( i2, idim, k, np, u )
c            write(i2,*)' v before smoothing'
c            call vecout( i2, idim, k, np, v )
             i = np(k)
             v(i) = (b(i)-a(i,1)*u(i+1))/a(i,0)
             v(i) = om*v(i)+(1d0-om)*u(i)
             i = np(k+1)-1
             v(i) = (b(i)-a(i,-1)*u(i-1))/a(i,0)
             v(i) = om*v(i)+(1d0-om)*u(i)
             do 65 i = np(k)+1, np(k+1) -2
               v(i) = (b(i)-a(i,-1)*u(i-1)-a(i,1)*u(i+1))/a(i,0)
               v(i) = om*v(i)+(1d0-om)*u(i)
c       write(i2,'(a,e11.3)')' v after',v(i)
   65       continue
c            write(i2,*)' v after smoothing',n
c            call vecout( i2, idim, k, np, v )
            call copy( idim, k, np, u, v )
   60    continue
      else
         write(i2,*)'smooth: ism out of range'
         write(*,*)' aborted'
         stop
      endif
c
c end of smooth
c
      end

      double precision function u(i2,ib,icase,x,f0,f1,aa,bb)
c
c  Exact solution as a function of x of convection-diffusion equation
c  Implemented only for icase = 2; c = 0.
c  u = c1 + c2 * exp( bb * x / a)
      implicit none
      integer icase, ib(13), i2
      double precision x, f0, f1, aa,bb, e, c1, c2
      if ( icase .eq. 2 ) then
         if ( ib(1).eq.1 .and. ib(2).eq.1 ) then
             e = exp(-bb/aa)
             c2 = (f1 - f0)/(1 - e)
             c1 = f0 - c2*e
             u = c1 + c2 * exp (bb*(x-1)/aa)
         else if ( ib(1).eq.1 .and. ib(2).eq.2) then
             c2 = f1*aa/bb
             c1 = f0 - c2*exp(-bb/aa)
             u = c1 + c2 * exp (bb*(x-1)/aa)
         else if ( ib(1).eq.2 .and. ib(2).eq.1 ) then
             c2 = aa*f0/bb
             c1 = f1 - c2*exp(bb/aa)
             u = c1 + c2 * exp (bb*x/aa)
         else
            u=0d0
            write(i2,*) ' Exact solution not implemented for this case'
         endif
      else
         u = 0d0
         write( i2,* )  ' Exact solution not implemented for this case'
      endif
      end
      subroutine uexact( i2, icase, idim, k, ib, np, uce,f0,f1,aa,bb)
c***********************************************************************
c
c  Programmer    F.N. van de Vosse
c  Version 1.0   date 18-01-92
c
c***********************************************************************
c
c     uce = exact continous solution on grid k
c
c ****************************************************************************
c
c                  INPUT  PARAMETERS
c
      implicit none
      integer i2, icase, idim, k, np(11), ib(13)
      double precision f0, f1, uce(idim),aa,bb
c ***************************************************************************
c
c                 LOCAL QUANTITIES
c
      integer i
      double precision h, x, u
c **************************************************************************
c
c            SUBROUTINES CALLED
c
c   u(icase,x) user provided function giving the exact solution in x
c
c ***************************************************************************
      h = 1d0/dble( np(k+1) - np(k) - 1 )
      if ( ib(4) .eq. 2 ) h = 1d0/dble( np(k+1) -np(k) )
      do 10 i = np(k), np(k+1)-1
         x = dble(i-np(k))*h
         if ( ib(4) .eq. 2 ) x = (dble(i-np(k))+.5d0)*h
         uce(i) = u(i2,ib,icase,x,f0,f1,aa,bb)
  10  continue
c
c end of uexact
c
      end
      subroutine vecout( i2, idim, k, np, v )
c **************************************************************************
c
c     prints vector v  on grid k
c
c ****************************************************************************
c
c                  INPUT  PARAMETERS
c
      implicit none
      integer idim, i2, k, np(11)
      double precision v(idim)
c ***************************************************************************
c
c                 LOCAL QUANTITIES
c
      integer i, n, m, ii
c **************************************************************************
c
c            SUBROUTINES CALLED
c
c                none
c ***************************************************************************
c
c                 I/O
c    output goes to logical unit number i2
c
c ***************************************************************************
      write( i2, * )' vector at grid', k
      n = np(k+1)-np(k)
      m = n/4
      do 10 i= 1, 4*m, 4
         ii = i+np(k)-1
         write( i2, '(4(i3,e11.3))' )
     $          i, v(ii), i+1, v(ii+1), i+2, v(ii+2), i+3, v(ii+3)
  10  continue
c
      write( i2, '(4(i3,e11.3))' ) (i,v(i+np(k)-1),i=4*m+1,n)
c
c end of vecout
c
      end