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