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