      SUBROUTINE DGESL(A,LDA,N,IPVT,B,JOB)                              
C***BEGIN PROLOGUE  DGESL                                               
C***DATE WRITTEN   780814   (YYMMDD)                                    
C***REVISION DATE  861211   (YYMMDD)                                    
C***CATEGORY NO.  D2A1                                                  
C***KEYWORDS  LIBRARY=SLATEC(LINPACK),                                  
C             TYPE=DOUBLE PRECISION(SGESL-S DGESL-D CGESL-C),           
C             LINEAR ALGEBRA,MATRIX,SOLVE                               
C***AUTHOR  MOLER, C. B., (U. OF NEW MEXICO)                            
C***PURPOSE  Solves the double precision system  A*X=B or  TRANS(A)*X=B 
C            using the factors computed by DGECO or DGEFA.              
C***DESCRIPTION                                                         
C                                                                       
C     DGESL solves the double precision system                          
C     A * X = B  or  TRANS(A) * X = B                                   
C     using the factors computed by DGECO or DGEFA.                     
C                                                                       
C     On Entry                                                          
C                                                                       
C        A       DOUBLE PRECISION(LDA, N)                               
C                the output from DGECO or DGEFA.                        
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        IPVT    INTEGER(N)                                             
C                the pivot vector from DGECO or DGEFA.                  
C                                                                       
C        B       DOUBLE PRECISION(N)                                    
C                the right hand side vector.                            
C                                                                       
C        JOB     INTEGER                                                
C                = 0         to solve  A*X = B ,                        
C                = nonzero   to solve  TRANS(A)*X = B  where            
C                            TRANS(A)  is the transpose.                
C                                                                       
C     On Return                                                         
C                                                                       
C        B       the solution vector  X .                               
C                                                                       
C     Error Condition                                                   
C                                                                       
C        A division by zero will occur if the input factor contains a   
C        zero on the diagonal.  Technically this indicates singularity  
C        but it is often caused by improper arguments or improper       
C        setting of LDA .  It will not occur if the subroutines are     
C        called correctly and if DGECO has set RCOND .GT. 0.0           
C        or DGEFA has set INFO .EQ. 0 .                                 
C                                                                       
C     To compute  INVERSE(A) * C  where  C  is a matrix                 
C     with  P  columns                                                  
C           CALL DGECO(A,LDA,N,IPVT,RCOND,Z)                            
C           IF (RCOND is too small) GO TO ...                           
C           DO 10 J = 1, P                                              
C              CALL DGESL(A,LDA,N,IPVT,C(1,J),0)                        
C        10 CONTINUE                                                    
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,DDOT                                                   
C***REFERENCES  DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W.,    
C                 *LINPACK USERS  GUIDE*, SIAM, 1979.                   
C***ROUTINES CALLED  DAXPY,DDOT                                         
C***END PROLOGUE  DGESL                                                 
      INTEGER LDA,N,IPVT(1),JOB                                         
      DOUBLE PRECISION A(LDA,1),B(1)                                    
C                                                                       
      DOUBLE PRECISION DDOT,T                                           
      INTEGER K,KB,L,NM1                                                
C***FIRST EXECUTABLE STATEMENT  DGESL                                   
      NM1 = N - 1                                                       
      IF (JOB .NE. 0) GO TO 50                                          
C                                                                       
C        JOB = 0 , SOLVE  A * X = B                                     
C        FIRST SOLVE  L*Y = B                                           
C                                                                       
         IF (NM1 .LT. 1) GO TO 30                                       
         DO 20 K = 1, NM1                                               
            L = IPVT(K)                                                 
            T = B(L)                                                    
            IF (L .EQ. K) GO TO 10                                      
               B(L) = B(K)                                              
               B(K) = T                                                 
   10       CONTINUE                                                    
            CALL DAXPY(N-K,T,A(K+1,K),1,B(K+1),1)                       
   20    CONTINUE                                                       
   30    CONTINUE                                                       
C                                                                       
C        NOW SOLVE  U*X = Y                                             
C                                                                       
         DO 40 KB = 1, N                                                
            K = N + 1 - KB                                              
            B(K) = B(K)/A(K,K)                                          
            T = -B(K)                                                   
            CALL DAXPY(K-1,T,A(1,K),1,B(1),1)                           
   40    CONTINUE                                                       
      GO TO 100                                                         
   50 CONTINUE                                                          
C                                                                       
C        JOB = NONZERO, SOLVE  TRANS(A) * X = B                         
C        FIRST SOLVE  TRANS(U)*Y = B                                    
C                                                                       
         DO 60 K = 1, N                                                 
            T = DDOT(K-1,A(1,K),1,B(1),1)                               
            B(K) = (B(K) - T)/A(K,K)                                    
   60    CONTINUE                                                       
C                                                                       
C        NOW SOLVE TRANS(L)*X = Y                                       
C                                                                       
         IF (NM1 .LT. 1) GO TO 90                                       
         DO 80 KB = 1, NM1                                              
            K = N - KB                                                  
            B(K) = B(K) + DDOT(N-K,A(K+1,K),1,B(K+1),1)                 
            L = IPVT(K)                                                 
            IF (L .EQ. K) GO TO 70                                      
               T = B(L)                                                 
               B(L) = B(K)                                              
               B(K) = T                                                 
   70       CONTINUE                                                    
   80    CONTINUE                                                       
   90    CONTINUE                                                       
  100 CONTINUE                                                          
      RETURN                                                            
      END                                                               

