BLAS/SRC/stbmv.f

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