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