![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
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