```

SUBROUTINE SGBSL( ABD, LDA, N, ML, MU, IPVT, B, JOB ) 2,4

c         Solves the real band system
c            A * X = B  or  transpose(A) * X = B
c         using the factors computed by SBGCO or SGBFA.

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

c     Input:

C        ABD     REAL(LDA, N)
c                the output from SBGCO or SGBFA.

C        LDA     INTEGER
c                the leading dimension of the array  ABD .

C        N       INTEGER
c                the order of the original matrix.

C        ML      INTEGER
c                number of diagonals below the main diagonal.

C        MU      INTEGER
c                number of diagonals above the main diagonal.

C        IPVT    INTEGER(N)
c                the pivot vector from SBGCO or SGBFA.

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

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

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

c     .. Scalar Arguments ..

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

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

INTEGER   K, KB, L, LA, LB, LM, M, NM1
REAL      T
c     ..
c     .. External Functions ..

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

EXTERNAL  SAXPY
c     ..
c     .. Intrinsic Functions ..

INTRINSIC MIN
c     ..

M   = MU + ML + 1
NM1 = N - 1

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

c                               ** first solve L*Y = B
IF( ML.NE.0 ) THEN

DO 10 K = 1, NM1

LM = MIN( ML, N - K )
L  = IPVT( K )
T  = B( L )

IF( L.NE.K ) THEN

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

END IF

CALL SAXPY( LM, T, ABD( M + 1,K ), 1, B( K + 1 ), 1 )

10       CONTINUE

END IF

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

K      = N + 1 - KB
B( K ) = B( K ) / ABD( M, K )
LM     = MIN( K, M ) - 1
LA     = M - LM
LB     = K - LM
T      = -B( K )

CALL SAXPY( LM, T, ABD( LA,K ), 1, B( LB ), 1 )

20    CONTINUE

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

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

LM     = MIN( K, M ) - 1
LA     = M - LM
LB     = K - LM
T      = SDOT( LM, ABD( LA,K ), 1, B( LB ), 1 )
B( K ) = ( B( K ) - T ) / ABD( M, K )

30    CONTINUE

c                                  ** now solve trans(L)*X = Y
IF( ML.NE.0 ) THEN

DO 40 KB = 1, NM1

K      = N - KB
LM     = MIN( ML, N - K )
B( K ) = B( K ) + SDOT( LM, ABD( M+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 IF

END
```