      SUBROUTINE DGECO(A,LDA,N,IPVT,RCOND,Z)                           
C***BEGIN PROLOGUE  DGECO                                              
C***DATE WRITTEN   780814   (YYMMDD)                                   
C***REVISION DATE  861211   (YYMMDD)                                   
C***CATEGORY NO.  D2A1                                                 
C***KEYWORDS  LIBRARY=SLATEC(LINPACK),                                 
C             TYPE=DOUBLE PRECISION(SGECO-S DGECO-D CGECO-C),          
C             CONDITION NUMBER,GENERAL MATRIX,LINEAR ALGEBRA,MATRIX,   
C             MATRIX FACTORIZATION                                     
C***AUTHOR  MOLER, C. B., (U. OF NEW MEXICO)                           
C***PURPOSE  Factors a double precision matrix by Gaussian elimination 
C            and estimates the condition of the matrix.                
C***DESCRIPTION                                                        
C                                                                      
C     DGECO factors a double precision matrix by Gaussian elimination  
C     and estimates the condition of the matrix.                       
C                                                                      
C     If  RCOND  is not needed, DGEFA is slightly faster.              
C     To solve  A*X = B , follow DGECO by DGESL.                       
C     To compute  INVERSE(A)*C , follow DGECO by DGESL.                
C     To compute  DETERMINANT(A) , follow DGECO by DGEDI.              
C     To compute  INVERSE(A) , follow DGECO by DGEDI.                  
C                                                                      
C     On Entry                                                         
C                                                                      
C        A       DOUBLE PRECISION(LDA, N)                              
C                the matrix to be factored.                            
C                                                                      
C        LDA     INTEGER                                               
C                the leading dimension of the array  A .               
C                                                                      
C        N       INTEGER                                               
C                the order of the matrix  A .                          
C                                                                      
C     On Return                                                        
C                                                                      
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                                                                      
C        IPVT    INTEGER(N)                                            
C                an INTEGER vector of pivot indices.                   
C                                                                      
C        RCOND   DOUBLE PRECISION                                      
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                                                                      
C        Z       DOUBLE PRECISION(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     LINPACK.  This version dated 08/14/78 .                          
C     Cleve Moler, University of New Mexico, Argonne National Lab.     
C                                                                      
C     Subroutines and Functions                                        
C                                                                      
C     LINPACK DGEFA                                                    
C     BLAS DAXPY,DDOT,DSCAL,DASUM                                      
C     Fortran DABS,DMAX1,DSIGN                                         
C***REFERENCES  DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W.,   
C                 *LINPACK USERS  GUIDE*, SIAM, 1979.                  
C***ROUTINES CALLED  DASUM,DAXPY,DDOT,DGEFA,DSCAL                      
C***END PROLOGUE  DGECO                                                
      INTEGER LDA,N,IPVT(1)                                            
      DOUBLE PRECISION A(LDA,1),Z(1)                                   
      DOUBLE PRECISION RCOND                                           
C                                                                      
      DOUBLE PRECISION DDOT,EK,T,WK,WKM                                
      DOUBLE PRECISION ANORM,S,DASUM,SM,YNORM                          
      INTEGER INFO,J,K,KB,KP1,L                                        
C                                                                      
C     COMPUTE 1-NORM OF A                                              
C                                                                      
C***FIRST EXECUTABLE STATEMENT  DGECO                                  
      ANORM = 0.0D0                                                    
      DO 10 J = 1, N                                                   
         ANORM = DMAX1(ANORM,DASUM(N,A(1,J),1))                        
   10 CONTINUE                                                         
C                                                                      
C     FACTOR                                                           
C                                                                      
      CALL DGEFA(A,LDA,N,IPVT,INFO)                                    
C                                                                      
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                                                                      
C     SOLVE TRANS(U)*W = E                                             
C                                                                      
      EK = 1.0D0                                                       
      DO 20 J = 1, N                                                   
         Z(J) = 0.0D0                                                  
   20 CONTINUE                                                         
      DO 100 K = 1, N                                                  
         IF (Z(K) .NE. 0.0D0) EK = DSIGN(EK,-Z(K))                     
         IF (DABS(EK-Z(K)) .LE. DABS(A(K,K))) GO TO 30                 
            S = DABS(A(K,K))/DABS(EK-Z(K))                             
            CALL DSCAL(N,S,Z,1)                                        
            EK = S*EK                                                  
   30    CONTINUE                                                      
         WK = EK - Z(K)                                                
         WKM = -EK - Z(K)                                              
         S = DABS(WK)                                                  
         SM = DABS(WKM)                                                
         IF (A(K,K) .EQ. 0.0D0) GO TO 40                               
            WK = WK/A(K,K)                                             
            WKM = WKM/A(K,K)                                           
         GO TO 50                                                      
   40    CONTINUE                                                      
            WK = 1.0D0                                                 
            WKM = 1.0D0                                                
   50    CONTINUE                                                      
         KP1 = K + 1                                                   
         IF (KP1 .GT. N) GO TO 90                                      
            DO 60 J = KP1, N                                           
               SM = SM + DABS(Z(J)+WKM*A(K,J))                         
               Z(J) = Z(J) + WK*A(K,J)                                 
               S = S + DABS(Z(J))                                      
   60       CONTINUE                                                   
            IF (S .GE. SM) GO TO 80                                    
               T = WKM - WK                                            
               WK = WKM                                                
               DO 70 J = KP1, N                                        
                  Z(J) = Z(J) + T*A(K,J)                               
   70          CONTINUE                                                
   80       CONTINUE                                                   
   90    CONTINUE                                                      
         Z(K) = WK                                                     
  100 CONTINUE                                                         
      S = 1.0D0/DASUM(N,Z,1)                                           
      CALL DSCAL(N,S,Z,1)                                              
C                                                                      
C     SOLVE TRANS(L)*Y = W                                             
C                                                                      
      DO 120 KB = 1, N                                                 
         K = N + 1 - KB                                                
         IF (K .LT. N) Z(K) = Z(K) + DDOT(N-K,A(K+1,K),1,Z(K+1),1)     
         IF (DABS(Z(K)) .LE. 1.0D0) GO TO 110                          
            S = 1.0D0/DABS(Z(K))                                       
            CALL DSCAL(N,S,Z,1)                                        
  110    CONTINUE                                                      
         L = IPVT(K)                                                   
         T = Z(L)                                                      
         Z(L) = Z(K)                                                   
         Z(K) = T                                                      
  120 CONTINUE                                                         
      S = 1.0D0/DASUM(N,Z,1)                                           
      CALL DSCAL(N,S,Z,1)                                              
C                                                                      
      YNORM = 1.0D0                                                    
C                                                                      
C     SOLVE L*V = Y                                                    
C                                                                      
      DO 140 K = 1, N                                                  
         L = IPVT(K)                                                   
         T = Z(L)                                                      
         Z(L) = Z(K)                                                   
         Z(K) = T                                                      
         IF (K .LT. N) CALL DAXPY(N-K,T,A(K+1,K),1,Z(K+1),1)           
         IF (DABS(Z(K)) .LE. 1.0D0) GO TO 130                          
            S = 1.0D0/DABS(Z(K))                                       
            CALL DSCAL(N,S,Z,1)                                        
            YNORM = S*YNORM                                            
  130    CONTINUE                                                      
  140 CONTINUE                                                         
      S = 1.0D0/DASUM(N,Z,1)                                           
      CALL DSCAL(N,S,Z,1)                                              
      YNORM = S*YNORM                                                  
C                                                                      
C     SOLVE  U*Z = V                                                   
C                                                                      
      DO 160 KB = 1, N                                                 
         K = N + 1 - KB                                                
         IF (DABS(Z(K)) .LE. DABS(A(K,K))) GO TO 150                   
            S = DABS(A(K,K))/DABS(Z(K))                                
            CALL DSCAL(N,S,Z,1)                                        
            YNORM = S*YNORM                                            
  150    CONTINUE                                                      
         IF (A(K,K) .NE. 0.0D0) Z(K) = Z(K)/A(K,K)                     
         IF (A(K,K) .EQ. 0.0D0) Z(K) = 1.0D0                           
         T = -Z(K)                                                     
         CALL DAXPY(K-1,T,A(1,K),1,Z(1),1)                             
  160 CONTINUE                                                         
C     MAKE ZNORM = 1.0                                                 
      S = 1.0D0/DASUM(N,Z,1)                                           
      CALL DSCAL(N,S,Z,1)                                              
      YNORM = S*YNORM                                                  
C                                                                      
      IF (ANORM .NE. 0.0D0) RCOND = YNORM/ANORM                        
      IF (ANORM .EQ. 0.0D0) RCOND = 0.0D0                              
      RETURN                                                           
      END                                                              

