      SUBROUTINE DGEFA(A,LDA,N,IPVT,INFO)                                  
C***BEGIN PROLOGUE  DGEFA                                                  
C***DATE WRITTEN   780814   (YYMMDD)                                       
C***REVISION DATE  861211   (YYMMDD)                                       
C***CATEGORY NO.  D2A1                                                     
C***KEYWORDS  LIBRARY=SLATEC(LINPACK),                                     
C             TYPE=DOUBLE PRECISION(SGEFA-S DGEFA-D CGEFA-C),              
C             GENERAL MATRIX,LINEAR ALGEBRA,MATRIX,MATRIX FACTORIZATION    
C***AUTHOR  MOLER, C. B., (U. OF NEW MEXICO)                               
C***PURPOSE  Factors a double precision matrix by Gaussian elimination.    
C***DESCRIPTION                                                            
C                                                                          
C     DGEFA factors a double precision matrix by Gaussian elimination.     
C                                                                          
C     DGEFA is usually called by DGECO, but it can be called               
C     directly with a saving in time if  RCOND  is not needed.             
C     (Time for DGECO) = (1 + 9/N)*(Time for DGEFA) .                      
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        INFO    INTEGER                                                   
C                = 0  normal value.                                        
C                = K  if  U(K,K) .EQ. 0.0 .  This is not an error          
C                     condition for this subroutine, but it does           
C                     indicate that DGESL or DGEDI will divide by zero     
C                     if called.  Use  RCOND  in DGECO for a reliable      
C                     indication of singularity.                           
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     BLAS DAXPY,DSCAL,IDAMAX                                              
C***REFERENCES  DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W.,       
C                 *LINPACK USERS  GUIDE*, SIAM, 1979.                      
C***ROUTINES CALLED  DAXPY,DSCAL,IDAMAX                                    
C***END PROLOGUE  DGEFA                                                    
      INTEGER LDA,N,IPVT(1),INFO                                           
      DOUBLE PRECISION A(LDA,1)                                            
C                                                                          
      DOUBLE PRECISION T                                                   
      INTEGER IDAMAX,J,K,KP1,L,NM1                                         
C                                                                          
C     GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING                           
C                                                                          
C***FIRST EXECUTABLE STATEMENT  DGEFA                                      
      INFO = 0                                                             
      NM1 = N - 1                                                          
      IF (NM1 .LT. 1) GO TO 70                                             
      DO 60 K = 1, NM1                                                     
         KP1 = K + 1                                                       
C                                                                          
C        FIND L = PIVOT INDEX                                              
C                                                                          
         L = IDAMAX(N-K+1,A(K,K),1) + K - 1                                
         IPVT(K) = L                                                       
C                                                                          
C        ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED             
C                                                                          
         IF (A(L,K) .EQ. 0.0D0) GO TO 40                                   
C                                                                          
C           INTERCHANGE IF NECESSARY                                       
C                                                                          
            IF (L .EQ. K) GO TO 10                                         
               T = A(L,K)                                                  
               A(L,K) = A(K,K)                                             
               A(K,K) = T                                                  
   10       CONTINUE                                                       
C                                                                          
C           COMPUTE MULTIPLIERS                                            
C                                                                          
            T = -1.0D0/A(K,K)                                              
            CALL DSCAL(N-K,T,A(K+1,K),1)                                   
C                                                                          
C           ROW ELIMINATION WITH COLUMN INDEXING                           
C                                                                          
            DO 30 J = KP1, N                                               
               T = A(L,J)                                                  
               IF (L .EQ. K) GO TO 20                                      
                  A(L,J) = A(K,J)                                          
                  A(K,J) = T                                               
   20          CONTINUE                                                    
               CALL DAXPY(N-K,T,A(K+1,K),1,A(K+1,J),1)                     
   30       CONTINUE                                                       
         GO TO 50                                                          
   40    CONTINUE                                                          
            INFO = K                                                       
   50    CONTINUE                                                          
   60 CONTINUE                                                             
   70 CONTINUE                                                             
      IPVT(N) = N                                                          
      IF (A(N,N) .EQ. 0.0D0) INFO = N                                      
      RETURN                                                               
      END                                                                  


