SUBROUTINE DTBSVF(UPLO,TRANS,DIAG,N,K,A,LDA,X,INCX)

*     .. Scalar Arguments ..

      INTEGER INCX,K,LDA,N

      CHARACTER DIAG,TRANS,UPLO

*     ..

*     .. Array Arguments ..

      DOUBLE PRECISION A(LDA,*),X(*)

*     ..

*

*  Purpose

*  =======

*

*  DTBSV  solves one of the systems of equations

*

*     A*x = b,   or   A'*x = b,

*

*  where b and x are n element vectors and A is an n by n unit, or

*  non-unit, upper or lower triangular band matrix, with ( k + 1 )

*  diagonals.

*

*  No test for singularity or near-singularity is included in this

*  routine. Such tests must be performed before calling this routine.

*

*  Arguments

*  ==========

*

*  UPLO   - CHARACTER*1.

*           On entry, UPLO specifies whether the matrix is an upper or

*           lower triangular matrix as follows:

*

*              UPLO = 'U' or 'u'   A is an upper triangular matrix.

*

*              UPLO = 'L' or 'l'   A is a lower triangular matrix.

*

*           Unchanged on exit.

*

*  TRANS  - CHARACTER*1.

*           On entry, TRANS specifies the equations to be solved as

*           follows:

*

*              TRANS = 'N' or 'n'   A*x = b.

*

*              TRANS = 'T' or 't'   A'*x = b.

*

*              TRANS = 'C' or 'c'   A'*x = b.

*

*           Unchanged on exit.

*

*  DIAG   - CHARACTER*1.

*           On entry, DIAG specifies whether or not A is unit

*           triangular as follows:

*

*              DIAG = 'U' or 'u'   A is assumed to be unit triangular.

*

*              DIAG = 'N' or 'n'   A is not assumed to be unit

*                                  triangular.

*

*           Unchanged on exit.

*

*  N      - INTEGER.

*           On entry, N specifies the order of the matrix A.

*           N must be at least zero.

*           Unchanged on exit.

*

*  K      - INTEGER.

*           On entry with UPLO = 'U' or 'u', K specifies the number of

*           super-diagonals of the matrix A.

*           On entry with UPLO = 'L' or 'l', K specifies the number of

*           sub-diagonals of the matrix A.

*           K must satisfy  0 .le. K.

*           Unchanged on exit.

*

*  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).

*           Before entry with UPLO = 'U' or 'u', the leading ( k + 1 )

*           by n part of the array A must contain the upper triangular

*           band part of the matrix of coefficients, supplied column by

*           column, with the leading diagonal of the matrix in row

*           ( k + 1 ) of the array, the first super-diagonal starting at

*           position 2 in row k, and so on. The top left k by k triangle

*           of the array A is not referenced.

*           The following program segment will transfer an upper

*           triangular band matrix from conventional full matrix storage

*           to band storage:

*

*                 DO 20, J = 1, N

*                    M = K + 1 - J

*                    DO 10, I = MAX( 1, J - K ), J

*                       A( M + I, J ) = matrix( I, J )

*              10    CONTINUE

*              20 CONTINUE

*

*           Before entry with UPLO = 'L' or 'l', the leading ( k + 1 )

*           by n part of the array A must contain the lower triangular

*           band part of the matrix of coefficients, supplied column by

*           column, with the leading diagonal of the matrix in row 1 of

*           the array, the first sub-diagonal starting at position 1 in

*           row 2, and so on. The bottom right k by k triangle of the

*           array A is not referenced.

*           The following program segment will transfer a lower

*           triangular band matrix from conventional full matrix storage

*           to band storage:

*

*                 DO 20, J = 1, N

*                    M = 1 - J

*                    DO 10, I = J, MIN( N, J + K )

*                       A( M + I, J ) = matrix( I, J )

*              10    CONTINUE

*              20 CONTINUE

*

*           Note that when DIAG = 'U' or 'u' the elements of the array A

*           corresponding to the diagonal elements of the matrix are not

*           referenced, but are assumed to be unity.

*           Unchanged on exit.

*

*  LDA    - INTEGER.

*           On entry, LDA specifies the first dimension of A as declared

*           in the calling (sub) program. LDA must be at least

*           ( k + 1 ).

*           Unchanged on exit.

*

*  X      - DOUBLE PRECISION array of dimension at least

*           ( 1 + ( n - 1 )*abs( INCX ) ).

*           Before entry, the incremented array X must contain the n

*           element right-hand side vector b. On exit, X is overwritten

*           with the solution vector x.

*

*  INCX   - INTEGER.

*           On entry, INCX specifies the increment for the elements of

*           X. INCX must not be zero.

*           Unchanged on exit.

*

*

*  Level 2 Blas routine.

*

*  -- Written on 22-October-1986.

*     Jack Dongarra, Argonne National Lab.

*     Jeremy Du Croz, Nag Central Office.

*     Sven Hammarling, Nag Central Office.

*     Richard Hanson, Sandia National Labs.

*

*

*     .. Parameters ..

      DOUBLE PRECISION ZERO

      PARAMETER (ZERO=0.0D+0)

*     ..

*     .. Local Scalars ..

      DOUBLE PRECISION TEMP

      INTEGER I,INFO,IX,J,JX,KPLUS1,KX,L

      LOGICAL NOUNIT

*     ..

*     .. External Functions ..

      LOGICAL LSAME

      EXTERNAL LSAME

*     ..

*     .. External Subroutines ..

      EXTERNAL XERBLA

*     ..

*     .. Intrinsic Functions ..

      INTRINSIC MAX,MIN

*     ..

*

*     Test the input parameters.

