1 | SUBROUTINE STBSV ( 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 | REAL A( LDA, * ), X( * ) |
---|
7 | * .. |
---|
8 | * |
---|
9 | * Purpose |
---|
10 | * ======= |
---|
11 | * |
---|
12 | * STBSV 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 - REAL 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 - REAL 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 | REAL ZERO |
---|
143 | PARAMETER ( ZERO = 0.0E+0 ) |
---|
144 | * .. Local Scalars .. |
---|
145 | REAL 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( 'STBSV ', 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 STBSV . |
---|
345 | * |
---|
346 | END |
---|