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 4.6
c
c     Example of Chebyshev acceleration on Jacobi method
c
c
c     file: ex6s46.f
c
      parameter (n=4,m=50)
      dimension u(n),v(n),g(n,n),c(n)
      data (g(1,j),j=1,n) /0.0,0.25,0.25,0.0/
      data (g(2,j),j=1,n) /0.25,0.0,0.0,0.25/
      data (g(3,j),j=1,n) /0.25,0.0,0.0,0.25/
      data (g(4,j),j=1,n) /0.0,0.25,0.25,0.0/           
      data c /-1.0,0.0,1.0,-1.0/
      data u /0.0,0.0,0.0,0.0/
      data delta /1.0e-4/
c
      print *
      print *,' Chebyshev acceleration example'
      print *,' Section 4.6, Kincaid-Cheney'
      print *
c
      b = 0.5
      a = -b
      gamma = 2.0/(2.0 - (b + a))
      alpha = (0.5*(b - a)/(2.0 - (b + a)))**2
      print 3,0,u
      call extrap(gamma,n,G,c,u,v)
      print 3,1,v
      rho = 1.0/(1.0 - 2.0*alpha)
      call cheb(rho,gamma,n,g,c,u,v)
      print 3,2,u
      do 2 k=3,m,2
         rho = 1.0/(1.0 - rho*alpha)
         call cheb(rho,gamma,n,G,c,v,u)
         print 3,k,v
         rho = 1.0/(1.0 - rho*alpha)
         call cheb(rho,gamma,n,G,c,u,v)
         print 3,k+1,u
         if (vnorm(n,u,v) .lt. delta) stop
 2    continue
c
 3    format (3x,i2,2x,4(e13.6,2x))
      stop
      end    
c
      subroutine extrap(gamma,n,G,c,u,v)
c
c     extrapolation procedure
c
      dimension u(n),v(n),G(n,n),c(n)
c
      do 2 i=1,n
         v(i) = gamma*c(i) + (1.0 - gamma)*u(i)
         do 2 j=1,n
            v(i) = v(i) + gamma*g(i,j)*u(j)
 2       continue
c
      return
      end
c
      subroutine cheb(rho,gamma,n,G,c,u,v)
c
c     chebyshev acceleration 
c
      dimension u(n),v(n),G(n,n),c(n)
c
      do 2 i=1,n
         u(i) = gamma*rho*c(i)+rho*(1.0-gamma)*v(i)+(1.0-rho)*u(i)
 2    continue
c
      do 3 i=1,n
          do 3 j=1,n
            u(i) = u(i) + rho*gamma*G(i,j)*v(j) 
 3    continue
c
      return
      end
c
      function vnorm(n,x,y) 
      real x(n),y(n)
c
c     compute max-norm of x - y
c
      temp = 0.0
      do 2 i=1,n
         temp = max(temp,abs(x(i) - y(i)))
 2    continue
      vnorm = temp
c
      return
      end