*

      INFO = 0

      IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN

          INFO = 1

      ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.

     +         .NOT.LSAME(TRANS,'C')) THEN

          INFO = 2

      ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN

          INFO = 3

      ELSE IF (N.LT.0) THEN

          INFO = 4

      ELSE IF (K.LT.0) THEN

          INFO = 5

      ELSE IF (LDA.LT. (K+1)) THEN

          INFO = 7

      ELSE IF (INCX.EQ.0) THEN

          INFO = 9

      END IF

      IF (INFO.NE.0) THEN

          CALL XERBLA('DTBSV ',INFO)

          RETURN

      END IF

*

*     Quick return if possible.

*

      IF (N.EQ.0) RETURN

*

      NOUNIT = LSAME(DIAG,'N')

*

*     Set up the start point in X if the increment is not unity. This

*     will be  ( N - 1 )*INCX  too small for descending loops.

*

      IF (INCX.LE.0) THEN

          KX = 1 - (N-1)*INCX

      ELSE IF (INCX.NE.1) THEN

          KX = 1

      END IF

*

*     Start the operations. In this version the elements of A are

*     accessed by sequentially with one pass through A.

*

      IF (LSAME(TRANS,'N')) THEN

*

*        Form  x := inv( A )*x.

*

          IF (LSAME(UPLO,'U')) THEN

              KPLUS1 = K + 1

              IF (INCX.EQ.1) THEN

                  DO 20 J = N,1,-1

                      IF (X(J).NE.ZERO) THEN

                          L = KPLUS1 - J

                          IF (NOUNIT) X(J) = X(J)/A(KPLUS1,J)

                          TEMP = X(J)

                          DO 10 I = J - 1,MAX(1,J-K),-1

                              X(I) = X(I) - TEMP*A(L+I,J)

   10                     CONTINUE

                      END IF

   20             CONTINUE

              ELSE

                  KX = KX + (N-1)*INCX

                  JX = KX

                  DO 40 J = N,1,-1

                      KX = KX - INCX

                      IF (X(JX).NE.ZERO) THEN

                          IX = KX

                          L = KPLUS1 - J

                          IF (NOUNIT) X(JX) = X(JX)/A(KPLUS1,J)

                          TEMP = X(JX)

                          DO 30 I = J - 1,MAX(1,J-K),-1

                              X(IX) = X(IX) - TEMP*A(L+I,J)

                              IX = IX - INCX

   30                     CONTINUE

                      END IF

                      JX = JX - INCX

   40             CONTINUE

              END IF

          ELSE

              IF (INCX.EQ.1) THEN

                  DO 60 J = 1,N

                      IF (X(J).NE.ZERO) THEN

                          L = 1 - J

                          IF (NOUNIT) X(J) = X(J)/A(1,J)

                          TEMP = X(J)

                          DO 50 I = J + 1,MIN(N,J+K)

                              X(I) = X(I) - TEMP*A(L+I,J)

   50                     CONTINUE

                      END IF

   60             CONTINUE

              ELSE

                  JX = KX

                  DO 80 J = 1,N

                      KX = KX + INCX

                      IF (X(JX).NE.ZERO) THEN

                          IX = KX

                          L = 1 - J

                          IF (NOUNIT) X(JX) = X(JX)/A(1,J)

                          TEMP = X(JX)

                          DO 70 I = J + 1,MIN(N,J+K)

                              X(IX) = X(IX) - TEMP*A(L+I,J)

                              IX = IX + INCX

   70                     CONTINUE

                      END IF

                      JX = JX + INCX

   80             CONTINUE

              END IF

          END IF

      ELSE

*

*        Form  x := inv( A')*x.

*

          IF (LSAME(UPLO,'U')) THEN

              KPLUS1 = K + 1

              IF (INCX.EQ.1) THEN

                  DO 100 J = 1,N

                      TEMP = X(J)

                      L = KPLUS1 - J

                      DO 90 I = MAX(1,J-K),J - 1

                          TEMP = TEMP - A(L+I,J)*X(I)

   90                 CONTINUE

                      IF (NOUNIT) TEMP = TEMP/A(KPLUS1,J)

                      X(J) = TEMP

  100             CONTINUE

              ELSE

                  JX = KX

                  DO 120 J = 1,N

                      TEMP = X(JX)

                      IX = KX

                      L = KPLUS1 - J

                      DO 110 I = MAX(1,J-K),J - 1

                          TEMP = TEMP - A(L+I,J)*X(IX)

                          IX = IX + INCX

  110                 CONTINUE

                      IF (NOUNIT) TEMP = TEMP/A(KPLUS1,J)

                      X(JX) = TEMP

                      JX = JX + INCX

                      IF (J.GT.K) KX = KX + INCX

  120             CONTINUE

              END IF

          ELSE

              IF (INCX.EQ.1) THEN

                  DO 140 J = N,1,-1

                      TEMP = X(J)

                      L = 1 - J

                      DO 130 I = MIN(N,J+K),J + 1,-1

                          TEMP = TEMP - A(L+I,J)*X(I)

  130                 CONTINUE

                      IF (NOUNIT) TEMP = TEMP/A(1,J)

                      X(J) = TEMP

  140             CONTINUE

              ELSE

                  KX = KX + (N-1)*INCX

                  JX = KX

                  DO 160 J = N,1,-1

                      TEMP = X(JX)

                      IX = KX

                      L = 1 - J

                      DO 150 I = MIN(N,J+K),J + 1,-1

                          TEMP = TEMP - A(L+I,J)*X(IX)

                          IX = IX - INCX

  150                 CONTINUE

                      IF (NOUNIT) TEMP = TEMP/A(1,J)

                      X(JX) = TEMP

                      JX = JX - INCX

                      IF ((N-J).GE.K) KX = KX - INCX

  160             CONTINUE

              END IF

          END IF

      END IF

*

      RETURN

*

*     End of DTBSV .

*

      END