```

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
```