root / src / blas / dtpmv.f @ 5
Historique | Voir | Annoter | Télécharger (8,96 ko)
1 | 1 | pfleura2 | SUBROUTINE DTPMV(UPLO,TRANS,DIAG,N,AP,X,INCX) |
---|---|---|---|
2 | 1 | pfleura2 | * .. Scalar Arguments .. |
3 | 1 | pfleura2 | INTEGER INCX,N |
4 | 1 | pfleura2 | CHARACTER DIAG,TRANS,UPLO |
5 | 1 | pfleura2 | * .. |
6 | 1 | pfleura2 | * .. Array Arguments .. |
7 | 1 | pfleura2 | DOUBLE PRECISION AP(*),X(*) |
8 | 1 | pfleura2 | * .. |
9 | 1 | pfleura2 | * |
10 | 1 | pfleura2 | * Purpose |
11 | 1 | pfleura2 | * ======= |
12 | 1 | pfleura2 | * |
13 | 1 | pfleura2 | * DTPMV performs one of the matrix-vector operations |
14 | 1 | pfleura2 | * |
15 | 1 | pfleura2 | * x := A*x, or x := A'*x, |
16 | 1 | pfleura2 | * |
17 | 1 | pfleura2 | * where x is an n element vector and A is an n by n unit, or non-unit, |
18 | 1 | pfleura2 | * upper or lower triangular matrix, supplied in packed form. |
19 | 1 | pfleura2 | * |
20 | 1 | pfleura2 | * Arguments |
21 | 1 | pfleura2 | * ========== |
22 | 1 | pfleura2 | * |
23 | 1 | pfleura2 | * UPLO - CHARACTER*1. |
24 | 1 | pfleura2 | * On entry, UPLO specifies whether the matrix is an upper or |
25 | 1 | pfleura2 | * lower triangular matrix as follows: |
26 | 1 | pfleura2 | * |
27 | 1 | pfleura2 | * UPLO = 'U' or 'u' A is an upper triangular matrix. |
28 | 1 | pfleura2 | * |
29 | 1 | pfleura2 | * UPLO = 'L' or 'l' A is a lower triangular matrix. |
30 | 1 | pfleura2 | * |
31 | 1 | pfleura2 | * Unchanged on exit. |
32 | 1 | pfleura2 | * |
33 | 1 | pfleura2 | * TRANS - CHARACTER*1. |
34 | 1 | pfleura2 | * On entry, TRANS specifies the operation to be performed as |
35 | 1 | pfleura2 | * follows: |
36 | 1 | pfleura2 | * |
37 | 1 | pfleura2 | * TRANS = 'N' or 'n' x := A*x. |
38 | 1 | pfleura2 | * |
39 | 1 | pfleura2 | * TRANS = 'T' or 't' x := A'*x. |
40 | 1 | pfleura2 | * |
41 | 1 | pfleura2 | * TRANS = 'C' or 'c' x := A'*x. |
42 | 1 | pfleura2 | * |
43 | 1 | pfleura2 | * Unchanged on exit. |
44 | 1 | pfleura2 | * |
45 | 1 | pfleura2 | * DIAG - CHARACTER*1. |
46 | 1 | pfleura2 | * On entry, DIAG specifies whether or not A is unit |
47 | 1 | pfleura2 | * triangular as follows: |
48 | 1 | pfleura2 | * |
49 | 1 | pfleura2 | * DIAG = 'U' or 'u' A is assumed to be unit triangular. |
50 | 1 | pfleura2 | * |
51 | 1 | pfleura2 | * DIAG = 'N' or 'n' A is not assumed to be unit |
52 | 1 | pfleura2 | * triangular. |
53 | 1 | pfleura2 | * |
54 | 1 | pfleura2 | * Unchanged on exit. |
55 | 1 | pfleura2 | * |
56 | 1 | pfleura2 | * N - INTEGER. |
57 | 1 | pfleura2 | * On entry, N specifies the order of the matrix A. |
58 | 1 | pfleura2 | * N must be at least zero. |
59 | 1 | pfleura2 | * Unchanged on exit. |
60 | 1 | pfleura2 | * |
61 | 1 | pfleura2 | * AP - DOUBLE PRECISION array of DIMENSION at least |
62 | 1 | pfleura2 | * ( ( n*( n + 1 ) )/2 ). |
63 | 1 | pfleura2 | * Before entry with UPLO = 'U' or 'u', the array AP must |
64 | 1 | pfleura2 | * contain the upper triangular matrix packed sequentially, |
65 | 1 | pfleura2 | * column by column, so that AP( 1 ) contains a( 1, 1 ), |
66 | 1 | pfleura2 | * AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 ) |
67 | 1 | pfleura2 | * respectively, and so on. |
68 | 1 | pfleura2 | * Before entry with UPLO = 'L' or 'l', the array AP must |
69 | 1 | pfleura2 | * contain the lower triangular matrix packed sequentially, |
70 | 1 | pfleura2 | * column by column, so that AP( 1 ) contains a( 1, 1 ), |
71 | 1 | pfleura2 | * AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 ) |
72 | 1 | pfleura2 | * respectively, and so on. |
73 | 1 | pfleura2 | * Note that when DIAG = 'U' or 'u', the diagonal elements of |
74 | 1 | pfleura2 | * A are not referenced, but are assumed to be unity. |
75 | 1 | pfleura2 | * Unchanged on exit. |
76 | 1 | pfleura2 | * |
77 | 1 | pfleura2 | * X - DOUBLE PRECISION array of dimension at least |
78 | 1 | pfleura2 | * ( 1 + ( n - 1 )*abs( INCX ) ). |
79 | 1 | pfleura2 | * Before entry, the incremented array X must contain the n |
80 | 1 | pfleura2 | * element vector x. On exit, X is overwritten with the |
81 | 1 | pfleura2 | * tranformed vector x. |
82 | 1 | pfleura2 | * |
83 | 1 | pfleura2 | * INCX - INTEGER. |
84 | 1 | pfleura2 | * On entry, INCX specifies the increment for the elements of |
85 | 1 | pfleura2 | * X. INCX must not be zero. |
86 | 1 | pfleura2 | * Unchanged on exit. |
87 | 1 | pfleura2 | * |
88 | 1 | pfleura2 | * |
89 | 1 | pfleura2 | * Level 2 Blas routine. |
90 | 1 | pfleura2 | * |
91 | 1 | pfleura2 | * -- Written on 22-October-1986. |
92 | 1 | pfleura2 | * Jack Dongarra, Argonne National Lab. |
93 | 1 | pfleura2 | * Jeremy Du Croz, Nag Central Office. |
94 | 1 | pfleura2 | * Sven Hammarling, Nag Central Office. |
95 | 1 | pfleura2 | * Richard Hanson, Sandia National Labs. |
96 | 1 | pfleura2 | * |
97 | 1 | pfleura2 | * |
98 | 1 | pfleura2 | * .. Parameters .. |
99 | 1 | pfleura2 | DOUBLE PRECISION ZERO |
100 | 1 | pfleura2 | PARAMETER (ZERO=0.0D+0) |
101 | 1 | pfleura2 | * .. |
102 | 1 | pfleura2 | * .. Local Scalars .. |
103 | 1 | pfleura2 | DOUBLE PRECISION TEMP |
104 | 1 | pfleura2 | INTEGER I,INFO,IX,J,JX,K,KK,KX |
105 | 1 | pfleura2 | LOGICAL NOUNIT |
106 | 1 | pfleura2 | * .. |
107 | 1 | pfleura2 | * .. External Functions .. |
108 | 1 | pfleura2 | LOGICAL LSAME |
109 | 1 | pfleura2 | EXTERNAL LSAME |
110 | 1 | pfleura2 | * .. |
111 | 1 | pfleura2 | * .. External Subroutines .. |
112 | 1 | pfleura2 | EXTERNAL XERBLA |
113 | 1 | pfleura2 | * .. |
114 | 1 | pfleura2 | * |
115 | 1 | pfleura2 | * Test the input parameters. |
116 | 1 | pfleura2 | * |
117 | 1 | pfleura2 | INFO = 0 |
118 | 1 | pfleura2 | IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN |
119 | 1 | pfleura2 | INFO = 1 |
120 | 1 | pfleura2 | ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND. |
121 | 1 | pfleura2 | + .NOT.LSAME(TRANS,'C')) THEN |
122 | 1 | pfleura2 | INFO = 2 |
123 | 1 | pfleura2 | ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN |
124 | 1 | pfleura2 | INFO = 3 |
125 | 1 | pfleura2 | ELSE IF (N.LT.0) THEN |
126 | 1 | pfleura2 | INFO = 4 |
127 | 1 | pfleura2 | ELSE IF (INCX.EQ.0) THEN |
128 | 1 | pfleura2 | INFO = 7 |
129 | 1 | pfleura2 | END IF |
130 | 1 | pfleura2 | IF (INFO.NE.0) THEN |
131 | 1 | pfleura2 | CALL XERBLA('DTPMV ',INFO) |
132 | 1 | pfleura2 | RETURN |
133 | 1 | pfleura2 | END IF |
134 | 1 | pfleura2 | * |
135 | 1 | pfleura2 | * Quick return if possible. |
136 | 1 | pfleura2 | * |
137 | 1 | pfleura2 | IF (N.EQ.0) RETURN |
138 | 1 | pfleura2 | * |
139 | 1 | pfleura2 | NOUNIT = LSAME(DIAG,'N') |
140 | 1 | pfleura2 | * |
141 | 1 | pfleura2 | * Set up the start point in X if the increment is not unity. This |
142 | 1 | pfleura2 | * will be ( N - 1 )*INCX too small for descending loops. |
143 | 1 | pfleura2 | * |
144 | 1 | pfleura2 | IF (INCX.LE.0) THEN |
145 | 1 | pfleura2 | KX = 1 - (N-1)*INCX |
146 | 1 | pfleura2 | ELSE IF (INCX.NE.1) THEN |
147 | 1 | pfleura2 | KX = 1 |
148 | 1 | pfleura2 | END IF |
149 | 1 | pfleura2 | * |
150 | 1 | pfleura2 | * Start the operations. In this version the elements of AP are |
151 | 1 | pfleura2 | * accessed sequentially with one pass through AP. |
152 | 1 | pfleura2 | * |
153 | 1 | pfleura2 | IF (LSAME(TRANS,'N')) THEN |
154 | 1 | pfleura2 | * |
155 | 1 | pfleura2 | * Form x:= A*x. |
156 | 1 | pfleura2 | * |
157 | 1 | pfleura2 | IF (LSAME(UPLO,'U')) THEN |
158 | 1 | pfleura2 | KK = 1 |
159 | 1 | pfleura2 | IF (INCX.EQ.1) THEN |
160 | 1 | pfleura2 | DO 20 J = 1,N |
161 | 1 | pfleura2 | IF (X(J).NE.ZERO) THEN |
162 | 1 | pfleura2 | TEMP = X(J) |
163 | 1 | pfleura2 | K = KK |
164 | 1 | pfleura2 | DO 10 I = 1,J - 1 |
165 | 1 | pfleura2 | X(I) = X(I) + TEMP*AP(K) |
166 | 1 | pfleura2 | K = K + 1 |
167 | 1 | pfleura2 | 10 CONTINUE |
168 | 1 | pfleura2 | IF (NOUNIT) X(J) = X(J)*AP(KK+J-1) |
169 | 1 | pfleura2 | END IF |
170 | 1 | pfleura2 | KK = KK + J |
171 | 1 | pfleura2 | 20 CONTINUE |
172 | 1 | pfleura2 | ELSE |
173 | 1 | pfleura2 | JX = KX |
174 | 1 | pfleura2 | DO 40 J = 1,N |
175 | 1 | pfleura2 | IF (X(JX).NE.ZERO) THEN |
176 | 1 | pfleura2 | TEMP = X(JX) |
177 | 1 | pfleura2 | IX = KX |
178 | 1 | pfleura2 | DO 30 K = KK,KK + J - 2 |
179 | 1 | pfleura2 | X(IX) = X(IX) + TEMP*AP(K) |
180 | 1 | pfleura2 | IX = IX + INCX |
181 | 1 | pfleura2 | 30 CONTINUE |
182 | 1 | pfleura2 | IF (NOUNIT) X(JX) = X(JX)*AP(KK+J-1) |
183 | 1 | pfleura2 | END IF |
184 | 1 | pfleura2 | JX = JX + INCX |
185 | 1 | pfleura2 | KK = KK + J |
186 | 1 | pfleura2 | 40 CONTINUE |
187 | 1 | pfleura2 | END IF |
188 | 1 | pfleura2 | ELSE |
189 | 1 | pfleura2 | KK = (N* (N+1))/2 |
190 | 1 | pfleura2 | IF (INCX.EQ.1) THEN |
191 | 1 | pfleura2 | DO 60 J = N,1,-1 |
192 | 1 | pfleura2 | IF (X(J).NE.ZERO) THEN |
193 | 1 | pfleura2 | TEMP = X(J) |
194 | 1 | pfleura2 | K = KK |
195 | 1 | pfleura2 | DO 50 I = N,J + 1,-1 |
196 | 1 | pfleura2 | X(I) = X(I) + TEMP*AP(K) |
197 | 1 | pfleura2 | K = K - 1 |
198 | 1 | pfleura2 | 50 CONTINUE |
199 | 1 | pfleura2 | IF (NOUNIT) X(J) = X(J)*AP(KK-N+J) |
200 | 1 | pfleura2 | END IF |
201 | 1 | pfleura2 | KK = KK - (N-J+1) |
202 | 1 | pfleura2 | 60 CONTINUE |
203 | 1 | pfleura2 | ELSE |
204 | 1 | pfleura2 | KX = KX + (N-1)*INCX |
205 | 1 | pfleura2 | JX = KX |
206 | 1 | pfleura2 | DO 80 J = N,1,-1 |
207 | 1 | pfleura2 | IF (X(JX).NE.ZERO) THEN |
208 | 1 | pfleura2 | TEMP = X(JX) |
209 | 1 | pfleura2 | IX = KX |
210 | 1 | pfleura2 | DO 70 K = KK,KK - (N- (J+1)),-1 |
211 | 1 | pfleura2 | X(IX) = X(IX) + TEMP*AP(K) |
212 | 1 | pfleura2 | IX = IX - INCX |
213 | 1 | pfleura2 | 70 CONTINUE |
214 | 1 | pfleura2 | IF (NOUNIT) X(JX) = X(JX)*AP(KK-N+J) |
215 | 1 | pfleura2 | END IF |
216 | 1 | pfleura2 | JX = JX - INCX |
217 | 1 | pfleura2 | KK = KK - (N-J+1) |
218 | 1 | pfleura2 | 80 CONTINUE |
219 | 1 | pfleura2 | END IF |
220 | 1 | pfleura2 | END IF |
221 | 1 | pfleura2 | ELSE |
222 | 1 | pfleura2 | * |
223 | 1 | pfleura2 | * Form x := A'*x. |
224 | 1 | pfleura2 | * |
225 | 1 | pfleura2 | IF (LSAME(UPLO,'U')) THEN |
226 | 1 | pfleura2 | KK = (N* (N+1))/2 |
227 | 1 | pfleura2 | IF (INCX.EQ.1) THEN |
228 | 1 | pfleura2 | DO 100 J = N,1,-1 |
229 | 1 | pfleura2 | TEMP = X(J) |
230 | 1 | pfleura2 | IF (NOUNIT) TEMP = TEMP*AP(KK) |
231 | 1 | pfleura2 | K = KK - 1 |
232 | 1 | pfleura2 | DO 90 I = J - 1,1,-1 |
233 | 1 | pfleura2 | TEMP = TEMP + AP(K)*X(I) |
234 | 1 | pfleura2 | K = K - 1 |
235 | 1 | pfleura2 | 90 CONTINUE |
236 | 1 | pfleura2 | X(J) = TEMP |
237 | 1 | pfleura2 | KK = KK - J |
238 | 1 | pfleura2 | 100 CONTINUE |
239 | 1 | pfleura2 | ELSE |
240 | 1 | pfleura2 | JX = KX + (N-1)*INCX |
241 | 1 | pfleura2 | DO 120 J = N,1,-1 |
242 | 1 | pfleura2 | TEMP = X(JX) |
243 | 1 | pfleura2 | IX = JX |
244 | 1 | pfleura2 | IF (NOUNIT) TEMP = TEMP*AP(KK) |
245 | 1 | pfleura2 | DO 110 K = KK - 1,KK - J + 1,-1 |
246 | 1 | pfleura2 | IX = IX - INCX |
247 | 1 | pfleura2 | TEMP = TEMP + AP(K)*X(IX) |
248 | 1 | pfleura2 | 110 CONTINUE |
249 | 1 | pfleura2 | X(JX) = TEMP |
250 | 1 | pfleura2 | JX = JX - INCX |
251 | 1 | pfleura2 | KK = KK - J |
252 | 1 | pfleura2 | 120 CONTINUE |
253 | 1 | pfleura2 | END IF |
254 | 1 | pfleura2 | ELSE |
255 | 1 | pfleura2 | KK = 1 |
256 | 1 | pfleura2 | IF (INCX.EQ.1) THEN |
257 | 1 | pfleura2 | DO 140 J = 1,N |
258 | 1 | pfleura2 | TEMP = X(J) |
259 | 1 | pfleura2 | IF (NOUNIT) TEMP = TEMP*AP(KK) |
260 | 1 | pfleura2 | K = KK + 1 |
261 | 1 | pfleura2 | DO 130 I = J + 1,N |
262 | 1 | pfleura2 | TEMP = TEMP + AP(K)*X(I) |
263 | 1 | pfleura2 | K = K + 1 |
264 | 1 | pfleura2 | 130 CONTINUE |
265 | 1 | pfleura2 | X(J) = TEMP |
266 | 1 | pfleura2 | KK = KK + (N-J+1) |
267 | 1 | pfleura2 | 140 CONTINUE |
268 | 1 | pfleura2 | ELSE |
269 | 1 | pfleura2 | JX = KX |
270 | 1 | pfleura2 | DO 160 J = 1,N |
271 | 1 | pfleura2 | TEMP = X(JX) |
272 | 1 | pfleura2 | IX = JX |
273 | 1 | pfleura2 | IF (NOUNIT) TEMP = TEMP*AP(KK) |
274 | 1 | pfleura2 | DO 150 K = KK + 1,KK + N - J |
275 | 1 | pfleura2 | IX = IX + INCX |
276 | 1 | pfleura2 | TEMP = TEMP + AP(K)*X(IX) |
277 | 1 | pfleura2 | 150 CONTINUE |
278 | 1 | pfleura2 | X(JX) = TEMP |
279 | 1 | pfleura2 | JX = JX + INCX |
280 | 1 | pfleura2 | KK = KK + (N-J+1) |
281 | 1 | pfleura2 | 160 CONTINUE |
282 | 1 | pfleura2 | END IF |
283 | 1 | pfleura2 | END IF |
284 | 1 | pfleura2 | END IF |
285 | 1 | pfleura2 | * |
286 | 1 | pfleura2 | RETURN |
287 | 1 | pfleura2 | * |
288 | 1 | pfleura2 | * End of DTPMV . |
289 | 1 | pfleura2 | * |
290 | 1 | pfleura2 | END |