SUBROUTINE SGECO( A, LDA, N, IPVT, RCOND, Z ) 2,17

c         Factors a real matrix by Gaussian elimination
c         and estimates the condition of the matrix.

c         Revision date:  8/1/82
c         Author:  Moler, C. B. (U. of New Mexico)

c         If  RCOND  is not needed, SGEFA is slightly faster.
c         To solve  A*X = B , follow SGECO by SGESL.

c     On entry

c        A       REAL(LDA, N)
c                the matrix to be factored.

c        LDA     INTEGER
c                the leading dimension of the array  A .

c        N       INTEGER
c                the order of the matrix  A .

c     On return

c        A       an upper triangular matrix and the multipliers
c                which were used to obtain it.
c                The factorization can be written  A = L*U , where
c                L  is a product of permutation and unit lower
c                triangular matrices and  U  is upper triangular.

c        IPVT    INTEGER(N)
c                an integer vector of pivot indices.

c        RCOND   REAL
c                an estimate of the reciprocal condition of  A .
c                For the system  A*X = B , relative perturbations
c                in  A  and  B  of size  epsilon  may cause
c                relative perturbations in  X  of size  epsilon/RCOND .
c                If  RCOND  is so small that the logical expression
c                           1.0 + RCOND .EQ. 1.0
c                is true, then  A  may be singular to working
c                precision.  In particular,  RCOND  is zero  if
c                exact singularity is detected or the estimate
c                underflows.

C        Z       REAL(N)
c                a work vector whose contents are usually unimportant.
c                If  A  is close to a singular matrix, then  Z  is
c                an approximate null vector in the sense that
c                norm(A*Z) = RCOND*norm(A)*norm(Z) .

c ------------------------------------------------------------------

c     .. Scalar Arguments ..

INTEGER   LDA, N
REAL      RCOND
c     ..
c     .. Array Arguments ..

INTEGER   IPVT( * )
REAL      A( LDA, * ), Z( * )
c     ..
c     .. Local Scalars ..

INTEGER   INFO, J, K, KB, KP1, L
REAL      ANORM, EK, S, SM, T, WK, WKM, YNORM
c     ..
c     .. External Functions ..

REAL      SASUM, SDOT
EXTERNAL  SASUM, SDOT
c     ..
c     .. External Subroutines ..

EXTERNAL  SAXPY, SGEFA, SSCAL
c     ..
c     .. Intrinsic Functions ..

INTRINSIC ABS, MAX, SIGN
c     ..

c                        ** compute 1-norm of A
ANORM  = 0.0E0
DO 10 J = 1, N
ANORM  = MAX( ANORM, SASUM( N,A( 1,J ),1 ) )
10 CONTINUE
c                                      ** factor

CALL SGEFA( A, LDA, N, IPVT, INFO )

c     RCOND = 1/(norm(A)*(estimate of norm(inverse(A)))) .
c     estimate = norm(Z)/norm(Y) where  A*Z = Y  and  trans(A)*Y = E .
c     trans(A) is the transpose of A.  The components of E  are
c     chosen to cause maximum local growth in the elements of W  where
c     trans(U)*W = E.  The vectors are frequently rescaled to avoid
c     overflow.

c                        ** solve trans(U)*W = E
EK = 1.0E0

DO 20 J = 1, N
Z( J ) = 0.0E0
20 CONTINUE

DO 50 K = 1, N

IF( Z( K ).NE.0.0E0 ) EK = SIGN( EK, -Z( K ) )

IF( ABS( EK - Z( K ) ).GT.ABS( A( K,K ) ) ) THEN

S  = ABS( A( K,K ) ) / ABS( EK - Z( K ) )

CALL SSCAL( N, S, Z, 1 )

EK = S*EK

END IF

WK   = EK - Z( K )
WKM  = -EK - Z( K )
S    = ABS( WK )
SM   = ABS( WKM )

IF( A( K,K ).NE.0.0E0 ) THEN

WK   = WK / A( K, K )
WKM  = WKM / A( K, K )

ELSE

WK   = 1.0E0
WKM  = 1.0E0

END IF

KP1  = K + 1

IF( KP1.LE.N ) THEN

DO 30 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 ) )
30       CONTINUE

IF( S.LT.SM ) THEN

T  = WKM - WK
WK = WKM

DO 40 J = KP1, N
Z( J ) = Z( J ) + T*A( K, J )
40          CONTINUE

END IF

END IF

Z( K ) = WK

50 CONTINUE

S  = 1.0E0 / SASUM( N, Z, 1 )

CALL SSCAL( N, S, Z, 1 )
c                                ** solve trans(L)*Y = W
DO 60 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 ) ).GT.1.0E0 ) THEN

S  = 1.0E0 / ABS( Z( K ) )

CALL SSCAL( N, S, Z, 1 )

END IF

L      = IPVT( K )
T      = Z( L )
Z( L ) = Z( K )
Z( K ) = T
60 CONTINUE

S  = 1.0E0 / SASUM( N, Z, 1 )

CALL SSCAL( N, S, Z, 1 )
c                                 ** solve L*V = Y
YNORM  = 1.0E0

DO 70 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 ) ).GT.1.0E0 ) THEN

S  = 1.0E0 / ABS( Z( K ) )

CALL SSCAL( N, S, Z, 1 )

YNORM  = S*YNORM
END IF

70 CONTINUE

S  = 1.0E0 / SASUM( N, Z, 1 )

CALL SSCAL( N, S, Z, 1 )
c                                  ** solve  U*Z = V
YNORM  = S*YNORM

DO 80 KB = 1, N

K  = N + 1 - KB

IF( ABS( Z( K ) ).GT.ABS( A( K,K ) ) ) THEN

S  = ABS( A( K,K ) ) / ABS( Z( K ) )

CALL SSCAL( N, S, Z, 1 )

YNORM  = S*YNORM

END IF

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 )

80 CONTINUE
c                                   ** 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

END