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