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