1 /* blas/dtrsv.f -- translated by f2c (version 20050501).
2 You must link the resulting object file with libf2c:
3 on Microsoft Windows system, link with libf2c.lib;
4 on Linux or Unix systems, link with .../path/to/libf2c.a -lm
5 or, if you install libf2c.a in a standard place, with -lf2c -lm
6 -- in that order, at the end of the command line, as in
7 cc *.o -lf2c -lm
8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
9
10 http://www.netlib.org/f2c/libf2c.zip
11 */
12
13 #ifdef __cplusplus
14 extern "C" {
15 #endif
16 #include "v3p_netlib.h"
17
18 /*< SUBROUTINE DTRSV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX ) >*/
dtrsv_(char * uplo,char * trans,char * diag,integer * n,doublereal * a,integer * lda,doublereal * x,integer * incx,ftnlen uplo_len,ftnlen trans_len,ftnlen diag_len)19 /* Subroutine */ int dtrsv_(char *uplo, char *trans, char *diag, integer *n,
20 doublereal *a, integer *lda, doublereal *x, integer *incx, ftnlen
21 uplo_len, ftnlen trans_len, ftnlen diag_len)
22 {
23 /* System generated locals */
24 integer a_dim1, a_offset, i__1, i__2;
25
26 /* Local variables */
27 integer i__, j, ix, jx, kx=0, info;
28 doublereal temp;
29 extern logical lsame_(const char *, const char *, ftnlen, ftnlen);
30 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
31 logical nounit;
32 (void)uplo_len;
33 (void)trans_len;
34 (void)diag_len;
35
36 /* .. Scalar Arguments .. */
37 /*< INTEGER INCX, LDA, N >*/
38 /*< CHARACTER*1 DIAG, TRANS, UPLO >*/
39 /* .. Array Arguments .. */
40 /*< DOUBLE PRECISION A( LDA, * ), X( * ) >*/
41 /* .. */
42
43 /* Purpose */
44 /* ======= */
45
46 /* DTRSV solves one of the systems of equations */
47
48 /* A*x = b, or A'*x = b, */
49
50 /* where b and x are n element vectors and A is an n by n unit, or */
51 /* non-unit, upper or lower triangular matrix. */
52
53 /* No test for singularity or near-singularity is included in this */
54 /* routine. Such tests must be performed before calling this routine. */
55
56 /* Parameters */
57 /* ========== */
58
59 /* UPLO - CHARACTER*1. */
60 /* On entry, UPLO specifies whether the matrix is an upper or */
61 /* lower triangular matrix as follows: */
62
63 /* UPLO = 'U' or 'u' A is an upper triangular matrix. */
64
65 /* UPLO = 'L' or 'l' A is a lower triangular matrix. */
66
67 /* Unchanged on exit. */
68
69 /* TRANS - CHARACTER*1. */
70 /* On entry, TRANS specifies the equations to be solved as */
71 /* follows: */
72
73 /* TRANS = 'N' or 'n' A*x = b. */
74
75 /* TRANS = 'T' or 't' A'*x = b. */
76
77 /* TRANS = 'C' or 'c' A'*x = b. */
78
79 /* Unchanged on exit. */
80
81 /* DIAG - CHARACTER*1. */
82 /* On entry, DIAG specifies whether or not A is unit */
83 /* triangular as follows: */
84
85 /* DIAG = 'U' or 'u' A is assumed to be unit triangular. */
86
87 /* DIAG = 'N' or 'n' A is not assumed to be unit */
88 /* triangular. */
89
90 /* Unchanged on exit. */
91
92 /* N - INTEGER. */
93 /* On entry, N specifies the order of the matrix A. */
94 /* N must be at least zero. */
95 /* Unchanged on exit. */
96
97 /* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). */
98 /* Before entry with UPLO = 'U' or 'u', the leading n by n */
99 /* upper triangular part of the array A must contain the upper */
100 /* triangular matrix and the strictly lower triangular part of */
101 /* A is not referenced. */
102 /* Before entry with UPLO = 'L' or 'l', the leading n by n */
103 /* lower triangular part of the array A must contain the lower */
104 /* triangular matrix and the strictly upper triangular part of */
105 /* A is not referenced. */
106 /* Note that when DIAG = 'U' or 'u', the diagonal elements of */
107 /* A are not referenced either, but are assumed to be unity. */
108 /* Unchanged on exit. */
109
110 /* LDA - INTEGER. */
111 /* On entry, LDA specifies the first dimension of A as declared */
112 /* in the calling (sub) program. LDA must be at least */
113 /* max( 1, n ). */
114 /* Unchanged on exit. */
115
116 /* X - DOUBLE PRECISION array of dimension at least */
117 /* ( 1 + ( n - 1 )*abs( INCX ) ). */
118 /* Before entry, the incremented array X must contain the n */
119 /* element right-hand side vector b. On exit, X is overwritten */
120 /* with the solution vector x. */
121
122 /* INCX - INTEGER. */
123 /* On entry, INCX specifies the increment for the elements of */
124 /* X. INCX must not be zero. */
125 /* Unchanged on exit. */
126
127
128 /* Level 2 Blas routine. */
129
130 /* -- Written on 22-October-1986. */
131 /* Jack Dongarra, Argonne National Lab. */
132 /* Jeremy Du Croz, Nag Central Office. */
133 /* Sven Hammarling, Nag Central Office. */
134 /* Richard Hanson, Sandia National Labs. */
135
136
137 /* .. Parameters .. */
138 /*< DOUBLE PRECISION ZERO >*/
139 /*< PARAMETER ( ZERO = 0.0D+0 ) >*/
140 /* .. Local Scalars .. */
141 /*< DOUBLE PRECISION TEMP >*/
142 /*< INTEGER I, INFO, IX, J, JX, KX >*/
143 /*< LOGICAL NOUNIT >*/
144 /* .. External Functions .. */
145 /*< LOGICAL LSAME >*/
146 /*< EXTERNAL LSAME >*/
147 /* .. External Subroutines .. */
148 /*< EXTERNAL XERBLA >*/
149 /* .. Intrinsic Functions .. */
150 /*< INTRINSIC MAX >*/
151 /* .. */
152 /* .. Executable Statements .. */
153
154 /* Test the input parameters. */
155
156 /*< INFO = 0 >*/
157 /* Parameter adjustments */
158 a_dim1 = *lda;
159 a_offset = 1 + a_dim1;
160 a -= a_offset;
161 --x;
162
163 /* Function Body */
164 info = 0;
165 /*< >*/
166 if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo, "L", (
167 ftnlen)1, (ftnlen)1)) {
168 /*< INFO = 1 >*/
169 info = 1;
170 /*< >*/
171 } else if (! lsame_(trans, "N", (ftnlen)1, (ftnlen)1) && ! lsame_(trans,
172 "T", (ftnlen)1, (ftnlen)1) && ! lsame_(trans, "C", (ftnlen)1, (
173 ftnlen)1)) {
174 /*< INFO = 2 >*/
175 info = 2;
176 /*< >*/
177 } else if (! lsame_(diag, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(diag,
178 "N", (ftnlen)1, (ftnlen)1)) {
179 /*< INFO = 3 >*/
180 info = 3;
181 /*< ELSE IF( N.LT.0 )THEN >*/
182 } else if (*n < 0) {
183 /*< INFO = 4 >*/
184 info = 4;
185 /*< ELSE IF( LDA.LT.MAX( 1, N ) )THEN >*/
186 } else if (*lda < max(1,*n)) {
187 /*< INFO = 6 >*/
188 info = 6;
189 /*< ELSE IF( INCX.EQ.0 )THEN >*/
190 } else if (*incx == 0) {
191 /*< INFO = 8 >*/
192 info = 8;
193 /*< END IF >*/
194 }
195 /*< IF( INFO.NE.0 )THEN >*/
196 if (info != 0) {
197 /*< CALL XERBLA( 'DTRSV ', INFO ) >*/
198 xerbla_("DTRSV ", &info, (ftnlen)6);
199 /*< RETURN >*/
200 return 0;
201 /*< END IF >*/
202 }
203
204 /* Quick return if possible. */
205
206 /*< >*/
207 if (*n == 0) {
208 return 0;
209 }
210
211 /*< NOUNIT = LSAME( DIAG, 'N' ) >*/
212 nounit = lsame_(diag, "N", (ftnlen)1, (ftnlen)1);
213
214 /* Set up the start point in X if the increment is not unity. This */
215 /* will be ( N - 1 )*INCX too small for descending loops. */
216
217 /*< IF( INCX.LE.0 )THEN >*/
218 if (*incx <= 0) {
219 /*< KX = 1 - ( N - 1 )*INCX >*/
220 kx = 1 - (*n - 1) * *incx;
221 /*< ELSE IF( INCX.NE.1 )THEN >*/
222 } else if (*incx != 1) {
223 /*< KX = 1 >*/
224 kx = 1;
225 /*< END IF >*/
226 }
227
228 /* Start the operations. In this version the elements of A are */
229 /* accessed sequentially with one pass through A. */
230
231 /*< IF( LSAME( TRANS, 'N' ) )THEN >*/
232 if (lsame_(trans, "N", (ftnlen)1, (ftnlen)1)) {
233
234 /* Form x := inv( A )*x. */
235
236 /*< IF( LSAME( UPLO, 'U' ) )THEN >*/
237 if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
238 /*< IF( INCX.EQ.1 )THEN >*/
239 if (*incx == 1) {
240 /*< DO 20, J = N, 1, -1 >*/
241 for (j = *n; j >= 1; --j) {
242 /*< IF( X( J ).NE.ZERO )THEN >*/
243 if (x[j] != 0.) {
244 /*< >*/
245 if (nounit) {
246 x[j] /= a[j + j * a_dim1];
247 }
248 /*< TEMP = X( J ) >*/
249 temp = x[j];
250 /*< DO 10, I = J - 1, 1, -1 >*/
251 for (i__ = j - 1; i__ >= 1; --i__) {
252 /*< X( I ) = X( I ) - TEMP*A( I, J ) >*/
253 x[i__] -= temp * a[i__ + j * a_dim1];
254 /*< 10 CONTINUE >*/
255 /* L10: */
256 }
257 /*< END IF >*/
258 }
259 /*< 20 CONTINUE >*/
260 /* L20: */
261 }
262 /*< ELSE >*/
263 } else {
264 /*< JX = KX + ( N - 1 )*INCX >*/
265 jx = kx + (*n - 1) * *incx;
266 /*< DO 40, J = N, 1, -1 >*/
267 for (j = *n; j >= 1; --j) {
268 /*< IF( X( JX ).NE.ZERO )THEN >*/
269 if (x[jx] != 0.) {
270 /*< >*/
271 if (nounit) {
272 x[jx] /= a[j + j * a_dim1];
273 }
274 /*< TEMP = X( JX ) >*/
275 temp = x[jx];
276 /*< IX = JX >*/
277 ix = jx;
278 /*< DO 30, I = J - 1, 1, -1 >*/
279 for (i__ = j - 1; i__ >= 1; --i__) {
280 /*< IX = IX - INCX >*/
281 ix -= *incx;
282 /*< X( IX ) = X( IX ) - TEMP*A( I, J ) >*/
283 x[ix] -= temp * a[i__ + j * a_dim1];
284 /*< 30 CONTINUE >*/
285 /* L30: */
286 }
287 /*< END IF >*/
288 }
289 /*< JX = JX - INCX >*/
290 jx -= *incx;
291 /*< 40 CONTINUE >*/
292 /* L40: */
293 }
294 /*< END IF >*/
295 }
296 /*< ELSE >*/
297 } else {
298 /*< IF( INCX.EQ.1 )THEN >*/
299 if (*incx == 1) {
300 /*< DO 60, J = 1, N >*/
301 i__1 = *n;
302 for (j = 1; j <= i__1; ++j) {
303 /*< IF( X( J ).NE.ZERO )THEN >*/
304 if (x[j] != 0.) {
305 /*< >*/
306 if (nounit) {
307 x[j] /= a[j + j * a_dim1];
308 }
309 /*< TEMP = X( J ) >*/
310 temp = x[j];
311 /*< DO 50, I = J + 1, N >*/
312 i__2 = *n;
313 for (i__ = j + 1; i__ <= i__2; ++i__) {
314 /*< X( I ) = X( I ) - TEMP*A( I, J ) >*/
315 x[i__] -= temp * a[i__ + j * a_dim1];
316 /*< 50 CONTINUE >*/
317 /* L50: */
318 }
319 /*< END IF >*/
320 }
321 /*< 60 CONTINUE >*/
322 /* L60: */
323 }
324 /*< ELSE >*/
325 } else {
326 /*< JX = KX >*/
327 jx = kx;
328 /*< DO 80, J = 1, N >*/
329 i__1 = *n;
330 for (j = 1; j <= i__1; ++j) {
331 /*< IF( X( JX ).NE.ZERO )THEN >*/
332 if (x[jx] != 0.) {
333 /*< >*/
334 if (nounit) {
335 x[jx] /= a[j + j * a_dim1];
336 }
337 /*< TEMP = X( JX ) >*/
338 temp = x[jx];
339 /*< IX = JX >*/
340 ix = jx;
341 /*< DO 70, I = J + 1, N >*/
342 i__2 = *n;
343 for (i__ = j + 1; i__ <= i__2; ++i__) {
344 /*< IX = IX + INCX >*/
345 ix += *incx;
346 /*< X( IX ) = X( IX ) - TEMP*A( I, J ) >*/
347 x[ix] -= temp * a[i__ + j * a_dim1];
348 /*< 70 CONTINUE >*/
349 /* L70: */
350 }
351 /*< END IF >*/
352 }
353 /*< JX = JX + INCX >*/
354 jx += *incx;
355 /*< 80 CONTINUE >*/
356 /* L80: */
357 }
358 /*< END IF >*/
359 }
360 /*< END IF >*/
361 }
362 /*< ELSE >*/
363 } else {
364
365 /* Form x := inv( A' )*x. */
366
367 /*< IF( LSAME( UPLO, 'U' ) )THEN >*/
368 if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
369 /*< IF( INCX.EQ.1 )THEN >*/
370 if (*incx == 1) {
371 /*< DO 100, J = 1, N >*/
372 i__1 = *n;
373 for (j = 1; j <= i__1; ++j) {
374 /*< TEMP = X( J ) >*/
375 temp = x[j];
376 /*< DO 90, I = 1, J - 1 >*/
377 i__2 = j - 1;
378 for (i__ = 1; i__ <= i__2; ++i__) {
379 /*< TEMP = TEMP - A( I, J )*X( I ) >*/
380 temp -= a[i__ + j * a_dim1] * x[i__];
381 /*< 90 CONTINUE >*/
382 /* L90: */
383 }
384 /*< >*/
385 if (nounit) {
386 temp /= a[j + j * a_dim1];
387 }
388 /*< X( J ) = TEMP >*/
389 x[j] = temp;
390 /*< 100 CONTINUE >*/
391 /* L100: */
392 }
393 /*< ELSE >*/
394 } else {
395 /*< JX = KX >*/
396 jx = kx;
397 /*< DO 120, J = 1, N >*/
398 i__1 = *n;
399 for (j = 1; j <= i__1; ++j) {
400 /*< TEMP = X( JX ) >*/
401 temp = x[jx];
402 /*< IX = KX >*/
403 ix = kx;
404 /*< DO 110, I = 1, J - 1 >*/
405 i__2 = j - 1;
406 for (i__ = 1; i__ <= i__2; ++i__) {
407 /*< TEMP = TEMP - A( I, J )*X( IX ) >*/
408 temp -= a[i__ + j * a_dim1] * x[ix];
409 /*< IX = IX + INCX >*/
410 ix += *incx;
411 /*< 110 CONTINUE >*/
412 /* L110: */
413 }
414 /*< >*/
415 if (nounit) {
416 temp /= a[j + j * a_dim1];
417 }
418 /*< X( JX ) = TEMP >*/
419 x[jx] = temp;
420 /*< JX = JX + INCX >*/
421 jx += *incx;
422 /*< 120 CONTINUE >*/
423 /* L120: */
424 }
425 /*< END IF >*/
426 }
427 /*< ELSE >*/
428 } else {
429 /*< IF( INCX.EQ.1 )THEN >*/
430 if (*incx == 1) {
431 /*< DO 140, J = N, 1, -1 >*/
432 for (j = *n; j >= 1; --j) {
433 /*< TEMP = X( J ) >*/
434 temp = x[j];
435 /*< DO 130, I = N, J + 1, -1 >*/
436 i__1 = j + 1;
437 for (i__ = *n; i__ >= i__1; --i__) {
438 /*< TEMP = TEMP - A( I, J )*X( I ) >*/
439 temp -= a[i__ + j * a_dim1] * x[i__];
440 /*< 130 CONTINUE >*/
441 /* L130: */
442 }
443 /*< >*/
444 if (nounit) {
445 temp /= a[j + j * a_dim1];
446 }
447 /*< X( J ) = TEMP >*/
448 x[j] = temp;
449 /*< 140 CONTINUE >*/
450 /* L140: */
451 }
452 /*< ELSE >*/
453 } else {
454 /*< KX = KX + ( N - 1 )*INCX >*/
455 kx += (*n - 1) * *incx;
456 /*< JX = KX >*/
457 jx = kx;
458 /*< DO 160, J = N, 1, -1 >*/
459 for (j = *n; j >= 1; --j) {
460 /*< TEMP = X( JX ) >*/
461 temp = x[jx];
462 /*< IX = KX >*/
463 ix = kx;
464 /*< DO 150, I = N, J + 1, -1 >*/
465 i__1 = j + 1;
466 for (i__ = *n; i__ >= i__1; --i__) {
467 /*< TEMP = TEMP - A( I, J )*X( IX ) >*/
468 temp -= a[i__ + j * a_dim1] * x[ix];
469 /*< IX = IX - INCX >*/
470 ix -= *incx;
471 /*< 150 CONTINUE >*/
472 /* L150: */
473 }
474 /*< >*/
475 if (nounit) {
476 temp /= a[j + j * a_dim1];
477 }
478 /*< X( JX ) = TEMP >*/
479 x[jx] = temp;
480 /*< JX = JX - INCX >*/
481 jx -= *incx;
482 /*< 160 CONTINUE >*/
483 /* L160: */
484 }
485 /*< END IF >*/
486 }
487 /*< END IF >*/
488 }
489 /*< END IF >*/
490 }
491
492 /*< RETURN >*/
493 return 0;
494
495 /* End of DTRSV . */
496
497 /*< END >*/
498 } /* dtrsv_ */
499
500 #ifdef __cplusplus
501 }
502 #endif
503