source: pacpussensors/trunk/Vislab/lib3dv-1.2.0/lib3dv/eigen/blas/dsbmv.f

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

Doc

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