c
c     Numerical Analysis:
c     The Mathematics of Scientific Computing
c     D.R. Kincaid & E.W. Cheney
c     Brooks/Cole Publ., 1990
c
c     Section 9.8
c
c     multigrid program used on the equation u'' = f, in one variable.
c
c     parameter m - the number of mesh refinements used in the method 
c               n - 2**m
c               k - number of iterations in the Gauss-Seidel method 
c                   to be used on each grid.
c
c 
c     file: mgrid1.f
c
      parameter (m=6, k=3, np1=65)
      dimension v(0:np1),w(0:np1)
c
c     the function f is the right-hand side of the differential equation. 
c     the function g is the solution function.
c
      real f,g
      f(x) = cos(x)
      g(x) = -cos(x) + x*(cos(1.0) - 1.0) + 1.0 
c
      print * 
      print *,' Multigrid method example'
      print *,' Section 9.8, Kincaid-Cheney'
      print *
c
c     first we define the solution of the difference equations on the
c     coarsest grid. this grid has three points, two of which are on the 
c     boundary.
c
      h = 0.5
      n = 1
      v(0) = 0.0
      v(2) = 0.0
      v(1) = -f(0.5)/8.0
      print *,' n =',n,' h =',h
      print *
c
c     now the main iteration (outer loop) begins.
c
      do 8 i=1,m-1
         h = 0.5*h
         n = 2*n + 1
         print *,' n =',n,' h =',h
c
c     the next part is the interpolation phase in which an approximate
c     solution on a finer grid is manufactured from an approximation on
c     a coarse grid.
c
         do 2 j=0,n-1
            w(2*j) = v(j)
 2       continue
c
         do 3 j=1,n-1
            w(2*j-1) = 0.5*(v(j-1) + v(j))
 3       continue
c 
         do 4 j=0,2*n-2
            v(j) = w(j)
 4       continue
c
c     next is the Gauss-Seidel iteration on the finer grid (k steps)
c
         do 6 p=1,k
            do 5 j=1,n
               v(j) = 0.5*(v(j-1) + v(j+1) - h*h*f(real(j)*h)) 
 5          continue
 6       continue
c
         do 7 j=0,n+1
            w(j) = abs(g(real(j)*h) - v(j))
 7       continue
c
         print *,' error norm =',vnorm(n+1,w)
         print *
 8    continue
c
      stop
      end 
c
      function vnorm(n,v)
c
c     infinity norm of vector v
c
      dimension v(0:n)
      real max
c
      temp = 0.0
      do 2 i=0,n
         temp = max(temp,abs(v(i)))
 2    continue
      vnorm = temp
c
      return
      end