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