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