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