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