share/lapack/blas/fortran/dtbsv.f

   1       SUBROUTINE DTBSV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX )
   2 *     .. Scalar Arguments ..
   3       INTEGER            INCX, K, LDA, N
   4       CHARACTER*1        DIAG, TRANS, UPLO
   5 *     .. Array Arguments ..
   6       DOUBLE PRECISION   A( LDA, * ), X( * )
   7 *     ..
   8 *
   9 *  Purpose
  10 *  =======
  11 *
  12 *  DTBSV  solves one of the systems of equations
  13 *
  14 *     A*x = b,   or   A'*x = b,
  15 *
  16 *  where b and x are n element vectors and A is an n by n unit, or
  17 *  non-unit, upper or lower triangular band matrix, with ( k + 1 )
  18 *  diagonals.
  19 *
  20 *  No test for singularity or near-singularity is included in this
  21 *  routine. Such tests must be performed before calling this routine.
  22 *
  23 *  Parameters
  24 *  ==========
  25 *
  26 *  UPLO   - CHARACTER*1.
  27 *           On entry, UPLO specifies whether the matrix is an upper or
  28 *           lower triangular matrix as follows:
  29 *
  30 *              UPLO = 'U' or 'u'   A is an upper triangular matrix.
  31 *
  32 *              UPLO = 'L' or 'l'   A is a lower triangular matrix.
  33 *
  34 *           Unchanged on exit.
  35 *
  36 *  TRANS  - CHARACTER*1.
  37 *           On entry, TRANS specifies the equations to be solved as
  38 *           follows:
  39 *
  40 *              TRANS = 'N' or 'n'   A*x = b.
  41 *
  42 *              TRANS = 'T' or 't'   A'*x = b.
  43 *
  44 *              TRANS = 'C' or 'c'   A'*x = b.
  45 *
  46 *           Unchanged on exit.
  47 *
  48 *  DIAG   - CHARACTER*1.
  49 *           On entry, DIAG specifies whether or not A is unit
  50 *           triangular as follows:
  51 *
  52 *              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
  53 *
  54 *              DIAG = 'N' or 'n'   A is not assumed to be unit
  55 *                                  triangular.
  56 *
  57 *           Unchanged on exit.
  58 *
  59 *  N      - INTEGER.
  60 *           On entry, N specifies the order of the matrix A.
  61 *           N must be at least zero.
  62 *           Unchanged on exit.
  63 *
  64 *  K      - INTEGER.
  65 *           On entry with UPLO = 'U' or 'u', K specifies the number of
  66 *           super-diagonals of the matrix A.
  67 *           On entry with UPLO = 'L' or 'l', K specifies the number of
  68 *           sub-diagonals of the matrix A.
  69 *           K must satisfy  0 .le. K.
  70 *           Unchanged on exit.
  71 *
  72 *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
  73 *           Before entry with UPLO = 'U' or 'u', the leading ( k + 1 )
  74 *           by n part of the array A must contain the upper triangular
  75 *           band part of the matrix of coefficients, supplied column by
  76 *           column, with the leading diagonal of the matrix in row
  77 *           ( k + 1 ) of the array, the first super-diagonal starting at
  78 *           position 2 in row k, and so on. The top left k by k triangle
  79 *           of the array A is not referenced.
  80 *           The following program segment will transfer an upper
  81 *           triangular band matrix from conventional full matrix storage
  82 *           to band storage:
  83 *
  84 *                 DO 20, J = 1, N
  85 *                    M = K + 1 - J
  86 *                    DO 10, I = MAX( 1, J - K ), J
  87 *                       A( M + I, J ) = matrix( I, J )
  88 *              10    CONTINUE
  89 *              20 CONTINUE
  90 *
  91 *           Before entry with UPLO = 'L' or 'l', the leading ( k + 1 )
  92 *           by n part of the array A must contain the lower triangular
  93 *           band part of the matrix of coefficients, supplied column by
  94 *           column, with the leading diagonal of the matrix in row 1 of
  95 *           the array, the first sub-diagonal starting at position 1 in
  96 *           row 2, and so on. The bottom right k by k triangle of the
  97 *           array A is not referenced.
  98 *           The following program segment will transfer a lower
  99 *           triangular band matrix from conventional full matrix storage
 100 *           to band storage:
 101 *
 102 *                 DO 20, J = 1, N
 103 *                    M = 1 - J
 104 *                    DO 10, I = J, MIN( N, J + K )
 105 *                       A( M + I, J ) = matrix( I, J )
 106 *              10    CONTINUE
 107 *              20 CONTINUE
 108 *
 109 *           Note that when DIAG = 'U' or 'u' the elements of the array A
 110 *           corresponding to the diagonal elements of the matrix are not
 111 *           referenced, but are assumed to be unity.
 112 *           Unchanged on exit.
 113 *
 114 *  LDA    - INTEGER.
 115 *           On entry, LDA specifies the first dimension of A as declared
 116 *           in the calling (sub) program. LDA must be at least
 117 *           ( k + 1 ).
 118 *           Unchanged on exit.
 119 *
 120 *  X      - DOUBLE PRECISION array of dimension at least
 121 *           ( 1 + ( n - 1 )*abs( INCX ) ).
 122 *           Before entry, the incremented array X must contain the n
 123 *           element right-hand side vector b. On exit, X is overwritten
 124 *           with the solution vector x.
 125 *
 126 *  INCX   - INTEGER.
 127 *           On entry, INCX specifies the increment for the elements of
 128 *           X. INCX must not be zero.
 129 *           Unchanged on exit.
 130 *
 131 *
 132 *  Level 2 Blas routine.
 133 *
 134 *  -- Written on 22-October-1986.
 135 *     Jack Dongarra, Argonne National Lab.
 136 *     Jeremy Du Croz, Nag Central Office.
 137 *     Sven Hammarling, Nag Central Office.
 138 *     Richard Hanson, Sandia National Labs.
 139 *
 140 *
 141 *     .. Parameters ..
 142       DOUBLE PRECISION   ZERO
 143       PARAMETER        ( ZERO = 0.0D+0 )
 144 *     .. Local Scalars ..
 145       DOUBLE PRECISION   TEMP
 146       INTEGER            I, INFO, IX, J, JX, KPLUS1, KX, L
 147       LOGICAL            NOUNIT
 148 *     .. External Functions ..
 149       LOGICAL            LSAME
 150       EXTERNAL           LSAME
 151 *     .. External Subroutines ..
 152       EXTERNAL           XERBLA
 153 *     .. Intrinsic Functions ..
 154       INTRINSIC          MAX, MIN
 155 *     ..
 156 *     .. Executable Statements ..
 157 *
 158 *     Test the input parameters.
 159 *
 160       INFO = 0
 161       IF     ( .NOT.LSAME( UPLO , 'U' ).AND.
 162      $         .NOT.LSAME( UPLO , 'L' )      )THEN
 163          INFO = 1
 164       ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND.
 165      $         .NOT.LSAME( TRANS, 'T' ).AND.
 166      $         .NOT.LSAME( TRANS, 'C' )      )THEN
 167          INFO = 2
 168       ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND.
 169      $         .NOT.LSAME( DIAG , 'N' )      )THEN
 170          INFO = 3
 171       ELSE IF( N.LT.0 )THEN
 172          INFO = 4
 173       ELSE IF( K.LT.0 )THEN
 174          INFO = 5
 175       ELSE IF( LDA.LT.( K + 1 ) )THEN
 176          INFO = 7
 177       ELSE IF( INCX.EQ.0 )THEN
 178          INFO = 9
 179       END IF
 180       IF( INFO.NE.0 )THEN
 181          CALL XERBLA( 'DTBSV ', INFO )
 182          RETURN
 183       END IF
 184 *
 185 *     Quick return if possible.
 186 *
 187       IF( N.EQ.0 )
 188      $   RETURN
 189 *
 190       NOUNIT = LSAME( DIAG, 'N' )
 191 *
 192 *     Set up the start point in X if the increment is not unity. This
 193 *     will be  ( N - 1 )*INCX  too small for descending loops.
 194 *
 195       IF( INCX.LE.0 )THEN
 196          KX = 1 - ( N - 1 )*INCX
 197       ELSE IF( INCX.NE.1 )THEN
 198          KX = 1
 199       END IF
 200 *
 201 *     Start the operations. In this version the elements of A are
 202 *     accessed by sequentially with one pass through A.
 203 *
 204       IF( LSAME( TRANS, 'N' ) )THEN
 205 *
 206 *        Form  x := inv( A )*x.
 207 *
 208          IF( LSAME( UPLO, 'U' ) )THEN
 209             KPLUS1 = K + 1
 210             IF( INCX.EQ.1 )THEN
 211                DO 20, J = N, 1, -1
 212                   IF( X( J ).NE.ZERO )THEN
 213                      L = KPLUS1 - J
 214                      IF( NOUNIT )
 215      $                  X( J ) = X( J )/A( KPLUS1, J )
 216                      TEMP = X( J )
 217                      DO 10, I = J - 1, MAX( 1, J - K ), -1
 218                         X( I ) = X( I ) - TEMP*A( L + I, J )
 219    10                CONTINUE
 220                   END IF
 221    20          CONTINUE
 222             ELSE
 223                KX = KX + ( N - 1 )*INCX
 224                JX = KX
 225                DO 40, J = N, 1, -1
 226                   KX = KX - INCX
 227                   IF( X( JX ).NE.ZERO )THEN
 228                      IX = KX
 229                      L  = KPLUS1 - J
 230                      IF( NOUNIT )
 231      $                  X( JX ) = X( JX )/A( KPLUS1, J )
 232                      TEMP = X( JX )
 233                      DO 30, I = J - 1, MAX( 1, J - K ), -1
 234                         X( IX ) = X( IX ) - TEMP*A( L + I, J )
 235                         IX      = IX      - INCX
 236    30                CONTINUE
 237                   END IF
 238                   JX = JX - INCX
 239    40          CONTINUE
 240             END IF
 241          ELSE
 242             IF( INCX.EQ.1 )THEN
 243                DO 60, J = 1, N
 244                   IF( X( J ).NE.ZERO )THEN
 245                      L = 1 - J
 246                      IF( NOUNIT )
 247      $                  X( J ) = X( J )/A( 1, J )
 248                      TEMP = X( J )
 249                      DO 50, I = J + 1, MIN( N, J + K )
 250                         X( I ) = X( I ) - TEMP*A( L + I, J )
 251    50                CONTINUE
 252                   END IF
 253    60          CONTINUE
 254             ELSE
 255                JX = KX
 256                DO 80, J = 1, N
 257                   KX = KX + INCX
 258                   IF( X( JX ).NE.ZERO )THEN
 259                      IX = KX
 260                      L  = 1  - J
 261                      IF( NOUNIT )
 262      $                  X( JX ) = X( JX )/A( 1, J )
 263                      TEMP = X( JX )
 264                      DO 70, I = J + 1, MIN( N, J + K )
 265                         X( IX ) = X( IX ) - TEMP*A( L + I, J )
 266                         IX      = IX      + INCX
 267    70                CONTINUE
 268                   END IF
 269                   JX = JX + INCX
 270    80          CONTINUE
 271             END IF
 272          END IF
 273       ELSE
 274 *
 275 *        Form  x := inv( A')*x.
 276 *
 277          IF( LSAME( UPLO, 'U' ) )THEN
 278             KPLUS1 = K + 1
 279             IF( INCX.EQ.1 )THEN
 280                DO 100, J = 1, N
 281                   TEMP = X( J )
 282                   L    = KPLUS1 - J
 283                   DO 90, I = MAX( 1, J - K ), J - 1
 284                      TEMP = TEMP - A( L + I, J )*X( I )
 285    90             CONTINUE
 286                   IF( NOUNIT )
 287      $               TEMP = TEMP/A( KPLUS1, J )
 288                   X( J ) = TEMP
 289   100          CONTINUE
 290             ELSE
 291                JX = KX
 292                DO 120, J = 1, N
 293                   TEMP = X( JX )
 294                   IX   = KX
 295                   L    = KPLUS1  - J
 296                   DO 110, I = MAX( 1, J - K ), J - 1
 297                      TEMP = TEMP - A( L + I, J )*X( IX )
 298                      IX   = IX   + INCX
 299   110             CONTINUE
 300                   IF( NOUNIT )
 301      $               TEMP = TEMP/A( KPLUS1, J )
 302                   X( JX ) = TEMP
 303                   JX      = JX   + INCX
 304                   IF( J.GT.K )
 305      $               KX = KX + INCX
 306   120          CONTINUE
 307             END IF
 308          ELSE
 309             IF( INCX.EQ.1 )THEN
 310                DO 140, J = N, 1, -1
 311                   TEMP = X( J )
 312                   L    = 1      - J
 313                   DO 130, I = MIN( N, J + K ), J + 1, -1
 314                      TEMP = TEMP - A( L + I, J )*X( I )
 315   130             CONTINUE
 316                   IF( NOUNIT )
 317      $               TEMP = TEMP/A( 1, J )
 318                   X( J ) = TEMP
 319   140          CONTINUE
 320             ELSE
 321                KX = KX + ( N - 1 )*INCX
 322                JX = KX
 323                DO 160, J = N, 1, -1
 324                   TEMP = X( JX )
 325                   IX   = KX
 326                   L    = 1       - J
 327                   DO 150, I = MIN( N, J + K ), J + 1, -1
 328                      TEMP = TEMP - A( L + I, J )*X( IX )
 329                      IX   = IX   - INCX
 330   150             CONTINUE
 331                   IF( NOUNIT )
 332      $               TEMP = TEMP/A( 1, J )
 333                   X( JX ) = TEMP
 334                   JX      = JX   - INCX
 335                   IF( ( N - J ).GE.K )
 336      $               KX = KX - INCX
 337   160          CONTINUE
 338             END IF
 339          END IF
 340       END IF
 341 *
 342       RETURN
 343 *
 344 *     End of DTBSV .
 345 *
 346       END