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 5.1 
c
c     Example to test power method for eigenvalues 
c
c
c     file: poweracc.f
c
      parameter (n=3,m=30) 
      dimension a(n,n),x(n),y(n),at(m),r(m)
      data (a(1,j),j=1,n)/6.,5.,-5./
      data (a(2,j),j=1,n)/2.,6.,-2./
      data (a(3,j),j=1,n)/2.,5.,-1./
      data (x(j),j=1,n)/-1.,1.,1./ 
c
      print *
      print *,' Power method with Aitken Acceleration'
      print *,' Section 5.1, Kincaid-Cheney'
      print *
      print 3
      print 4,0,x
c
      do 2 k=1,m
         call prod(n,a,x,y)
         r(k) = y(2)/x(2)
         call dot(n,x,y,t) 
         call dot(n,x,x,s) 
         if ((k .ge. 3) .and. (k .le. m-2)) then
            at(k) = (r(k-2)*r(k) - r(k-1)**2)/(r(k)-2.*r(k-1)+r(k-2))
         endif
         s = t/s 
         call normal(n,y) 
         call store(n,y,x)
         print 4,k,x,r(k),s,at(k)
 2    continue
c
 3    format (1x,'k',7x,'x(1)',10x,'x(2)',10x,'x(3)',9x,'ratio',
     +        5x,'rayleigh quot.',3x,'acc. quot.')
 4    format (1x,i2,1x,6(e13.6,1x))
      stop
      end 
c
      subroutine dot(n,x,y,t) 
c
c     compute dot product t=<x,y>
c
      dimension x(n),y(n)
c
      t = 0.0
      do 2 j=1,n
        t = t + x(j)*y(j) 
 2    continue
c   
      return
      end 
c 
      subroutine prod(n,a,x,y)
c
c     compute y=ax
c
      dimension a(n,n),x(n),y(n)
c
      do 2 i=1,n
         y(i) = 0.0
         do 2 j=1,n
           y(i)=y(i)+a(i,j)*x(j) 
 2    continue
c      
      return
      end 
c 
      subroutine store(n,x,y)
c
c     put x into y
c
      dimension x(n),y(n)
c
      do 2 j=1,n
         y(j)=x(j)
 2    continue
c    
      return
      end 
c 
      subroutine norm(n,x,t) 
c
c     compute t= max-norm of x
c
      real x(n),max
c
      t = 0.0
      do 2 j=1,n
         t = max(t,abs(x(j)))
 2    continue
c   
      return
      end 
c 
      subroutine normal(n,x)
c
c     normalize x with max norm
c
      dimension x(n)
c
      call norm (n,x,t)
      do 2 j=1,n
         x(j) = x(j)/t
 2    continue
c   
      return
      end