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.5
c
c     Example of Gaussian elimination followed by iterative improvement
c
c
c     file: ex2s45.f
c
      parameter (n=4,ia=4)
      dimension a(ia,n),b(n),r(n),l(n),s(n),x(n)
      dimension aa(ia,n),e(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./
c
      print *  
      print *,' Gaussian elimination followed by iterative improvement'
      print *,' Section 3.5, Kincaid-Cheney'
      print *
c
c     save a
c
      do 2 i=1,n
         do 2 j=1,n
            aa(i,j) = a(i,j)
 2    continue
c
      call gauss(a,ia,n,l,s)
      call solve(a,ia,n,l,b,x)
      print *,' Initial solution is:'
      print 5,x
      print *
c
      do 4 k=1,5
         call residual(aa,ia,n,b,x,r)
         print *,' At iteration',k
         print *,' Residual vector is:'
         print 5,r
         call solve(a,ia,n,l,r,e)
         print *,' Error vector is:'
         print 5,e
c
         do 3 i=1,n
            x(i) = x(i) + e(i)
 3       continue
c
         print *,' Solution is:'
         print 5,x
         print *
 4    continue
c
 5    format (3x,4(e13.6,2x))
      stop
      end
c
      subroutine residual(a,ia,n,b,x,r)
c
c     computes residual vector r = b-Ax
c
      dimension a(ia,n),b(n),r(n),x(n)
      double precision sum
      do 3 i=1,n
         sum = 0.0
         do 2 j=1,n
            sum = sum + a(i,j)*x(j)
 2       continue
         r(i) = b(i) - sum
 3    continue
c
      return
      end
c
      subroutine gauss(a,ia,n,l,s)    
c
c     gaussian elimination with scaled row pivoting
c
      dimension a(ia,n),l(n),s(n)    
      real max
      do 2 i=1,n
         l(i) = i    
         s(i) = 0.0  
         do 2 j=1,n
            s(i) = max(s(i),abs(a(i,j)))  
 2    continue
      do 4 k=1,n-1
         rmax = 0.0  
         do 3 i=k,n
            r = abs(a(l(i),k))/s(l(i))      
            if (r .gt. rmax) then
               j = i       
               rmax = r    
            endif
 3       continue    
c
         lk = l(j)   
         l(j) = l(k) 
         l(k) = lk   
         do 4 i = k+1,n
            xmult = a(l(i),k)/a(lk,k)       
            a(l(i),k) = xmult     
            do 4 j = k+1,n
               a(l(i),j) = a(l(i),j) - xmult*a(lk,j)     
 4    continue
c
      return      
      end 
c
      subroutine solve(a,ia,n,l,b,x)  
c
c     back substitution LUx = b
c
      dimension a(ia,n),l(n),b(n),x(n),z(4) 
c
      do 3 i=1,n
         sum = b(l(i))
         do 2 j=1,i-1
            sum = sum - a(l(i),j)*z(j)
 2       continue
         z(i) = sum
 3    continue
      do 5 i=n,1,-1
         sum = z(i)       
         do 4 j=i+1,n      
            sum = sum - a(l(i),j)*x(j)    
 4       continue
         x(i) = sum/a(l(i),i)      
 5    continue
c
      return      
      end