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 Conjugate Gradient method
c
c
c     file: cg.f
c
      parameter (n=4)
      dimension a(n,n),r(n),b(n),v(n),x(n),z(n)
      data (a(1,i),i=1,n) /420.,210.,140.,105./
      data (a(2,i),i=1,n) /210.,140.,105.,84./
      data (a(3,i),i=1,n) /140.,105.,84.,70./
      data (a(4,i),i=1,n) /105.,84.,70.,60./
      data (b(i),i=1,n) /875.,539.,399.,319./
      data m/15/
      data (x(i),i=1,n) /0.,0.,0.,0./
      data delta /1.0e-6/
c
      print *
      print *,' Iterative refinement example'
      print *,' Section 4.5, Kincaid-Cheney'
      print *
c
      call residual(n,a,x,b,r)
      print *,' The residual vector corresponding to x(0) is:'
      print 7,r
      print *
c
      do 2 i=1,n
         v(i) = r(i)
 2    continue
c
      c = prod(n,r,r)
      do 6 k=1,M
         if (sqrt(prod(n,v,v)) .lt. delta) stop 
         call mult(n,a,v,z)
         t = c/prod(n,v,z)
         do 3 i=1,n
            x(i) = x(i) + t*v(i)
 3       continue
         do 4 i=1,n
            r(i) = r(i) - t*z(i)
 4       continue
c
         d = prod(n,r,r)
c
         do 5 i=1,n
            v(i) = r(i) + (d/c)*v(i)
 5       continue
c
         c = d
c
         print *,' At iteration',k
         print *,' Solution is:'
         print 7,x
         print *,' Residual vector is:'
         print 7,r
         print *
 6    continue
c
 7    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
c
      return
      end
c
      subroutine residual(n,a,x,b,r)
c
c     compute the residual r = b-Ax
c
      dimension a(n,n),x(n),b(n),r(n)
      call mult(n,a,x,r)
c
      do 2 i=1,n
         r(i) = b(i) - r(i)
 2    continue
c
      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)
c
      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
c
      return
      end