```

SUBROUTINE SGBFA( ABD, LDA, N, ML, MU, IPVT, INFO ) 1,3

c         Factors a real band matrix by elimination.

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

c     SGBFA is usually called by SBGCO, but it can be called
c     directly with a saving in time if  RCOND  is not needed.

c     Input:  same as SGBCO

c     On return:

c        ABD,IPVT    same as SGBCO

c        INFO    INTEGER
c                = 0  normal value.
c                = k  if  u(k,k) .eq. 0.0 .  This is not an error
c                     condition for this subroutine, but it does
c                     indicate that SGBSL will divide by zero if
c                     called.  Use  RCOND  in SBGCO for a reliable
c                     indication of singularity.

c     (see SGBCO for description of band storage mode)

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

c     .. Scalar Arguments ..

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

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

INTEGER   I, I0, J, J0, J1, JU, JZ, K, KP1, L, LM, M, MM, NM1
REAL      T
c     ..
c     .. External Functions ..

INTEGER   ISAMAX
EXTERNAL  ISAMAX
c     ..
c     .. External Subroutines ..

EXTERNAL  SAXPY, SSCAL
c     ..
c     .. Intrinsic Functions ..

INTRINSIC MAX, MIN
c     ..

M    = ML + MU + 1
INFO = 0
c                        ** zero initial fill-in columns
J0 = MU + 2
J1 = MIN( N, M ) - 1

DO 20 JZ = J0, J1

I0 = M + 1 - JZ

DO 10 I = I0, ML
ABD( I, JZ ) = 0.0E0
10    CONTINUE

20 CONTINUE

JZ = J1
JU = 0
c                       ** Gaussian elimination with partial pivoting
NM1  = N - 1

DO 50 K = 1, NM1

KP1 = K + 1
c                                  ** zero next fill-in column
JZ = JZ + 1

IF( JZ.LE.N ) THEN

DO 30 I = 1, ML
ABD( I, JZ ) = 0.0E0
30       CONTINUE

END IF
c                                  ** find L = pivot index
LM  = MIN( ML, N - K )
L   = ISAMAX( LM + 1, ABD( M, K ), 1 ) + M - 1
IPVT( K ) = L + K - M

IF( ABD( L,K ).EQ.0.0E0 ) THEN
c                                      ** zero pivot implies this column
INFO = K

ELSE
c                                ** interchange if necessary
IF( L.NE.M ) THEN

T           = ABD( L, K )
ABD( L, K ) = ABD( M, K )
ABD( M, K ) = T
END IF
c                                      ** compute multipliers
T  = - 1.0E0 / ABD( M, K )

CALL SSCAL( LM, T, ABD( M + 1,K ), 1 )

c                               ** row elimination with column indexing

JU = MIN( MAX( JU,MU + IPVT( K ) ), N )
MM = M

DO 40 J = KP1, JU

L  = L - 1
MM = MM - 1
T  = ABD( L, J )

IF( L.NE.MM ) THEN

ABD( L, J ) = ABD( MM, J )
ABD( MM, J ) = T

END IF

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

40       CONTINUE

END IF

50 CONTINUE

IPVT( N ) = N
IF( ABD( M,N ).EQ.0.0E0 ) INFO = N

END
```