Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
ElmerCSC
GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/mathlibs/src/blas/chbmv.f
5198 views
1
SUBROUTINE CHBMV ( UPLO, N, K, ALPHA, A, LDA, X, INCX,
2
$ BETA, Y, INCY )
3
* .. Scalar Arguments ..
4
COMPLEX ALPHA, BETA
5
INTEGER INCX, INCY, K, LDA, N
6
CHARACTER*1 UPLO
7
* .. Array Arguments ..
8
COMPLEX A( LDA, * ), X( * ), Y( * )
9
* ..
10
*
11
* Purpose
12
* =======
13
*
14
* CHBMV performs the matrix-vector operation
15
*
16
* y := alpha*A*x + beta*y,
17
*
18
* where alpha and beta are scalars, x and y are n element vectors and
19
* A is an n by n hermitian band matrix, with k super-diagonals.
20
*
21
* Parameters
22
* ==========
23
*
24
* UPLO - CHARACTER*1.
25
* On entry, UPLO specifies whether the upper or lower
26
* triangular part of the band matrix A is being supplied as
27
* follows:
28
*
29
* UPLO = 'U' or 'u' The upper triangular part of A is
30
* being supplied.
31
*
32
* UPLO = 'L' or 'l' The lower triangular part of A is
33
* being supplied.
34
*
35
* Unchanged on exit.
36
*
37
* N - INTEGER.
38
* On entry, N specifies the order of the matrix A.
39
* N must be at least zero.
40
* Unchanged on exit.
41
*
42
* K - INTEGER.
43
* On entry, K specifies the number of super-diagonals of the
44
* matrix A. K must satisfy 0 .le. K.
45
* Unchanged on exit.
46
*
47
* ALPHA - COMPLEX .
48
* On entry, ALPHA specifies the scalar alpha.
49
* Unchanged on exit.
50
*
51
* A - COMPLEX array of DIMENSION ( LDA, n ).
52
* Before entry with UPLO = 'U' or 'u', the leading ( k + 1 )
53
* by n part of the array A must contain the upper triangular
54
* band part of the hermitian matrix, supplied column by
55
* column, with the leading diagonal of the matrix in row
56
* ( k + 1 ) of the array, the first super-diagonal starting at
57
* position 2 in row k, and so on. The top left k by k triangle
58
* of the array A is not referenced.
59
* The following program segment will transfer the upper
60
* triangular part of a hermitian band matrix from conventional
61
* full matrix storage to band storage:
62
*
63
* DO 20, J = 1, N
64
* M = K + 1 - J
65
* DO 10, I = MAX( 1, J - K ), J
66
* A( M + I, J ) = matrix( I, J )
67
* 10 CONTINUE
68
* 20 CONTINUE
69
*
70
* Before entry with UPLO = 'L' or 'l', the leading ( k + 1 )
71
* by n part of the array A must contain the lower triangular
72
* band part of the hermitian matrix, supplied column by
73
* column, with the leading diagonal of the matrix in row 1 of
74
* the array, the first sub-diagonal starting at position 1 in
75
* row 2, and so on. The bottom right k by k triangle of the
76
* array A is not referenced.
77
* The following program segment will transfer the lower
78
* triangular part of a hermitian band matrix from conventional
79
* full matrix storage to band storage:
80
*
81
* DO 20, J = 1, N
82
* M = 1 - J
83
* DO 10, I = J, MIN( N, J + K )
84
* A( M + I, J ) = matrix( I, J )
85
* 10 CONTINUE
86
* 20 CONTINUE
87
*
88
* Note that the imaginary parts of the diagonal elements need
89
* not be set and are assumed to be zero.
90
* Unchanged on exit.
91
*
92
* LDA - INTEGER.
93
* On entry, LDA specifies the first dimension of A as declared
94
* in the calling (sub) program. LDA must be at least
95
* ( k + 1 ).
96
* Unchanged on exit.
97
*
98
* X - COMPLEX array of DIMENSION at least
99
* ( 1 + ( n - 1 )*abs( INCX ) ).
100
* Before entry, the incremented array X must contain the
101
* vector x.
102
* Unchanged on exit.
103
*
104
* INCX - INTEGER.
105
* On entry, INCX specifies the increment for the elements of
106
* X. INCX must not be zero.
107
* Unchanged on exit.
108
*
109
* BETA - COMPLEX .
110
* On entry, BETA specifies the scalar beta.
111
* Unchanged on exit.
112
*
113
* Y - COMPLEX array of DIMENSION at least
114
* ( 1 + ( n - 1 )*abs( INCY ) ).
115
* Before entry, the incremented array Y must contain the
116
* vector y. On exit, Y is overwritten by the updated vector y.
117
*
118
* INCY - INTEGER.
119
* On entry, INCY specifies the increment for the elements of
120
* Y. INCY must not be zero.
121
* Unchanged on exit.
122
*
123
*
124
* Level 2 Blas routine.
125
*
126
* -- Written on 22-October-1986.
127
* Jack Dongarra, Argonne National Lab.
128
* Jeremy Du Croz, Nag Central Office.
129
* Sven Hammarling, Nag Central Office.
130
* Richard Hanson, Sandia National Labs.
131
*
132
*
133
* .. Parameters ..
134
COMPLEX ONE
135
PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) )
136
COMPLEX ZERO
137
PARAMETER ( ZERO = ( 0.0E+0, 0.0E+0 ) )
138
* .. Local Scalars ..
139
COMPLEX TEMP1, TEMP2
140
INTEGER I, INFO, IX, IY, J, JX, JY, KPLUS1, KX, KY, L
141
* .. External Functions ..
142
LOGICAL LSAME
143
EXTERNAL LSAME
144
* .. External Subroutines ..
145
EXTERNAL XERBLA
146
* .. Intrinsic Functions ..
147
INTRINSIC CONJG, MAX, MIN, REAL
148
* ..
149
* .. Executable Statements ..
150
*
151
* Test the input parameters.
152
*
153
INFO = 0
154
IF ( .NOT.LSAME( UPLO, 'U' ).AND.
155
$ .NOT.LSAME( UPLO, 'L' ) )THEN
156
INFO = 1
157
ELSE IF( N.LT.0 )THEN
158
INFO = 2
159
ELSE IF( K.LT.0 )THEN
160
INFO = 3
161
ELSE IF( LDA.LT.( K + 1 ) )THEN
162
INFO = 6
163
ELSE IF( INCX.EQ.0 )THEN
164
INFO = 8
165
ELSE IF( INCY.EQ.0 )THEN
166
INFO = 11
167
END IF
168
IF( INFO.NE.0 )THEN
169
CALL XERBLA( 'CHBMV ', INFO )
170
RETURN
171
END IF
172
*
173
* Quick return if possible.
174
*
175
IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
176
$ RETURN
177
*
178
* Set up the start points in X and Y.
179
*
180
IF( INCX.GT.0 )THEN
181
KX = 1
182
ELSE
183
KX = 1 - ( N - 1 )*INCX
184
END IF
185
IF( INCY.GT.0 )THEN
186
KY = 1
187
ELSE
188
KY = 1 - ( N - 1 )*INCY
189
END IF
190
*
191
* Start the operations. In this version the elements of the array A
192
* are accessed sequentially with one pass through A.
193
*
194
* First form y := beta*y.
195
*
196
IF( BETA.NE.ONE )THEN
197
IF( INCY.EQ.1 )THEN
198
IF( BETA.EQ.ZERO )THEN
199
DO 10, I = 1, N
200
Y( I ) = ZERO
201
10 CONTINUE
202
ELSE
203
DO 20, I = 1, N
204
Y( I ) = BETA*Y( I )
205
20 CONTINUE
206
END IF
207
ELSE
208
IY = KY
209
IF( BETA.EQ.ZERO )THEN
210
DO 30, I = 1, N
211
Y( IY ) = ZERO
212
IY = IY + INCY
213
30 CONTINUE
214
ELSE
215
DO 40, I = 1, N
216
Y( IY ) = BETA*Y( IY )
217
IY = IY + INCY
218
40 CONTINUE
219
END IF
220
END IF
221
END IF
222
IF( ALPHA.EQ.ZERO )
223
$ RETURN
224
IF( LSAME( UPLO, 'U' ) )THEN
225
*
226
* Form y when upper triangle of A is stored.
227
*
228
KPLUS1 = K + 1
229
IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
230
DO 60, J = 1, N
231
TEMP1 = ALPHA*X( J )
232
TEMP2 = ZERO
233
L = KPLUS1 - J
234
DO 50, I = MAX( 1, J - K ), J - 1
235
Y( I ) = Y( I ) + TEMP1*A( L + I, J )
236
TEMP2 = TEMP2 + CONJG( A( L + I, J ) )*X( I )
237
50 CONTINUE
238
Y( J ) = Y( J ) + TEMP1*REAL( A( KPLUS1, J ) )
239
$ + ALPHA*TEMP2
240
60 CONTINUE
241
ELSE
242
JX = KX
243
JY = KY
244
DO 80, J = 1, N
245
TEMP1 = ALPHA*X( JX )
246
TEMP2 = ZERO
247
IX = KX
248
IY = KY
249
L = KPLUS1 - J
250
DO 70, I = MAX( 1, J - K ), J - 1
251
Y( IY ) = Y( IY ) + TEMP1*A( L + I, J )
252
TEMP2 = TEMP2 + CONJG( A( L + I, J ) )*X( IX )
253
IX = IX + INCX
254
IY = IY + INCY
255
70 CONTINUE
256
Y( JY ) = Y( JY ) + TEMP1*REAL( A( KPLUS1, J ) )
257
$ + ALPHA*TEMP2
258
JX = JX + INCX
259
JY = JY + INCY
260
IF( J.GT.K )THEN
261
KX = KX + INCX
262
KY = KY + INCY
263
END IF
264
80 CONTINUE
265
END IF
266
ELSE
267
*
268
* Form y when lower triangle of A is stored.
269
*
270
IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
271
DO 100, J = 1, N
272
TEMP1 = ALPHA*X( J )
273
TEMP2 = ZERO
274
Y( J ) = Y( J ) + TEMP1*REAL( A( 1, J ) )
275
L = 1 - J
276
DO 90, I = J + 1, MIN( N, J + K )
277
Y( I ) = Y( I ) + TEMP1*A( L + I, J )
278
TEMP2 = TEMP2 + CONJG( A( L + I, J ) )*X( I )
279
90 CONTINUE
280
Y( J ) = Y( J ) + ALPHA*TEMP2
281
100 CONTINUE
282
ELSE
283
JX = KX
284
JY = KY
285
DO 120, J = 1, N
286
TEMP1 = ALPHA*X( JX )
287
TEMP2 = ZERO
288
Y( JY ) = Y( JY ) + TEMP1*REAL( A( 1, J ) )
289
L = 1 - J
290
IX = JX
291
IY = JY
292
DO 110, I = J + 1, MIN( N, J + K )
293
IX = IX + INCX
294
IY = IY + INCY
295
Y( IY ) = Y( IY ) + TEMP1*A( L + I, J )
296
TEMP2 = TEMP2 + CONJG( A( L + I, J ) )*X( IX )
297
110 CONTINUE
298
Y( JY ) = Y( JY ) + ALPHA*TEMP2
299
JX = JX + INCX
300
JY = JY + INCY
301
120 CONTINUE
302
END IF
303
END IF
304
*
305
RETURN
306
*
307
* End of CHBMV .
308
*
309
END
310
311