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.3 
c
c     Example of Gaussian elimination with scaled row pivoting 
c     Solving Lz=Pb and Ux=z 
c
c
c     file: paxeb.f
c
      parameter (n = 3)
      dimension a(n,n),s(n),b(n),x(n),z(n)
      integer p(n)
      data (a(1,j),j=1,n) /2.0,3.0,-6.0/
      data (a(2,j),j=1,n) /1.0,-6.0,8.0/
      data (a(3,j),j=1,n) /3.0,-2.0,1.0/
      data (b(i),i=1,n) /-1.0,3.0,2.0/
c
      print *
      print *,' Gaussian elimination with scaled row pivoting'
      print *,' Section 4.3, Kincaid-Cheney'
      print *
c
      call gauss(n,a,s,p)
c
      do 3 i=1,n
         sum = b(p(i))
         do 2 j=1,i-1
            sum = sum - a(p(i),j)*z(j)
 2       continue
         z(i) = sum
 3    continue
c
      do 5 i=n,1,-1
         sum = z(i)
         do 4 j=i+1,n
            sum = sum - a(p(i),j)*x(j)
 4       continue
         x(i) = sum/a(p(i),i)
 5    continue
c
      do 6 i=1,n
         print 7,i,x(i)
 6    continue
c
 7    format (3x,'x(',i2,') =',e13.6)
      stop
      end
c
      subroutine gauss(n,a,s,p)
c 
c     gaussian elimination
c
      dimension a(n,n),s(n)   
      integer p(n)
      real max
c
      do 9 i=1,n
         p(i) = i
         smax = 0.0 
         do 8 j=1,n
            smax = max(smax,abs(a(i,j)))
 8       continue
         s(i) = smax
 9    continue
c
      do 13 k=1,n-1
         rmax = 0.0
         do 10 i=k,n
            r = abs(a(p(i),k))/s(p(i))
            if (r .gt. rmax) then
               j = i
               rmax = r
            endif
 10      continue
c
         pk = p(j)
         p(j) = p(k)
         p(k) = pk
c
         do 12 i=k+1,n      
            z = a(p(i),k)/a(p(k),k)       
            a(p(i),k) = z
            do 11 j=k+1,n    
               a(p(i),j) = a(p(i),j) - z*a(p(k),j)      
 11         continue
 12      continue  
 13   continue    
c
      return
      end