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