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