      SUBROUTINE SNBFA(ABE,LDA,N,ML,MU,IPVT,INFO)
C
C***DESCRIPTION
C      SNBFA FACTORS A REAL BAND MATRIX BY ELIMINATION.
C
C     SNBFA IS USUALLY CALLED BY SNBCO, BUT IT CAN BE CALLED
C     DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED.
C
C     ON ENTRY
C
C        ABE     REAL(LDA, NC)
C                CONTAINS THE MATRIX IN BAND STORAGE.  THE ROWS
C                OF THE ORIGINAL MATRIX ARE STORED IN THE ROWS
C                OF ABE AND THE DIAGONALS OF THE ORIGINAL MATRIX
C                ARE STORED IN COLUMNS 1 THROUGH ML+MU+1 OF ABE.
C                NC MUST BE .GE. 2*ML+MU+1 .
C                SEE THE COMMENTS BELOW FOR DETAILS.
C
C        LDA     INTEGER
C                THE LEADING DIMENSION OF THE ARRAY ABE.
C                LDA MUST BE .GE. N .
C
C        N       INTEGER
C                THE ORDER OF THE ORIGINAL MATRIX.
C
C        ML      INTEGER
C                NUMBER OF DIAGONALS BELOW THE MAIN DIAGONAL.
C                0 .LE. ML .LT. N .
C
C        MU      INTEGER
C                NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL.
C                0 .LE. MU .LT. N .
C                MORE EFFICIENT IF ML .LE. MU .
C
C     ON RETURN
C
C        ABE     AN UPPER TRIANGULAR MATRIX IN BAND STORAGE
C                AND THE MULTIPLIERS WHICH WERE USED TO OBTAIN IT.
C                THE FACTORIZATION CAN BE WRITTEN  A = L*U  WHERE
C                L IS A PRODUCT OF PERMUTATION AND UNIT LOWER
C                TRIANGULAR MATRICES AND  U  IS UPPER TRIANGULAR.
C
C        IPVT    INTEGER(N)
C                AN INTEGER VECTOR OF PIVOT INDICES.
C
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 SNBSL WILL DIVIDE BY ZERO IF
C                CALLED.  USE RCOND IN SNBCO FOR A RELIABLE
C                INDICATION OF SINGULARITY.
C
C     BAND STORAGE
C
C           IF  A  IS A BAND MATRIX, THE FOLLOWING PROGRAM SEGMENT
C           WILL SET UP THE INPUT.
C
C                   ML = (BAND WIDTH BELOW THE DIAGONAL)
C                   MU = (BAND WIDTH ABOVE THE DIAGONAL)
C                   DO 20 I = 1, N
C                      J1 = MAX0(1, I-ML)
C                      J2 = MIN0(N, I+MU)
C                      DO 10 J = J1, J2
C                         K = J - I + ML + 1
C                         ABE(I,K) = A(I,J)
C                10    CONTINUE
C                20 CONTINUE
C
C           THIS USES COLUMNS  1  THROUGH  ML+MU+1  OF ABE .
C           FURTHERMORE,  ML  ADDITIONAL COLUMNS ARE NEEDED IN
C           ABE  STARTING WITH COLUMN  ML+MU+2  FOR ELEMENTS
C           GENERATED DURING THE TRIANGULARIZATION.  THE TOTAL
C           NUMBER OF COLUMNS NEEDED IN  ABE  IS  2*ML+MU+1 .
C
C     EXAMPLE..  IF THE ORIGINAL MATRIX IS
C
C           11 12 13  0  0  0
C           21 22 23 24  0  0
C            0 32 33 34 35  0
C            0  0 43 44 45 46
C            0  0  0 54 55 56
C            0  0  0  0 65 66
C
C      THEN  N = 6, ML = 1, MU = 2, LDA .GE. 6  AND ABE SHOULD CONTAIN
C
C            * 11 12 13  +     , * = NOT USED
C           21 22 23 24  +     , + = USED FOR PIVOTING
C           32 33 34 35  +
C           43 44 45 46  +
C           54 55 56  *  +
C           65 66  *  *  +
C
C     WRITTEN BY E.A.VOORHEES, LOS ALAMOS SCIENTIFIC LABORATORY.
C     MODIFIED BY A.H.MORRIS, NAVAL SURFACE WEAPONS CENTER.
C
C     SUBROUTINES AND FUNCTIONS
C
C      ISAMAX,SAXPY,SSCAL,SSWAP
C     FORTRAN  MIN0
C
      INTEGER LDA,N,ML,MU,IPVT(1),INFO
      REAL ABE(LDA,1)
      INTEGER ML1,MB,M,N1,LDB,I,J,K,L,LM,LM1,LM2,MP,ISAMAX
      INTEGER MJ,LMK,LL,JJ,J1
      REAL T
C
      ML1=ML+1
      MB=ML+MU
      M=ML+MU+1
      N1=N-1
      LDB=LDA-1
      INFO=0
C
C     SET FILL-IN COLUMNS TO ZERO
C
      IF(N.LE.1)GO TO 50
      IF(ML.EQ.0)GO TO 7
      DO 6 J=1,ML
        MJ=M+J
        DO 5 I=1,N
          ABE(I,MJ)=0.0E0
    5   CONTINUE
    6 CONTINUE
    7 CONTINUE
C
C     GAUSSIAN ELIMINATION WITH PARTIAL ELIMINATION
C
      DO 40 K=1,N1
        LM=MIN0(N-K,ML)
        LMK=LM+K
        LM1=LM+1
        LM2=ML1-LM
C
C     SEARCH FOR PIVOT INDEX
C
        L=-ISAMAX(LM1,ABE(LMK,LM2),LDB)+LM1+K
        IPVT(K)=L
        MP=MIN0(MB,N-K)
C
C     SWAP ROWS IF NECESSARY
C
        LL = ML1+K-L
        IF(L.NE.K)CALL SSWAP(MP+1,ABE(K,ML1),LDA,ABE(L,LL),LDA)
C
C     SKIP COLUMN REDUCTION IF PIVOT IS ZERO
C
        IF(ABE(K,ML1).EQ.0.0E0) GO TO 20
C
C     COMPUTE MULTIPLIERS
C
        IF(ML.EQ.0)GO TO 30
        T=-1.0/ABE(K,ML1)
        CALL SSCAL(LM,T,ABE(LMK,LM2),LDB)
C
C     ROW ELIMINATION WITH COLUMN INDEXING
C
        DO 10 J=1,MP
          JJ=ML1+J
          J1=LM2+J
          CALL SAXPY(LM,ABE(K,JJ),ABE(LMK,LM2),LDB,ABE(LMK,J1),LDB)
   10   CONTINUE
        GO TO 30
   20   INFO=K
   30   CONTINUE
   40 CONTINUE
   50 CONTINUE
      IPVT(N)=N
      IF(ABE(N,ML1).EQ.0.0E0) INFO=N
      RETURN
      END
