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