|
5164
|
1 SUBROUTINE ZTBSV ( 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 COMPLEX*16 A( LDA, * ), X( * ) |
|
|
7 * .. |
|
|
8 * |
|
|
9 * Purpose |
|
|
10 * ======= |
|
|
11 * |
|
|
12 * ZTBSV solves one of the systems of equations |
|
|
13 * |
|
|
14 * A*x = b, or A'*x = b, or conjg( 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' conjg( 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 - COMPLEX*16 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 - COMPLEX*16 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 COMPLEX*16 ZERO |
|
|
143 PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) |
|
|
144 * .. Local Scalars .. |
|
|
145 COMPLEX*16 TEMP |
|
|
146 INTEGER I, INFO, IX, J, JX, KPLUS1, KX, L |
|
|
147 LOGICAL NOCONJ, NOUNIT |
|
|
148 * .. External Functions .. |
|
|
149 LOGICAL LSAME |
|
|
150 EXTERNAL LSAME |
|
|
151 * .. External Subroutines .. |
|
|
152 EXTERNAL XERBLA |
|
|
153 * .. Intrinsic Functions .. |
|
|
154 INTRINSIC DCONJG, 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( 'ZTBSV ', INFO ) |
|
|
182 RETURN |
|
|
183 END IF |
|
|
184 * |
|
|
185 * Quick return if possible. |
|
|
186 * |
|
|
187 IF( N.EQ.0 ) |
|
|
188 $ RETURN |
|
|
189 * |
|
|
190 NOCONJ = LSAME( TRANS, 'T' ) |
|
|
191 NOUNIT = LSAME( DIAG , 'N' ) |
|
|
192 * |
|
|
193 * Set up the start point in X if the increment is not unity. This |
|
|
194 * will be ( N - 1 )*INCX too small for descending loops. |
|
|
195 * |
|
|
196 IF( INCX.LE.0 )THEN |
|
|
197 KX = 1 - ( N - 1 )*INCX |
|
|
198 ELSE IF( INCX.NE.1 )THEN |
|
|
199 KX = 1 |
|
|
200 END IF |
|
|
201 * |
|
|
202 * Start the operations. In this version the elements of A are |
|
|
203 * accessed by sequentially with one pass through A. |
|
|
204 * |
|
|
205 IF( LSAME( TRANS, 'N' ) )THEN |
|
|
206 * |
|
|
207 * Form x := inv( A )*x. |
|
|
208 * |
|
|
209 IF( LSAME( UPLO, 'U' ) )THEN |
|
|
210 KPLUS1 = K + 1 |
|
|
211 IF( INCX.EQ.1 )THEN |
|
|
212 DO 20, J = N, 1, -1 |
|
|
213 IF( X( J ).NE.ZERO )THEN |
|
|
214 L = KPLUS1 - J |
|
|
215 IF( NOUNIT ) |
|
|
216 $ X( J ) = X( J )/A( KPLUS1, J ) |
|
|
217 TEMP = X( J ) |
|
|
218 DO 10, I = J - 1, MAX( 1, J - K ), -1 |
|
|
219 X( I ) = X( I ) - TEMP*A( L + I, J ) |
|
|
220 10 CONTINUE |
|
|
221 END IF |
|
|
222 20 CONTINUE |
|
|
223 ELSE |
|
|
224 KX = KX + ( N - 1 )*INCX |
|
|
225 JX = KX |
|
|
226 DO 40, J = N, 1, -1 |
|
|
227 KX = KX - INCX |
|
|
228 IF( X( JX ).NE.ZERO )THEN |
|
|
229 IX = KX |
|
|
230 L = KPLUS1 - J |
|
|
231 IF( NOUNIT ) |
|
|
232 $ X( JX ) = X( JX )/A( KPLUS1, J ) |
|
|
233 TEMP = X( JX ) |
|
|
234 DO 30, I = J - 1, MAX( 1, J - K ), -1 |
|
|
235 X( IX ) = X( IX ) - TEMP*A( L + I, J ) |
|
|
236 IX = IX - INCX |
|
|
237 30 CONTINUE |
|
|
238 END IF |
|
|
239 JX = JX - INCX |
|
|
240 40 CONTINUE |
|
|
241 END IF |
|
|
242 ELSE |
|
|
243 IF( INCX.EQ.1 )THEN |
|
|
244 DO 60, J = 1, N |
|
|
245 IF( X( J ).NE.ZERO )THEN |
|
|
246 L = 1 - J |
|
|
247 IF( NOUNIT ) |
|
|
248 $ X( J ) = X( J )/A( 1, J ) |
|
|
249 TEMP = X( J ) |
|
|
250 DO 50, I = J + 1, MIN( N, J + K ) |
|
|
251 X( I ) = X( I ) - TEMP*A( L + I, J ) |
|
|
252 50 CONTINUE |
|
|
253 END IF |
|
|
254 60 CONTINUE |
|
|
255 ELSE |
|
|
256 JX = KX |
|
|
257 DO 80, J = 1, N |
|
|
258 KX = KX + INCX |
|
|
259 IF( X( JX ).NE.ZERO )THEN |
|
|
260 IX = KX |
|
|
261 L = 1 - J |
|
|
262 IF( NOUNIT ) |
|
|
263 $ X( JX ) = X( JX )/A( 1, J ) |
|
|
264 TEMP = X( JX ) |
|
|
265 DO 70, I = J + 1, MIN( N, J + K ) |
|
|
266 X( IX ) = X( IX ) - TEMP*A( L + I, J ) |
|
|
267 IX = IX + INCX |
|
|
268 70 CONTINUE |
|
|
269 END IF |
|
|
270 JX = JX + INCX |
|
|
271 80 CONTINUE |
|
|
272 END IF |
|
|
273 END IF |
|
|
274 ELSE |
|
|
275 * |
|
|
276 * Form x := inv( A' )*x or x := inv( conjg( A') )*x. |
|
|
277 * |
|
|
278 IF( LSAME( UPLO, 'U' ) )THEN |
|
|
279 KPLUS1 = K + 1 |
|
|
280 IF( INCX.EQ.1 )THEN |
|
|
281 DO 110, J = 1, N |
|
|
282 TEMP = X( J ) |
|
|
283 L = KPLUS1 - J |
|
|
284 IF( NOCONJ )THEN |
|
|
285 DO 90, I = MAX( 1, J - K ), J - 1 |
|
|
286 TEMP = TEMP - A( L + I, J )*X( I ) |
|
|
287 90 CONTINUE |
|
|
288 IF( NOUNIT ) |
|
|
289 $ TEMP = TEMP/A( KPLUS1, J ) |
|
|
290 ELSE |
|
|
291 DO 100, I = MAX( 1, J - K ), J - 1 |
|
|
292 TEMP = TEMP - DCONJG( A( L + I, J ) )*X( I ) |
|
|
293 100 CONTINUE |
|
|
294 IF( NOUNIT ) |
|
|
295 $ TEMP = TEMP/DCONJG( A( KPLUS1, J ) ) |
|
|
296 END IF |
|
|
297 X( J ) = TEMP |
|
|
298 110 CONTINUE |
|
|
299 ELSE |
|
|
300 JX = KX |
|
|
301 DO 140, J = 1, N |
|
|
302 TEMP = X( JX ) |
|
|
303 IX = KX |
|
|
304 L = KPLUS1 - J |
|
|
305 IF( NOCONJ )THEN |
|
|
306 DO 120, I = MAX( 1, J - K ), J - 1 |
|
|
307 TEMP = TEMP - A( L + I, J )*X( IX ) |
|
|
308 IX = IX + INCX |
|
|
309 120 CONTINUE |
|
|
310 IF( NOUNIT ) |
|
|
311 $ TEMP = TEMP/A( KPLUS1, J ) |
|
|
312 ELSE |
|
|
313 DO 130, I = MAX( 1, J - K ), J - 1 |
|
|
314 TEMP = TEMP - DCONJG( A( L + I, J ) )*X( IX ) |
|
|
315 IX = IX + INCX |
|
|
316 130 CONTINUE |
|
|
317 IF( NOUNIT ) |
|
|
318 $ TEMP = TEMP/DCONJG( A( KPLUS1, J ) ) |
|
|
319 END IF |
|
|
320 X( JX ) = TEMP |
|
|
321 JX = JX + INCX |
|
|
322 IF( J.GT.K ) |
|
|
323 $ KX = KX + INCX |
|
|
324 140 CONTINUE |
|
|
325 END IF |
|
|
326 ELSE |
|
|
327 IF( INCX.EQ.1 )THEN |
|
|
328 DO 170, J = N, 1, -1 |
|
|
329 TEMP = X( J ) |
|
|
330 L = 1 - J |
|
|
331 IF( NOCONJ )THEN |
|
|
332 DO 150, I = MIN( N, J + K ), J + 1, -1 |
|
|
333 TEMP = TEMP - A( L + I, J )*X( I ) |
|
|
334 150 CONTINUE |
|
|
335 IF( NOUNIT ) |
|
|
336 $ TEMP = TEMP/A( 1, J ) |
|
|
337 ELSE |
|
|
338 DO 160, I = MIN( N, J + K ), J + 1, -1 |
|
|
339 TEMP = TEMP - DCONJG( A( L + I, J ) )*X( I ) |
|
|
340 160 CONTINUE |
|
|
341 IF( NOUNIT ) |
|
|
342 $ TEMP = TEMP/DCONJG( A( 1, J ) ) |
|
|
343 END IF |
|
|
344 X( J ) = TEMP |
|
|
345 170 CONTINUE |
|
|
346 ELSE |
|
|
347 KX = KX + ( N - 1 )*INCX |
|
|
348 JX = KX |
|
|
349 DO 200, J = N, 1, -1 |
|
|
350 TEMP = X( JX ) |
|
|
351 IX = KX |
|
|
352 L = 1 - J |
|
|
353 IF( NOCONJ )THEN |
|
|
354 DO 180, I = MIN( N, J + K ), J + 1, -1 |
|
|
355 TEMP = TEMP - A( L + I, J )*X( IX ) |
|
|
356 IX = IX - INCX |
|
|
357 180 CONTINUE |
|
|
358 IF( NOUNIT ) |
|
|
359 $ TEMP = TEMP/A( 1, J ) |
|
|
360 ELSE |
|
|
361 DO 190, I = MIN( N, J + K ), J + 1, -1 |
|
|
362 TEMP = TEMP - DCONJG( A( L + I, J ) )*X( IX ) |
|
|
363 IX = IX - INCX |
|
|
364 190 CONTINUE |
|
|
365 IF( NOUNIT ) |
|
|
366 $ TEMP = TEMP/DCONJG( A( 1, J ) ) |
|
|
367 END IF |
|
|
368 X( JX ) = TEMP |
|
|
369 JX = JX - INCX |
|
|
370 IF( ( N - J ).GE.K ) |
|
|
371 $ KX = KX - INCX |
|
|
372 200 CONTINUE |
|
|
373 END IF |
|
|
374 END IF |
|
|
375 END IF |
|
|
376 * |
|
|
377 RETURN |
|
|
378 * |
|
|
379 * End of ZTBSV . |
|
|
380 * |
|
|
381 END |
