SRC/zlatbs.f

   1 *> \brief \b ZLATBS solves a triangular banded system of equations.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZLATBS + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlatbs.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlatbs.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlatbs.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X,
  22 *                          SCALE, CNORM, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          DIAG, NORMIN, TRANS, UPLO
  26 *       INTEGER            INFO, KD, LDAB, N
  27 *       DOUBLE PRECISION   SCALE
  28 *       ..
  29 *       .. Array Arguments ..
  30 *       DOUBLE PRECISION   CNORM( * )
  31 *       COMPLEX*16         AB( LDAB, * ), X( * )
  32 *       ..
  33 *
  34 *
  35 *> \par Purpose:
  36 *  =============
  37 *>
  38 *> \verbatim
  39 *>
  40 *> ZLATBS solves one of the triangular systems
  41 *>
  42 *>    A * x = s*b,  A**T * x = s*b,  or  A**H * x = s*b,
  43 *>
  44 *> with scaling to prevent overflow, where A is an upper or lower
  45 *> triangular band matrix.  Here A**T denotes the transpose of A, x and b
  46 *> are n-element vectors, and s is a scaling factor, usually less than
  47 *> or equal to 1, chosen so that the components of x will be less than
  48 *> the overflow threshold.  If the unscaled problem will not cause
  49 *> overflow, the Level 2 BLAS routine ZTBSV is called.  If the matrix A
  50 *> is singular (A(j,j) = 0 for some j), then s is set to 0 and a
  51 *> non-trivial solution to A*x = 0 is returned.
  52 *> \endverbatim
  53 *
  54 *  Arguments:
  55 *  ==========
  56 *
  57 *> \param[in] UPLO
  58 *> \verbatim
  59 *>          UPLO is CHARACTER*1
  60 *>          Specifies whether the matrix A is upper or lower triangular.
  61 *>          = 'U':  Upper triangular
  62 *>          = 'L':  Lower triangular
  63 *> \endverbatim
  64 *>
  65 *> \param[in] TRANS
  66 *> \verbatim
  67 *>          TRANS is CHARACTER*1
  68 *>          Specifies the operation applied to A.
  69 *>          = 'N':  Solve A * x = s*b     (No transpose)
  70 *>          = 'T':  Solve A**T * x = s*b  (Transpose)
  71 *>          = 'C':  Solve A**H * x = s*b  (Conjugate transpose)
  72 *> \endverbatim
  73 *>
  74 *> \param[in] DIAG
  75 *> \verbatim
  76 *>          DIAG is CHARACTER*1
  77 *>          Specifies whether or not the matrix A is unit triangular.
  78 *>          = 'N':  Non-unit triangular
  79 *>          = 'U':  Unit triangular
  80 *> \endverbatim
  81 *>
  82 *> \param[in] NORMIN
  83 *> \verbatim
  84 *>          NORMIN is CHARACTER*1
  85 *>          Specifies whether CNORM has been set or not.
  86 *>          = 'Y':  CNORM contains the column norms on entry
  87 *>          = 'N':  CNORM is not set on entry.  On exit, the norms will
  88 *>                  be computed and stored in CNORM.
  89 *> \endverbatim
  90 *>
  91 *> \param[in] N
  92 *> \verbatim
  93 *>          N is INTEGER
  94 *>          The order of the matrix A.  N >= 0.
  95 *> \endverbatim
  96 *>
  97 *> \param[in] KD
  98 *> \verbatim
  99 *>          KD is INTEGER
 100 *>          The number of subdiagonals or superdiagonals in the
 101 *>          triangular matrix A.  KD >= 0.
 102 *> \endverbatim
 103 *>
 104 *> \param[in] AB
 105 *> \verbatim
 106 *>          AB is COMPLEX*16 array, dimension (LDAB,N)
 107 *>          The upper or lower triangular band matrix A, stored in the
 108 *>          first KD+1 rows of the array. The j-th column of A is stored
 109 *>          in the j-th column of the array AB as follows:
 110 *>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
 111 *>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
 112 *> \endverbatim
 113 *>
 114 *> \param[in] LDAB
 115 *> \verbatim
 116 *>          LDAB is INTEGER
 117 *>          The leading dimension of the array AB.  LDAB >= KD+1.
 118 *> \endverbatim
 119 *>
 120 *> \param[in,out] X
 121 *> \verbatim
 122 *>          X is COMPLEX*16 array, dimension (N)
 123 *>          On entry, the right hand side b of the triangular system.
 124 *>          On exit, X is overwritten by the solution vector x.
 125 *> \endverbatim
 126 *>
 127 *> \param[out] SCALE
 128 *> \verbatim
 129 *>          SCALE is DOUBLE PRECISION
 130 *>          The scaling factor s for the triangular system
 131 *>             A * x = s*b,  A**T * x = s*b,  or  A**H * x = s*b.
 132 *>          If SCALE = 0, the matrix A is singular or badly scaled, and
 133 *>          the vector x is an exact or approximate solution to A*x = 0.
 134 *> \endverbatim
 135 *>
 136 *> \param[in,out] CNORM
 137 *> \verbatim
 138 *>          CNORM is DOUBLE PRECISION array, dimension (N)
 139 *>
 140 *>          If NORMIN = 'Y', CNORM is an input argument and CNORM(j)
 141 *>          contains the norm of the off-diagonal part of the j-th column
 142 *>          of A.  If TRANS = 'N', CNORM(j) must be greater than or equal
 143 *>          to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)
 144 *>          must be greater than or equal to the 1-norm.
 145 *>
 146 *>          If NORMIN = 'N', CNORM is an output argument and CNORM(j)
 147 *>          returns the 1-norm of the offdiagonal part of the j-th column
 148 *>          of A.
 149 *> \endverbatim
 150 *>
 151 *> \param[out] INFO
 152 *> \verbatim
 153 *>          INFO is INTEGER
 154 *>          = 0:  successful exit
 155 *>          < 0:  if INFO = -k, the k-th argument had an illegal value
 156 *> \endverbatim
 157 *
 158 *  Authors:
 159 *  ========
 160 *
 161 *> \author Univ. of Tennessee
 162 *> \author Univ. of California Berkeley
 163 *> \author Univ. of Colorado Denver
 164 *> \author NAG Ltd.
 165 *
 166 *> \date September 2012
 167 *
 168 *> \ingroup complex16OTHERauxiliary
 169 *
 170 *> \par Further Details:
 171 *  =====================
 172 *>
 173 *> \verbatim
 174 *>
 175 *>  A rough bound on x is computed; if that is less than overflow, ZTBSV
 176 *>  is called, otherwise, specific code is used which checks for possible
 177 *>  overflow or divide-by-zero at every operation.
 178 *>
 179 *>  A columnwise scheme is used for solving A*x = b.  The basic algorithm
 180 *>  if A is lower triangular is
 181 *>
 182 *>       x[1:n] := b[1:n]
 183 *>       for j = 1, ..., n
 184 *>            x(j) := x(j) / A(j,j)
 185 *>            x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
 186 *>       end
 187 *>
 188 *>  Define bounds on the components of x after j iterations of the loop:
 189 *>     M(j) = bound on x[1:j]
 190 *>     G(j) = bound on x[j+1:n]
 191 *>  Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
 192 *>
 193 *>  Then for iteration j+1 we have
 194 *>     M(j+1) <= G(j) / | A(j+1,j+1) |
 195 *>     G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
 196 *>            <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
 197 *>
 198 *>  where CNORM(j+1) is greater than or equal to the infinity-norm of
 199 *>  column j+1 of A, not counting the diagonal.  Hence
 200 *>
 201 *>     G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
 202 *>                  1<=i<=j
 203 *>  and
 204 *>
 205 *>     |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
 206 *>                                   1<=i< j
 207 *>
 208 *>  Since |x(j)| <= M(j), we use the Level 2 BLAS routine ZTBSV if the
 209 *>  reciprocal of the largest M(j), j=1,..,n, is larger than
 210 *>  max(underflow, 1/overflow).
 211 *>
 212 *>  The bound on x(j) is also used to determine when a step in the
 213 *>  columnwise method can be performed without fear of overflow.  If
 214 *>  the computed bound is greater than a large constant, x is scaled to
 215 *>  prevent overflow, but if the bound overflows, x is set to 0, x(j) to
 216 *>  1, and scale to 0, and a non-trivial solution to A*x = 0 is found.
 217 *>
 218 *>  Similarly, a row-wise scheme is used to solve A**T *x = b  or
 219 *>  A**H *x = b.  The basic algorithm for A upper triangular is
 220 *>
 221 *>       for j = 1, ..., n
 222 *>            x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)
 223 *>       end
 224 *>
 225 *>  We simultaneously compute two bounds
 226 *>       G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j
 227 *>       M(j) = bound on x(i), 1<=i<=j
 228 *>
 229 *>  The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we
 230 *>  add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.
 231 *>  Then the bound on x(j) is
 232 *>
 233 *>       M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
 234 *>
 235 *>            <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
 236 *>                      1<=i<=j
 237 *>
 238 *>  and we can safely call ZTBSV if 1/M(n) and 1/G(n) are both greater
 239 *>  than max(underflow, 1/overflow).
 240 *> \endverbatim
 241 *>
 242 *  =====================================================================
 243       SUBROUTINE ZLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X,
 244      $                   SCALE, CNORM, INFO )
 245 *
 246 *  -- LAPACK auxiliary routine (version 3.4.2) --
 247 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 248 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 249 *     September 2012
 250 *
 251 *     .. Scalar Arguments ..
 252       CHARACTER          DIAG, NORMIN, TRANS, UPLO
 253       INTEGER            INFO, KD, LDAB, N
 254       DOUBLE PRECISION   SCALE
 255 *     ..
 256 *     .. Array Arguments ..
 257       DOUBLE PRECISION   CNORM( * )
 258       COMPLEX*16         AB( LDAB, * ), X( * )
 259 *     ..
 260 *
 261 *  =====================================================================
 262 *
 263 *     .. Parameters ..
 264       DOUBLE PRECISION   ZERO, HALF, ONE, TWO
 265       PARAMETER          ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0,
 266      $                   TWO = 2.0D+0 )
 267 *     ..
 268 *     .. Local Scalars ..
 269       LOGICAL            NOTRAN, NOUNIT, UPPER
 270       INTEGER            I, IMAX, J, JFIRST, JINC, JLAST, JLEN, MAIND
 271       DOUBLE PRECISION   BIGNUM, GROW, REC, SMLNUM, TJJ, TMAX, TSCAL,
 272      $                   XBND, XJ, XMAX
 273       COMPLEX*16         CSUMJ, TJJS, USCAL, ZDUM
 274 *     ..
 275 *     .. External Functions ..
 276       LOGICAL            LSAME
 277       INTEGER            IDAMAX, IZAMAX
 278       DOUBLE PRECISION   DLAMCH, DZASUM
 279       COMPLEX*16         ZDOTC, ZDOTU, ZLADIV
 280       EXTERNAL           LSAME, IDAMAX, IZAMAX, DLAMCH, DZASUM, ZDOTC,
 281      $                   ZDOTU, ZLADIV
 282 *     ..
 283 *     .. External Subroutines ..
 284       EXTERNAL           DSCAL, XERBLA, ZAXPY, ZDSCAL, ZTBSV
 285 *     ..
 286 *     .. Intrinsic Functions ..
 287       INTRINSIC          ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN
 288 *     ..
 289 *     .. Statement Functions ..
 290       DOUBLE PRECISION   CABS1, CABS2
 291 *     ..
 292 *     .. Statement Function definitions ..
 293       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
 294       CABS2( ZDUM ) = ABS( DBLE( ZDUM ) / 2.D0 ) +
 295      $                ABS( DIMAG( ZDUM ) / 2.D0 )
 296 *     ..
 297 *     .. Executable Statements ..
 298 *
 299       INFO = 0
 300       UPPER = LSAME( UPLO, 'U' )
 301       NOTRAN = LSAME( TRANS, 'N' )
 302       NOUNIT = LSAME( DIAG, 'N' )
 303 *
 304 *     Test the input parameters.
 305 *
 306       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 307          INFO = -1
 308       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
 309      $         LSAME( TRANS, 'C' ) ) THEN
 310          INFO = -2
 311       ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
 312          INFO = -3
 313       ELSE IF( .NOT.LSAME( NORMIN, 'Y' ) .AND. .NOT.
 314      $         LSAME( NORMIN, 'N' ) ) THEN
 315          INFO = -4
 316       ELSE IF( N.LT.0 ) THEN
 317          INFO = -5
 318       ELSE IF( KD.LT.0 ) THEN
 319          INFO = -6
 320       ELSE IF( LDAB.LT.KD+1 ) THEN
 321          INFO = -8
 322       END IF
 323       IF( INFO.NE.0 ) THEN
 324          CALL XERBLA( 'ZLATBS', -INFO )
 325          RETURN
 326       END IF
 327 *
 328 *     Quick return if possible
 329 *
 330       IF( N.EQ.0 )
 331      $   RETURN
 332 *
 333 *     Determine machine dependent parameters to control overflow.
 334 *
 335       SMLNUM = DLAMCH( 'Safe minimum' )
 336       BIGNUM = ONE / SMLNUM
 337       CALL DLABAD( SMLNUM, BIGNUM )
 338       SMLNUM = SMLNUM / DLAMCH( 'Precision' )
 339       BIGNUM = ONE / SMLNUM
 340       SCALE = ONE
 341 *
 342       IF( LSAME( NORMIN, 'N' ) ) THEN
 343 *
 344 *        Compute the 1-norm of each column, not including the diagonal.
 345 *
 346          IF( UPPER ) THEN
 347 *
 348 *           A is upper triangular.
 349 *
 350             DO 10 J = 1, N
 351                JLEN = MIN( KD, J-1 )
 352                CNORM( J ) = DZASUM( JLEN, AB( KD+1-JLEN, J ), 1 )
 353    10       CONTINUE
 354          ELSE
 355 *
 356 *           A is lower triangular.
 357 *
 358             DO 20 J = 1, N
 359                JLEN = MIN( KD, N-J )
 360                IF( JLEN.GT.0 ) THEN
 361                   CNORM( J ) = DZASUM( JLEN, AB( 2, J ), 1 )
 362                ELSE
 363                   CNORM( J ) = ZERO
 364                END IF
 365    20       CONTINUE
 366          END IF
 367       END IF
 368 *
 369 *     Scale the column norms by TSCAL if the maximum element in CNORM is
 370 *     greater than BIGNUM/2.
 371 *
 372       IMAX = IDAMAX( N, CNORM, 1 )
 373       TMAX = CNORM( IMAX )
 374       IF( TMAX.LE.BIGNUM*HALF ) THEN
 375          TSCAL = ONE
 376       ELSE
 377          TSCAL = HALF / ( SMLNUM*TMAX )
 378          CALL DSCAL( N, TSCAL, CNORM, 1 )
 379       END IF
 380 *
 381 *     Compute a bound on the computed solution vector to see if the
 382 *     Level 2 BLAS routine ZTBSV can be used.
 383 *
 384       XMAX = ZERO
 385       DO 30 J = 1, N
 386          XMAX = MAX( XMAX, CABS2( X( J ) ) )
 387    30 CONTINUE
 388       XBND = XMAX
 389       IF( NOTRAN ) THEN
 390 *
 391 *        Compute the growth in A * x = b.
 392 *
 393          IF( UPPER ) THEN
 394             JFIRST = N
 395             JLAST = 1
 396             JINC = -1
 397             MAIND = KD + 1
 398          ELSE
 399             JFIRST = 1
 400             JLAST = N
 401             JINC = 1
 402             MAIND = 1
 403          END IF
 404 *
 405          IF( TSCAL.NE.ONE ) THEN
 406             GROW = ZERO
 407             GO TO 60
 408          END IF
 409 *
 410          IF( NOUNIT ) THEN
 411 *
 412 *           A is non-unit triangular.
 413 *
 414 *           Compute GROW = 1/G(j) and XBND = 1/M(j).
 415 *           Initially, G(0) = max{x(i), i=1,...,n}.
 416 *
 417             GROW = HALF / MAX( XBND, SMLNUM )
 418             XBND = GROW
 419             DO 40 J = JFIRST, JLAST, JINC
 420 *
 421 *              Exit the loop if the growth factor is too small.
 422 *
 423                IF( GROW.LE.SMLNUM )
 424      $            GO TO 60
 425 *
 426                TJJS = AB( MAIND, J )
 427                TJJ = CABS1( TJJS )
 428 *
 429                IF( TJJ.GE.SMLNUM ) THEN
 430 *
 431 *                 M(j) = G(j-1) / abs(A(j,j))
 432 *
 433                   XBND = MIN( XBND, MIN( ONE, TJJ )*GROW )
 434                ELSE
 435 *
 436 *                 M(j) could overflow, set XBND to 0.
 437 *
 438                   XBND = ZERO
 439                END IF
 440 *
 441                IF( TJJ+CNORM( J ).GE.SMLNUM ) THEN
 442 *
 443 *                 G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
 444 *
 445                   GROW = GROW*( TJJ / ( TJJ+CNORM( J ) ) )
 446                ELSE
 447 *
 448 *                 G(j) could overflow, set GROW to 0.
 449 *
 450                   GROW = ZERO
 451                END IF
 452    40       CONTINUE
 453             GROW = XBND
 454          ELSE
 455 *
 456 *           A is unit triangular.
 457 *
 458 *           Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
 459 *
 460             GROW = MIN( ONE, HALF / MAX( XBND, SMLNUM ) )
 461             DO 50 J = JFIRST, JLAST, JINC
 462 *
 463 *              Exit the loop if the growth factor is too small.
 464 *
 465                IF( GROW.LE.SMLNUM )
 466      $            GO TO 60
 467 *
 468 *              G(j) = G(j-1)*( 1 + CNORM(j) )
 469 *
 470                GROW = GROW*( ONE / ( ONE+CNORM( J ) ) )
 471    50       CONTINUE
 472          END IF
 473    60    CONTINUE
 474 *
 475       ELSE
 476 *
 477 *        Compute the growth in A**T * x = b  or  A**H * x = b.
 478 *
 479          IF( UPPER ) THEN
 480             JFIRST = 1
 481             JLAST = N
 482             JINC = 1
 483             MAIND = KD + 1
 484          ELSE
 485             JFIRST = N
 486             JLAST = 1
 487             JINC = -1
 488             MAIND = 1
 489          END IF
 490 *
 491          IF( TSCAL.NE.ONE ) THEN
 492             GROW = ZERO
 493             GO TO 90
 494          END IF
 495 *
 496          IF( NOUNIT ) THEN
 497 *
 498 *           A is non-unit triangular.
 499 *
 500 *           Compute GROW = 1/G(j) and XBND = 1/M(j).
 501 *           Initially, M(0) = max{x(i), i=1,...,n}.
 502 *
 503             GROW = HALF / MAX( XBND, SMLNUM )
 504             XBND = GROW
 505             DO 70 J = JFIRST, JLAST, JINC
 506 *
 507 *              Exit the loop if the growth factor is too small.
 508 *
 509                IF( GROW.LE.SMLNUM )
 510      $            GO TO 90
 511 *
 512 *              G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
 513 *
 514                XJ = ONE + CNORM( J )
 515                GROW = MIN( GROW, XBND / XJ )
 516 *
 517                TJJS = AB( MAIND, J )
 518                TJJ = CABS1( TJJS )
 519 *
 520                IF( TJJ.GE.SMLNUM ) THEN
 521 *
 522 *                 M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
 523 *
 524                   IF( XJ.GT.TJJ )
 525      $               XBND = XBND*( TJJ / XJ )
 526                ELSE
 527 *
 528 *                 M(j) could overflow, set XBND to 0.
 529 *
 530                   XBND = ZERO
 531                END IF
 532    70       CONTINUE
 533             GROW = MIN( GROW, XBND )
 534          ELSE
 535 *
 536 *           A is unit triangular.
 537 *
 538 *           Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
 539 *
 540             GROW = MIN( ONE, HALF / MAX( XBND, SMLNUM ) )
 541             DO 80 J = JFIRST, JLAST, JINC
 542 *
 543 *              Exit the loop if the growth factor is too small.
 544 *
 545                IF( GROW.LE.SMLNUM )
 546      $            GO TO 90
 547 *
 548 *              G(j) = ( 1 + CNORM(j) )*G(j-1)
 549 *
 550                XJ = ONE + CNORM( J )
 551                GROW = GROW / XJ
 552    80       CONTINUE
 553          END IF
 554    90    CONTINUE
 555       END IF
 556 *
 557       IF( ( GROW*TSCAL ).GT.SMLNUM ) THEN
 558 *
 559 *        Use the Level 2 BLAS solve if the reciprocal of the bound on
 560 *        elements of X is not too small.
 561 *
 562          CALL ZTBSV( UPLO, TRANS, DIAG, N, KD, AB, LDAB, X, 1 )
 563       ELSE
 564 *
 565 *        Use a Level 1 BLAS solve, scaling intermediate results.
 566 *
 567          IF( XMAX.GT.BIGNUM*HALF ) THEN
 568 *
 569 *           Scale X so that its components are less than or equal to
 570 *           BIGNUM in absolute value.
 571 *
 572             SCALE = ( BIGNUM*HALF ) / XMAX
 573             CALL ZDSCAL( N, SCALE, X, 1 )
 574             XMAX = BIGNUM
 575          ELSE
 576             XMAX = XMAX*TWO
 577          END IF
 578 *
 579          IF( NOTRAN ) THEN
 580 *
 581 *           Solve A * x = b
 582 *
 583             DO 120 J = JFIRST, JLAST, JINC
 584 *
 585 *              Compute x(j) = b(j) / A(j,j), scaling x if necessary.
 586 *
 587                XJ = CABS1( X( J ) )
 588                IF( NOUNIT ) THEN
 589                   TJJS = AB( MAIND, J )*TSCAL
 590                ELSE
 591                   TJJS = TSCAL
 592                   IF( TSCAL.EQ.ONE )
 593      $               GO TO 110
 594                END IF
 595                TJJ = CABS1( TJJS )
 596                IF( TJJ.GT.SMLNUM ) THEN
 597 *
 598 *                    abs(A(j,j)) > SMLNUM:
 599 *
 600                   IF( TJJ.LT.ONE ) THEN
 601                      IF( XJ.GT.TJJ*BIGNUM ) THEN
 602 *
 603 *                          Scale x by 1/b(j).
 604 *
 605                         REC = ONE / XJ
 606                         CALL ZDSCAL( N, REC, X, 1 )
 607                         SCALE = SCALE*REC
 608                         XMAX = XMAX*REC
 609                      END IF
 610                   END IF
 611                   X( J ) = ZLADIV( X( J ), TJJS )
 612                   XJ = CABS1( X( J ) )
 613                ELSE IF( TJJ.GT.ZERO ) THEN
 614 *
 615 *                    0 < abs(A(j,j)) <= SMLNUM:
 616 *
 617                   IF( XJ.GT.TJJ*BIGNUM ) THEN
 618 *
 619 *                       Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
 620 *                       to avoid overflow when dividing by A(j,j).
 621 *
 622                      REC = ( TJJ*BIGNUM ) / XJ
 623                      IF( CNORM( J ).GT.ONE ) THEN
 624 *
 625 *                          Scale by 1/CNORM(j) to avoid overflow when
 626 *                          multiplying x(j) times column j.
 627 *
 628                         REC = REC / CNORM( J )
 629                      END IF
 630                      CALL ZDSCAL( N, REC, X, 1 )
 631                      SCALE = SCALE*REC
 632                      XMAX = XMAX*REC
 633                   END IF
 634                   X( J ) = ZLADIV( X( J ), TJJS )
 635                   XJ = CABS1( X( J ) )
 636                ELSE
 637 *
 638 *                    A(j,j) = 0:  Set x(1:n) = 0, x(j) = 1, and
 639 *                    scale = 0, and compute a solution to A*x = 0.
 640 *
 641                   DO 100 I = 1, N
 642                      X( I ) = ZERO
 643   100             CONTINUE
 644                   X( J ) = ONE
 645                   XJ = ONE
 646                   SCALE = ZERO
 647                   XMAX = ZERO
 648                END IF
 649   110          CONTINUE
 650 *
 651 *              Scale x if necessary to avoid overflow when adding a
 652 *              multiple of column j of A.
 653 *
 654                IF( XJ.GT.ONE ) THEN
 655                   REC = ONE / XJ
 656                   IF( CNORM( J ).GT.( BIGNUM-XMAX )*REC ) THEN
 657 *
 658 *                    Scale x by 1/(2*abs(x(j))).
 659 *
 660                      REC = REC*HALF
 661                      CALL ZDSCAL( N, REC, X, 1 )
 662                      SCALE = SCALE*REC
 663                   END IF
 664                ELSE IF( XJ*CNORM( J ).GT.( BIGNUM-XMAX ) ) THEN
 665 *
 666 *                 Scale x by 1/2.
 667 *
 668                   CALL ZDSCAL( N, HALF, X, 1 )
 669                   SCALE = SCALE*HALF
 670                END IF
 671 *
 672                IF( UPPER ) THEN
 673                   IF( J.GT.1 ) THEN
 674 *
 675 *                    Compute the update
 676 *                       x(max(1,j-kd):j-1) := x(max(1,j-kd):j-1) -
 677 *                                             x(j)* A(max(1,j-kd):j-1,j)
 678 *
 679                      JLEN = MIN( KD, J-1 )
 680                      CALL ZAXPY( JLEN, -X( J )*TSCAL,
 681      $                           AB( KD+1-JLEN, J ), 1, X( J-JLEN ), 1 )
 682                      I = IZAMAX( J-1, X, 1 )
 683                      XMAX = CABS1( X( I ) )
 684                   END IF
 685                ELSE IF( J.LT.N ) THEN
 686 *
 687 *                 Compute the update
 688 *                    x(j+1:min(j+kd,n)) := x(j+1:min(j+kd,n)) -
 689 *                                          x(j) * A(j+1:min(j+kd,n),j)
 690 *
 691                   JLEN = MIN( KD, N-J )
 692                   IF( JLEN.GT.0 )
 693      $               CALL ZAXPY( JLEN, -X( J )*TSCAL, AB( 2, J ), 1,
 694      $                           X( J+1 ), 1 )
 695                   I = J + IZAMAX( N-J, X( J+1 ), 1 )
 696                   XMAX = CABS1( X( I ) )
 697                END IF
 698   120       CONTINUE
 699 *
 700          ELSE IF( LSAME( TRANS, 'T' ) ) THEN
 701 *
 702 *           Solve A**T * x = b
 703 *
 704             DO 170 J = JFIRST, JLAST, JINC
 705 *
 706 *              Compute x(j) = b(j) - sum A(k,j)*x(k).
 707 *                                    k<>j
 708 *
 709                XJ = CABS1( X( J ) )
 710                USCAL = TSCAL
 711                REC = ONE / MAX( XMAX, ONE )
 712                IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN
 713 *
 714 *                 If x(j) could overflow, scale x by 1/(2*XMAX).
 715 *
 716                   REC = REC*HALF
 717                   IF( NOUNIT ) THEN
 718                      TJJS = AB( MAIND, J )*TSCAL
 719                   ELSE
 720                      TJJS = TSCAL
 721                   END IF
 722                   TJJ = CABS1( TJJS )
 723                   IF( TJJ.GT.ONE ) THEN
 724 *
 725 *                       Divide by A(j,j) when scaling x if A(j,j) > 1.
 726 *
 727                      REC = MIN( ONE, REC*TJJ )
 728                      USCAL = ZLADIV( USCAL, TJJS )
 729                   END IF
 730                   IF( REC.LT.ONE ) THEN
 731                      CALL ZDSCAL( N, REC, X, 1 )
 732                      SCALE = SCALE*REC
 733                      XMAX = XMAX*REC
 734                   END IF
 735                END IF
 736 *
 737                CSUMJ = ZERO
 738                IF( USCAL.EQ.DCMPLX( ONE ) ) THEN
 739 *
 740 *                 If the scaling needed for A in the dot product is 1,
 741 *                 call ZDOTU to perform the dot product.
 742 *
 743                   IF( UPPER ) THEN
 744                      JLEN = MIN( KD, J-1 )
 745                      CSUMJ = ZDOTU( JLEN, AB( KD+1-JLEN, J ), 1,
 746      $                       X( J-JLEN ), 1 )
 747                   ELSE
 748                      JLEN = MIN( KD, N-J )
 749                      IF( JLEN.GT.1 )
 750      $                  CSUMJ = ZDOTU( JLEN, AB( 2, J ), 1, X( J+1 ),
 751      $                          1 )
 752                   END IF
 753                ELSE
 754 *
 755 *                 Otherwise, use in-line code for the dot product.
 756 *
 757                   IF( UPPER ) THEN
 758                      JLEN = MIN( KD, J-1 )
 759                      DO 130 I = 1, JLEN
 760                         CSUMJ = CSUMJ + ( AB( KD+I-JLEN, J )*USCAL )*
 761      $                          X( J-JLEN-1+I )
 762   130                CONTINUE
 763                   ELSE
 764                      JLEN = MIN( KD, N-J )
 765                      DO 140 I = 1, JLEN
 766                         CSUMJ = CSUMJ + ( AB( I+1, J )*USCAL )*X( J+I )
 767   140                CONTINUE
 768                   END IF
 769                END IF
 770 *
 771                IF( USCAL.EQ.DCMPLX( TSCAL ) ) THEN
 772 *
 773 *                 Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
 774 *                 was not used to scale the dotproduct.
 775 *
 776                   X( J ) = X( J ) - CSUMJ
 777                   XJ = CABS1( X( J ) )
 778                   IF( NOUNIT ) THEN
 779 *
 780 *                    Compute x(j) = x(j) / A(j,j), scaling if necessary.
 781 *
 782                      TJJS = AB( MAIND, J )*TSCAL
 783                   ELSE
 784                      TJJS = TSCAL
 785                      IF( TSCAL.EQ.ONE )
 786      $                  GO TO 160
 787                   END IF
 788                   TJJ = CABS1( TJJS )
 789                   IF( TJJ.GT.SMLNUM ) THEN
 790 *
 791 *                       abs(A(j,j)) > SMLNUM:
 792 *
 793                      IF( TJJ.LT.ONE ) THEN
 794                         IF( XJ.GT.TJJ*BIGNUM ) THEN
 795 *
 796 *                             Scale X by 1/abs(x(j)).
 797 *
 798                            REC = ONE / XJ
 799                            CALL ZDSCAL( N, REC, X, 1 )
 800                            SCALE = SCALE*REC
 801                            XMAX = XMAX*REC
 802                         END IF
 803                      END IF
 804                      X( J ) = ZLADIV( X( J ), TJJS )
 805                   ELSE IF( TJJ.GT.ZERO ) THEN
 806 *
 807 *                       0 < abs(A(j,j)) <= SMLNUM:
 808 *
 809                      IF( XJ.GT.TJJ*BIGNUM ) THEN
 810 *
 811 *                          Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
 812 *
 813                         REC = ( TJJ*BIGNUM ) / XJ
 814                         CALL ZDSCAL( N, REC, X, 1 )
 815                         SCALE = SCALE*REC
 816                         XMAX = XMAX*REC
 817                      END IF
 818                      X( J ) = ZLADIV( X( J ), TJJS )
 819                   ELSE
 820 *
 821 *                       A(j,j) = 0:  Set x(1:n) = 0, x(j) = 1, and
 822 *                       scale = 0 and compute a solution to A**T *x = 0.
 823 *
 824                      DO 150 I = 1, N
 825                         X( I ) = ZERO
 826   150                CONTINUE
 827                      X( J ) = ONE
 828                      SCALE = ZERO
 829                      XMAX = ZERO
 830                   END IF
 831   160             CONTINUE
 832                ELSE
 833 *
 834 *                 Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
 835 *                 product has already been divided by 1/A(j,j).
 836 *
 837                   X( J ) = ZLADIV( X( J ), TJJS ) - CSUMJ
 838                END IF
 839                XMAX = MAX( XMAX, CABS1( X( J ) ) )
 840   170       CONTINUE
 841 *
 842          ELSE
 843 *
 844 *           Solve A**H * x = b
 845 *
 846             DO 220 J = JFIRST, JLAST, JINC
 847 *
 848 *              Compute x(j) = b(j) - sum A(k,j)*x(k).
 849 *                                    k<>j
 850 *
 851                XJ = CABS1( X( J ) )
 852                USCAL = TSCAL
 853                REC = ONE / MAX( XMAX, ONE )
 854                IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN
 855 *
 856 *                 If x(j) could overflow, scale x by 1/(2*XMAX).
 857 *
 858                   REC = REC*HALF
 859                   IF( NOUNIT ) THEN
 860                      TJJS = DCONJG( AB( MAIND, J ) )*TSCAL
 861                   ELSE
 862                      TJJS = TSCAL
 863                   END IF
 864                   TJJ = CABS1( TJJS )
 865                   IF( TJJ.GT.ONE ) THEN
 866 *
 867 *                       Divide by A(j,j) when scaling x if A(j,j) > 1.
 868 *
 869                      REC = MIN( ONE, REC*TJJ )
 870                      USCAL = ZLADIV( USCAL, TJJS )
 871                   END IF
 872                   IF( REC.LT.ONE ) THEN
 873                      CALL ZDSCAL( N, REC, X, 1 )
 874                      SCALE = SCALE*REC
 875                      XMAX = XMAX*REC
 876                   END IF
 877                END IF
 878 *
 879                CSUMJ = ZERO
 880                IF( USCAL.EQ.DCMPLX( ONE ) ) THEN
 881 *
 882 *                 If the scaling needed for A in the dot product is 1,
 883 *                 call ZDOTC to perform the dot product.
 884 *
 885                   IF( UPPER ) THEN
 886                      JLEN = MIN( KD, J-1 )
 887                      CSUMJ = ZDOTC( JLEN, AB( KD+1-JLEN, J ), 1,
 888      $                       X( J-JLEN ), 1 )
 889                   ELSE
 890                      JLEN = MIN( KD, N-J )
 891                      IF( JLEN.GT.1 )
 892      $                  CSUMJ = ZDOTC( JLEN, AB( 2, J ), 1, X( J+1 ),
 893      $                          1 )
 894                   END IF
 895                ELSE
 896 *
 897 *                 Otherwise, use in-line code for the dot product.
 898 *
 899                   IF( UPPER ) THEN
 900                      JLEN = MIN( KD, J-1 )
 901                      DO 180 I = 1, JLEN
 902                         CSUMJ = CSUMJ + ( DCONJG( AB( KD+I-JLEN, J ) )*
 903      $                          USCAL )*X( J-JLEN-1+I )
 904   180                CONTINUE
 905                   ELSE
 906                      JLEN = MIN( KD, N-J )
 907                      DO 190 I = 1, JLEN
 908                         CSUMJ = CSUMJ + ( DCONJG( AB( I+1, J ) )*USCAL )
 909      $                          *X( J+I )
 910   190                CONTINUE
 911                   END IF
 912                END IF
 913 *
 914                IF( USCAL.EQ.DCMPLX( TSCAL ) ) THEN
 915 *
 916 *                 Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
 917 *                 was not used to scale the dotproduct.
 918 *
 919                   X( J ) = X( J ) - CSUMJ
 920                   XJ = CABS1( X( J ) )
 921                   IF( NOUNIT ) THEN
 922 *
 923 *                    Compute x(j) = x(j) / A(j,j), scaling if necessary.
 924 *
 925                      TJJS = DCONJG( AB( MAIND, J ) )*TSCAL
 926                   ELSE
 927                      TJJS = TSCAL
 928                      IF( TSCAL.EQ.ONE )
 929      $                  GO TO 210
 930                   END IF
 931                   TJJ = CABS1( TJJS )
 932                   IF( TJJ.GT.SMLNUM ) THEN
 933 *
 934 *                       abs(A(j,j)) > SMLNUM:
 935 *
 936                      IF( TJJ.LT.ONE ) THEN
 937                         IF( XJ.GT.TJJ*BIGNUM ) THEN
 938 *
 939 *                             Scale X by 1/abs(x(j)).
 940 *
 941                            REC = ONE / XJ
 942                            CALL ZDSCAL( N, REC, X, 1 )
 943                            SCALE = SCALE*REC
 944                            XMAX = XMAX*REC
 945                         END IF
 946                      END IF
 947                      X( J ) = ZLADIV( X( J ), TJJS )
 948                   ELSE IF( TJJ.GT.ZERO ) THEN
 949 *
 950 *                       0 < abs(A(j,j)) <= SMLNUM:
 951 *
 952                      IF( XJ.GT.TJJ*BIGNUM ) THEN
 953 *
 954 *                          Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
 955 *
 956                         REC = ( TJJ*BIGNUM ) / XJ
 957                         CALL ZDSCAL( N, REC, X, 1 )
 958                         SCALE = SCALE*REC
 959                         XMAX = XMAX*REC
 960                      END IF
 961                      X( J ) = ZLADIV( X( J ), TJJS )
 962                   ELSE
 963 *
 964 *                       A(j,j) = 0:  Set x(1:n) = 0, x(j) = 1, and
 965 *                       scale = 0 and compute a solution to A**H *x = 0.
 966 *
 967                      DO 200 I = 1, N
 968                         X( I ) = ZERO
 969   200                CONTINUE
 970                      X( J ) = ONE
 971                      SCALE = ZERO
 972                      XMAX = ZERO
 973                   END IF
 974   210             CONTINUE
 975                ELSE
 976 *
 977 *                 Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
 978 *                 product has already been divided by 1/A(j,j).
 979 *
 980                   X( J ) = ZLADIV( X( J ), TJJS ) - CSUMJ
 981                END IF
 982                XMAX = MAX( XMAX, CABS1( X( J ) ) )
 983   220       CONTINUE
 984          END IF
 985          SCALE = SCALE / TSCAL
 986       END IF
 987 *
 988 *     Scale the column norms by 1/TSCAL for return.
 989 *
 990       IF( TSCAL.NE.ONE ) THEN
 991          CALL DSCAL( N, ONE / TSCAL, CNORM, 1 )
 992       END IF
 993 *
 994       RETURN
 995 *
 996 *     End of ZLATBS
 997 *
 998       END