SUBROUTINE sgeco (a, lda, n, ipvt, rcond, z)  
!     sgeco factors a real matrix by gaussian elimination
!     and estimates the condition of the matrix.
!
!     if  rcond  is not needed, sgefa is slightly faster.
!     to solve  a*x = b , follow sgeco by sgesl.
!     to compute  inverse(a)*c , follow sgeco by sgesl.
!     to compute  determinant(a) , follow sgeco by sgedi.
!     to compute  inverse(a) , follow sgeco by sgedi.
!
!     on entry
!        a       real(lda, n)
!                the matrix to be factored.
!        lda     integer
!                the leading dimension of the array  a .
!        n       integer
!                the order of the matrix  a .
!     on return
!        a       an upper triangular matrix and the multipliers
!                which were used to obtain it.
!                the factorization can be written  a = l*u  where
!                l  is a product of permutation and unit lower
!                triangular matrices and  u  is upper triangular.
!        ipvt    integer(n)
!                an integer vector of pivot indices.
!        rcond   real
!                an estimate of the reciprocal condition of  a .
!                for the system  a*x = b , relative perturbations
!                in  a  and  b  of size  epsilon  may cause
!                relative perturbations in  x  of size  epsilon/rcond .
!                if  rcond  is so small that the logical expression
!                           1.0 + rcond .eq. 1.0
!                is true, then  a  may be singular to working
!                precision.  in particular,  rcond  is zero  if
!                exact singularity is detected or the estimate
!                underflows.
!        z       real(n)
!                a work vector whose contents are usually unimportant.
!                if  a  is close to a singular matrix, then  z  is
!                an approximate null vector in the sense that
!                norm(a*z) = rcond*norm(a)*norm(z) .
!     linpack. this version dated 08/14/78 .
!     cleve moler, university of new mexico, argonne national lab.
!
!     subroutines and functions
!     linpack sgefa
!     blas saxpy,sdot,sscal,sasum
!     fortran abs,amax1,sign
implicit none

INTEGER, INTENT(IN) :: n, lda, ipvt(1)
REAL   , INTENT(OUT):: rcond, z(1)
REAL   , INTENT(IN) :: a(lda,1)
!     internal variables
real :: sdot, ek, t, wk, wkm
real :: anorm, s, sasum, sm, ynorm
integer :: info, j, k, kb, kp1, l
!     compute 1-norm of a
anorm = 0.0e0
do 10 j = 1, n
   anorm = amax1 (anorm, sasum (n, a (1, j), 1) )
10 end do  
!     factor
call sgefa (a, lda, n, ipvt, info)  
!     rcond = 1/(norm(a)*(estimate of norm(inverse(a)))) .
!     estimate = norm(z)/norm(y) where  a*z = y  and  trans(a)*y = e .
!     trans(a)  is the transpose of a .  the components of  e  are
!     chosen to cause maximum local growth in the elements of w  where
!     trans(u)*w = e .  the vectors are frequently rescaled to avoid
!     overflow.
!     solve trans(u)*w = e
ek = 1.0e0  
do 20 j = 1, n  
   z (j) = 0.0e0  
20 end do  
do 100 k = 1, n  
   if (z (k) .ne.0.0e0) ek = sign (ek, - z (k) )  
   if (abs (ek - z (k) ) .le.abs (a (k, k) ) ) goto 30  
   s = abs (a (k, k) ) / abs (ek - z (k) )  
   call sscal (n, s, z, 1)  
   ek = s * ek  
   30    continue  
   wk = ek - z (k)  
   wkm = - ek - z (k)  
   s = abs (wk)  
   sm = abs (wkm)  
   if (a (k, k) .eq.0.0e0) goto 40  
   wk = wk / a (k, k)  
   wkm = wkm / a (k, k)  
   goto 50  
   40    continue  
   wk = 1.0e0  
   wkm = 1.0e0  
   50    continue  
   kp1 = k + 1  
   if (kp1.gt.n) goto 90  
   do 60 j = kp1, n  
      sm = sm + abs (z (j) + wkm * a (k, j) )  
      z (j) = z (j) + wk * a (k, j)  
      s = s + abs (z (j) )  
   60    end do  
   if (s.ge.sm) goto 80  
   t = wkm - wk  
   wk = wkm  
   do 70 j = kp1, n  
      z (j) = z (j) + t * a (k, j)  
   70    end do  
   80    continue  
   90    continue  
   z (k) = wk  
  100 end do  
s = 1.0e0 / sasum (n, z, 1)  
call sscal (n, s, z, 1)  
!
!     solve trans(l)*y = w
do 120 kb = 1, n  
   k = n + 1 - kb  
   if (k.lt.n) z (k) = z (k) + sdot (n - k, a (k + 1, k), 1, z (k &
    + 1), 1)
   if (abs (z (k) ) .le.1.0e0) goto 110  
   s = 1.0e0 / abs (z (k) )  
   call sscal (n, s, z, 1)  
  110    continue  
   l = ipvt (k)  
   t = z (l)  
   z (l) = z (k)  
   z (k) = t  
  120 end do  
s = 1.0e0 / sasum (n, z, 1)  
call sscal (n, s, z, 1)  

ynorm = 1.0e0  

!     solve l*v = y
do 140 k = 1, n  
   l = ipvt (k)  
   t = z (l)  
   z (l) = z (k)  
   z (k) = t  
   if (k.lt.n) call saxpy (n - k, t, a (k + 1, k), 1, z (k + 1), &
    1)
   if (abs (z (k) ) .le.1.0e0) goto 130  
   s = 1.0e0 / abs (z (k) )  
   call sscal (n, s, z, 1)  
   ynorm = s * ynorm  
  130    continue  
  140 end do  
s = 1.0e0 / sasum (n, z, 1)  
call sscal (n, s, z, 1)  
ynorm = s * ynorm  

!     solve  u*z = v
do 160 kb = 1, n  
   k = n + 1 - kb  
   if (abs (z (k) ) .le.abs (a (k, k) ) ) goto 150  
   s = abs (a (k, k) ) / abs (z (k) )  
   call sscal (n, s, z, 1)  
   ynorm = s * ynorm  
  150    continue  
   if (a (k, k) .ne.0.0e0) z (k) = z (k) / a (k, k)  
   if (a (k, k) .eq.0.0e0) z (k) = 1.0e0  
   t = - z (k)  
   call saxpy (k - 1, t, a (1, k), 1, z (1), 1)  
  160 end do  
!     make znorm = 1.0
s = 1.0e0 / sasum (n, z, 1)  
call sscal (n, s, z, 1)  
ynorm = s * ynorm  

if (anorm.ne.0.0e0) rcond = ynorm / anorm  
if (anorm.eq.0.0e0) rcond = 0.0e0

return  
END SUBROUTINE sgeco
