SUBROUTINE DGESLM (A,LDA,N,IPVT,B)
      INTEGER LDA,N,IPVT(N)
      DOUBLE PRECISION A(LDA,N),B(N)
C
C     DGESLM SOLVES THE (REAL) LINEAR SYSTEM
C            A * X = B
C     USING THE FACTORS COMPUTED BY DGECOM OR DGEFAM.
C
C     ON ENTRY:
C
C        A       DOUBLE PRECISION(LDA,N)
C                THE OUTPUT FROM DGECOM OR DGEFAM.
C
C        LDA     INTEGER
C                THE LEADING DIMENSION OF THE ARRAY A AS DECLARED
C                IN THE MAIN CALLING PROGRAM.
C
C        N       INTEGER
C                THE ORDER OF THE MATRIX A.
C
C        IPVT    INTEGER(N)
C                THE PIVOT VECTOR FROM DGECOM OR DGEFAM.
C
C        B       DOUBLE PRECISION(N)
C                THE RIGHT HAND SIDE VECTOR.
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 DGECOM HAS SET RCOND .GT. 0.0
C        OR DGEFAM HAS SET INFO .EQ. 0.
C
C     TO COMPUTE  INVERSE(A) * B  WHERE B IS A MATRIX WITH M COLUMNS
C             CALL DGECOM (A,LDA,N,IPVT,RCOND,Z)
C             IF (RCOND IS TOO SMALL) GO TO ...
C             DO 10 J=1,M
C                CALL DGESLM (A,LDA,N,IPVT,B(1,J))
C          10 CONTINUE
C
C     THIS VERSION IS ADAPTED FROM THE LINPACK SUBROUTINE DGESL BY
C     REPLACEMENT OF CALLS TO THE BLAS WITH IN-LINE CODE.  MODIFICA-
C     TION DONE BY ALAN J. LAUB, UNIVERSITY OF SOUTHERN CALIFORNIA,
C     APRIL 1980.
C
C     FORTRAN FUNCTIONS CALLED:
C
C     NONE
C
C     INTERNAL VARIABLES:
C
      INTEGER I,K,KB,KM1,KP1,L,NM1
      DOUBLE PRECISION T
C
      IF (N .EQ. 1) GO TO 60
      NM1 = N-1
C
C     FIRST SOLVE L*Y = B
C
      DO 30 K=1,NM1
         KP1 = K+1
         L = IPVT(K)
         T = B(L)
         IF (L .EQ. K) GO TO 10
         B(L) = B(K)
         B(K) = T
   10    CONTINUE
         DO 20 I=KP1,N
            B(I) = B(I)-T*A(I,K)
   20    CONTINUE
   30 CONTINUE
C
C     NOW SOLVE  U*X = Y
C
      DO 50 KB=1,NM1
         K = N+1-KB
         KM1 = K-1
         B(K) = B(K)/A(K,K)
         T = -B(K)
         DO 40 I=1,KM1
            B(I) = B(I)+T*A(I,K)
   40    CONTINUE
   50 CONTINUE
   60 CONTINUE
      B(1) = B(1)/A(1,1)
      RETURN
      END