      SUBROUTINE BANFAC ( W, NROWW, NROW, NBANDL, NBANDU, IFLAG )
C  FROM  * A PRACTICAL GUIDE TO SPLINES *  BY C. DE BOOR
C  RETURNS IN  W  THE LU-FACTORIZATION (WITHOUT PIVOTING) OF THE BANDED
C  MATRIX  A  OF ORDER  NROW  WITH  (NBANDL + 1 + NBANDU) BANDS OR DIAG-
C  ONALS IN THE WORK ARRAY  W .
C
C******  I N P U T  ******
C  W.....WORK ARRAY OF SIZE  (NROWW,NROW)  CONTAINING THE INTERESTING
C        PART OF A BANDED MATRIX  A , WITH THE DIAGONALS OR BANDS OF  A
C        STORED IN THE ROWS OF  W , WHILE COLUMNS OF  A  CORRESPOND TO
C        COLUMNS OF  W . THIS IS THE STORAGE MODE USED IN  LINPACK  AND
C        RESULTS IN EFFICIENT INNERMOST LOOPS.
C           EXPLICITLY,  A  HAS  NBANDL  BANDS BELOW THE DIAGONAL
C                            +     1     (MAIN) DIAGONAL
C                            +   NBANDU  BANDS ABOVE THE DIAGONAL
C        AND THUS, WITH    MIDDLE = NBANDU + 1,
C          A(I+J,J)  IS IN  W(I+MIDDLE,J)  FOR I=-NBANDU,...,NBANDL
C                                              J=1,...,NROW .
C        FOR EXAMPLE, THE INTERESTING ENTRIES OF A (1,2)-BANDED MATRIX
C        OF ORDER  9  WOULD APPEAR IN THE FIRST  1+1+2 = 4  ROWS OF  W
C        AS FOLLOWS.
C                          13 24 35 46 57 68 79
C                       12 23 34 45 56 67 78 89
C                    11 22 33 44 55 66 77 88 99
C                    21 32 43 54 65 76 87 98
C
C        ALL OTHER ENTRIES OF  W  NOT IDENTIFIED IN THIS WAY WITH AN EN-
C        TRY OF  A  ARE NEVER REFERENCED .
C  NROWW.....ROW DIMENSION OF THE WORK ARRAY  W .
C        MUST BE  .GE.  NBANDL + 1 + NBANDU  .
C  NBANDL.....NUMBER OF BANDS OF  A  BELOW THE MAIN DIAGONAL
C  NBANDU.....NUMBER OF BANDS OF  A  ABOVE THE MAIN DIAGONAL .
C
C******  O U T P U T  ******
C  IFLAG.....INTEGER INDICATING SUCCESS( = 1) OR FAILURE ( = 2) .
C     IF  IFLAG = 1, THEN
C  W.....CONTAINS THE LU-FACTORIZATION OF  A  INTO A UNIT LOWER TRIANGU-
C        LAR MATRIX  L  AND AN UPPER TRIANGULAR MATRIX  U (BOTH BANDED)
C        AND STORED IN CUSTOMARY FASHION OVER THE CORRESPONDING ENTRIES
C        OF  A . THIS MAKES IT POSSIBLE TO SOLVE ANY PARTICULAR LINEAR
C        SYSTEM  A*X = B  FOR  X  BY A
C              CALL BANSLV ( W, NROWW, NROW, NBANDL, NBANDU, B )
C        WITH THE SOLUTION X  CONTAINED IN  B  ON RETURN .
C     IF  IFLAG = 2, THEN
C        ONE OF  NROW-1, NBANDL,NBANDU FAILED TO BE NONNEGATIVE, OR ELSE
C        ONE OF THE POTENTIAL PIVOTS WAS FOUND TO BE ZERO INDICATING
C        THAT  A  DOES NOT HAVE AN LU-FACTORIZATION. THIS IMPLIES THAT
C        A  IS SINGULAR IN CASE IT IS TOTALLY POSITIVE .
C
C******  M E T H O D  ******
C     GAUSS ELIMINATION  W I T H O U T  PIVOTING IS USED. THE ROUTINE IS
C  INTENDED FOR USE WITH MATRICES  A  WHICH DO NOT REQUIRE ROW INTER-
C  CHANGES DURING FACTORIZATION, ESPECIALLY FOR THE  T O T A L L Y
C  P O S I T I V E  MATRICES WHICH OCCUR IN SPLINE CALCULATIONS.
C     THE ROUTINE SHOULD NOT BE USED FOR AN ARBITRARY BANDED MATRIX.
C
      REAL W(NROWW,NROW),   FACTOR,PIVOT
