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.7
c
c     Example of Steepest Descent method
c
c
c     file: steepd.f
c
      parameter (n=4)
      dimension a(n,n),b(n),x(n),v(n),y(n)
      data (a(1,j),j=1,n) /420.,210.,140.,105./
      data (a(2,j),j=1,n) /210.,140.,105.,84./
      data (a(3,j),j=1,n) /140.,105.,84.,70./
      data (a(4,j),j=1,n) /105.,84.,70.,60./
      data (b(i),i=1,n) /875.,539.,399.,319./
      data m/200/
      data (x(i),i=1,n) /0.,0.,0.,0./
c
      print *
      print *,' Steepest Descent method example'
      print *,' Section 4.7, Kincaid-Cheney'
      print *
      print *,' At iteration ',0,' initial x is:'
      print 5,x
      print *
c
      do 4 k=1,m
         call mult(n,a,x,y)
         do 2 i=1,n 
            v(i) = b(i) - y(i)
 2       continue
         c = prod(n,v,v)
         call mult(n,a,v,y)
         d = prod(n,v,y)
         t = c/d
         do 3 i=1,n
            x(i) = x(i) + t*v(i)
 3       continue
         if (0 .eq. mod(k,10)) then
            print *,' At iteration ',k,' solution is:'
            print 5,x
            print *
         endif
 4    continue
c
 5    format (3x,4(e13.6,2x))
      stop
      end
c 
      function prod(n,x,y)
c
c     compute the vector product
c
      dimension x(n),y(n)
      sum = 0.0
      do 2 i=1,n
         sum = sum + x(i)*y(i)
 2    continue
      prod = sum
      return
      end
c
      subroutine mult(n,a,x,y)
c
c     compute the matrix-vector product
c
      dimension a(n,n),x(n),y(n)
      do 3 i=1,n
         sum = 0.0
         do 2 j=1,n
         sum = sum + a(i,j)*x(j)
 2       continue
         y(i) = sum
 3    continue
      return
      end