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