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