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 Inverse Power method
c
c
c     file: ex2s51.f
c
      parameter (n=3,m=25)
      dimension a(n,n),x(n),xx(n),y(n),s(n)
      integer p(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) /3.,7.,-13./
c
      print *
      print *,' Inverse Power method example'
      print *,' Section 5.1, Kincaid-Cheney'
      print *
      print 3
      print 4,0,x
c
      call gauss (n,n,a,p,s)
c
      do 2 k=1,m
         call store (n,x,xx)
         call solve (n,n,a,p,xx,y)
         r = y(1)/x(1)
         call dot(n,x,y,t)
         call dot(n,x,x,ss)
         ss = t/ss
         call normal (n,y)
         call store (n,y,x)
         print 4,k,x,r,ss
 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 gauss(n,ia,a,p,s)
c
c     gaussian elimination
c
      dimension a(ia,n),s(n)
      integer p(n)
      real max
c
      do 3 i=1,n
         p(i) = i
         smax = 0.0
         do 2 j=1,n
            smax = max(smax,abs(a(i,j)))
 2       continue
         s(i) = smax
 3    continue
c
      do 7 k=1,n-1
         rmax = 0.0
         do 4 i=k,n
            r = abs(a(p(i),k))/s(p(i))
            if (r .gt. rmax) then
               j = i
               rmax = r
            endif
 4       continue
         pk = p(j)
         p(j) = p(k)
         p(k) = pk
         do 6 i=k+1,n
            z = a(p(i),k)/a(p(k),k)
            a(p(i),k) = z
            do 5 j=k+1,n
               a(p(i),j) = a(p(i),j) - z*a(p(k),j)
 5          continue
 6       continue
 7    continue
c
      return
      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 
c 
      subroutine solve(ia,n,a,p,b,x)
c
c     solve for x
c
      dimension a(ia,n),b(n),x(n)
      integer p(n)
c
      do 2 k=1,n-1
         do 2 i=k+1,n
            b(p(i)) = b(p(i)) - a(p(i),k)*b(p(k))
 2    continue
c
      do 4 i=n,1,-1
         sum = b(p(i))
         do 3 j=i+1,n
            sum = sum - a(p(i),j)*x(j)
 3       continue
         x(i) = sum/a(p(i),i)
 4    continue
c
      return
      end