1 SUBROUTINE CHBMV(UPLO,N,K,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
2 * .. Scalar Arguments ..
3 COMPLEX ALPHA,BETA
4 INTEGER INCX,INCY,K,LDA,N
5 CHARACTER UPLO
6 * ..
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 * Arguments
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 * Further Details
124 * ===============
125 *
126 * Level 2 Blas routine.
127 * The vector and matrix arguments are not referenced when N = 0, or M = 0
128 *
129 * -- Written on 22-October-1986.
130 * Jack Dongarra, Argonne National Lab.
131 * Jeremy Du Croz, Nag Central Office.
132 * Sven Hammarling, Nag Central Office.
133 * Richard Hanson, Sandia National Labs.
134 *
135 * =====================================================================
136 *
137 * .. Parameters ..
138 COMPLEX ONE
139 PARAMETER (ONE= (1.0E+0,0.0E+0))
140 COMPLEX ZERO
141 PARAMETER (ZERO= (0.0E+0,0.0E+0))
142 * ..
143 * .. Local Scalars ..
144 COMPLEX TEMP1,TEMP2
145 INTEGER I,INFO,IX,IY,J,JX,JY,KPLUS1,KX,KY,L
146 * ..
147 * .. External Functions ..
148 LOGICAL LSAME
149 EXTERNAL LSAME
150 * ..
151 * .. External Subroutines ..
152 EXTERNAL XERBLA
153 * ..
154 * .. Intrinsic Functions ..
155 INTRINSIC CONJG,MAX,MIN,REAL
156 * ..
157 *
158 * Test the input parameters.
159 *
160 INFO = 0
161 IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
162 INFO = 1
163 ELSE IF (N.LT.0) THEN
164 INFO = 2
165 ELSE IF (K.LT.0) THEN
166 INFO = 3
167 ELSE IF (LDA.LT. (K+1)) THEN
168 INFO = 6
169 ELSE IF (INCX.EQ.0) THEN
170 INFO = 8
171 ELSE IF (INCY.EQ.0) THEN
172 INFO = 11
173 END IF
174 IF (INFO.NE.0) THEN
175 CALL XERBLA('CHBMV ',INFO)
176 RETURN
177 END IF
178 *
179 * Quick return if possible.
180 *
181 IF ((N.EQ.0) .OR. ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
182 *
183 * Set up the start points in X and Y.
184 *
185 IF (INCX.GT.0) THEN
186 KX = 1
187 ELSE
188 KX = 1 - (N-1)*INCX
189 END IF
190 IF (INCY.GT.0) THEN
191 KY = 1
192 ELSE
193 KY = 1 - (N-1)*INCY
194 END IF
195 *
196 * Start the operations. In this version the elements of the array A
197 * are accessed sequentially with one pass through A.
198 *
199 * First form y := beta*y.
200 *
201 IF (BETA.NE.ONE) THEN
202 IF (INCY.EQ.1) THEN
203 IF (BETA.EQ.ZERO) THEN
204 DO 10 I = 1,N
205 Y(I) = ZERO
206 10 CONTINUE
207 ELSE
208 DO 20 I = 1,N
209 Y(I) = BETA*Y(I)
210 20 CONTINUE
211 END IF
212 ELSE
213 IY = KY
214 IF (BETA.EQ.ZERO) THEN
215 DO 30 I = 1,N
216 Y(IY) = ZERO
217 IY = IY + INCY
218 30 CONTINUE
219 ELSE
220 DO 40 I = 1,N
221 Y(IY) = BETA*Y(IY)
222 IY = IY + INCY
223 40 CONTINUE
224 END IF
225 END IF
226 END IF
227 IF (ALPHA.EQ.ZERO) RETURN
228 IF (LSAME(UPLO,'U')) THEN
229 *
230 * Form y when upper triangle of A is stored.
231 *
232 KPLUS1 = K + 1
233 IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
234 DO 60 J = 1,N
235 TEMP1 = ALPHA*X(J)
236 TEMP2 = ZERO
237 L = KPLUS1 - J
238 DO 50 I = MAX(1,J-K),J - 1
239 Y(I) = Y(I) + TEMP1*A(L+I,J)
240 TEMP2 = TEMP2 + CONJG(A(L+I,J))*X(I)
241 50 CONTINUE
242 Y(J) = Y(J) + TEMP1*REAL(A(KPLUS1,J)) + ALPHA*TEMP2
243 60 CONTINUE
244 ELSE
245 JX = KX
246 JY = KY
247 DO 80 J = 1,N
248 TEMP1 = ALPHA*X(JX)
249 TEMP2 = ZERO
250 IX = KX
251 IY = KY
252 L = KPLUS1 - J
253 DO 70 I = MAX(1,J-K),J - 1
254 Y(IY) = Y(IY) + TEMP1*A(L+I,J)
255 TEMP2 = TEMP2 + CONJG(A(L+I,J))*X(IX)
256 IX = IX + INCX
257 IY = IY + INCY
258 70 CONTINUE
259 Y(JY) = Y(JY) + TEMP1*REAL(A(KPLUS1,J)) + ALPHA*TEMP2
260 JX = JX + INCX
261 JY = JY + INCY
262 IF (J.GT.K) THEN
263 KX = KX + INCX
264 KY = KY + INCY
265 END IF
266 80 CONTINUE
267 END IF
268 ELSE
269 *
270 * Form y when lower triangle of A is stored.
271 *
272 IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
273 DO 100 J = 1,N
274 TEMP1 = ALPHA*X(J)
275 TEMP2 = ZERO
276 Y(J) = Y(J) + TEMP1*REAL(A(1,J))
277 L = 1 - J
278 DO 90 I = J + 1,MIN(N,J+K)
279 Y(I) = Y(I) + TEMP1*A(L+I,J)
280 TEMP2 = TEMP2 + CONJG(A(L+I,J))*X(I)
281 90 CONTINUE
282 Y(J) = Y(J) + ALPHA*TEMP2
283 100 CONTINUE
284 ELSE
285 JX = KX
286 JY = KY
287 DO 120 J = 1,N
288 TEMP1 = ALPHA*X(JX)
289 TEMP2 = ZERO
290 Y(JY) = Y(JY) + TEMP1*REAL(A(1,J))
291 L = 1 - J
292 IX = JX
293 IY = JY
294 DO 110 I = J + 1,MIN(N,J+K)
295 IX = IX + INCX
296 IY = IY + INCY
297 Y(IY) = Y(IY) + TEMP1*A(L+I,J)
298 TEMP2 = TEMP2 + CONJG(A(L+I,J))*X(IX)
299 110 CONTINUE
300 Y(JY) = Y(JY) + ALPHA*TEMP2
301 JX = JX + INCX
302 JY = JY + INCY
303 120 CONTINUE
304 END IF
305 END IF
306 *
307 RETURN
308 *
309 * End of CHBMV .
310 *
311 END
2 * .. Scalar Arguments ..
3 COMPLEX ALPHA,BETA
4 INTEGER INCX,INCY,K,LDA,N
5 CHARACTER UPLO
6 * ..
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 * Arguments
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 * Further Details
124 * ===============
125 *
126 * Level 2 Blas routine.
127 * The vector and matrix arguments are not referenced when N = 0, or M = 0
128 *
129 * -- Written on 22-October-1986.
130 * Jack Dongarra, Argonne National Lab.
131 * Jeremy Du Croz, Nag Central Office.
132 * Sven Hammarling, Nag Central Office.
133 * Richard Hanson, Sandia National Labs.
134 *
135 * =====================================================================
136 *
137 * .. Parameters ..
138 COMPLEX ONE
139 PARAMETER (ONE= (1.0E+0,0.0E+0))
140 COMPLEX ZERO
141 PARAMETER (ZERO= (0.0E+0,0.0E+0))
142 * ..
143 * .. Local Scalars ..
144 COMPLEX TEMP1,TEMP2
145 INTEGER I,INFO,IX,IY,J,JX,JY,KPLUS1,KX,KY,L
146 * ..
147 * .. External Functions ..
148 LOGICAL LSAME
149 EXTERNAL LSAME
150 * ..
151 * .. External Subroutines ..
152 EXTERNAL XERBLA
153 * ..
154 * .. Intrinsic Functions ..
155 INTRINSIC CONJG,MAX,MIN,REAL
156 * ..
157 *
158 * Test the input parameters.
159 *
160 INFO = 0
161 IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
162 INFO = 1
163 ELSE IF (N.LT.0) THEN
164 INFO = 2
165 ELSE IF (K.LT.0) THEN
166 INFO = 3
167 ELSE IF (LDA.LT. (K+1)) THEN
168 INFO = 6
169 ELSE IF (INCX.EQ.0) THEN
170 INFO = 8
171 ELSE IF (INCY.EQ.0) THEN
172 INFO = 11
173 END IF
174 IF (INFO.NE.0) THEN
175 CALL XERBLA('CHBMV ',INFO)
176 RETURN
177 END IF
178 *
179 * Quick return if possible.
180 *
181 IF ((N.EQ.0) .OR. ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
182 *
183 * Set up the start points in X and Y.
184 *
185 IF (INCX.GT.0) THEN
186 KX = 1
187 ELSE
188 KX = 1 - (N-1)*INCX
189 END IF
190 IF (INCY.GT.0) THEN
191 KY = 1
192 ELSE
193 KY = 1 - (N-1)*INCY
194 END IF
195 *
196 * Start the operations. In this version the elements of the array A
197 * are accessed sequentially with one pass through A.
198 *
199 * First form y := beta*y.
200 *
201 IF (BETA.NE.ONE) THEN
202 IF (INCY.EQ.1) THEN
203 IF (BETA.EQ.ZERO) THEN
204 DO 10 I = 1,N
205 Y(I) = ZERO
206 10 CONTINUE
207 ELSE
208 DO 20 I = 1,N
209 Y(I) = BETA*Y(I)
210 20 CONTINUE
211 END IF
212 ELSE
213 IY = KY
214 IF (BETA.EQ.ZERO) THEN
215 DO 30 I = 1,N
216 Y(IY) = ZERO
217 IY = IY + INCY
218 30 CONTINUE
219 ELSE
220 DO 40 I = 1,N
221 Y(IY) = BETA*Y(IY)
222 IY = IY + INCY
223 40 CONTINUE
224 END IF
225 END IF
226 END IF
227 IF (ALPHA.EQ.ZERO) RETURN
228 IF (LSAME(UPLO,'U')) THEN
229 *
230 * Form y when upper triangle of A is stored.
231 *
232 KPLUS1 = K + 1
233 IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
234 DO 60 J = 1,N
235 TEMP1 = ALPHA*X(J)
236 TEMP2 = ZERO
237 L = KPLUS1 - J
238 DO 50 I = MAX(1,J-K),J - 1
239 Y(I) = Y(I) + TEMP1*A(L+I,J)
240 TEMP2 = TEMP2 + CONJG(A(L+I,J))*X(I)
241 50 CONTINUE
242 Y(J) = Y(J) + TEMP1*REAL(A(KPLUS1,J)) + ALPHA*TEMP2
243 60 CONTINUE
244 ELSE
245 JX = KX
246 JY = KY
247 DO 80 J = 1,N
248 TEMP1 = ALPHA*X(JX)
249 TEMP2 = ZERO
250 IX = KX
251 IY = KY
252 L = KPLUS1 - J
253 DO 70 I = MAX(1,J-K),J - 1
254 Y(IY) = Y(IY) + TEMP1*A(L+I,J)
255 TEMP2 = TEMP2 + CONJG(A(L+I,J))*X(IX)
256 IX = IX + INCX
257 IY = IY + INCY
258 70 CONTINUE
259 Y(JY) = Y(JY) + TEMP1*REAL(A(KPLUS1,J)) + ALPHA*TEMP2
260 JX = JX + INCX
261 JY = JY + INCY
262 IF (J.GT.K) THEN
263 KX = KX + INCX
264 KY = KY + INCY
265 END IF
266 80 CONTINUE
267 END IF
268 ELSE
269 *
270 * Form y when lower triangle of A is stored.
271 *
272 IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
273 DO 100 J = 1,N
274 TEMP1 = ALPHA*X(J)
275 TEMP2 = ZERO
276 Y(J) = Y(J) + TEMP1*REAL(A(1,J))
277 L = 1 - J
278 DO 90 I = J + 1,MIN(N,J+K)
279 Y(I) = Y(I) + TEMP1*A(L+I,J)
280 TEMP2 = TEMP2 + CONJG(A(L+I,J))*X(I)
281 90 CONTINUE
282 Y(J) = Y(J) + ALPHA*TEMP2
283 100 CONTINUE
284 ELSE
285 JX = KX
286 JY = KY
287 DO 120 J = 1,N
288 TEMP1 = ALPHA*X(JX)
289 TEMP2 = ZERO
290 Y(JY) = Y(JY) + TEMP1*REAL(A(1,J))
291 L = 1 - J
292 IX = JX
293 IY = JY
294 DO 110 I = J + 1,MIN(N,J+K)
295 IX = IX + INCX
296 IY = IY + INCY
297 Y(IY) = Y(IY) + TEMP1*A(L+I,J)
298 TEMP2 = TEMP2 + CONJG(A(L+I,J))*X(IX)
299 110 CONTINUE
300 Y(JY) = Y(JY) + ALPHA*TEMP2
301 JX = JX + INCX
302 JY = JY + INCY
303 120 CONTINUE
304 END IF
305 END IF
306 *
307 RETURN
308 *
309 * End of CHBMV .
310 *
311 END