source: pacpussensors/trunk/Vislab/lib3dv/eigen/blas/zhbmv.f@ 138

Last change on this file since 138 was 136, checked in by ldecherf, 8 years ago

Doc

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