SUBROUTINE SGESL( A, LDA, N, IPVT, B, JOB ) 2,4

c         Solves the real system
c            A * X = B  or  transpose(A) * X = B
c         using the factors computed by SGECO or SGEFA.

c         Revision date:  8/1/82
c         Author:  Moler, C. B. (U. of New Mexico)

c     On entry

c        A       REAL(LDA, N)
c                the output from SGECO or SGEFA.

c        LDA     INTEGER
c                the leading dimension of the array  A

c        N       INTEGER
c                the order of the matrix  A

c        IPVT    INTEGER(N)
c                the pivot vector from SGECO or SGEFA.

c        B       REAL(N)
c                the right hand side vector.

c        JOB     INTEGER
c                = 0         to solve  A*X = B ,
c                = nonzero   to solve  transpose(A)*X = B

c     On return

c        B       the solution vector  X

c     Error condition

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 SGECO has set RCOND .GT. 0.0
c        or SGEFA has set INFO .EQ. 0 .

c     To compute  inverse(a) * c  where  c  is a matrix
c     with  p  columns
c           call sgeco(a,lda,n,ipvt,rcond,z)
c           if (rcond is too small) go to ...
c           do 10 j = 1, p
c              call sgesl(a,lda,n,ipvt,c(1,j),0)
c        10 continue

c ---------------------------------------------------------------------

c     .. Scalar Arguments ..

      INTEGER   JOB, LDA, N
c     ..
c     .. Array Arguments ..

      INTEGER   IPVT( * )
      REAL      A( LDA, * ), B( * )
c     ..
c     .. Local Scalars ..

      INTEGER   K, KB, L, NM1
      REAL      T
c     ..
c     .. External Functions ..

      REAL      SDOT
      EXTERNAL  SDOT
c     ..
c     .. External Subroutines ..

      EXTERNAL  SAXPY
c     ..


      NM1  = N - 1

      IF( JOB.EQ.0 ) THEN
c                                 ** solve  A * X = B

c                                     ** first solve  L*Y = B
         DO 10 K = 1, NM1

            L  = IPVT( K )
            T  = B( L )

            IF( L.NE.K ) THEN

               B( L ) = B( K )
               B( K ) = T

            END IF

            CALL SAXPY( N - K, T, A( K+1, K ), 1, B( K+1 ), 1 )

   10    CONTINUE
c                                    ** now solve  U*X = Y
         DO 20 KB = 1, N

            K      = N + 1 - KB
            B( K ) = B( K ) / A( K, K )
            T      = - B( K )

            CALL SAXPY( K-1, T, A( 1, K ), 1, B(1), 1 )

   20    CONTINUE


      ELSE
c                         ** solve  trans(A) * X = B

c                                    ** first solve  trans(U)*Y = B
         DO 30 K = 1, N

            T      = SDOT( K - 1, A( 1,K ), 1, B( 1 ), 1 )
            B( K ) = ( B( K ) - T ) / A( K, K )

   30    CONTINUE

c                                    ** now solve  trans(l)*x = y
         DO 40 KB = 1, NM1

            K      = N - KB
            B( K ) = B( K ) + SDOT( N - K, A( K+1, K ), 1, B( K+1 ), 1)
            L      = IPVT( K )

            IF( L.NE.K ) THEN

               T      = B( L )
               B( L ) = B( K )
               B( K ) = T

            END IF

   40    CONTINUE

      END IF

      END