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     complete V-cycle in the multigrid algorithm to solve
c     uxx = g(x), with u(0)=u(l)=0.
c     parameter  m  gives the number of grids to be employed. 
c     parameter  k  gives the number of iterations to be performed 
c                   in each gauss-seidel iteration. 
c     parameter  np1 is simply n+1
c
c
c     file: mgrid2.f
c
      parameter (m=7, np1=128)
      dimension v(m,0:np1), f(m,0:np1), w(0:np1), error(15)
      integer p
      real g,true
      g(x) = cos(x)
      true(x) = -g(x) + x*(g(1.0) - 1.0) + 1.0
c
      print *   
      print *,' V-cycle in the multigrid method'
      print *,' Section 9.8, Kincaid-Cheney'
      print *
c
      error(1) = 1.0
      do 12 k=2,15
c
c     initialize arrays and
c     put best available guess into v**m
c
         n = 2**m - 1
         h = 1.0/real(n+1)
         print *,' n =',n,' h =',h
c
         do 1 i=1,m
            do 1 j=0,np1
               f(i,j) = 0.0
               v(i,j) = 0.0
 1       continue
c
c        put data from the bvp into f**m
c
         do 2 j=1,n
            f(m,j) = g(real(j)*h) 
 2       continue
c
         do 6 i=m,2,-1 
c
c        gauss-seidel iteration (k steps)
c
            do 3 p=1,k
               do 3 j=1,n
                  v(i,j) = 0.5*(v(i,j-1) + v(i,j+1) - (h**2)*f(i,j)) 
 3          continue
c
c           residual vector
c
            do 4 j=1,n
               w(j) = f(i,j) - (v(i,j-1) - 2.0*v(i,j) + v(i,j+1))/h**2 
 4          continue
c
c           apply restriction operator
c
            do 5 j=1,(n-1)/2
               f(i-1,j) = w(2*j)
 5          continue
c
            h = 2.0*h
            n = (n-1)/2
            print *,' n =',n,' h =',h
 6       continue
c
c        solve coarsest (smallest) system exactly
c
         v(1,1) = -g(0.5)/8.0
         print *,' Bottom of V-cycle'
c
c        end downward phase of V-cycle
c        start upward phase of V-cycle
c
         do 10 i=2,m
            h = 0.5*h
            n = 2*n + 1
            print *,' n =',n,' h =',h
c
c           apply extension operation, i.e.
c           interpolation up from coarse grid to fine grid
c
            do 7 j=0,(n+1)/2
               w(2*j) = v(i-1,j)
 7          continue
c
            do 8 j=1,(n+1)/2
               w(2*j-1) = 0.5*(v(i-1,j-1) + v(i-1,j))
 8          continue
c
c           add correction to extension operator
c
            do 9 j=0,n+1
               v(i,j) = v(i,j) + w(j)
 9          continue
c
c           gauss-seidel iteration (k steps)
c
            do 10 p=1,k
               do 10 j=1,n
                  v(i,j) = 0.5*(v(i,j-1) + v(i,j+1) - (h**2)*f(i,j)) 
 10      continue
c
c        output phase
c
         print *,' k =',k,' m =',m
         do 11 j=0,np1
            w(j) = true(real(j)*h) - v(m,j)
 11      continue
         error(k) = vnorm(np1,w)
         print *,' maximun error =', error(k),
     +           ' ratio =', error(k)/error(k-1)
c
 12   continue
c
      stop
      end 
c
      function vnorm(n,v)
c
c     infinity norm of vectorv 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