C
      IFLAG = 1
      MIDDLE = NBANDU + 1
C                         W(MIDDLE,.) CONTAINS THE MAIN DIAGONAL OF  A .
      NROWM1 = NROW - 1
      IF (NROWM1)                       999,900,1
    1 IF (NBANDL .GT. 0)                GO TO 10
C                A IS UPPER TRIANGULAR. CHECK THAT DIAGONAL IS NONZERO .
      DO 5 I=1,NROWM1
         IF (W(MIDDLE,I) .EQ. 0.)       GO TO 999
    5    CONTINUE
                                        GO TO 900
   10 IF (NBANDU .GT. 0)                GO TO 20
C              A IS LOWER TRIANGULAR. CHECK THAT DIAGONAL IS NONZERO AND
C                 DIVIDE EACH COLUMN BY ITS DIAGONAL .
      DO 15 I=1,NROWM1
         PIVOT = W(MIDDLE,I)
         IF(PIVOT .EQ. 0.)              GO TO 999
         JMAX = MIN0(NBANDL, NROW - I)
         JBEG = MIDDLE + 1
         JEND = MIDDLE + JMAX
         DO 15 J=JBEG,JEND
   15       W(J,I) = W(J,I)/PIVOT
                                        GO TO 900
C
C        A  IS NOT JUST A TRIANGULAR MATRIX. CONSTRUCT LU FACTORIZATION
   20 DO 50 I=1,NROWM1
C                                  W(MIDDLE,I)  IS PIVOT FOR I-TH STEP .
         PIVOT = W(MIDDLE,I)
         IF (PIVOT .EQ. 0.)             GO TO 999
C                 JMAX  IS THE NUMBER OF (NONZERO) ENTRIES IN COLUMN  I
C                     BELOW THE DIAGONAL .
         JMAX = MIN0(NBANDL,NROW - I)
C              DIVIDE EACH ENTRY IN COLUMN  I  BELOW DIAGONAL BY PIVOT .
         JBEG = MIDDLE + 1
         JEND = MIDDLE + JMAX
         DO 32 J=JBEG,JEND
   32       W(J,I) = W(J,I)/PIVOT
C                 KMAX  IS THE NUMBER OF (NONZERO) ENTRIES IN ROW  I  TO
C                     THE RIGHT OF THE DIAGONAL .
         KMAX = MIN0(NBANDU,NROW - I)
C                  SUBTRACT  A(I,I+K)*(I-TH COLUMN) FROM (I+K)-TH COLUMN
C                  (BELOW ROW  I ) .
         DO 40 K=1,KMAX
            IPK = I + K
            MIDMK = MIDDLE - K
            FACTOR = W(MIDMK,IPK)
            DO 40 J=1,JMAX
               MJ = MIDDLE + J
               MDJ = MIDMK + J
   40          W(MDJ,IPK) = W(MDJ,IPK) - W(MJ,I)*FACTOR
   50    CONTINUE
C                                       CHECK THE LAST DIAGONAL ENTRY .
  900 IF (W(MIDDLE,NROW) .NE. 0.)       RETURN
  999 IFLAG = 2
                                        RETURN
      END
