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 U(trans)z=c and L(trans)Py=z 
c
c
c     file: yaec.f
c
      parameter (n = 3)
      dimension a(n,n),s(n),c(n),y(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 (c(i),i=1,n) /6.0,-5.0,3.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 j=1,n
         sum = c(j)
         do 2 i=1,j-1
            sum = sum - a(p(i),j)*z(i)
 2       continue
         z(j) = sum/a(p(j),j)
 3    continue
      do 5 j=n,1,-1
         sum = z(j)
         do 4 i=j+1,n
            sum = sum - a(p(i),j)*y(p(i))
 4       continue
         y(p(j)) = sum
 5    continue
c
      do 6 i=1,n
         print 7,i,y(i)
 6    continue
c
 7    format (3x,'y(',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