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