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