1 /* ./src_f77/zlatps.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
8 /* Table of constant values */
9
10 static integer c__1 = 1;
11 static doublereal c_b36 = .5;
12
zlatps_(char * uplo,char * trans,char * diag,char * normin,integer * n,doublecomplex * ap,doublecomplex * x,doublereal * scale,doublereal * cnorm,integer * info,ftnlen uplo_len,ftnlen trans_len,ftnlen diag_len,ftnlen normin_len)13 /* Subroutine */ int zlatps_(char *uplo, char *trans, char *diag, char *
14 normin, integer *n, doublecomplex *ap, doublecomplex *x, doublereal *
15 scale, doublereal *cnorm, integer *info, ftnlen uplo_len, ftnlen
16 trans_len, ftnlen diag_len, ftnlen normin_len)
17 {
18 /* System generated locals */
19 integer i__1, i__2, i__3, i__4, i__5;
20 doublereal d__1, d__2, d__3, d__4;
21 doublecomplex z__1, z__2, z__3, z__4;
22
23 /* Builtin functions */
24 double d_imag(doublecomplex *);
25 void d_cnjg(doublecomplex *, doublecomplex *);
26
27 /* Local variables */
28 static integer i__, j, ip;
29 static doublereal xj, rec, tjj;
30 static integer jinc, jlen;
31 static doublereal xbnd;
32 static integer imax;
33 static doublereal tmax;
34 static doublecomplex tjjs;
35 static doublereal xmax, grow;
36 extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
37 integer *);
38 extern logical lsame_(char *, char *, ftnlen, ftnlen);
39 static doublereal tscal;
40 static doublecomplex uscal;
41 static integer jlast;
42 static doublecomplex csumj;
43 extern /* Double Complex */ VOID zdotc_(doublecomplex *, integer *,
44 doublecomplex *, integer *, doublecomplex *, integer *);
45 static logical upper;
46 extern /* Double Complex */ VOID zdotu_(doublecomplex *, integer *,
47 doublecomplex *, integer *, doublecomplex *, integer *);
48 extern /* Subroutine */ int zaxpy_(integer *, doublecomplex *,
49 doublecomplex *, integer *, doublecomplex *, integer *), ztpsv_(
50 char *, char *, char *, integer *, doublecomplex *, doublecomplex
51 *, integer *, ftnlen, ftnlen, ftnlen), dlabad_(doublereal *,
52 doublereal *);
53 extern doublereal dlamch_(char *, ftnlen);
54 extern integer idamax_(integer *, doublereal *, integer *);
55 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen), zdscal_(
56 integer *, doublereal *, doublecomplex *, integer *);
57 static doublereal bignum;
58 extern integer izamax_(integer *, doublecomplex *, integer *);
59 extern /* Double Complex */ VOID zladiv_(doublecomplex *, doublecomplex *,
60 doublecomplex *);
61 static logical notran;
62 static integer jfirst;
63 extern doublereal dzasum_(integer *, doublecomplex *, integer *);
64 static doublereal smlnum;
65 static logical nounit;
66
67
68 /* -- LAPACK auxiliary routine (version 3.0) -- */
69 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
70 /* Courant Institute, Argonne National Lab, and Rice University */
71 /* October 31, 1992 */
72
73 /* .. Scalar Arguments .. */
74 /* .. */
75 /* .. Array Arguments .. */
76 /* .. */
77
78 /* Purpose */
79 /* ======= */
80
81 /* ZLATPS solves one of the triangular systems */
82
83 /* A * x = s*b, A**T * x = s*b, or A**H * x = s*b, */
84
85 /* with scaling to prevent overflow, where A is an upper or lower */
86 /* triangular matrix stored in packed form. Here A**T denotes the */
87 /* transpose of A, A**H denotes the conjugate transpose of A, x and b */
88 /* are n-element vectors, and s is a scaling factor, usually less than */
89 /* or equal to 1, chosen so that the components of x will be less than */
90 /* the overflow threshold. If the unscaled problem will not cause */
91 /* overflow, the Level 2 BLAS routine ZTPSV is called. If the matrix A */
92 /* is singular (A(j,j) = 0 for some j), then s is set to 0 and a */
93 /* non-trivial solution to A*x = 0 is returned. */
94
95 /* Arguments */
96 /* ========= */
97
98 /* UPLO (input) CHARACTER*1 */
99 /* Specifies whether the matrix A is upper or lower triangular. */
100 /* = 'U': Upper triangular */
101 /* = 'L': Lower triangular */
102
103 /* TRANS (input) CHARACTER*1 */
104 /* Specifies the operation applied to A. */
105 /* = 'N': Solve A * x = s*b (No transpose) */
106 /* = 'T': Solve A**T * x = s*b (Transpose) */
107 /* = 'C': Solve A**H * x = s*b (Conjugate transpose) */
108
109 /* DIAG (input) CHARACTER*1 */
110 /* Specifies whether or not the matrix A is unit triangular. */
111 /* = 'N': Non-unit triangular */
112 /* = 'U': Unit triangular */
113
114 /* NORMIN (input) CHARACTER*1 */
115 /* Specifies whether CNORM has been set or not. */
116 /* = 'Y': CNORM contains the column norms on entry */
117 /* = 'N': CNORM is not set on entry. On exit, the norms will */
118 /* be computed and stored in CNORM. */
119
120 /* N (input) INTEGER */
121 /* The order of the matrix A. N >= 0. */
122
123 /* AP (input) COMPLEX*16 array, dimension (N*(N+1)/2) */
124 /* The upper or lower triangular matrix A, packed columnwise in */
125 /* a linear array. The j-th column of A is stored in the array */
126 /* AP as follows: */
127 /* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; */
128 /* if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. */
129
130 /* X (input/output) COMPLEX*16 array, dimension (N) */
131 /* On entry, the right hand side b of the triangular system. */
132 /* On exit, X is overwritten by the solution vector x. */
133
134 /* SCALE (output) DOUBLE PRECISION */
135 /* The scaling factor s for the triangular system */
136 /* A * x = s*b, A**T * x = s*b, or A**H * x = s*b. */
137 /* If SCALE = 0, the matrix A is singular or badly scaled, and */
138 /* the vector x is an exact or approximate solution to A*x = 0. */
139
140 /* CNORM (input or output) DOUBLE PRECISION array, dimension (N) */
141
142 /* If NORMIN = 'Y', CNORM is an input argument and CNORM(j) */
143 /* contains the norm of the off-diagonal part of the j-th column */
144 /* of A. If TRANS = 'N', CNORM(j) must be greater than or equal */
145 /* to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) */
146 /* must be greater than or equal to the 1-norm. */
147
148 /* If NORMIN = 'N', CNORM is an output argument and CNORM(j) */
149 /* returns the 1-norm of the offdiagonal part of the j-th column */
150 /* of A. */
151
152 /* INFO (output) INTEGER */
153 /* = 0: successful exit */
154 /* < 0: if INFO = -k, the k-th argument had an illegal value */
155
156 /* Further Details */
157 /* ======= ======= */
158
159 /* A rough bound on x is computed; if that is less than overflow, ZTPSV */
160 /* is called, otherwise, specific code is used which checks for possible */
161 /* overflow or divide-by-zero at every operation. */
162
163 /* A columnwise scheme is used for solving A*x = b. The basic algorithm */
164 /* if A is lower triangular is */
165
166 /* x[1:n] := b[1:n] */
167 /* for j = 1, ..., n */
168 /* x(j) := x(j) / A(j,j) */
169 /* x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j] */
170 /* end */
171
172 /* Define bounds on the components of x after j iterations of the loop: */
173 /* M(j) = bound on x[1:j] */
174 /* G(j) = bound on x[j+1:n] */
175 /* Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}. */
176
177 /* Then for iteration j+1 we have */
178 /* M(j+1) <= G(j) / | A(j+1,j+1) | */
179 /* G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] | */
180 /* <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | ) */
181
182 /* where CNORM(j+1) is greater than or equal to the infinity-norm of */
183 /* column j+1 of A, not counting the diagonal. Hence */
184
185 /* G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | ) */
186 /* 1<=i<=j */
187 /* and */
188
189 /* |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| ) */
190 /* 1<=i< j */
191
192 /* Since |x(j)| <= M(j), we use the Level 2 BLAS routine ZTPSV if the */
193 /* reciprocal of the largest M(j), j=1,..,n, is larger than */
194 /* max(underflow, 1/overflow). */
195
196 /* The bound on x(j) is also used to determine when a step in the */
197 /* columnwise method can be performed without fear of overflow. If */
198 /* the computed bound is greater than a large constant, x is scaled to */
199 /* prevent overflow, but if the bound overflows, x is set to 0, x(j) to */
200 /* 1, and scale to 0, and a non-trivial solution to A*x = 0 is found. */
201
202 /* Similarly, a row-wise scheme is used to solve A**T *x = b or */
203 /* A**H *x = b. The basic algorithm for A upper triangular is */
204
205 /* for j = 1, ..., n */
206 /* x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j) */
207 /* end */
208
209 /* We simultaneously compute two bounds */
210 /* G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j */
211 /* M(j) = bound on x(i), 1<=i<=j */
212
213 /* The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we */
214 /* add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1. */
215 /* Then the bound on x(j) is */
216
217 /* M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) | */
218
219 /* <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| ) */
220 /* 1<=i<=j */
221
222 /* and we can safely call ZTPSV if 1/M(n) and 1/G(n) are both greater */
223 /* than max(underflow, 1/overflow). */
224
225 /* ===================================================================== */
226
227 /* .. Parameters .. */
228 /* .. */
229 /* .. Local Scalars .. */
230 /* .. */
231 /* .. External Functions .. */
232 /* .. */
233 /* .. External Subroutines .. */
234 /* .. */
235 /* .. Intrinsic Functions .. */
236 /* .. */
237 /* .. Statement Functions .. */
238 /* .. */
239 /* .. Statement Function definitions .. */
240 /* .. */
241 /* .. Executable Statements .. */
242
243 /* Parameter adjustments */
244 --cnorm;
245 --x;
246 --ap;
247
248 /* Function Body */
249 *info = 0;
250 upper = lsame_(uplo, "U", (ftnlen)1, (ftnlen)1);
251 notran = lsame_(trans, "N", (ftnlen)1, (ftnlen)1);
252 nounit = lsame_(diag, "N", (ftnlen)1, (ftnlen)1);
253
254 /* Test the input parameters. */
255
256 if (! upper && ! lsame_(uplo, "L", (ftnlen)1, (ftnlen)1)) {
257 *info = -1;
258 } else if (! notran && ! lsame_(trans, "T", (ftnlen)1, (ftnlen)1) && !
259 lsame_(trans, "C", (ftnlen)1, (ftnlen)1)) {
260 *info = -2;
261 } else if (! nounit && ! lsame_(diag, "U", (ftnlen)1, (ftnlen)1)) {
262 *info = -3;
263 } else if (! lsame_(normin, "Y", (ftnlen)1, (ftnlen)1) && ! lsame_(normin,
264 "N", (ftnlen)1, (ftnlen)1)) {
265 *info = -4;
266 } else if (*n < 0) {
267 *info = -5;
268 }
269 if (*info != 0) {
270 i__1 = -(*info);
271 xerbla_("ZLATPS", &i__1, (ftnlen)6);
272 return 0;
273 }
274
275 /* Quick return if possible */
276
277 if (*n == 0) {
278 return 0;
279 }
280
281 /* Determine machine dependent parameters to control overflow. */
282
283 smlnum = dlamch_("Safe minimum", (ftnlen)12);
284 bignum = 1. / smlnum;
285 dlabad_(&smlnum, &bignum);
286 smlnum /= dlamch_("Precision", (ftnlen)9);
287 bignum = 1. / smlnum;
288 *scale = 1.;
289
290 if (lsame_(normin, "N", (ftnlen)1, (ftnlen)1)) {
291
292 /* Compute the 1-norm of each column, not including the diagonal. */
293
294 if (upper) {
295
296 /* A is upper triangular. */
297
298 ip = 1;
299 i__1 = *n;
300 for (j = 1; j <= i__1; ++j) {
301 i__2 = j - 1;
302 cnorm[j] = dzasum_(&i__2, &ap[ip], &c__1);
303 ip += j;
304 /* L10: */
305 }
306 } else {
307
308 /* A is lower triangular. */
309
310 ip = 1;
311 i__1 = *n - 1;
312 for (j = 1; j <= i__1; ++j) {
313 i__2 = *n - j;
314 cnorm[j] = dzasum_(&i__2, &ap[ip + 1], &c__1);
315 ip = ip + *n - j + 1;
316 /* L20: */
317 }
318 cnorm[*n] = 0.;
319 }
320 }
321
322 /* Scale the column norms by TSCAL if the maximum element in CNORM is */
323 /* greater than BIGNUM/2. */
324
325 imax = idamax_(n, &cnorm[1], &c__1);
326 tmax = cnorm[imax];
327 if (tmax <= bignum * .5) {
328 tscal = 1.;
329 } else {
330 tscal = .5 / (smlnum * tmax);
331 dscal_(n, &tscal, &cnorm[1], &c__1);
332 }
333
334 /* Compute a bound on the computed solution vector to see if the */
335 /* Level 2 BLAS routine ZTPSV can be used. */
336
337 xmax = 0.;
338 i__1 = *n;
339 for (j = 1; j <= i__1; ++j) {
340 /* Computing MAX */
341 i__2 = j;
342 d__3 = xmax, d__4 = (d__1 = x[i__2].r / 2., abs(d__1)) + (d__2 =
343 d_imag(&x[j]) / 2., abs(d__2));
344 xmax = max(d__3,d__4);
345 /* L30: */
346 }
347 xbnd = xmax;
348 if (notran) {
349
350 /* Compute the growth in A * x = b. */
351
352 if (upper) {
353 jfirst = *n;
354 jlast = 1;
355 jinc = -1;
356 } else {
357 jfirst = 1;
358 jlast = *n;
359 jinc = 1;
360 }
361
362 if (tscal != 1.) {
363 grow = 0.;
364 goto L60;
365 }
366
367 if (nounit) {
368
369 /* A is non-unit triangular. */
370
371 /* Compute GROW = 1/G(j) and XBND = 1/M(j). */
372 /* Initially, G(0) = max{x(i), i=1,...,n}. */
373
374 grow = .5 / max(xbnd,smlnum);
375 xbnd = grow;
376 ip = jfirst * (jfirst + 1) / 2;
377 jlen = *n;
378 i__1 = jlast;
379 i__2 = jinc;
380 for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
381
382 /* Exit the loop if the growth factor is too small. */
383
384 if (grow <= smlnum) {
385 goto L60;
386 }
387
388 i__3 = ip;
389 tjjs.r = ap[i__3].r, tjjs.i = ap[i__3].i;
390 tjj = (d__1 = tjjs.r, abs(d__1)) + (d__2 = d_imag(&tjjs), abs(
391 d__2));
392
393 if (tjj >= smlnum) {
394
395 /* M(j) = G(j-1) / abs(A(j,j)) */
396
397 /* Computing MIN */
398 d__1 = xbnd, d__2 = min(1.,tjj) * grow;
399 xbnd = min(d__1,d__2);
400 } else {
401
402 /* M(j) could overflow, set XBND to 0. */
403
404 xbnd = 0.;
405 }
406
407 if (tjj + cnorm[j] >= smlnum) {
408
409 /* G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) ) */
410
411 grow *= tjj / (tjj + cnorm[j]);
412 } else {
413
414 /* G(j) could overflow, set GROW to 0. */
415
416 grow = 0.;
417 }
418 ip += jinc * jlen;
419 --jlen;
420 /* L40: */
421 }
422 grow = xbnd;
423 } else {
424
425 /* A is unit triangular. */
426
427 /* Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}. */
428
429 /* Computing MIN */
430 d__1 = 1., d__2 = .5 / max(xbnd,smlnum);
431 grow = min(d__1,d__2);
432 i__2 = jlast;
433 i__1 = jinc;
434 for (j = jfirst; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
435
436 /* Exit the loop if the growth factor is too small. */
437
438 if (grow <= smlnum) {
439 goto L60;
440 }
441
442 /* G(j) = G(j-1)*( 1 + CNORM(j) ) */
443
444 grow *= 1. / (cnorm[j] + 1.);
445 /* L50: */
446 }
447 }
448 L60:
449
450 ;
451 } else {
452
453 /* Compute the growth in A**T * x = b or A**H * x = b. */
454
455 if (upper) {
456 jfirst = 1;
457 jlast = *n;
458 jinc = 1;
459 } else {
460 jfirst = *n;
461 jlast = 1;
462 jinc = -1;
463 }
464
465 if (tscal != 1.) {
466 grow = 0.;
467 goto L90;
468 }
469
470 if (nounit) {
471
472 /* A is non-unit triangular. */
473
474 /* Compute GROW = 1/G(j) and XBND = 1/M(j). */
475 /* Initially, M(0) = max{x(i), i=1,...,n}. */
476
477 grow = .5 / max(xbnd,smlnum);
478 xbnd = grow;
479 ip = jfirst * (jfirst + 1) / 2;
480 jlen = 1;
481 i__1 = jlast;
482 i__2 = jinc;
483 for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
484
485 /* Exit the loop if the growth factor is too small. */
486
487 if (grow <= smlnum) {
488 goto L90;
489 }
490
491 /* G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) ) */
492
493 xj = cnorm[j] + 1.;
494 /* Computing MIN */
495 d__1 = grow, d__2 = xbnd / xj;
496 grow = min(d__1,d__2);
497
498 i__3 = ip;
499 tjjs.r = ap[i__3].r, tjjs.i = ap[i__3].i;
500 tjj = (d__1 = tjjs.r, abs(d__1)) + (d__2 = d_imag(&tjjs), abs(
501 d__2));
502
503 if (tjj >= smlnum) {
504
505 /* M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j)) */
506
507 if (xj > tjj) {
508 xbnd *= tjj / xj;
509 }
510 } else {
511
512 /* M(j) could overflow, set XBND to 0. */
513
514 xbnd = 0.;
515 }
516 ++jlen;
517 ip += jinc * jlen;
518 /* L70: */
519 }
520 grow = min(grow,xbnd);
521 } else {
522
523 /* A is unit triangular. */
524
525 /* Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}. */
526
527 /* Computing MIN */
528 d__1 = 1., d__2 = .5 / max(xbnd,smlnum);
529 grow = min(d__1,d__2);
530 i__2 = jlast;
531 i__1 = jinc;
532 for (j = jfirst; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
533
534 /* Exit the loop if the growth factor is too small. */
535
536 if (grow <= smlnum) {
537 goto L90;
538 }
539
540 /* G(j) = ( 1 + CNORM(j) )*G(j-1) */
541
542 xj = cnorm[j] + 1.;
543 grow /= xj;
544 /* L80: */
545 }
546 }
547 L90:
548 ;
549 }
550
551 if (grow * tscal > smlnum) {
552
553 /* Use the Level 2 BLAS solve if the reciprocal of the bound on */
554 /* elements of X is not too small. */
555
556 ztpsv_(uplo, trans, diag, n, &ap[1], &x[1], &c__1, (ftnlen)1, (ftnlen)
557 1, (ftnlen)1);
558 } else {
559
560 /* Use a Level 1 BLAS solve, scaling intermediate results. */
561
562 if (xmax > bignum * .5) {
563
564 /* Scale X so that its components are less than or equal to */
565 /* BIGNUM in absolute value. */
566
567 *scale = bignum * .5 / xmax;
568 zdscal_(n, scale, &x[1], &c__1);
569 xmax = bignum;
570 } else {
571 xmax *= 2.;
572 }
573
574 if (notran) {
575
576 /* Solve A * x = b */
577
578 ip = jfirst * (jfirst + 1) / 2;
579 i__1 = jlast;
580 i__2 = jinc;
581 for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
582
583 /* Compute x(j) = b(j) / A(j,j), scaling x if necessary. */
584
585 i__3 = j;
586 xj = (d__1 = x[i__3].r, abs(d__1)) + (d__2 = d_imag(&x[j]),
587 abs(d__2));
588 if (nounit) {
589 i__3 = ip;
590 z__1.r = tscal * ap[i__3].r, z__1.i = tscal * ap[i__3].i;
591 tjjs.r = z__1.r, tjjs.i = z__1.i;
592 } else {
593 tjjs.r = tscal, tjjs.i = 0.;
594 if (tscal == 1.) {
595 goto L110;
596 }
597 }
598 tjj = (d__1 = tjjs.r, abs(d__1)) + (d__2 = d_imag(&tjjs), abs(
599 d__2));
600 if (tjj > smlnum) {
601
602 /* abs(A(j,j)) > SMLNUM: */
603
604 if (tjj < 1.) {
605 if (xj > tjj * bignum) {
606
607 /* Scale x by 1/b(j). */
608
609 rec = 1. / xj;
610 zdscal_(n, &rec, &x[1], &c__1);
611 *scale *= rec;
612 xmax *= rec;
613 }
614 }
615 i__3 = j;
616 zladiv_(&z__1, &x[j], &tjjs);
617 x[i__3].r = z__1.r, x[i__3].i = z__1.i;
618 i__3 = j;
619 xj = (d__1 = x[i__3].r, abs(d__1)) + (d__2 = d_imag(&x[j])
620 , abs(d__2));
621 } else if (tjj > 0.) {
622
623 /* 0 < abs(A(j,j)) <= SMLNUM: */
624
625 if (xj > tjj * bignum) {
626
627 /* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM */
628 /* to avoid overflow when dividing by A(j,j). */
629
630 rec = tjj * bignum / xj;
631 if (cnorm[j] > 1.) {
632
633 /* Scale by 1/CNORM(j) to avoid overflow when */
634 /* multiplying x(j) times column j. */
635
636 rec /= cnorm[j];
637 }
638 zdscal_(n, &rec, &x[1], &c__1);
639 *scale *= rec;
640 xmax *= rec;
641 }
642 i__3 = j;
643 zladiv_(&z__1, &x[j], &tjjs);
644 x[i__3].r = z__1.r, x[i__3].i = z__1.i;
645 i__3 = j;
646 xj = (d__1 = x[i__3].r, abs(d__1)) + (d__2 = d_imag(&x[j])
647 , abs(d__2));
648 } else {
649
650 /* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and */
651 /* scale = 0, and compute a solution to A*x = 0. */
652
653 i__3 = *n;
654 for (i__ = 1; i__ <= i__3; ++i__) {
655 i__4 = i__;
656 x[i__4].r = 0., x[i__4].i = 0.;
657 /* L100: */
658 }
659 i__3 = j;
660 x[i__3].r = 1., x[i__3].i = 0.;
661 xj = 1.;
662 *scale = 0.;
663 xmax = 0.;
664 }
665 L110:
666
667 /* Scale x if necessary to avoid overflow when adding a */
668 /* multiple of column j of A. */
669
670 if (xj > 1.) {
671 rec = 1. / xj;
672 if (cnorm[j] > (bignum - xmax) * rec) {
673
674 /* Scale x by 1/(2*abs(x(j))). */
675
676 rec *= .5;
677 zdscal_(n, &rec, &x[1], &c__1);
678 *scale *= rec;
679 }
680 } else if (xj * cnorm[j] > bignum - xmax) {
681
682 /* Scale x by 1/2. */
683
684 zdscal_(n, &c_b36, &x[1], &c__1);
685 *scale *= .5;
686 }
687
688 if (upper) {
689 if (j > 1) {
690
691 /* Compute the update */
692 /* x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j) */
693
694 i__3 = j - 1;
695 i__4 = j;
696 z__2.r = -x[i__4].r, z__2.i = -x[i__4].i;
697 z__1.r = tscal * z__2.r, z__1.i = tscal * z__2.i;
698 zaxpy_(&i__3, &z__1, &ap[ip - j + 1], &c__1, &x[1], &
699 c__1);
700 i__3 = j - 1;
701 i__ = izamax_(&i__3, &x[1], &c__1);
702 i__3 = i__;
703 xmax = (d__1 = x[i__3].r, abs(d__1)) + (d__2 = d_imag(
704 &x[i__]), abs(d__2));
705 }
706 ip -= j;
707 } else {
708 if (j < *n) {
709
710 /* Compute the update */
711 /* x(j+1:n) := x(j+1:n) - x(j) * A(j+1:n,j) */
712
713 i__3 = *n - j;
714 i__4 = j;
715 z__2.r = -x[i__4].r, z__2.i = -x[i__4].i;
716 z__1.r = tscal * z__2.r, z__1.i = tscal * z__2.i;
717 zaxpy_(&i__3, &z__1, &ap[ip + 1], &c__1, &x[j + 1], &
718 c__1);
719 i__3 = *n - j;
720 i__ = j + izamax_(&i__3, &x[j + 1], &c__1);
721 i__3 = i__;
722 xmax = (d__1 = x[i__3].r, abs(d__1)) + (d__2 = d_imag(
723 &x[i__]), abs(d__2));
724 }
725 ip = ip + *n - j + 1;
726 }
727 /* L120: */
728 }
729
730 } else if (lsame_(trans, "T", (ftnlen)1, (ftnlen)1)) {
731
732 /* Solve A**T * x = b */
733
734 ip = jfirst * (jfirst + 1) / 2;
735 jlen = 1;
736 i__2 = jlast;
737 i__1 = jinc;
738 for (j = jfirst; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
739
740 /* Compute x(j) = b(j) - sum A(k,j)*x(k). */
741 /* k<>j */
742
743 i__3 = j;
744 xj = (d__1 = x[i__3].r, abs(d__1)) + (d__2 = d_imag(&x[j]),
745 abs(d__2));
746 uscal.r = tscal, uscal.i = 0.;
747 rec = 1. / max(xmax,1.);
748 if (cnorm[j] > (bignum - xj) * rec) {
749
750 /* If x(j) could overflow, scale x by 1/(2*XMAX). */
751
752 rec *= .5;
753 if (nounit) {
754 i__3 = ip;
755 z__1.r = tscal * ap[i__3].r, z__1.i = tscal * ap[i__3]
756 .i;
757 tjjs.r = z__1.r, tjjs.i = z__1.i;
758 } else {
759 tjjs.r = tscal, tjjs.i = 0.;
760 }
761 tjj = (d__1 = tjjs.r, abs(d__1)) + (d__2 = d_imag(&tjjs),
762 abs(d__2));
763 if (tjj > 1.) {
764
765 /* Divide by A(j,j) when scaling x if A(j,j) > 1. */
766
767 /* Computing MIN */
768 d__1 = 1., d__2 = rec * tjj;
769 rec = min(d__1,d__2);
770 zladiv_(&z__1, &uscal, &tjjs);
771 uscal.r = z__1.r, uscal.i = z__1.i;
772 }
773 if (rec < 1.) {
774 zdscal_(n, &rec, &x[1], &c__1);
775 *scale *= rec;
776 xmax *= rec;
777 }
778 }
779
780 csumj.r = 0., csumj.i = 0.;
781 if (uscal.r == 1. && uscal.i == 0.) {
782
783 /* If the scaling needed for A in the dot product is 1, */
784 /* call ZDOTU to perform the dot product. */
785
786 if (upper) {
787 i__3 = j - 1;
788 zdotu_(&z__1, &i__3, &ap[ip - j + 1], &c__1, &x[1], &
789 c__1);
790 csumj.r = z__1.r, csumj.i = z__1.i;
791 } else if (j < *n) {
792 i__3 = *n - j;
793 zdotu_(&z__1, &i__3, &ap[ip + 1], &c__1, &x[j + 1], &
794 c__1);
795 csumj.r = z__1.r, csumj.i = z__1.i;
796 }
797 } else {
798
799 /* Otherwise, use in-line code for the dot product. */
800
801 if (upper) {
802 i__3 = j - 1;
803 for (i__ = 1; i__ <= i__3; ++i__) {
804 i__4 = ip - j + i__;
805 z__3.r = ap[i__4].r * uscal.r - ap[i__4].i *
806 uscal.i, z__3.i = ap[i__4].r * uscal.i +
807 ap[i__4].i * uscal.r;
808 i__5 = i__;
809 z__2.r = z__3.r * x[i__5].r - z__3.i * x[i__5].i,
810 z__2.i = z__3.r * x[i__5].i + z__3.i * x[
811 i__5].r;
812 z__1.r = csumj.r + z__2.r, z__1.i = csumj.i +
813 z__2.i;
814 csumj.r = z__1.r, csumj.i = z__1.i;
815 /* L130: */
816 }
817 } else if (j < *n) {
818 i__3 = *n - j;
819 for (i__ = 1; i__ <= i__3; ++i__) {
820 i__4 = ip + i__;
821 z__3.r = ap[i__4].r * uscal.r - ap[i__4].i *
822 uscal.i, z__3.i = ap[i__4].r * uscal.i +
823 ap[i__4].i * uscal.r;
824 i__5 = j + i__;
825 z__2.r = z__3.r * x[i__5].r - z__3.i * x[i__5].i,
826 z__2.i = z__3.r * x[i__5].i + z__3.i * x[
827 i__5].r;
828 z__1.r = csumj.r + z__2.r, z__1.i = csumj.i +
829 z__2.i;
830 csumj.r = z__1.r, csumj.i = z__1.i;
831 /* L140: */
832 }
833 }
834 }
835
836 z__1.r = tscal, z__1.i = 0.;
837 if (uscal.r == z__1.r && uscal.i == z__1.i) {
838
839 /* Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j) */
840 /* was not used to scale the dotproduct. */
841
842 i__3 = j;
843 i__4 = j;
844 z__1.r = x[i__4].r - csumj.r, z__1.i = x[i__4].i -
845 csumj.i;
846 x[i__3].r = z__1.r, x[i__3].i = z__1.i;
847 i__3 = j;
848 xj = (d__1 = x[i__3].r, abs(d__1)) + (d__2 = d_imag(&x[j])
849 , abs(d__2));
850 if (nounit) {
851
852 /* Compute x(j) = x(j) / A(j,j), scaling if necessary. */
853
854 i__3 = ip;
855 z__1.r = tscal * ap[i__3].r, z__1.i = tscal * ap[i__3]
856 .i;
857 tjjs.r = z__1.r, tjjs.i = z__1.i;
858 } else {
859 tjjs.r = tscal, tjjs.i = 0.;
860 if (tscal == 1.) {
861 goto L160;
862 }
863 }
864 tjj = (d__1 = tjjs.r, abs(d__1)) + (d__2 = d_imag(&tjjs),
865 abs(d__2));
866 if (tjj > smlnum) {
867
868 /* abs(A(j,j)) > SMLNUM: */
869
870 if (tjj < 1.) {
871 if (xj > tjj * bignum) {
872
873 /* Scale X by 1/abs(x(j)). */
874
875 rec = 1. / xj;
876 zdscal_(n, &rec, &x[1], &c__1);
877 *scale *= rec;
878 xmax *= rec;
879 }
880 }
881 i__3 = j;
882 zladiv_(&z__1, &x[j], &tjjs);
883 x[i__3].r = z__1.r, x[i__3].i = z__1.i;
884 } else if (tjj > 0.) {
885
886 /* 0 < abs(A(j,j)) <= SMLNUM: */
887
888 if (xj > tjj * bignum) {
889
890 /* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM. */
891
892 rec = tjj * bignum / xj;
893 zdscal_(n, &rec, &x[1], &c__1);
894 *scale *= rec;
895 xmax *= rec;
896 }
897 i__3 = j;
898 zladiv_(&z__1, &x[j], &tjjs);
899 x[i__3].r = z__1.r, x[i__3].i = z__1.i;
900 } else {
901
902 /* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and */
903 /* scale = 0 and compute a solution to A**T *x = 0. */
904
905 i__3 = *n;
906 for (i__ = 1; i__ <= i__3; ++i__) {
907 i__4 = i__;
908 x[i__4].r = 0., x[i__4].i = 0.;
909 /* L150: */
910 }
911 i__3 = j;
912 x[i__3].r = 1., x[i__3].i = 0.;
913 *scale = 0.;
914 xmax = 0.;
915 }
916 L160:
917 ;
918 } else {
919
920 /* Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot */
921 /* product has already been divided by 1/A(j,j). */
922
923 i__3 = j;
924 zladiv_(&z__2, &x[j], &tjjs);
925 z__1.r = z__2.r - csumj.r, z__1.i = z__2.i - csumj.i;
926 x[i__3].r = z__1.r, x[i__3].i = z__1.i;
927 }
928 /* Computing MAX */
929 i__3 = j;
930 d__3 = xmax, d__4 = (d__1 = x[i__3].r, abs(d__1)) + (d__2 =
931 d_imag(&x[j]), abs(d__2));
932 xmax = max(d__3,d__4);
933 ++jlen;
934 ip += jinc * jlen;
935 /* L170: */
936 }
937
938 } else {
939
940 /* Solve A**H * x = b */
941
942 ip = jfirst * (jfirst + 1) / 2;
943 jlen = 1;
944 i__1 = jlast;
945 i__2 = jinc;
946 for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
947
948 /* Compute x(j) = b(j) - sum A(k,j)*x(k). */
949 /* k<>j */
950
951 i__3 = j;
952 xj = (d__1 = x[i__3].r, abs(d__1)) + (d__2 = d_imag(&x[j]),
953 abs(d__2));
954 uscal.r = tscal, uscal.i = 0.;
955 rec = 1. / max(xmax,1.);
956 if (cnorm[j] > (bignum - xj) * rec) {
957
958 /* If x(j) could overflow, scale x by 1/(2*XMAX). */
959
960 rec *= .5;
961 if (nounit) {
962 d_cnjg(&z__2, &ap[ip]);
963 z__1.r = tscal * z__2.r, z__1.i = tscal * z__2.i;
964 tjjs.r = z__1.r, tjjs.i = z__1.i;
965 } else {
966 tjjs.r = tscal, tjjs.i = 0.;
967 }
968 tjj = (d__1 = tjjs.r, abs(d__1)) + (d__2 = d_imag(&tjjs),
969 abs(d__2));
970 if (tjj > 1.) {
971
972 /* Divide by A(j,j) when scaling x if A(j,j) > 1. */
973
974 /* Computing MIN */
975 d__1 = 1., d__2 = rec * tjj;
976 rec = min(d__1,d__2);
977 zladiv_(&z__1, &uscal, &tjjs);
978 uscal.r = z__1.r, uscal.i = z__1.i;
979 }
980 if (rec < 1.) {
981 zdscal_(n, &rec, &x[1], &c__1);
982 *scale *= rec;
983 xmax *= rec;
984 }
985 }
986
987 csumj.r = 0., csumj.i = 0.;
988 if (uscal.r == 1. && uscal.i == 0.) {
989
990 /* If the scaling needed for A in the dot product is 1, */
991 /* call ZDOTC to perform the dot product. */
992
993 if (upper) {
994 i__3 = j - 1;
995 zdotc_(&z__1, &i__3, &ap[ip - j + 1], &c__1, &x[1], &
996 c__1);
997 csumj.r = z__1.r, csumj.i = z__1.i;
998 } else if (j < *n) {
999 i__3 = *n - j;
1000 zdotc_(&z__1, &i__3, &ap[ip + 1], &c__1, &x[j + 1], &
1001 c__1);
1002 csumj.r = z__1.r, csumj.i = z__1.i;
1003 }
1004 } else {
1005
1006 /* Otherwise, use in-line code for the dot product. */
1007
1008 if (upper) {
1009 i__3 = j - 1;
1010 for (i__ = 1; i__ <= i__3; ++i__) {
1011 d_cnjg(&z__4, &ap[ip - j + i__]);
1012 z__3.r = z__4.r * uscal.r - z__4.i * uscal.i,
1013 z__3.i = z__4.r * uscal.i + z__4.i *
1014 uscal.r;
1015 i__4 = i__;
1016 z__2.r = z__3.r * x[i__4].r - z__3.i * x[i__4].i,
1017 z__2.i = z__3.r * x[i__4].i + z__3.i * x[
1018 i__4].r;
1019 z__1.r = csumj.r + z__2.r, z__1.i = csumj.i +
1020 z__2.i;
1021 csumj.r = z__1.r, csumj.i = z__1.i;
1022 /* L180: */
1023 }
1024 } else if (j < *n) {
1025 i__3 = *n - j;
1026 for (i__ = 1; i__ <= i__3; ++i__) {
1027 d_cnjg(&z__4, &ap[ip + i__]);
1028 z__3.r = z__4.r * uscal.r - z__4.i * uscal.i,
1029 z__3.i = z__4.r * uscal.i + z__4.i *
1030 uscal.r;
1031 i__4 = j + i__;
1032 z__2.r = z__3.r * x[i__4].r - z__3.i * x[i__4].i,
1033 z__2.i = z__3.r * x[i__4].i + z__3.i * x[
1034 i__4].r;
1035 z__1.r = csumj.r + z__2.r, z__1.i = csumj.i +
1036 z__2.i;
1037 csumj.r = z__1.r, csumj.i = z__1.i;
1038 /* L190: */
1039 }
1040 }
1041 }
1042
1043 z__1.r = tscal, z__1.i = 0.;
1044 if (uscal.r == z__1.r && uscal.i == z__1.i) {
1045
1046 /* Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j) */
1047 /* was not used to scale the dotproduct. */
1048
1049 i__3 = j;
1050 i__4 = j;
1051 z__1.r = x[i__4].r - csumj.r, z__1.i = x[i__4].i -
1052 csumj.i;
1053 x[i__3].r = z__1.r, x[i__3].i = z__1.i;
1054 i__3 = j;
1055 xj = (d__1 = x[i__3].r, abs(d__1)) + (d__2 = d_imag(&x[j])
1056 , abs(d__2));
1057 if (nounit) {
1058
1059 /* Compute x(j) = x(j) / A(j,j), scaling if necessary. */
1060
1061 d_cnjg(&z__2, &ap[ip]);
1062 z__1.r = tscal * z__2.r, z__1.i = tscal * z__2.i;
1063 tjjs.r = z__1.r, tjjs.i = z__1.i;
1064 } else {
1065 tjjs.r = tscal, tjjs.i = 0.;
1066 if (tscal == 1.) {
1067 goto L210;
1068 }
1069 }
1070 tjj = (d__1 = tjjs.r, abs(d__1)) + (d__2 = d_imag(&tjjs),
1071 abs(d__2));
1072 if (tjj > smlnum) {
1073
1074 /* abs(A(j,j)) > SMLNUM: */
1075
1076 if (tjj < 1.) {
1077 if (xj > tjj * bignum) {
1078
1079 /* Scale X by 1/abs(x(j)). */
1080
1081 rec = 1. / xj;
1082 zdscal_(n, &rec, &x[1], &c__1);
1083 *scale *= rec;
1084 xmax *= rec;
1085 }
1086 }
1087 i__3 = j;
1088 zladiv_(&z__1, &x[j], &tjjs);
1089 x[i__3].r = z__1.r, x[i__3].i = z__1.i;
1090 } else if (tjj > 0.) {
1091
1092 /* 0 < abs(A(j,j)) <= SMLNUM: */
1093
1094 if (xj > tjj * bignum) {
1095
1096 /* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM. */
1097
1098 rec = tjj * bignum / xj;
1099 zdscal_(n, &rec, &x[1], &c__1);
1100 *scale *= rec;
1101 xmax *= rec;
1102 }
1103 i__3 = j;
1104 zladiv_(&z__1, &x[j], &tjjs);
1105 x[i__3].r = z__1.r, x[i__3].i = z__1.i;
1106 } else {
1107
1108 /* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and */
1109 /* scale = 0 and compute a solution to A**H *x = 0. */
1110
1111 i__3 = *n;
1112 for (i__ = 1; i__ <= i__3; ++i__) {
1113 i__4 = i__;
1114 x[i__4].r = 0., x[i__4].i = 0.;
1115 /* L200: */
1116 }
1117 i__3 = j;
1118 x[i__3].r = 1., x[i__3].i = 0.;
1119 *scale = 0.;
1120 xmax = 0.;
1121 }
1122 L210:
1123 ;
1124 } else {
1125
1126 /* Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot */
1127 /* product has already been divided by 1/A(j,j). */
1128
1129 i__3 = j;
1130 zladiv_(&z__2, &x[j], &tjjs);
1131 z__1.r = z__2.r - csumj.r, z__1.i = z__2.i - csumj.i;
1132 x[i__3].r = z__1.r, x[i__3].i = z__1.i;
1133 }
1134 /* Computing MAX */
1135 i__3 = j;
1136 d__3 = xmax, d__4 = (d__1 = x[i__3].r, abs(d__1)) + (d__2 =
1137 d_imag(&x[j]), abs(d__2));
1138 xmax = max(d__3,d__4);
1139 ++jlen;
1140 ip += jinc * jlen;
1141 /* L220: */
1142 }
1143 }
1144 *scale /= tscal;
1145 }
1146
1147 /* Scale the column norms by 1/TSCAL for return. */
1148
1149 if (tscal != 1.) {
1150 d__1 = 1. / tscal;
1151 dscal_(n, &d__1, &cnorm[1], &c__1);
1152 }
1153
1154 return 0;
1155
1156 /* End of ZLATPS */
1157
1158 } /* zlatps_ */
1159
1160