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