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