source: pacpussensors/trunk/Vislab/lib3dv/eigen/blas/dspmv.f@ 140

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

Doc
File size: 7.8 KB
Line 
1 SUBROUTINE DSPMV(UPLO,N,ALPHA,AP,X,INCX,BETA,Y,INCY)
2* .. Scalar Arguments ..
3 DOUBLE PRECISION ALPHA,BETA
4 INTEGER INCX,INCY,N
5 CHARACTER UPLO
6* ..
7* .. Array Arguments ..
8 DOUBLE PRECISION AP(*),X(*),Y(*)
9* ..
10*
11* Purpose
12* =======
13*
14* DSPMV 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 matrix, supplied in packed form.
20*
21* Arguments
22* ==========
23*
24* UPLO - CHARACTER*1.
25* On entry, UPLO specifies whether the upper or lower
26* triangular part of the matrix A is supplied in the packed
27* array AP as follows:
28*
29* UPLO = 'U' or 'u' The upper triangular part of A is
30* supplied in AP.
31*
32* UPLO = 'L' or 'l' The lower triangular part of A is
33* supplied in AP.
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* ALPHA - DOUBLE PRECISION.
43* On entry, ALPHA specifies the scalar alpha.
44* Unchanged on exit.
45*
46* AP - DOUBLE PRECISION array of DIMENSION at least
47* ( ( n*( n + 1 ) )/2 ).
48* Before entry with UPLO = 'U' or 'u', the array AP must
49* contain the upper triangular part of the symmetric matrix
50* packed sequentially, column by column, so that AP( 1 )
51* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
52* and a( 2, 2 ) respectively, and so on.
53* Before entry with UPLO = 'L' or 'l', the array AP must
54* contain the lower triangular part of the symmetric matrix
55* packed sequentially, column by column, so that AP( 1 )
56* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
57* and a( 3, 1 ) respectively, and so on.
58* Unchanged on exit.
59*
60* X - DOUBLE PRECISION array of dimension at least
61* ( 1 + ( n - 1 )*abs( INCX ) ).
62* Before entry, the incremented array X must contain the n
63* element vector x.
64* Unchanged on exit.
65*
66* INCX - INTEGER.
67* On entry, INCX specifies the increment for the elements of
68* X. INCX must not be zero.
69* Unchanged on exit.
70*
71* BETA - DOUBLE PRECISION.
72* On entry, BETA specifies the scalar beta. When BETA is
73* supplied as zero then Y need not be set on input.
74* Unchanged on exit.
75*
76* Y - DOUBLE PRECISION array of dimension at least
77* ( 1 + ( n - 1 )*abs( INCY ) ).
78* Before entry, the incremented array Y must contain the n
79* element vector y. On exit, Y is overwritten by the updated
80* vector y.
81*
82* INCY - INTEGER.
83* On entry, INCY specifies the increment for the elements of
84* Y. INCY must not be zero.
85* Unchanged on exit.
86*
87* Further Details
88* ===============
89*
90* Level 2 Blas routine.
91*
92* -- Written on 22-October-1986.
93* Jack Dongarra, Argonne National Lab.
94* Jeremy Du Croz, Nag Central Office.
95* Sven Hammarling, Nag Central Office.
96* Richard Hanson, Sandia National Labs.
97*
98* =====================================================================
99*
100* .. Parameters ..
101 DOUBLE PRECISION ONE,ZERO
102 PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
103* ..
104* .. Local Scalars ..
105 DOUBLE PRECISION TEMP1,TEMP2
106 INTEGER I,INFO,IX,IY,J,JX,JY,K,KK,KX,KY
107* ..
108* .. External Functions ..
109 LOGICAL LSAME
110 EXTERNAL LSAME
111* ..
112* .. External Subroutines ..
113 EXTERNAL XERBLA
114* ..
115*
116* Test the input parameters.
117*
118 INFO = 0
119 IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
120 INFO = 1
121 ELSE IF (N.LT.0) THEN
122 INFO = 2
123 ELSE IF (INCX.EQ.0) THEN
124 INFO = 6
125 ELSE IF (INCY.EQ.0) THEN
126 INFO = 9
127 END IF
128 IF (INFO.NE.0) THEN
129 CALL XERBLA('DSPMV ',INFO)
130 RETURN
131 END IF
132*
133* Quick return if possible.
134*
135 IF ((N.EQ.0) .OR. ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
136*
137* Set up the start points in X and Y.
138*
139 IF (INCX.GT.0) THEN
140 KX = 1
141 ELSE
142 KX = 1 - (N-1)*INCX
143 END IF
144 IF (INCY.GT.0) THEN
145 KY = 1
146 ELSE
147 KY = 1 - (N-1)*INCY
148 END IF
149*
150* Start the operations. In this version the elements of the array AP
151* are accessed sequentially with one pass through AP.
152*
153* First form y := beta*y.
154*
155 IF (BETA.NE.ONE) THEN
156 IF (INCY.EQ.1) THEN
157 IF (BETA.EQ.ZERO) THEN
158 DO 10 I = 1,N
159 Y(I) = ZERO
160 10 CONTINUE
161 ELSE
162 DO 20 I = 1,N
163 Y(I) = BETA*Y(I)
164 20 CONTINUE
165 END IF
166 ELSE
167 IY = KY
168 IF (BETA.EQ.ZERO) THEN
169 DO 30 I = 1,N
170 Y(IY) = ZERO
171 IY = IY + INCY
172 30 CONTINUE
173 ELSE
174 DO 40 I = 1,N
175 Y(IY) = BETA*Y(IY)
176 IY = IY + INCY
177 40 CONTINUE
178 END IF
179 END IF
180 END IF
181 IF (ALPHA.EQ.ZERO) RETURN
182 KK = 1
183 IF (LSAME(UPLO,'U')) THEN
184*
185* Form y when AP contains the upper triangle.
186*
187 IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
188 DO 60 J = 1,N
189 TEMP1 = ALPHA*X(J)
190 TEMP2 = ZERO
191 K = KK
192 DO 50 I = 1,J - 1
193 Y(I) = Y(I) + TEMP1*AP(K)
194 TEMP2 = TEMP2 + AP(K)*X(I)
195 K = K + 1
196 50 CONTINUE
197 Y(J) = Y(J) + TEMP1*AP(KK+J-1) + ALPHA*TEMP2
198 KK = KK + J
199 60 CONTINUE
200 ELSE
201 JX = KX
202 JY = KY
203 DO 80 J = 1,N
204 TEMP1 = ALPHA*X(JX)
205 TEMP2 = ZERO
206 IX = KX
207 IY = KY
208 DO 70 K = KK,KK + J - 2
209 Y(IY) = Y(IY) + TEMP1*AP(K)
210 TEMP2 = TEMP2 + AP(K)*X(IX)
211 IX = IX + INCX
212 IY = IY + INCY
213 70 CONTINUE
214 Y(JY) = Y(JY) + TEMP1*AP(KK+J-1) + ALPHA*TEMP2
215 JX = JX + INCX
216 JY = JY + INCY
217 KK = KK + J
218 80 CONTINUE
219 END IF
220 ELSE
221*
222* Form y when AP contains the lower triangle.
223*
224 IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
225 DO 100 J = 1,N
226 TEMP1 = ALPHA*X(J)
227 TEMP2 = ZERO
228 Y(J) = Y(J) + TEMP1*AP(KK)
229 K = KK + 1
230 DO 90 I = J + 1,N
231 Y(I) = Y(I) + TEMP1*AP(K)
232 TEMP2 = TEMP2 + AP(K)*X(I)
233 K = K + 1
234 90 CONTINUE
235 Y(J) = Y(J) + ALPHA*TEMP2
236 KK = KK + (N-J+1)
237 100 CONTINUE
238 ELSE
239 JX = KX
240 JY = KY
241 DO 120 J = 1,N
242 TEMP1 = ALPHA*X(JX)
243 TEMP2 = ZERO
244 Y(JY) = Y(JY) + TEMP1*AP(KK)
245 IX = JX
246 IY = JY
247 DO 110 K = KK + 1,KK + N - J
248 IX = IX + INCX
249 IY = IY + INCY
250 Y(IY) = Y(IY) + TEMP1*AP(K)
251 TEMP2 = TEMP2 + AP(K)*X(IX)
252 110 CONTINUE
253 Y(JY) = Y(JY) + ALPHA*TEMP2
254 JX = JX + INCX
255 JY = JY + INCY
256 KK = KK + (N-J+1)
257 120 CONTINUE
258 END IF
259 END IF
260*
261 RETURN
262*
263* End of DSPMV .
264*
265 END
Note: See TracBrowser for help on using the repository browser.