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