      SUBROUTINE DGESL(A,LDA,N,IPVT,B,JOB)                                      
      INTEGER LDA,N,IPVT(1),JOB                                                 
      DOUBLE PRECISION A(LDA,1),B(1)                                            
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                                                                               
C     INTERNAL VARIABLES                                                        
C                                                                               
      DOUBLE PRECISION DDOT,T                                                   
      INTEGER K,KB,L,NM1                                                        
C                                                                               
      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                                                                       
