001:       SUBROUTINE CLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X,
002:      $                   SCALE, CNORM, INFO )
003: *
004: *  -- LAPACK auxiliary routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       CHARACTER          DIAG, NORMIN, TRANS, UPLO
010:       INTEGER            INFO, KD, LDAB, N
011:       REAL               SCALE
012: *     ..
013: *     .. Array Arguments ..
014:       REAL               CNORM( * )
015:       COMPLEX            AB( LDAB, * ), X( * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  CLATBS solves one of the triangular systems
022: *
023: *     A * x = s*b,  A**T * x = s*b,  or  A**H * x = s*b,
024: *
025: *  with scaling to prevent overflow, where A is an upper or lower
026: *  triangular band matrix.  Here A' denotes the transpose of A, x and b
027: *  are n-element vectors, and s is a scaling factor, usually less than
028: *  or equal to 1, chosen so that the components of x will be less than
029: *  the overflow threshold.  If the unscaled problem will not cause
030: *  overflow, the Level 2 BLAS routine CTBSV is called.  If the matrix A
031: *  is singular (A(j,j) = 0 for some j), then s is set to 0 and a
032: *  non-trivial solution to A*x = 0 is returned.
033: *
034: *  Arguments
035: *  =========
036: *
037: *  UPLO    (input) CHARACTER*1
038: *          Specifies whether the matrix A is upper or lower triangular.
039: *          = 'U':  Upper triangular
040: *          = 'L':  Lower triangular
041: *
042: *  TRANS   (input) CHARACTER*1
043: *          Specifies the operation applied to A.
044: *          = 'N':  Solve A * x = s*b     (No transpose)
045: *          = 'T':  Solve A**T * x = s*b  (Transpose)
046: *          = 'C':  Solve A**H * x = s*b  (Conjugate transpose)
047: *
048: *  DIAG    (input) CHARACTER*1
049: *          Specifies whether or not the matrix A is unit triangular.
050: *          = 'N':  Non-unit triangular
051: *          = 'U':  Unit triangular
052: *
053: *  NORMIN  (input) CHARACTER*1
054: *          Specifies whether CNORM has been set or not.
055: *          = 'Y':  CNORM contains the column norms on entry
056: *          = 'N':  CNORM is not set on entry.  On exit, the norms will
057: *                  be computed and stored in CNORM.
058: *
059: *  N       (input) INTEGER
060: *          The order of the matrix A.  N >= 0.
061: *
062: *  KD      (input) INTEGER
063: *          The number of subdiagonals or superdiagonals in the
064: *          triangular matrix A.  KD >= 0.
065: *
066: *  AB      (input) COMPLEX array, dimension (LDAB,N)
067: *          The upper or lower triangular band matrix A, stored in the
068: *          first KD+1 rows of the array. The j-th column of A is stored
069: *          in the j-th column of the array AB as follows:
070: *          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
071: *          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
072: *
073: *  LDAB    (input) INTEGER
074: *          The leading dimension of the array AB.  LDAB >= KD+1.
075: *
076: *  X       (input/output) COMPLEX array, dimension (N)
077: *          On entry, the right hand side b of the triangular system.
078: *          On exit, X is overwritten by the solution vector x.
079: *
080: *  SCALE   (output) REAL
081: *          The scaling factor s for the triangular system
082: *             A * x = s*b,  A**T * x = s*b,  or  A**H * x = s*b.
083: *          If SCALE = 0, the matrix A is singular or badly scaled, and
084: *          the vector x is an exact or approximate solution to A*x = 0.
085: *
086: *  CNORM   (input or output) REAL array, dimension (N)
087: *
088: *          If NORMIN = 'Y', CNORM is an input argument and CNORM(j)
089: *          contains the norm of the off-diagonal part of the j-th column
090: *          of A.  If TRANS = 'N', CNORM(j) must be greater than or equal
091: *          to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)
092: *          must be greater than or equal to the 1-norm.
093: *
094: *          If NORMIN = 'N', CNORM is an output argument and CNORM(j)
095: *          returns the 1-norm of the offdiagonal part of the j-th column
096: *          of A.
097: *
098: *  INFO    (output) INTEGER
099: *          = 0:  successful exit
100: *          < 0:  if INFO = -k, the k-th argument had an illegal value
101: *
102: *  Further Details
103: *  ======= =======
104: *
105: *  A rough bound on x is computed; if that is less than overflow, CTBSV
106: *  is called, otherwise, specific code is used which checks for possible
107: *  overflow or divide-by-zero at every operation.
108: *
109: *  A columnwise scheme is used for solving A*x = b.  The basic algorithm
110: *  if A is lower triangular is
111: *
112: *       x[1:n] := b[1:n]
113: *       for j = 1, ..., n
114: *            x(j) := x(j) / A(j,j)
115: *            x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
116: *       end
117: *
118: *  Define bounds on the components of x after j iterations of the loop:
119: *     M(j) = bound on x[1:j]
120: *     G(j) = bound on x[j+1:n]
121: *  Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
122: *
123: *  Then for iteration j+1 we have
124: *     M(j+1) <= G(j) / | A(j+1,j+1) |
125: *     G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
126: *            <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
127: *
128: *  where CNORM(j+1) is greater than or equal to the infinity-norm of
129: *  column j+1 of A, not counting the diagonal.  Hence
130: *
131: *     G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
132: *                  1<=i<=j
133: *  and
134: *
135: *     |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
136: *                                   1<=i< j
137: *
138: *  Since |x(j)| <= M(j), we use the Level 2 BLAS routine CTBSV if the
139: *  reciprocal of the largest M(j), j=1,..,n, is larger than
140: *  max(underflow, 1/overflow).
141: *
142: *  The bound on x(j) is also used to determine when a step in the
143: *  columnwise method can be performed without fear of overflow.  If
144: *  the computed bound is greater than a large constant, x is scaled to
145: *  prevent overflow, but if the bound overflows, x is set to 0, x(j) to
146: *  1, and scale to 0, and a non-trivial solution to A*x = 0 is found.
147: *
148: *  Similarly, a row-wise scheme is used to solve A**T *x = b  or
149: *  A**H *x = b.  The basic algorithm for A upper triangular is
150: *
151: *       for j = 1, ..., n
152: *            x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)
153: *       end
154: *
155: *  We simultaneously compute two bounds
156: *       G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j
157: *       M(j) = bound on x(i), 1<=i<=j
158: *
159: *  The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we
160: *  add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.
161: *  Then the bound on x(j) is
162: *
163: *       M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
164: *
165: *            <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
166: *                      1<=i<=j
167: *
168: *  and we can safely call CTBSV if 1/M(n) and 1/G(n) are both greater
169: *  than max(underflow, 1/overflow).
170: *
171: *  =====================================================================
172: *
173: *     .. Parameters ..
174:       REAL               ZERO, HALF, ONE, TWO
175:       PARAMETER          ( ZERO = 0.0E+0, HALF = 0.5E+0, ONE = 1.0E+0,
176:      $                   TWO = 2.0E+0 )
177: *     ..
178: *     .. Local Scalars ..
179:       LOGICAL            NOTRAN, NOUNIT, UPPER
180:       INTEGER            I, IMAX, J, JFIRST, JINC, JLAST, JLEN, MAIND
181:       REAL               BIGNUM, GROW, REC, SMLNUM, TJJ, TMAX, TSCAL,
182:      $                   XBND, XJ, XMAX
183:       COMPLEX            CSUMJ, TJJS, USCAL, ZDUM
184: *     ..
185: *     .. External Functions ..
186:       LOGICAL            LSAME
187:       INTEGER            ICAMAX, ISAMAX
188:       REAL               SCASUM, SLAMCH
189:       COMPLEX            CDOTC, CDOTU, CLADIV
190:       EXTERNAL           LSAME, ICAMAX, ISAMAX, SCASUM, SLAMCH, CDOTC,
191:      $                   CDOTU, CLADIV
192: *     ..
193: *     .. External Subroutines ..
194:       EXTERNAL           CAXPY, CSSCAL, CTBSV, SLABAD, SSCAL, XERBLA
195: *     ..
196: *     .. Intrinsic Functions ..
197:       INTRINSIC          ABS, AIMAG, CMPLX, CONJG, MAX, MIN, REAL
198: *     ..
199: *     .. Statement Functions ..
200:       REAL               CABS1, CABS2
201: *     ..
202: *     .. Statement Function definitions ..
203:       CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
204:       CABS2( ZDUM ) = ABS( REAL( ZDUM ) / 2. ) +
205:      $                ABS( AIMAG( ZDUM ) / 2. )
206: *     ..
207: *     .. Executable Statements ..
208: *
209:       INFO = 0
210:       UPPER = LSAME( UPLO, 'U' )
211:       NOTRAN = LSAME( TRANS, 'N' )
212:       NOUNIT = LSAME( DIAG, 'N' )
213: *
214: *     Test the input parameters.
215: *
216:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
217:          INFO = -1
218:       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
219:      $         LSAME( TRANS, 'C' ) ) THEN
220:          INFO = -2
221:       ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
222:          INFO = -3
223:       ELSE IF( .NOT.LSAME( NORMIN, 'Y' ) .AND. .NOT.
224:      $         LSAME( NORMIN, 'N' ) ) THEN
225:          INFO = -4
226:       ELSE IF( N.LT.0 ) THEN
227:          INFO = -5
228:       ELSE IF( KD.LT.0 ) THEN
229:          INFO = -6
230:       ELSE IF( LDAB.LT.KD+1 ) THEN
231:          INFO = -8
232:       END IF
233:       IF( INFO.NE.0 ) THEN
234:          CALL XERBLA( 'CLATBS', -INFO )
235:          RETURN
236:       END IF
237: *
238: *     Quick return if possible
239: *
240:       IF( N.EQ.0 )
241:      $   RETURN
242: *
243: *     Determine machine dependent parameters to control overflow.
244: *
245:       SMLNUM = SLAMCH( 'Safe minimum' )
246:       BIGNUM = ONE / SMLNUM
247:       CALL SLABAD( SMLNUM, BIGNUM )
248:       SMLNUM = SMLNUM / SLAMCH( 'Precision' )
249:       BIGNUM = ONE / SMLNUM
250:       SCALE = ONE
251: *
252:       IF( LSAME( NORMIN, 'N' ) ) THEN
253: *
254: *        Compute the 1-norm of each column, not including the diagonal.
255: *
256:          IF( UPPER ) THEN
257: *
258: *           A is upper triangular.
259: *
260:             DO 10 J = 1, N
261:                JLEN = MIN( KD, J-1 )
262:                CNORM( J ) = SCASUM( JLEN, AB( KD+1-JLEN, J ), 1 )
263:    10       CONTINUE
264:          ELSE
265: *
266: *           A is lower triangular.
267: *
268:             DO 20 J = 1, N
269:                JLEN = MIN( KD, N-J )
270:                IF( JLEN.GT.0 ) THEN
271:                   CNORM( J ) = SCASUM( JLEN, AB( 2, J ), 1 )
272:                ELSE
273:                   CNORM( J ) = ZERO
274:                END IF
275:    20       CONTINUE
276:          END IF
277:       END IF
278: *
279: *     Scale the column norms by TSCAL if the maximum element in CNORM is
280: *     greater than BIGNUM/2.
281: *
282:       IMAX = ISAMAX( N, CNORM, 1 )
283:       TMAX = CNORM( IMAX )
284:       IF( TMAX.LE.BIGNUM*HALF ) THEN
285:          TSCAL = ONE
286:       ELSE
287:          TSCAL = HALF / ( SMLNUM*TMAX )
288:          CALL SSCAL( N, TSCAL, CNORM, 1 )
289:       END IF
290: *
291: *     Compute a bound on the computed solution vector to see if the
292: *     Level 2 BLAS routine CTBSV can be used.
293: *
294:       XMAX = ZERO
295:       DO 30 J = 1, N
296:          XMAX = MAX( XMAX, CABS2( X( J ) ) )
297:    30 CONTINUE
298:       XBND = XMAX
299:       IF( NOTRAN ) THEN
300: *
301: *        Compute the growth in A * x = b.
302: *
303:          IF( UPPER ) THEN
304:             JFIRST = N
305:             JLAST = 1
306:             JINC = -1
307:             MAIND = KD + 1
308:          ELSE
309:             JFIRST = 1
310:             JLAST = N
311:             JINC = 1
312:             MAIND = 1
313:          END IF
314: *
315:          IF( TSCAL.NE.ONE ) THEN
316:             GROW = ZERO
317:             GO TO 60
318:          END IF
319: *
320:          IF( NOUNIT ) THEN
321: *
322: *           A is non-unit triangular.
323: *
324: *           Compute GROW = 1/G(j) and XBND = 1/M(j).
325: *           Initially, G(0) = max{x(i), i=1,...,n}.
326: *
327:             GROW = HALF / MAX( XBND, SMLNUM )
328:             XBND = GROW
329:             DO 40 J = JFIRST, JLAST, JINC
330: *
331: *              Exit the loop if the growth factor is too small.
332: *
333:                IF( GROW.LE.SMLNUM )
334:      $            GO TO 60
335: *
336:                TJJS = AB( MAIND, J )
337:                TJJ = CABS1( TJJS )
338: *
339:                IF( TJJ.GE.SMLNUM ) THEN
340: *
341: *                 M(j) = G(j-1) / abs(A(j,j))
342: *
343:                   XBND = MIN( XBND, MIN( ONE, TJJ )*GROW )
344:                ELSE
345: *
346: *                 M(j) could overflow, set XBND to 0.
347: *
348:                   XBND = ZERO
349:                END IF
350: *
351:                IF( TJJ+CNORM( J ).GE.SMLNUM ) THEN
352: *
353: *                 G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
354: *
355:                   GROW = GROW*( TJJ / ( TJJ+CNORM( J ) ) )
356:                ELSE
357: *
358: *                 G(j) could overflow, set GROW to 0.
359: *
360:                   GROW = ZERO
361:                END IF
362:    40       CONTINUE
363:             GROW = XBND
364:          ELSE
365: *
366: *           A is unit triangular.
367: *
368: *           Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
369: *
370:             GROW = MIN( ONE, HALF / MAX( XBND, SMLNUM ) )
371:             DO 50 J = JFIRST, JLAST, JINC
372: *
373: *              Exit the loop if the growth factor is too small.
374: *
375:                IF( GROW.LE.SMLNUM )
376:      $            GO TO 60
377: *
378: *              G(j) = G(j-1)*( 1 + CNORM(j) )
379: *
380:                GROW = GROW*( ONE / ( ONE+CNORM( J ) ) )
381:    50       CONTINUE
382:          END IF
383:    60    CONTINUE
384: *
385:       ELSE
386: *
387: *        Compute the growth in A**T * x = b  or  A**H * x = b.
388: *
389:          IF( UPPER ) THEN
390:             JFIRST = 1
391:             JLAST = N
392:             JINC = 1
393:             MAIND = KD + 1
394:          ELSE
395:             JFIRST = N
396:             JLAST = 1
397:             JINC = -1
398:             MAIND = 1
399:          END IF
400: *
401:          IF( TSCAL.NE.ONE ) THEN
402:             GROW = ZERO
403:             GO TO 90
404:          END IF
405: *
406:          IF( NOUNIT ) THEN
407: *
408: *           A is non-unit triangular.
409: *
410: *           Compute GROW = 1/G(j) and XBND = 1/M(j).
411: *           Initially, M(0) = max{x(i), i=1,...,n}.
412: *
413:             GROW = HALF / MAX( XBND, SMLNUM )
414:             XBND = GROW
415:             DO 70 J = JFIRST, JLAST, JINC
416: *
417: *              Exit the loop if the growth factor is too small.
418: *
419:                IF( GROW.LE.SMLNUM )
420:      $            GO TO 90
421: *
422: *              G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
423: *
424:                XJ = ONE + CNORM( J )
425:                GROW = MIN( GROW, XBND / XJ )
426: *
427:                TJJS = AB( MAIND, J )
428:                TJJ = CABS1( TJJS )
429: *
430:                IF( TJJ.GE.SMLNUM ) THEN
431: *
432: *                 M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
433: *
434:                   IF( XJ.GT.TJJ )
435:      $               XBND = XBND*( TJJ / XJ )
436:                ELSE
437: *
438: *                 M(j) could overflow, set XBND to 0.
439: *
440:                   XBND = ZERO
441:                END IF
442:    70       CONTINUE
443:             GROW = MIN( GROW, XBND )
444:          ELSE
445: *
446: *           A is unit triangular.
447: *
448: *           Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
449: *
450:             GROW = MIN( ONE, HALF / MAX( XBND, SMLNUM ) )
451:             DO 80 J = JFIRST, JLAST, JINC
452: *
453: *              Exit the loop if the growth factor is too small.
454: *
455:                IF( GROW.LE.SMLNUM )
456:      $            GO TO 90
457: *
458: *              G(j) = ( 1 + CNORM(j) )*G(j-1)
459: *
460:                XJ = ONE + CNORM( J )
461:                GROW = GROW / XJ
462:    80       CONTINUE
463:          END IF
464:    90    CONTINUE
465:       END IF
466: *
467:       IF( ( GROW*TSCAL ).GT.SMLNUM ) THEN
468: *
469: *        Use the Level 2 BLAS solve if the reciprocal of the bound on
470: *        elements of X is not too small.
471: *
472:          CALL CTBSV( UPLO, TRANS, DIAG, N, KD, AB, LDAB, X, 1 )
473:       ELSE
474: *
475: *        Use a Level 1 BLAS solve, scaling intermediate results.
476: *
477:          IF( XMAX.GT.BIGNUM*HALF ) THEN
478: *
479: *           Scale X so that its components are less than or equal to
480: *           BIGNUM in absolute value.
481: *
482:             SCALE = ( BIGNUM*HALF ) / XMAX
483:             CALL CSSCAL( N, SCALE, X, 1 )
484:             XMAX = BIGNUM
485:          ELSE
486:             XMAX = XMAX*TWO
487:          END IF
488: *
489:          IF( NOTRAN ) THEN
490: *
491: *           Solve A * x = b
492: *
493:             DO 110 J = JFIRST, JLAST, JINC
494: *
495: *              Compute x(j) = b(j) / A(j,j), scaling x if necessary.
496: *
497:                XJ = CABS1( X( J ) )
498:                IF( NOUNIT ) THEN
499:                   TJJS = AB( MAIND, J )*TSCAL
500:                ELSE
501:                   TJJS = TSCAL
502:                   IF( TSCAL.EQ.ONE )
503:      $               GO TO 105
504:                END IF
505:                   TJJ = CABS1( TJJS )
506:                   IF( TJJ.GT.SMLNUM ) THEN
507: *
508: *                    abs(A(j,j)) > SMLNUM:
509: *
510:                      IF( TJJ.LT.ONE ) THEN
511:                         IF( XJ.GT.TJJ*BIGNUM ) THEN
512: *
513: *                          Scale x by 1/b(j).
514: *
515:                            REC = ONE / XJ
516:                            CALL CSSCAL( N, REC, X, 1 )
517:                            SCALE = SCALE*REC
518:                            XMAX = XMAX*REC
519:                         END IF
520:                      END IF
521:                      X( J ) = CLADIV( X( J ), TJJS )
522:                      XJ = CABS1( X( J ) )
523:                   ELSE IF( TJJ.GT.ZERO ) THEN
524: *
525: *                    0 < abs(A(j,j)) <= SMLNUM:
526: *
527:                      IF( XJ.GT.TJJ*BIGNUM ) THEN
528: *
529: *                       Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
530: *                       to avoid overflow when dividing by A(j,j).
531: *
532:                         REC = ( TJJ*BIGNUM ) / XJ
533:                         IF( CNORM( J ).GT.ONE ) THEN
534: *
535: *                          Scale by 1/CNORM(j) to avoid overflow when
536: *                          multiplying x(j) times column j.
537: *
538:                            REC = REC / CNORM( J )
539:                         END IF
540:                         CALL CSSCAL( N, REC, X, 1 )
541:                         SCALE = SCALE*REC
542:                         XMAX = XMAX*REC
543:                      END IF
544:                      X( J ) = CLADIV( X( J ), TJJS )
545:                      XJ = CABS1( X( J ) )
546:                   ELSE
547: *
548: *                    A(j,j) = 0:  Set x(1:n) = 0, x(j) = 1, and
549: *                    scale = 0, and compute a solution to A*x = 0.
550: *
551:                      DO 100 I = 1, N
552:                         X( I ) = ZERO
553:   100                CONTINUE
554:                      X( J ) = ONE
555:                      XJ = ONE
556:                      SCALE = ZERO
557:                      XMAX = ZERO
558:                   END IF
559:   105          CONTINUE
560: *
561: *              Scale x if necessary to avoid overflow when adding a
562: *              multiple of column j of A.
563: *
564:                IF( XJ.GT.ONE ) THEN
565:                   REC = ONE / XJ
566:                   IF( CNORM( J ).GT.( BIGNUM-XMAX )*REC ) THEN
567: *
568: *                    Scale x by 1/(2*abs(x(j))).
569: *
570:                      REC = REC*HALF
571:                      CALL CSSCAL( N, REC, X, 1 )
572:                      SCALE = SCALE*REC
573:                   END IF
574:                ELSE IF( XJ*CNORM( J ).GT.( BIGNUM-XMAX ) ) THEN
575: *
576: *                 Scale x by 1/2.
577: *
578:                   CALL CSSCAL( N, HALF, X, 1 )
579:                   SCALE = SCALE*HALF
580:                END IF
581: *
582:                IF( UPPER ) THEN
583:                   IF( J.GT.1 ) THEN
584: *
585: *                    Compute the update
586: *                       x(max(1,j-kd):j-1) := x(max(1,j-kd):j-1) -
587: *                                             x(j)* A(max(1,j-kd):j-1,j)
588: *
589:                      JLEN = MIN( KD, J-1 )
590:                      CALL CAXPY( JLEN, -X( J )*TSCAL,
591:      $                           AB( KD+1-JLEN, J ), 1, X( J-JLEN ), 1 )
592:                      I = ICAMAX( J-1, X, 1 )
593:                      XMAX = CABS1( X( I ) )
594:                   END IF
595:                ELSE IF( J.LT.N ) THEN
596: *
597: *                 Compute the update
598: *                    x(j+1:min(j+kd,n)) := x(j+1:min(j+kd,n)) -
599: *                                          x(j) * A(j+1:min(j+kd,n),j)
600: *
601:                   JLEN = MIN( KD, N-J )
602:                   IF( JLEN.GT.0 )
603:      $               CALL CAXPY( JLEN, -X( J )*TSCAL, AB( 2, J ), 1,
604:      $                           X( J+1 ), 1 )
605:                   I = J + ICAMAX( N-J, X( J+1 ), 1 )
606:                   XMAX = CABS1( X( I ) )
607:                END IF
608:   110       CONTINUE
609: *
610:          ELSE IF( LSAME( TRANS, 'T' ) ) THEN
611: *
612: *           Solve A**T * x = b
613: *
614:             DO 150 J = JFIRST, JLAST, JINC
615: *
616: *              Compute x(j) = b(j) - sum A(k,j)*x(k).
617: *                                    k<>j
618: *
619:                XJ = CABS1( X( J ) )
620:                USCAL = TSCAL
621:                REC = ONE / MAX( XMAX, ONE )
622:                IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN
623: *
624: *                 If x(j) could overflow, scale x by 1/(2*XMAX).
625: *
626:                   REC = REC*HALF
627:                   IF( NOUNIT ) THEN
628:                      TJJS = AB( MAIND, J )*TSCAL
629:                   ELSE
630:                      TJJS = TSCAL
631:                   END IF
632:                      TJJ = CABS1( TJJS )
633:                      IF( TJJ.GT.ONE ) THEN
634: *
635: *                       Divide by A(j,j) when scaling x if A(j,j) > 1.
636: *
637:                         REC = MIN( ONE, REC*TJJ )
638:                         USCAL = CLADIV( USCAL, TJJS )
639:                      END IF
640:                   IF( REC.LT.ONE ) THEN
641:                      CALL CSSCAL( N, REC, X, 1 )
642:                      SCALE = SCALE*REC
643:                      XMAX = XMAX*REC
644:                   END IF
645:                END IF
646: *
647:                CSUMJ = ZERO
648:                IF( USCAL.EQ.CMPLX( ONE ) ) THEN
649: *
650: *                 If the scaling needed for A in the dot product is 1,
651: *                 call CDOTU to perform the dot product.
652: *
653:                   IF( UPPER ) THEN
654:                      JLEN = MIN( KD, J-1 )
655:                      CSUMJ = CDOTU( JLEN, AB( KD+1-JLEN, J ), 1,
656:      $                       X( J-JLEN ), 1 )
657:                   ELSE
658:                      JLEN = MIN( KD, N-J )
659:                      IF( JLEN.GT.1 )
660:      $                  CSUMJ = CDOTU( JLEN, AB( 2, J ), 1, X( J+1 ),
661:      $                          1 )
662:                   END IF
663:                ELSE
664: *
665: *                 Otherwise, use in-line code for the dot product.
666: *
667:                   IF( UPPER ) THEN
668:                      JLEN = MIN( KD, J-1 )
669:                      DO 120 I = 1, JLEN
670:                         CSUMJ = CSUMJ + ( AB( KD+I-JLEN, J )*USCAL )*
671:      $                          X( J-JLEN-1+I )
672:   120                CONTINUE
673:                   ELSE
674:                      JLEN = MIN( KD, N-J )
675:                      DO 130 I = 1, JLEN
676:                         CSUMJ = CSUMJ + ( AB( I+1, J )*USCAL )*X( J+I )
677:   130                CONTINUE
678:                   END IF
679:                END IF
680: *
681:                IF( USCAL.EQ.CMPLX( TSCAL ) ) THEN
682: *
683: *                 Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
684: *                 was not used to scale the dotproduct.
685: *
686:                   X( J ) = X( J ) - CSUMJ
687:                   XJ = CABS1( X( J ) )
688:                   IF( NOUNIT ) THEN
689: *
690: *                    Compute x(j) = x(j) / A(j,j), scaling if necessary.
691: *
692:                      TJJS = AB( MAIND, J )*TSCAL
693:                   ELSE
694:                      TJJS = TSCAL
695:                      IF( TSCAL.EQ.ONE )
696:      $                  GO TO 145
697:                   END IF
698:                      TJJ = CABS1( TJJS )
699:                      IF( TJJ.GT.SMLNUM ) THEN
700: *
701: *                       abs(A(j,j)) > SMLNUM:
702: *
703:                         IF( TJJ.LT.ONE ) THEN
704:                            IF( XJ.GT.TJJ*BIGNUM ) THEN
705: *
706: *                             Scale X by 1/abs(x(j)).
707: *
708:                               REC = ONE / XJ
709:                               CALL CSSCAL( N, REC, X, 1 )
710:                               SCALE = SCALE*REC
711:                               XMAX = XMAX*REC
712:                            END IF
713:                         END IF
714:                         X( J ) = CLADIV( X( J ), TJJS )
715:                      ELSE IF( TJJ.GT.ZERO ) THEN
716: *
717: *                       0 < abs(A(j,j)) <= SMLNUM:
718: *
719:                         IF( XJ.GT.TJJ*BIGNUM ) THEN
720: *
721: *                          Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
722: *
723:                            REC = ( TJJ*BIGNUM ) / XJ
724:                            CALL CSSCAL( N, REC, X, 1 )
725:                            SCALE = SCALE*REC
726:                            XMAX = XMAX*REC
727:                         END IF
728:                         X( J ) = CLADIV( X( J ), TJJS )
729:                      ELSE
730: *
731: *                       A(j,j) = 0:  Set x(1:n) = 0, x(j) = 1, and
732: *                       scale = 0 and compute a solution to A**T *x = 0.
733: *
734:                         DO 140 I = 1, N
735:                            X( I ) = ZERO
736:   140                   CONTINUE
737:                         X( J ) = ONE
738:                         SCALE = ZERO
739:                         XMAX = ZERO
740:                      END IF
741:   145             CONTINUE
742:                ELSE
743: *
744: *                 Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
745: *                 product has already been divided by 1/A(j,j).
746: *
747:                   X( J ) = CLADIV( X( J ), TJJS ) - CSUMJ
748:                END IF
749:                XMAX = MAX( XMAX, CABS1( X( J ) ) )
750:   150       CONTINUE
751: *
752:          ELSE
753: *
754: *           Solve A**H * x = b
755: *
756:             DO 190 J = JFIRST, JLAST, JINC
757: *
758: *              Compute x(j) = b(j) - sum A(k,j)*x(k).
759: *                                    k<>j
760: *
761:                XJ = CABS1( X( J ) )
762:                USCAL = TSCAL
763:                REC = ONE / MAX( XMAX, ONE )
764:                IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN
765: *
766: *                 If x(j) could overflow, scale x by 1/(2*XMAX).
767: *
768:                   REC = REC*HALF
769:                   IF( NOUNIT ) THEN
770:                      TJJS = CONJG( AB( MAIND, J ) )*TSCAL
771:                   ELSE
772:                      TJJS = TSCAL
773:                   END IF
774:                      TJJ = CABS1( TJJS )
775:                      IF( TJJ.GT.ONE ) THEN
776: *
777: *                       Divide by A(j,j) when scaling x if A(j,j) > 1.
778: *
779:                         REC = MIN( ONE, REC*TJJ )
780:                         USCAL = CLADIV( USCAL, TJJS )
781:                      END IF
782:                   IF( REC.LT.ONE ) THEN
783:                      CALL CSSCAL( N, REC, X, 1 )
784:                      SCALE = SCALE*REC
785:                      XMAX = XMAX*REC
786:                   END IF
787:                END IF
788: *
789:                CSUMJ = ZERO
790:                IF( USCAL.EQ.CMPLX( ONE ) ) THEN
791: *
792: *                 If the scaling needed for A in the dot product is 1,
793: *                 call CDOTC to perform the dot product.
794: *
795:                   IF( UPPER ) THEN
796:                      JLEN = MIN( KD, J-1 )
797:                      CSUMJ = CDOTC( JLEN, AB( KD+1-JLEN, J ), 1,
798:      $                       X( J-JLEN ), 1 )
799:                   ELSE
800:                      JLEN = MIN( KD, N-J )
801:                      IF( JLEN.GT.1 )
802:      $                  CSUMJ = CDOTC( JLEN, AB( 2, J ), 1, X( J+1 ),
803:      $                          1 )
804:                   END IF
805:                ELSE
806: *
807: *                 Otherwise, use in-line code for the dot product.
808: *
809:                   IF( UPPER ) THEN
810:                      JLEN = MIN( KD, J-1 )
811:                      DO 160 I = 1, JLEN
812:                         CSUMJ = CSUMJ + ( CONJG( AB( KD+I-JLEN, J ) )*
813:      $                          USCAL )*X( J-JLEN-1+I )
814:   160                CONTINUE
815:                   ELSE
816:                      JLEN = MIN( KD, N-J )
817:                      DO 170 I = 1, JLEN
818:                         CSUMJ = CSUMJ + ( CONJG( AB( I+1, J ) )*USCAL )*
819:      $                          X( J+I )
820:   170                CONTINUE
821:                   END IF
822:                END IF
823: *
824:                IF( USCAL.EQ.CMPLX( TSCAL ) ) THEN
825: *
826: *                 Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
827: *                 was not used to scale the dotproduct.
828: *
829:                   X( J ) = X( J ) - CSUMJ
830:                   XJ = CABS1( X( J ) )
831:                   IF( NOUNIT ) THEN
832: *
833: *                    Compute x(j) = x(j) / A(j,j), scaling if necessary.
834: *
835:                      TJJS = CONJG( AB( MAIND, J ) )*TSCAL
836:                   ELSE
837:                      TJJS = TSCAL
838:                      IF( TSCAL.EQ.ONE )
839:      $                  GO TO 185
840:                   END IF
841:                      TJJ = CABS1( TJJS )
842:                      IF( TJJ.GT.SMLNUM ) THEN
843: *
844: *                       abs(A(j,j)) > SMLNUM:
845: *
846:                         IF( TJJ.LT.ONE ) THEN
847:                            IF( XJ.GT.TJJ*BIGNUM ) THEN
848: *
849: *                             Scale X by 1/abs(x(j)).
850: *
851:                               REC = ONE / XJ
852:                               CALL CSSCAL( N, REC, X, 1 )
853:                               SCALE = SCALE*REC
854:                               XMAX = XMAX*REC
855:                            END IF
856:                         END IF
857:                         X( J ) = CLADIV( X( J ), TJJS )
858:                      ELSE IF( TJJ.GT.ZERO ) THEN
859: *
860: *                       0 < abs(A(j,j)) <= SMLNUM:
861: *
862:                         IF( XJ.GT.TJJ*BIGNUM ) THEN
863: *
864: *                          Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
865: *
866:                            REC = ( TJJ*BIGNUM ) / XJ
867:                            CALL CSSCAL( N, REC, X, 1 )
868:                            SCALE = SCALE*REC
869:                            XMAX = XMAX*REC
870:                         END IF
871:                         X( J ) = CLADIV( X( J ), TJJS )
872:                      ELSE
873: *
874: *                       A(j,j) = 0:  Set x(1:n) = 0, x(j) = 1, and
875: *                       scale = 0 and compute a solution to A**H *x = 0.
876: *
877:                         DO 180 I = 1, N
878:                            X( I ) = ZERO
879:   180                   CONTINUE
880:                         X( J ) = ONE
881:                         SCALE = ZERO
882:                         XMAX = ZERO
883:                      END IF
884:   185             CONTINUE
885:                ELSE
886: *
887: *                 Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
888: *                 product has already been divided by 1/A(j,j).
889: *
890:                   X( J ) = CLADIV( X( J ), TJJS ) - CSUMJ
891:                END IF
892:                XMAX = MAX( XMAX, CABS1( X( J ) ) )
893:   190       CONTINUE
894:          END IF
895:          SCALE = SCALE / TSCAL
896:       END IF
897: *
898: *     Scale the column norms by 1/TSCAL for return.
899: *
900:       IF( TSCAL.NE.ONE ) THEN
901:          CALL SSCAL( N, ONE / TSCAL, CNORM, 1 )
902:       END IF
903: *
904:       RETURN
905: *
906: *     End of CLATBS
907: *
908:       END
909: