LAPACK: SRC/slatbs.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X,
00002      $                   SCALE, CNORM, INFO )
00003 *
00004 *  -- LAPACK auxiliary routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       CHARACTER          DIAG, NORMIN, TRANS, UPLO
00011       INTEGER            INFO, KD, LDAB, N
00012       REAL               SCALE
00013 *     ..
00014 *     .. Array Arguments ..
00015       REAL               AB( LDAB, * ), CNORM( * ), X( * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  SLATBS solves one of the triangular systems
00022 *
00023 *     A *x = s*b  or  A'*x = s*b
00024 *
00025 *  with scaling to prevent overflow, where A is an upper or lower
00026 *  triangular band matrix.  Here A' denotes the transpose of A, x and b
00027 *  are n-element vectors, and s is a scaling factor, usually less than
00028 *  or equal to 1, chosen so that the components of x will be less than
00029 *  the overflow threshold.  If the unscaled problem will not cause
00030 *  overflow, the Level 2 BLAS routine STBSV is called.  If the matrix A
00031 *  is singular (A(j,j) = 0 for some j), then s is set to 0 and a
00032 *  non-trivial solution to A*x = 0 is returned.
00033 *
00034 *  Arguments
00035 *  =========
00036 *
00037 *  UPLO    (input) CHARACTER*1
00038 *          Specifies whether the matrix A is upper or lower triangular.
00039 *          = 'U':  Upper triangular
00040 *          = 'L':  Lower triangular
00041 *
00042 *  TRANS   (input) CHARACTER*1
00043 *          Specifies the operation applied to A.
00044 *          = 'N':  Solve A * x = s*b  (No transpose)
00045 *          = 'T':  Solve A'* x = s*b  (Transpose)
00046 *          = 'C':  Solve A'* x = s*b  (Conjugate transpose = Transpose)
00047 *
00048 *  DIAG    (input) CHARACTER*1
00049 *          Specifies whether or not the matrix A is unit triangular.
00050 *          = 'N':  Non-unit triangular
00051 *          = 'U':  Unit triangular
00052 *
00053 *  NORMIN  (input) CHARACTER*1
00054 *          Specifies whether CNORM has been set or not.
00055 *          = 'Y':  CNORM contains the column norms on entry
00056 *          = 'N':  CNORM is not set on entry.  On exit, the norms will
00057 *                  be computed and stored in CNORM.
00058 *
00059 *  N       (input) INTEGER
00060 *          The order of the matrix A.  N >= 0.
00061 *
00062 *  KD      (input) INTEGER
00063 *          The number of subdiagonals or superdiagonals in the
00064 *          triangular matrix A.  KD >= 0.
00065 *
00066 *  AB      (input) REAL array, dimension (LDAB,N)
00067 *          The upper or lower triangular band matrix A, stored in the
00068 *          first KD+1 rows of the array. The j-th column of A is stored
00069 *          in the j-th column of the array AB as follows:
00070 *          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
00071 *          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
00072 *
00073 *  LDAB    (input) INTEGER
00074 *          The leading dimension of the array AB.  LDAB >= KD+1.
00075 *
00076 *  X       (input/output) REAL array, dimension (N)
00077 *          On entry, the right hand side b of the triangular system.
00078 *          On exit, X is overwritten by the solution vector x.
00079 *
00080 *  SCALE   (output) REAL
00081 *          The scaling factor s for the triangular system
00082 *             A * x = s*b  or  A'* x = s*b.
00083 *          If SCALE = 0, the matrix A is singular or badly scaled, and
00084 *          the vector x is an exact or approximate solution to A*x = 0.
00085 *
00086 *  CNORM   (input or output) REAL array, dimension (N)
00087 *
00088 *          If NORMIN = 'Y', CNORM is an input argument and CNORM(j)
00089 *          contains the norm of the off-diagonal part of the j-th column
00090 *          of A.  If TRANS = 'N', CNORM(j) must be greater than or equal
00091 *          to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)
00092 *          must be greater than or equal to the 1-norm.
00093 *
00094 *          If NORMIN = 'N', CNORM is an output argument and CNORM(j)
00095 *          returns the 1-norm of the offdiagonal part of the j-th column
00096 *          of A.
00097 *
00098 *  INFO    (output) INTEGER
00099 *          = 0:  successful exit
00100 *          < 0:  if INFO = -k, the k-th argument had an illegal value
00101 *
00102 *  Further Details
00103 *  ======= =======
00104 *
00105 *  A rough bound on x is computed; if that is less than overflow, STBSV
00106 *  is called, otherwise, specific code is used which checks for possible
00107 *  overflow or divide-by-zero at every operation.
00108 *
00109 *  A columnwise scheme is used for solving A*x = b.  The basic algorithm
00110 *  if A is lower triangular is
00111 *
00112 *       x[1:n] := b[1:n]
00113 *       for j = 1, ..., n
00114 *            x(j) := x(j) / A(j,j)
00115 *            x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
00116 *       end
00117 *
00118 *  Define bounds on the components of x after j iterations of the loop:
00119 *     M(j) = bound on x[1:j]
00120 *     G(j) = bound on x[j+1:n]
00121 *  Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
00122 *
00123 *  Then for iteration j+1 we have
00124 *     M(j+1) <= G(j) / | A(j+1,j+1) |
00125 *     G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
00126 *            <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
00127 *
00128 *  where CNORM(j+1) is greater than or equal to the infinity-norm of
00129 *  column j+1 of A, not counting the diagonal.  Hence
00130 *
00131 *     G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
00132 *                  1<=i<=j
00133 *  and
00134 *
00135 *     |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
00136 *                                   1<=i< j
00137 *
00138 *  Since |x(j)| <= M(j), we use the Level 2 BLAS routine STBSV if the
00139 *  reciprocal of the largest M(j), j=1,..,n, is larger than
00140 *  max(underflow, 1/overflow).
00141 *
00142 *  The bound on x(j) is also used to determine when a step in the
00143 *  columnwise method can be performed without fear of overflow.  If
00144 *  the computed bound is greater than a large constant, x is scaled to
00145 *  prevent overflow, but if the bound overflows, x is set to 0, x(j) to
00146 *  1, and scale to 0, and a non-trivial solution to A*x = 0 is found.
00147 *
00148 *  Similarly, a row-wise scheme is used to solve A'*x = b.  The basic
00149 *  algorithm for A upper triangular is
00150 *
00151 *       for j = 1, ..., n
00152 *            x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)
00153 *       end
00154 *
00155 *  We simultaneously compute two bounds
00156 *       G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j
00157 *       M(j) = bound on x(i), 1<=i<=j
00158 *
00159 *  The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we
00160 *  add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.
00161 *  Then the bound on x(j) is
00162 *
00163 *       M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
00164 *
00165 *            <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
00166 *                      1<=i<=j
00167 *
00168 *  and we can safely call STBSV if 1/M(n) and 1/G(n) are both greater
00169 *  than max(underflow, 1/overflow).
00170 *
00171 *  =====================================================================
00172 *
00173 *     .. Parameters ..
00174       REAL               ZERO, HALF, ONE
00175       PARAMETER          ( ZERO = 0.0E+0, HALF = 0.5E+0, ONE = 1.0E+0 )
00176 *     ..
00177 *     .. Local Scalars ..
00178       LOGICAL            NOTRAN, NOUNIT, UPPER
00179       INTEGER            I, IMAX, J, JFIRST, JINC, JLAST, JLEN, MAIND
00180       REAL               BIGNUM, GROW, REC, SMLNUM, SUMJ, TJJ, TJJS,
00181      $                   TMAX, TSCAL, USCAL, XBND, XJ, XMAX
00182 *     ..
00183 *     .. External Functions ..
00184       LOGICAL            LSAME
00185       INTEGER            ISAMAX
00186       REAL               SASUM, SDOT, SLAMCH
00187       EXTERNAL           LSAME, ISAMAX, SASUM, SDOT, SLAMCH
00188 *     ..
00189 *     .. External Subroutines ..
00190       EXTERNAL           SAXPY, SSCAL, STBSV, XERBLA
00191 *     ..
00192 *     .. Intrinsic Functions ..
00193       INTRINSIC          ABS, MAX, MIN
00194 *     ..
00195 *     .. Executable Statements ..
00196 *
00197       INFO = 0
00198       UPPER = LSAME( UPLO, 'U' )
00199       NOTRAN = LSAME( TRANS, 'N' )
00200       NOUNIT = LSAME( DIAG, 'N' )
00201 *
00202 *     Test the input parameters.
00203 *
00204       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00205          INFO = -1
00206       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
00207      $         LSAME( TRANS, 'C' ) ) THEN
00208          INFO = -2
00209       ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
00210          INFO = -3
00211       ELSE IF( .NOT.LSAME( NORMIN, 'Y' ) .AND. .NOT.
00212      $         LSAME( NORMIN, 'N' ) ) THEN
00213          INFO = -4
00214       ELSE IF( N.LT.0 ) THEN
00215          INFO = -5
00216       ELSE IF( KD.LT.0 ) THEN
00217          INFO = -6
00218       ELSE IF( LDAB.LT.KD+1 ) THEN
00219          INFO = -8
00220       END IF
00221       IF( INFO.NE.0 ) THEN
00222          CALL XERBLA( 'SLATBS', -INFO )
00223          RETURN
00224       END IF
00225 *
00226 *     Quick return if possible
00227 *
00228       IF( N.EQ.0 )
00229      $   RETURN
00230 *
00231 *     Determine machine dependent parameters to control overflow.
00232 *
00233       SMLNUM = SLAMCH( 'Safe minimum' ) / SLAMCH( 'Precision' )
00234       BIGNUM = ONE / SMLNUM
00235       SCALE = ONE
00236 *
00237       IF( LSAME( NORMIN, 'N' ) ) THEN
00238 *
00239 *        Compute the 1-norm of each column, not including the diagonal.
00240 *
00241          IF( UPPER ) THEN
00242 *
00243 *           A is upper triangular.
00244 *
00245             DO 10 J = 1, N
00246                JLEN = MIN( KD, J-1 )
00247                CNORM( J ) = SASUM( JLEN, AB( KD+1-JLEN, J ), 1 )
00248    10       CONTINUE
00249          ELSE
00250 *
00251 *           A is lower triangular.
00252 *
00253             DO 20 J = 1, N
00254                JLEN = MIN( KD, N-J )
00255                IF( JLEN.GT.0 ) THEN
00256                   CNORM( J ) = SASUM( JLEN, AB( 2, J ), 1 )
00257                ELSE
00258                   CNORM( J ) = ZERO
00259                END IF
00260    20       CONTINUE
00261          END IF
00262       END IF
00263 *
00264 *     Scale the column norms by TSCAL if the maximum element in CNORM is
00265 *     greater than BIGNUM.
00266 *
00267       IMAX = ISAMAX( N, CNORM, 1 )
00268       TMAX = CNORM( IMAX )
00269       IF( TMAX.LE.BIGNUM ) THEN
00270          TSCAL = ONE
00271       ELSE
00272          TSCAL = ONE / ( SMLNUM*TMAX )
00273          CALL SSCAL( N, TSCAL, CNORM, 1 )
00274       END IF
00275 *
00276 *     Compute a bound on the computed solution vector to see if the
00277 *     Level 2 BLAS routine STBSV can be used.
00278 *
00279       J = ISAMAX( N, X, 1 )
00280       XMAX = ABS( X( J ) )
00281       XBND = XMAX
00282       IF( NOTRAN ) THEN
00283 *
00284 *        Compute the growth in A * x = b.
00285 *
00286          IF( UPPER ) THEN
00287             JFIRST = N
00288             JLAST = 1
00289             JINC = -1
00290             MAIND = KD + 1
00291          ELSE
00292             JFIRST = 1
00293             JLAST = N
00294             JINC = 1
00295             MAIND = 1
00296          END IF
00297 *
00298          IF( TSCAL.NE.ONE ) THEN
00299             GROW = ZERO
00300             GO TO 50
00301          END IF
00302 *
00303          IF( NOUNIT ) THEN
00304 *
00305 *           A is non-unit triangular.
00306 *
00307 *           Compute GROW = 1/G(j) and XBND = 1/M(j).
00308 *           Initially, G(0) = max{x(i), i=1,...,n}.
00309 *
00310             GROW = ONE / MAX( XBND, SMLNUM )
00311             XBND = GROW
00312             DO 30 J = JFIRST, JLAST, JINC
00313 *
00314 *              Exit the loop if the growth factor is too small.
00315 *
00316                IF( GROW.LE.SMLNUM )
00317      $            GO TO 50
00318 *
00319 *              M(j) = G(j-1) / abs(A(j,j))
00320 *
00321                TJJ = ABS( AB( MAIND, J ) )
00322                XBND = MIN( XBND, MIN( ONE, TJJ )*GROW )
00323                IF( TJJ+CNORM( J ).GE.SMLNUM ) THEN
00324 *
00325 *                 G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
00326 *
00327                   GROW = GROW*( TJJ / ( TJJ+CNORM( J ) ) )
00328                ELSE
00329 *
00330 *                 G(j) could overflow, set GROW to 0.
00331 *
00332                   GROW = ZERO
00333                END IF
00334    30       CONTINUE
00335             GROW = XBND
00336          ELSE
00337 *
00338 *           A is unit triangular.
00339 *
00340 *           Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
00341 *
00342             GROW = MIN( ONE, ONE / MAX( XBND, SMLNUM ) )
00343             DO 40 J = JFIRST, JLAST, JINC
00344 *
00345 *              Exit the loop if the growth factor is too small.
00346 *
00347                IF( GROW.LE.SMLNUM )
00348      $            GO TO 50
00349 *
00350 *              G(j) = G(j-1)*( 1 + CNORM(j) )
00351 *
00352                GROW = GROW*( ONE / ( ONE+CNORM( J ) ) )
00353    40       CONTINUE
00354          END IF
00355    50    CONTINUE
00356 *
00357       ELSE
00358 *
00359 *        Compute the growth in A' * x = b.
00360 *
00361          IF( UPPER ) THEN
00362             JFIRST = 1
00363             JLAST = N
00364             JINC = 1
00365             MAIND = KD + 1
00366          ELSE
00367             JFIRST = N
00368             JLAST = 1
00369             JINC = -1
00370             MAIND = 1
00371          END IF
00372 *
00373          IF( TSCAL.NE.ONE ) THEN
00374             GROW = ZERO
00375             GO TO 80
00376          END IF
00377 *
00378          IF( NOUNIT ) THEN
00379 *
00380 *           A is non-unit triangular.
00381 *
00382 *           Compute GROW = 1/G(j) and XBND = 1/M(j).
00383 *           Initially, M(0) = max{x(i), i=1,...,n}.
00384 *
00385             GROW = ONE / MAX( XBND, SMLNUM )
00386             XBND = GROW
00387             DO 60 J = JFIRST, JLAST, JINC
00388 *
00389 *              Exit the loop if the growth factor is too small.
00390 *
00391                IF( GROW.LE.SMLNUM )
00392      $            GO TO 80
00393 *
00394 *              G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
00395 *
00396                XJ = ONE + CNORM( J )
00397                GROW = MIN( GROW, XBND / XJ )
00398 *
00399 *              M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
00400 *
00401                TJJ = ABS( AB( MAIND, J ) )
00402                IF( XJ.GT.TJJ )
00403      $            XBND = XBND*( TJJ / XJ )
00404    60       CONTINUE
00405             GROW = MIN( GROW, XBND )
00406          ELSE
00407 *
00408 *           A is unit triangular.
00409 *
00410 *           Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
00411 *
00412             GROW = MIN( ONE, ONE / MAX( XBND, SMLNUM ) )
00413             DO 70 J = JFIRST, JLAST, JINC
00414 *
00415 *              Exit the loop if the growth factor is too small.
00416 *
00417                IF( GROW.LE.SMLNUM )
00418      $            GO TO 80
00419 *
00420 *              G(j) = ( 1 + CNORM(j) )*G(j-1)
00421 *
00422                XJ = ONE + CNORM( J )
00423                GROW = GROW / XJ
00424    70       CONTINUE
00425          END IF
00426    80    CONTINUE
00427       END IF
00428 *
00429       IF( ( GROW*TSCAL ).GT.SMLNUM ) THEN
00430 *
00431 *        Use the Level 2 BLAS solve if the reciprocal of the bound on
00432 *        elements of X is not too small.
00433 *
00434          CALL STBSV( UPLO, TRANS, DIAG, N, KD, AB, LDAB, X, 1 )
00435       ELSE
00436 *
00437 *        Use a Level 1 BLAS solve, scaling intermediate results.
00438 *
00439          IF( XMAX.GT.BIGNUM ) THEN
00440 *
00441 *           Scale X so that its components are less than or equal to
00442 *           BIGNUM in absolute value.
00443 *
00444             SCALE = BIGNUM / XMAX
00445             CALL SSCAL( N, SCALE, X, 1 )
00446             XMAX = BIGNUM
00447          END IF
00448 *
00449          IF( NOTRAN ) THEN
00450 *
00451 *           Solve A * x = b
00452 *
00453             DO 100 J = JFIRST, JLAST, JINC
00454 *
00455 *              Compute x(j) = b(j) / A(j,j), scaling x if necessary.
00456 *
00457                XJ = ABS( X( J ) )
00458                IF( NOUNIT ) THEN
00459                   TJJS = AB( MAIND, J )*TSCAL
00460                ELSE
00461                   TJJS = TSCAL
00462                   IF( TSCAL.EQ.ONE )
00463      $               GO TO 95
00464                END IF
00465                   TJJ = ABS( TJJS )
00466                   IF( TJJ.GT.SMLNUM ) THEN
00467 *
00468 *                    abs(A(j,j)) > SMLNUM:
00469 *
00470                      IF( TJJ.LT.ONE ) THEN
00471                         IF( XJ.GT.TJJ*BIGNUM ) THEN
00472 *
00473 *                          Scale x by 1/b(j).
00474 *
00475                            REC = ONE / XJ
00476                            CALL SSCAL( N, REC, X, 1 )
00477                            SCALE = SCALE*REC
00478                            XMAX = XMAX*REC
00479                         END IF
00480                      END IF
00481                      X( J ) = X( J ) / TJJS
00482                      XJ = ABS( X( J ) )
00483                   ELSE IF( TJJ.GT.ZERO ) THEN
00484 *
00485 *                    0 < abs(A(j,j)) <= SMLNUM:
00486 *
00487                      IF( XJ.GT.TJJ*BIGNUM ) THEN
00488 *
00489 *                       Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
00490 *                       to avoid overflow when dividing by A(j,j).
00491 *
00492                         REC = ( TJJ*BIGNUM ) / XJ
00493                         IF( CNORM( J ).GT.ONE ) THEN
00494 *
00495 *                          Scale by 1/CNORM(j) to avoid overflow when
00496 *                          multiplying x(j) times column j.
00497 *
00498                            REC = REC / CNORM( J )
00499                         END IF
00500                         CALL SSCAL( N, REC, X, 1 )
00501                         SCALE = SCALE*REC
00502                         XMAX = XMAX*REC
00503                      END IF
00504                      X( J ) = X( J ) / TJJS
00505                      XJ = ABS( X( J ) )
00506                   ELSE
00507 *
00508 *                    A(j,j) = 0:  Set x(1:n) = 0, x(j) = 1, and
00509 *                    scale = 0, and compute a solution to A*x = 0.
00510 *
00511                      DO 90 I = 1, N
00512                         X( I ) = ZERO
00513    90                CONTINUE
00514                      X( J ) = ONE
00515                      XJ = ONE
00516                      SCALE = ZERO
00517                      XMAX = ZERO
00518                   END IF
00519    95          CONTINUE
00520 *
00521 *              Scale x if necessary to avoid overflow when adding a
00522 *              multiple of column j of A.
00523 *
00524                IF( XJ.GT.ONE ) THEN
00525                   REC = ONE / XJ
00526                   IF( CNORM( J ).GT.( BIGNUM-XMAX )*REC ) THEN
00527 *
00528 *                    Scale x by 1/(2*abs(x(j))).
00529 *
00530                      REC = REC*HALF
00531                      CALL SSCAL( N, REC, X, 1 )
00532                      SCALE = SCALE*REC
00533                   END IF
00534                ELSE IF( XJ*CNORM( J ).GT.( BIGNUM-XMAX ) ) THEN
00535 *
00536 *                 Scale x by 1/2.
00537 *
00538                   CALL SSCAL( N, HALF, X, 1 )
00539                   SCALE = SCALE*HALF
00540                END IF
00541 *
00542                IF( UPPER ) THEN
00543                   IF( J.GT.1 ) THEN
00544 *
00545 *                    Compute the update
00546 *                       x(max(1,j-kd):j-1) := x(max(1,j-kd):j-1) -
00547 *                                             x(j)* A(max(1,j-kd):j-1,j)
00548 *
00549                      JLEN = MIN( KD, J-1 )
00550                      CALL SAXPY( JLEN, -X( J )*TSCAL,
00551      $                           AB( KD+1-JLEN, J ), 1, X( J-JLEN ), 1 )
00552                      I = ISAMAX( J-1, X, 1 )
00553                      XMAX = ABS( X( I ) )
00554                   END IF
00555                ELSE IF( J.LT.N ) THEN
00556 *
00557 *                 Compute the update
00558 *                    x(j+1:min(j+kd,n)) := x(j+1:min(j+kd,n)) -
00559 *                                          x(j) * A(j+1:min(j+kd,n),j)
00560 *
00561                   JLEN = MIN( KD, N-J )
00562                   IF( JLEN.GT.0 )
00563      $               CALL SAXPY( JLEN, -X( J )*TSCAL, AB( 2, J ), 1,
00564      $                           X( J+1 ), 1 )
00565                   I = J + ISAMAX( N-J, X( J+1 ), 1 )
00566                   XMAX = ABS( X( I ) )
00567                END IF
00568   100       CONTINUE
00569 *
00570          ELSE
00571 *
00572 *           Solve A' * x = b
00573 *
00574             DO 140 J = JFIRST, JLAST, JINC
00575 *
00576 *              Compute x(j) = b(j) - sum A(k,j)*x(k).
00577 *                                    k<>j
00578 *
00579                XJ = ABS( X( J ) )
00580                USCAL = TSCAL
00581                REC = ONE / MAX( XMAX, ONE )
00582                IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN
00583 *
00584 *                 If x(j) could overflow, scale x by 1/(2*XMAX).
00585 *
00586                   REC = REC*HALF
00587                   IF( NOUNIT ) THEN
00588                      TJJS = AB( MAIND, J )*TSCAL
00589                   ELSE
00590                      TJJS = TSCAL
00591                   END IF
00592                      TJJ = ABS( TJJS )
00593                      IF( TJJ.GT.ONE ) THEN
00594 *
00595 *                       Divide by A(j,j) when scaling x if A(j,j) > 1.
00596 *
00597                         REC = MIN( ONE, REC*TJJ )
00598                         USCAL = USCAL / TJJS
00599                      END IF
00600                   IF( REC.LT.ONE ) THEN
00601                      CALL SSCAL( N, REC, X, 1 )
00602                      SCALE = SCALE*REC
00603                      XMAX = XMAX*REC
00604                   END IF
00605                END IF
00606 *
00607                SUMJ = ZERO
00608                IF( USCAL.EQ.ONE ) THEN
00609 *
00610 *                 If the scaling needed for A in the dot product is 1,
00611 *                 call SDOT to perform the dot product.
00612 *
00613                   IF( UPPER ) THEN
00614                      JLEN = MIN( KD, J-1 )
00615                      SUMJ = SDOT( JLEN, AB( KD+1-JLEN, J ), 1,
00616      $                      X( J-JLEN ), 1 )
00617                   ELSE
00618                      JLEN = MIN( KD, N-J )
00619                      IF( JLEN.GT.0 )
00620      $                  SUMJ = SDOT( JLEN, AB( 2, J ), 1, X( J+1 ), 1 )
00621                   END IF
00622                ELSE
00623 *
00624 *                 Otherwise, use in-line code for the dot product.
00625 *
00626                   IF( UPPER ) THEN
00627                      JLEN = MIN( KD, J-1 )
00628                      DO 110 I = 1, JLEN
00629                         SUMJ = SUMJ + ( AB( KD+I-JLEN, J )*USCAL )*
00630      $                         X( J-JLEN-1+I )
00631   110                CONTINUE
00632                   ELSE
00633                      JLEN = MIN( KD, N-J )
00634                      DO 120 I = 1, JLEN
00635                         SUMJ = SUMJ + ( AB( I+1, J )*USCAL )*X( J+I )
00636   120                CONTINUE
00637                   END IF
00638                END IF
00639 *
00640                IF( USCAL.EQ.TSCAL ) THEN
00641 *
00642 *                 Compute x(j) := ( x(j) - sumj ) / A(j,j) if 1/A(j,j)
00643 *                 was not used to scale the dotproduct.
00644 *
00645                   X( J ) = X( J ) - SUMJ
00646                   XJ = ABS( X( J ) )
00647                   IF( NOUNIT ) THEN
00648 *
00649 *                    Compute x(j) = x(j) / A(j,j), scaling if necessary.
00650 *
00651                      TJJS = AB( MAIND, J )*TSCAL
00652                   ELSE
00653                      TJJS = TSCAL
00654                      IF( TSCAL.EQ.ONE )
00655      $                  GO TO 135
00656                   END IF
00657                      TJJ = ABS( TJJS )
00658                      IF( TJJ.GT.SMLNUM ) THEN
00659 *
00660 *                       abs(A(j,j)) > SMLNUM:
00661 *
00662                         IF( TJJ.LT.ONE ) THEN
00663                            IF( XJ.GT.TJJ*BIGNUM ) THEN
00664 *
00665 *                             Scale X by 1/abs(x(j)).
00666 *
00667                               REC = ONE / XJ
00668                               CALL SSCAL( N, REC, X, 1 )
00669                               SCALE = SCALE*REC
00670                               XMAX = XMAX*REC
00671                            END IF
00672                         END IF
00673                         X( J ) = X( J ) / TJJS
00674                      ELSE IF( TJJ.GT.ZERO ) THEN
00675 *
00676 *                       0 < abs(A(j,j)) <= SMLNUM:
00677 *
00678                         IF( XJ.GT.TJJ*BIGNUM ) THEN
00679 *
00680 *                          Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
00681 *
00682                            REC = ( TJJ*BIGNUM ) / XJ
00683                            CALL SSCAL( N, REC, X, 1 )
00684                            SCALE = SCALE*REC
00685                            XMAX = XMAX*REC
00686                         END IF
00687                         X( J ) = X( J ) / TJJS
00688                      ELSE
00689 *
00690 *                       A(j,j) = 0:  Set x(1:n) = 0, x(j) = 1, and
00691 *                       scale = 0, and compute a solution to A'*x = 0.
00692 *
00693                         DO 130 I = 1, N
00694                            X( I ) = ZERO
00695   130                   CONTINUE
00696                         X( J ) = ONE
00697                         SCALE = ZERO
00698                         XMAX = ZERO
00699                      END IF
00700   135             CONTINUE
00701                ELSE
00702 *
00703 *                 Compute x(j) := x(j) / A(j,j) - sumj if the dot
00704 *                 product has already been divided by 1/A(j,j).
00705 *
00706                   X( J ) = X( J ) / TJJS - SUMJ
00707                END IF
00708                XMAX = MAX( XMAX, ABS( X( J ) ) )
00709   140       CONTINUE
00710          END IF
00711          SCALE = SCALE / TSCAL
00712       END IF
00713 *
00714 *     Scale the column norms by 1/TSCAL for return.
00715 *
00716       IF( TSCAL.NE.ONE ) THEN
00717          CALL SSCAL( N, ONE / TSCAL, CNORM, 1 )
00718       END IF
00719 *
00720       RETURN
00721 *
00722 *     End of SLATBS
00723 *
00724       END