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