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 of power method
c
c
c     file: ex1s51.f
c
      parameter (n=3,m=30) 
      dimension a(n,n),x(n),y(n) 
      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 example'
      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 = y(2)/x(2)
         call dot(n,x,y,t) 
         call dot(n,x,x,s) 
         s = t/s 
         call normal(n,y) 
         call store(n,y,x)
         print 4,k,x,r,s
 2    continue   
c
 3    format (1x,'k',7x,'x(1)',10x,'x(2)',10x,'x(3)',9x,'ratio',
     +         5x,'rayleigh quot.')
 4    format (1x,i2,1x,5(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