SRC/dlatrs.f

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