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