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