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