1*> \brief \b STBSV 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 STBSV(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*> STBSV solves one of the systems of equations 28*> 29*> A*x = b, or A**T*x = b, 30*> 31*> where b and x are n element vectors and A is an n by n unit, or 32*> non-unit, upper or lower triangular band matrix, with ( k + 1 ) 33*> diagonals. 34*> 35*> No test for singularity or near-singularity is included in this 36*> routine. Such tests must be performed before calling this routine. 37*> \endverbatim 38* 39* Arguments: 40* ========== 41* 42*> \param[in] UPLO 43*> \verbatim 44*> UPLO is CHARACTER*1 45*> On entry, UPLO specifies whether the matrix is an upper or 46*> lower triangular matrix as follows: 47*> 48*> UPLO = 'U' or 'u' A is an upper triangular matrix. 49*> 50*> UPLO = 'L' or 'l' A is a lower triangular matrix. 51*> \endverbatim 52*> 53*> \param[in] TRANS 54*> \verbatim 55*> TRANS is CHARACTER*1 56*> On entry, TRANS specifies the equations to be solved as 57*> follows: 58*> 59*> TRANS = 'N' or 'n' A*x = b. 60*> 61*> TRANS = 'T' or 't' A**T*x = b. 62*> 63*> TRANS = 'C' or 'c' A**T*x = b. 64*> \endverbatim 65*> 66*> \param[in] DIAG 67*> \verbatim 68*> DIAG is CHARACTER*1 69*> On entry, DIAG specifies whether or not A is unit 70*> triangular as follows: 71*> 72*> DIAG = 'U' or 'u' A is assumed to be unit triangular. 73*> 74*> DIAG = 'N' or 'n' A is not assumed to be unit 75*> triangular. 76*> \endverbatim 77*> 78*> \param[in] N 79*> \verbatim 80*> N is INTEGER 81*> On entry, N specifies the order of the matrix A. 82*> N must be at least zero. 83*> \endverbatim 84*> 85*> \param[in] K 86*> \verbatim 87*> K is INTEGER 88*> On entry with UPLO = 'U' or 'u', K specifies the number of 89*> super-diagonals of the matrix A. 90*> On entry with UPLO = 'L' or 'l', K specifies the number of 91*> sub-diagonals of the matrix A. 92*> K must satisfy 0 .le. K. 93*> \endverbatim 94*> 95*> \param[in] A 96*> \verbatim 97*> A is REAL array, dimension ( LDA, N ) 98*> Before entry with UPLO = 'U' or 'u', the leading ( k + 1 ) 99*> by n part of the array A must contain the upper triangular 100*> band part of the matrix of coefficients, supplied column by 101*> column, with the leading diagonal of the matrix in row 102*> ( k + 1 ) of the array, the first super-diagonal starting at 103*> position 2 in row k, and so on. The top left k by k triangle 104*> of the array A is not referenced. 105*> The following program segment will transfer an upper 106*> triangular band matrix from conventional full matrix storage 107*> to band storage: 108*> 109*> DO 20, J = 1, N 110*> M = K + 1 - J 111*> DO 10, I = MAX( 1, J - K ), J 112*> A( M + I, J ) = matrix( I, J ) 113*> 10 CONTINUE 114*> 20 CONTINUE 115*> 116*> Before entry with UPLO = 'L' or 'l', the leading ( k + 1 ) 117*> by n part of the array A must contain the lower triangular 118*> band part of the matrix of coefficients, supplied column by 119*> column, with the leading diagonal of the matrix in row 1 of 120*> the array, the first sub-diagonal starting at position 1 in 121*> row 2, and so on. The bottom right k by k triangle of the 122*> array A is not referenced. 123*> The following program segment will transfer a lower 124*> triangular band matrix from conventional full matrix storage 125*> to band storage: 126*> 127*> DO 20, J = 1, N 128*> M = 1 - J 129*> DO 10, I = J, MIN( N, J + K ) 130*> A( M + I, J ) = matrix( I, J ) 131*> 10 CONTINUE 132*> 20 CONTINUE 133*> 134*> Note that when DIAG = 'U' or 'u' the elements of the array A 135*> corresponding to the diagonal elements of the matrix are not 136*> referenced, but are assumed to be unity. 137*> \endverbatim 138*> 139*> \param[in] LDA 140*> \verbatim 141*> LDA is INTEGER 142*> On entry, LDA specifies the first dimension of A as declared 143*> in the calling (sub) program. LDA must be at least 144*> ( k + 1 ). 145*> \endverbatim 146*> 147*> \param[in,out] X 148*> \verbatim 149*> X is REAL array, dimension at least 150*> ( 1 + ( n - 1 )*abs( INCX ) ). 151*> Before entry, the incremented array X must contain the n 152*> element right-hand side vector b. On exit, X is overwritten 153*> with the solution vector x. 154*> \endverbatim 155*> 156*> \param[in] INCX 157*> \verbatim 158*> INCX is INTEGER 159*> On entry, INCX specifies the increment for the elements of 160*> X. INCX must not be zero. 161*> \endverbatim 162* 163* Authors: 164* ======== 165* 166*> \author Univ. of Tennessee 167*> \author Univ. of California Berkeley 168*> \author Univ. of Colorado Denver 169*> \author NAG Ltd. 170* 171*> \ingroup single_blas_level2 172* 173*> \par Further Details: 174* ===================== 175*> 176*> \verbatim 177*> 178*> Level 2 Blas routine. 179*> 180*> -- Written on 22-October-1986. 181*> Jack Dongarra, Argonne National Lab. 182*> Jeremy Du Croz, Nag Central Office. 183*> Sven Hammarling, Nag Central Office. 184*> Richard Hanson, Sandia National Labs. 185*> \endverbatim 186*> 187* ===================================================================== 188 SUBROUTINE STBSV(UPLO,TRANS,DIAG,N,K,A,LDA,X,INCX) 189* 190* -- Reference BLAS level2 routine -- 191* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- 192* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 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('STBSV ',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 by sequentially with one pass through A. 265* 266 IF (LSAME(TRANS,'N')) THEN 267* 268* Form x := inv( A )*x. 269* 270 IF (LSAME(UPLO,'U')) THEN 271 KPLUS1 = K + 1 272 IF (INCX.EQ.1) THEN 273 DO 20 J = N,1,-1 274 IF (X(J).NE.ZERO) THEN 275 L = KPLUS1 - J 276 IF (NOUNIT) X(J) = X(J)/A(KPLUS1,J) 277 TEMP = X(J) 278 DO 10 I = J - 1,MAX(1,J-K),-1 279 X(I) = X(I) - TEMP*A(L+I,J) 280 10 CONTINUE 281 END IF 282 20 CONTINUE 283 ELSE 284 KX = KX + (N-1)*INCX 285 JX = KX 286 DO 40 J = N,1,-1 287 KX = KX - INCX 288 IF (X(JX).NE.ZERO) THEN 289 IX = KX 290 L = KPLUS1 - J 291 IF (NOUNIT) X(JX) = X(JX)/A(KPLUS1,J) 292 TEMP = X(JX) 293 DO 30 I = J - 1,MAX(1,J-K),-1 294 X(IX) = X(IX) - TEMP*A(L+I,J) 295 IX = IX - INCX 296 30 CONTINUE 297 END IF 298 JX = JX - INCX 299 40 CONTINUE 300 END IF 301 ELSE 302 IF (INCX.EQ.1) THEN 303 DO 60 J = 1,N 304 IF (X(J).NE.ZERO) THEN 305 L = 1 - J 306 IF (NOUNIT) X(J) = X(J)/A(1,J) 307 TEMP = X(J) 308 DO 50 I = J + 1,MIN(N,J+K) 309 X(I) = X(I) - TEMP*A(L+I,J) 310 50 CONTINUE 311 END IF 312 60 CONTINUE 313 ELSE 314 JX = KX 315 DO 80 J = 1,N 316 KX = KX + INCX 317 IF (X(JX).NE.ZERO) THEN 318 IX = KX 319 L = 1 - J 320 IF (NOUNIT) X(JX) = X(JX)/A(1,J) 321 TEMP = X(JX) 322 DO 70 I = J + 1,MIN(N,J+K) 323 X(IX) = X(IX) - TEMP*A(L+I,J) 324 IX = IX + INCX 325 70 CONTINUE 326 END IF 327 JX = JX + INCX 328 80 CONTINUE 329 END IF 330 END IF 331 ELSE 332* 333* Form x := inv( A**T)*x. 334* 335 IF (LSAME(UPLO,'U')) THEN 336 KPLUS1 = K + 1 337 IF (INCX.EQ.1) THEN 338 DO 100 J = 1,N 339 TEMP = X(J) 340 L = KPLUS1 - J 341 DO 90 I = MAX(1,J-K),J - 1 342 TEMP = TEMP - A(L+I,J)*X(I) 343 90 CONTINUE 344 IF (NOUNIT) TEMP = TEMP/A(KPLUS1,J) 345 X(J) = TEMP 346 100 CONTINUE 347 ELSE 348 JX = KX 349 DO 120 J = 1,N 350 TEMP = X(JX) 351 IX = KX 352 L = KPLUS1 - J 353 DO 110 I = MAX(1,J-K),J - 1 354 TEMP = TEMP - A(L+I,J)*X(IX) 355 IX = IX + INCX 356 110 CONTINUE 357 IF (NOUNIT) TEMP = TEMP/A(KPLUS1,J) 358 X(JX) = TEMP 359 JX = JX + INCX 360 IF (J.GT.K) KX = KX + INCX 361 120 CONTINUE 362 END IF 363 ELSE 364 IF (INCX.EQ.1) THEN 365 DO 140 J = N,1,-1 366 TEMP = X(J) 367 L = 1 - J 368 DO 130 I = MIN(N,J+K),J + 1,-1 369 TEMP = TEMP - A(L+I,J)*X(I) 370 130 CONTINUE 371 IF (NOUNIT) TEMP = TEMP/A(1,J) 372 X(J) = TEMP 373 140 CONTINUE 374 ELSE 375 KX = KX + (N-1)*INCX 376 JX = KX 377 DO 160 J = N,1,-1 378 TEMP = X(JX) 379 IX = KX 380 L = 1 - J 381 DO 150 I = MIN(N,J+K),J + 1,-1 382 TEMP = TEMP - A(L+I,J)*X(IX) 383 IX = IX - INCX 384 150 CONTINUE 385 IF (NOUNIT) TEMP = TEMP/A(1,J) 386 X(JX) = TEMP 387 JX = JX - INCX 388 IF ((N-J).GE.K) KX = KX - INCX 389 160 CONTINUE 390 END IF 391 END IF 392 END IF 393* 394 RETURN 395* 396* End of STBSV 397* 398 END 399