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