1*> \brief \b CLATBS solves a triangular banded system of equations. 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download CLATBS + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clatbs.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clatbs.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clatbs.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE CLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X, 22* SCALE, CNORM, INFO ) 23* 24* .. Scalar Arguments .. 25* CHARACTER DIAG, NORMIN, TRANS, UPLO 26* INTEGER INFO, KD, LDAB, N 27* REAL SCALE 28* .. 29* .. Array Arguments .. 30* REAL CNORM( * ) 31* COMPLEX AB( LDAB, * ), X( * ) 32* .. 33* 34* 35*> \par Purpose: 36* ============= 37*> 38*> \verbatim 39*> 40*> CLATBS solves one of the triangular systems 41*> 42*> A * x = s*b, A**T * x = s*b, or A**H * x = s*b, 43*> 44*> with scaling to prevent overflow, where A is an upper or lower 45*> triangular band matrix. Here A**T denotes the transpose of A, x and b 46*> are n-element vectors, and s is a scaling factor, usually less than 47*> or equal to 1, chosen so that the components of x will be less than 48*> the overflow threshold. If the unscaled problem will not cause 49*> overflow, the Level 2 BLAS routine CTBSV is called. If the matrix A 50*> is singular (A(j,j) = 0 for some j), then s is set to 0 and a 51*> non-trivial solution to A*x = 0 is returned. 52*> \endverbatim 53* 54* Arguments: 55* ========== 56* 57*> \param[in] UPLO 58*> \verbatim 59*> UPLO is CHARACTER*1 60*> Specifies whether the matrix A is upper or lower triangular. 61*> = 'U': Upper triangular 62*> = 'L': Lower triangular 63*> \endverbatim 64*> 65*> \param[in] TRANS 66*> \verbatim 67*> TRANS is CHARACTER*1 68*> Specifies the operation applied to A. 69*> = 'N': Solve A * x = s*b (No transpose) 70*> = 'T': Solve A**T * x = s*b (Transpose) 71*> = 'C': Solve A**H * x = s*b (Conjugate transpose) 72*> \endverbatim 73*> 74*> \param[in] DIAG 75*> \verbatim 76*> DIAG is CHARACTER*1 77*> Specifies whether or not the matrix A is unit triangular. 78*> = 'N': Non-unit triangular 79*> = 'U': Unit triangular 80*> \endverbatim 81*> 82*> \param[in] NORMIN 83*> \verbatim 84*> NORMIN is CHARACTER*1 85*> Specifies whether CNORM has been set or not. 86*> = 'Y': CNORM contains the column norms on entry 87*> = 'N': CNORM is not set on entry. On exit, the norms will 88*> be computed and stored in CNORM. 89*> \endverbatim 90*> 91*> \param[in] N 92*> \verbatim 93*> N is INTEGER 94*> The order of the matrix A. N >= 0. 95*> \endverbatim 96*> 97*> \param[in] KD 98*> \verbatim 99*> KD is INTEGER 100*> The number of subdiagonals or superdiagonals in the 101*> triangular matrix A. KD >= 0. 102*> \endverbatim 103*> 104*> \param[in] AB 105*> \verbatim 106*> AB is COMPLEX array, dimension (LDAB,N) 107*> The upper or lower triangular band matrix A, stored in the 108*> first KD+1 rows of the array. The j-th column of A is stored 109*> in the j-th column of the array AB as follows: 110*> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; 111*> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). 112*> \endverbatim 113*> 114*> \param[in] LDAB 115*> \verbatim 116*> LDAB is INTEGER 117*> The leading dimension of the array AB. LDAB >= KD+1. 118*> \endverbatim 119*> 120*> \param[in,out] X 121*> \verbatim 122*> X is COMPLEX array, dimension (N) 123*> On entry, the right hand side b of the triangular system. 124*> On exit, X is overwritten by the solution vector x. 125*> \endverbatim 126*> 127*> \param[out] SCALE 128*> \verbatim 129*> SCALE is REAL 130*> The scaling factor s for the triangular system 131*> A * x = s*b, A**T * x = s*b, or A**H * x = s*b. 132*> If SCALE = 0, the matrix A is singular or badly scaled, and 133*> the vector x is an exact or approximate solution to A*x = 0. 134*> \endverbatim 135*> 136*> \param[in,out] CNORM 137*> \verbatim 138*> CNORM is REAL array, dimension (N) 139*> 140*> If NORMIN = 'Y', CNORM is an input argument and CNORM(j) 141*> contains the norm of the off-diagonal part of the j-th column 142*> of A. If TRANS = 'N', CNORM(j) must be greater than or equal 143*> to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) 144*> must be greater than or equal to the 1-norm. 145*> 146*> If NORMIN = 'N', CNORM is an output argument and CNORM(j) 147*> returns the 1-norm of the offdiagonal part of the j-th column 148*> of A. 149*> \endverbatim 150*> 151*> \param[out] INFO 152*> \verbatim 153*> INFO is INTEGER 154*> = 0: successful exit 155*> < 0: if INFO = -k, the k-th argument had an illegal value 156*> \endverbatim 157* 158* Authors: 159* ======== 160* 161*> \author Univ. of Tennessee 162*> \author Univ. of California Berkeley 163*> \author Univ. of Colorado Denver 164*> \author NAG Ltd. 165* 166*> \ingroup complexOTHERauxiliary 167* 168*> \par Further Details: 169* ===================== 170*> 171*> \verbatim 172*> 173*> A rough bound on x is computed; if that is less than overflow, CTBSV 174*> is called, otherwise, specific code is used which checks for possible 175*> overflow or divide-by-zero at every operation. 176*> 177*> A columnwise scheme is used for solving A*x = b. The basic algorithm 178*> if A is lower triangular is 179*> 180*> x[1:n] := b[1:n] 181*> for j = 1, ..., n 182*> x(j) := x(j) / A(j,j) 183*> x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j] 184*> end 185*> 186*> Define bounds on the components of x after j iterations of the loop: 187*> M(j) = bound on x[1:j] 188*> G(j) = bound on x[j+1:n] 189*> Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}. 190*> 191*> Then for iteration j+1 we have 192*> M(j+1) <= G(j) / | A(j+1,j+1) | 193*> G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] | 194*> <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | ) 195*> 196*> where CNORM(j+1) is greater than or equal to the infinity-norm of 197*> column j+1 of A, not counting the diagonal. Hence 198*> 199*> G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | ) 200*> 1<=i<=j 201*> and 202*> 203*> |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| ) 204*> 1<=i< j 205*> 206*> Since |x(j)| <= M(j), we use the Level 2 BLAS routine CTBSV if the 207*> reciprocal of the largest M(j), j=1,..,n, is larger than 208*> max(underflow, 1/overflow). 209*> 210*> The bound on x(j) is also used to determine when a step in the 211*> columnwise method can be performed without fear of overflow. If 212*> the computed bound is greater than a large constant, x is scaled to 213*> prevent overflow, but if the bound overflows, x is set to 0, x(j) to 214*> 1, and scale to 0, and a non-trivial solution to A*x = 0 is found. 215*> 216*> Similarly, a row-wise scheme is used to solve A**T *x = b or 217*> A**H *x = b. The basic algorithm for A upper triangular is 218*> 219*> for j = 1, ..., n 220*> x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j) 221*> end 222*> 223*> We simultaneously compute two bounds 224*> G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j 225*> M(j) = bound on x(i), 1<=i<=j 226*> 227*> The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we 228*> add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1. 229*> Then the bound on x(j) is 230*> 231*> M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) | 232*> 233*> <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| ) 234*> 1<=i<=j 235*> 236*> and we can safely call CTBSV if 1/M(n) and 1/G(n) are both greater 237*> than max(underflow, 1/overflow). 238*> \endverbatim 239*> 240* ===================================================================== 241 SUBROUTINE CLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X, 242 $ SCALE, CNORM, INFO ) 243* 244* -- LAPACK auxiliary routine -- 245* -- LAPACK is a software package provided by Univ. of Tennessee, -- 246* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 247* 248* .. Scalar Arguments .. 249 CHARACTER DIAG, NORMIN, TRANS, UPLO 250 INTEGER INFO, KD, LDAB, N 251 REAL SCALE 252* .. 253* .. Array Arguments .. 254 REAL CNORM( * ) 255 COMPLEX AB( LDAB, * ), X( * ) 256* .. 257* 258* ===================================================================== 259* 260* .. Parameters .. 261 REAL ZERO, HALF, ONE, TWO 262 PARAMETER ( ZERO = 0.0E+0, HALF = 0.5E+0, ONE = 1.0E+0, 263 $ TWO = 2.0E+0 ) 264* .. 265* .. Local Scalars .. 266 LOGICAL NOTRAN, NOUNIT, UPPER 267 INTEGER I, IMAX, J, JFIRST, JINC, JLAST, JLEN, MAIND 268 REAL BIGNUM, GROW, REC, SMLNUM, TJJ, TMAX, TSCAL, 269 $ XBND, XJ, XMAX 270 COMPLEX CSUMJ, TJJS, USCAL, ZDUM 271* .. 272* .. External Functions .. 273 LOGICAL LSAME 274 INTEGER ICAMAX, ISAMAX 275 REAL SCASUM, SLAMCH 276 COMPLEX CDOTC, CDOTU, CLADIV 277 EXTERNAL LSAME, ICAMAX, ISAMAX, SCASUM, SLAMCH, CDOTC, 278 $ CDOTU, CLADIV 279* .. 280* .. External Subroutines .. 281 EXTERNAL CAXPY, CSSCAL, CTBSV, SLABAD, SSCAL, XERBLA 282* .. 283* .. Intrinsic Functions .. 284 INTRINSIC ABS, AIMAG, CMPLX, CONJG, MAX, MIN, REAL 285* .. 286* .. Statement Functions .. 287 REAL CABS1, CABS2 288* .. 289* .. Statement Function definitions .. 290 CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) ) 291 CABS2( ZDUM ) = ABS( REAL( ZDUM ) / 2. ) + 292 $ ABS( AIMAG( ZDUM ) / 2. ) 293* .. 294* .. Executable Statements .. 295* 296 INFO = 0 297 UPPER = LSAME( UPLO, 'U' ) 298 NOTRAN = LSAME( TRANS, 'N' ) 299 NOUNIT = LSAME( DIAG, 'N' ) 300* 301* Test the input parameters. 302* 303 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 304 INFO = -1 305 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. 306 $ LSAME( TRANS, 'C' ) ) THEN 307 INFO = -2 308 ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN 309 INFO = -3 310 ELSE IF( .NOT.LSAME( NORMIN, 'Y' ) .AND. .NOT. 311 $ LSAME( NORMIN, 'N' ) ) THEN 312 INFO = -4 313 ELSE IF( N.LT.0 ) THEN 314 INFO = -5 315 ELSE IF( KD.LT.0 ) THEN 316 INFO = -6 317 ELSE IF( LDAB.LT.KD+1 ) THEN 318 INFO = -8 319 END IF 320 IF( INFO.NE.0 ) THEN 321 CALL XERBLA( 'CLATBS', -INFO ) 322 RETURN 323 END IF 324* 325* Quick return if possible 326* 327 IF( N.EQ.0 ) 328 $ RETURN 329* 330* Determine machine dependent parameters to control overflow. 331* 332 SMLNUM = SLAMCH( 'Safe minimum' ) 333 BIGNUM = ONE / SMLNUM 334 CALL SLABAD( SMLNUM, BIGNUM ) 335 SMLNUM = SMLNUM / SLAMCH( 'Precision' ) 336 BIGNUM = ONE / SMLNUM 337 SCALE = ONE 338* 339 IF( LSAME( NORMIN, 'N' ) ) THEN 340* 341* Compute the 1-norm of each column, not including the diagonal. 342* 343 IF( UPPER ) THEN 344* 345* A is upper triangular. 346* 347 DO 10 J = 1, N 348 JLEN = MIN( KD, J-1 ) 349 CNORM( J ) = SCASUM( JLEN, AB( KD+1-JLEN, J ), 1 ) 350 10 CONTINUE 351 ELSE 352* 353* A is lower triangular. 354* 355 DO 20 J = 1, N 356 JLEN = MIN( KD, N-J ) 357 IF( JLEN.GT.0 ) THEN 358 CNORM( J ) = SCASUM( JLEN, AB( 2, J ), 1 ) 359 ELSE 360 CNORM( J ) = ZERO 361 END IF 362 20 CONTINUE 363 END IF 364 END IF 365* 366* Scale the column norms by TSCAL if the maximum element in CNORM is 367* greater than BIGNUM/2. 368* 369 IMAX = ISAMAX( N, CNORM, 1 ) 370 TMAX = CNORM( IMAX ) 371 IF( TMAX.LE.BIGNUM*HALF ) THEN 372 TSCAL = ONE 373 ELSE 374 TSCAL = HALF / ( SMLNUM*TMAX ) 375 CALL SSCAL( N, TSCAL, CNORM, 1 ) 376 END IF 377* 378* Compute a bound on the computed solution vector to see if the 379* Level 2 BLAS routine CTBSV can be used. 380* 381 XMAX = ZERO 382 DO 30 J = 1, N 383 XMAX = MAX( XMAX, CABS2( X( J ) ) ) 384 30 CONTINUE 385 XBND = XMAX 386 IF( NOTRAN ) THEN 387* 388* Compute the growth in A * x = b. 389* 390 IF( UPPER ) THEN 391 JFIRST = N 392 JLAST = 1 393 JINC = -1 394 MAIND = KD + 1 395 ELSE 396 JFIRST = 1 397 JLAST = N 398 JINC = 1 399 MAIND = 1 400 END IF 401* 402 IF( TSCAL.NE.ONE ) THEN 403 GROW = ZERO 404 GO TO 60 405 END IF 406* 407 IF( NOUNIT ) THEN 408* 409* A is non-unit triangular. 410* 411* Compute GROW = 1/G(j) and XBND = 1/M(j). 412* Initially, G(0) = max{x(i), i=1,...,n}. 413* 414 GROW = HALF / MAX( XBND, SMLNUM ) 415 XBND = GROW 416 DO 40 J = JFIRST, JLAST, JINC 417* 418* Exit the loop if the growth factor is too small. 419* 420 IF( GROW.LE.SMLNUM ) 421 $ GO TO 60 422* 423 TJJS = AB( MAIND, J ) 424 TJJ = CABS1( TJJS ) 425* 426 IF( TJJ.GE.SMLNUM ) THEN 427* 428* M(j) = G(j-1) / abs(A(j,j)) 429* 430 XBND = MIN( XBND, MIN( ONE, TJJ )*GROW ) 431 ELSE 432* 433* M(j) could overflow, set XBND to 0. 434* 435 XBND = ZERO 436 END IF 437* 438 IF( TJJ+CNORM( J ).GE.SMLNUM ) THEN 439* 440* G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) ) 441* 442 GROW = GROW*( TJJ / ( TJJ+CNORM( J ) ) ) 443 ELSE 444* 445* G(j) could overflow, set GROW to 0. 446* 447 GROW = ZERO 448 END IF 449 40 CONTINUE 450 GROW = XBND 451 ELSE 452* 453* A is unit triangular. 454* 455* Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}. 456* 457 GROW = MIN( ONE, HALF / MAX( XBND, SMLNUM ) ) 458 DO 50 J = JFIRST, JLAST, JINC 459* 460* Exit the loop if the growth factor is too small. 461* 462 IF( GROW.LE.SMLNUM ) 463 $ GO TO 60 464* 465* G(j) = G(j-1)*( 1 + CNORM(j) ) 466* 467 GROW = GROW*( ONE / ( ONE+CNORM( J ) ) ) 468 50 CONTINUE 469 END IF 470 60 CONTINUE 471* 472 ELSE 473* 474* Compute the growth in A**T * x = b or A**H * x = b. 475* 476 IF( UPPER ) THEN 477 JFIRST = 1 478 JLAST = N 479 JINC = 1 480 MAIND = KD + 1 481 ELSE 482 JFIRST = N 483 JLAST = 1 484 JINC = -1 485 MAIND = 1 486 END IF 487* 488 IF( TSCAL.NE.ONE ) THEN 489 GROW = ZERO 490 GO TO 90 491 END IF 492* 493 IF( NOUNIT ) THEN 494* 495* A is non-unit triangular. 496* 497* Compute GROW = 1/G(j) and XBND = 1/M(j). 498* Initially, M(0) = max{x(i), i=1,...,n}. 499* 500 GROW = HALF / MAX( XBND, SMLNUM ) 501 XBND = GROW 502 DO 70 J = JFIRST, JLAST, JINC 503* 504* Exit the loop if the growth factor is too small. 505* 506 IF( GROW.LE.SMLNUM ) 507 $ GO TO 90 508* 509* G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) ) 510* 511 XJ = ONE + CNORM( J ) 512 GROW = MIN( GROW, XBND / XJ ) 513* 514 TJJS = AB( MAIND, J ) 515 TJJ = CABS1( TJJS ) 516* 517 IF( TJJ.GE.SMLNUM ) THEN 518* 519* M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j)) 520* 521 IF( XJ.GT.TJJ ) 522 $ XBND = XBND*( TJJ / XJ ) 523 ELSE 524* 525* M(j) could overflow, set XBND to 0. 526* 527 XBND = ZERO 528 END IF 529 70 CONTINUE 530 GROW = MIN( GROW, XBND ) 531 ELSE 532* 533* A is unit triangular. 534* 535* Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}. 536* 537 GROW = MIN( ONE, HALF / MAX( XBND, SMLNUM ) ) 538 DO 80 J = JFIRST, JLAST, JINC 539* 540* Exit the loop if the growth factor is too small. 541* 542 IF( GROW.LE.SMLNUM ) 543 $ GO TO 90 544* 545* G(j) = ( 1 + CNORM(j) )*G(j-1) 546* 547 XJ = ONE + CNORM( J ) 548 GROW = GROW / XJ 549 80 CONTINUE 550 END IF 551 90 CONTINUE 552 END IF 553* 554 IF( ( GROW*TSCAL ).GT.SMLNUM ) THEN 555* 556* Use the Level 2 BLAS solve if the reciprocal of the bound on 557* elements of X is not too small. 558* 559 CALL CTBSV( UPLO, TRANS, DIAG, N, KD, AB, LDAB, X, 1 ) 560 ELSE 561* 562* Use a Level 1 BLAS solve, scaling intermediate results. 563* 564 IF( XMAX.GT.BIGNUM*HALF ) THEN 565* 566* Scale X so that its components are less than or equal to 567* BIGNUM in absolute value. 568* 569 SCALE = ( BIGNUM*HALF ) / XMAX 570 CALL CSSCAL( N, SCALE, X, 1 ) 571 XMAX = BIGNUM 572 ELSE 573 XMAX = XMAX*TWO 574 END IF 575* 576 IF( NOTRAN ) THEN 577* 578* Solve A * x = b 579* 580 DO 110 J = JFIRST, JLAST, JINC 581* 582* Compute x(j) = b(j) / A(j,j), scaling x if necessary. 583* 584 XJ = CABS1( X( J ) ) 585 IF( NOUNIT ) THEN 586 TJJS = AB( MAIND, J )*TSCAL 587 ELSE 588 TJJS = TSCAL 589 IF( TSCAL.EQ.ONE ) 590 $ GO TO 105 591 END IF 592 TJJ = CABS1( TJJS ) 593 IF( TJJ.GT.SMLNUM ) THEN 594* 595* abs(A(j,j)) > SMLNUM: 596* 597 IF( TJJ.LT.ONE ) THEN 598 IF( XJ.GT.TJJ*BIGNUM ) THEN 599* 600* Scale x by 1/b(j). 601* 602 REC = ONE / XJ 603 CALL CSSCAL( N, REC, X, 1 ) 604 SCALE = SCALE*REC 605 XMAX = XMAX*REC 606 END IF 607 END IF 608 X( J ) = CLADIV( X( J ), TJJS ) 609 XJ = CABS1( X( J ) ) 610 ELSE IF( TJJ.GT.ZERO ) THEN 611* 612* 0 < abs(A(j,j)) <= SMLNUM: 613* 614 IF( XJ.GT.TJJ*BIGNUM ) THEN 615* 616* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM 617* to avoid overflow when dividing by A(j,j). 618* 619 REC = ( TJJ*BIGNUM ) / XJ 620 IF( CNORM( J ).GT.ONE ) THEN 621* 622* Scale by 1/CNORM(j) to avoid overflow when 623* multiplying x(j) times column j. 624* 625 REC = REC / CNORM( J ) 626 END IF 627 CALL CSSCAL( N, REC, X, 1 ) 628 SCALE = SCALE*REC 629 XMAX = XMAX*REC 630 END IF 631 X( J ) = CLADIV( X( J ), TJJS ) 632 XJ = CABS1( X( J ) ) 633 ELSE 634* 635* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and 636* scale = 0, and compute a solution to A*x = 0. 637* 638 DO 100 I = 1, N 639 X( I ) = ZERO 640 100 CONTINUE 641 X( J ) = ONE 642 XJ = ONE 643 SCALE = ZERO 644 XMAX = ZERO 645 END IF 646 105 CONTINUE 647* 648* Scale x if necessary to avoid overflow when adding a 649* multiple of column j of A. 650* 651 IF( XJ.GT.ONE ) THEN 652 REC = ONE / XJ 653 IF( CNORM( J ).GT.( BIGNUM-XMAX )*REC ) THEN 654* 655* Scale x by 1/(2*abs(x(j))). 656* 657 REC = REC*HALF 658 CALL CSSCAL( N, REC, X, 1 ) 659 SCALE = SCALE*REC 660 END IF 661 ELSE IF( XJ*CNORM( J ).GT.( BIGNUM-XMAX ) ) THEN 662* 663* Scale x by 1/2. 664* 665 CALL CSSCAL( N, HALF, X, 1 ) 666 SCALE = SCALE*HALF 667 END IF 668* 669 IF( UPPER ) THEN 670 IF( J.GT.1 ) THEN 671* 672* Compute the update 673* x(max(1,j-kd):j-1) := x(max(1,j-kd):j-1) - 674* x(j)* A(max(1,j-kd):j-1,j) 675* 676 JLEN = MIN( KD, J-1 ) 677 CALL CAXPY( JLEN, -X( J )*TSCAL, 678 $ AB( KD+1-JLEN, J ), 1, X( J-JLEN ), 1 ) 679 I = ICAMAX( J-1, X, 1 ) 680 XMAX = CABS1( X( I ) ) 681 END IF 682 ELSE IF( J.LT.N ) THEN 683* 684* Compute the update 685* x(j+1:min(j+kd,n)) := x(j+1:min(j+kd,n)) - 686* x(j) * A(j+1:min(j+kd,n),j) 687* 688 JLEN = MIN( KD, N-J ) 689 IF( JLEN.GT.0 ) 690 $ CALL CAXPY( JLEN, -X( J )*TSCAL, AB( 2, J ), 1, 691 $ X( J+1 ), 1 ) 692 I = J + ICAMAX( N-J, X( J+1 ), 1 ) 693 XMAX = CABS1( X( I ) ) 694 END IF 695 110 CONTINUE 696* 697 ELSE IF( LSAME( TRANS, 'T' ) ) THEN 698* 699* Solve A**T * x = b 700* 701 DO 150 J = JFIRST, JLAST, JINC 702* 703* Compute x(j) = b(j) - sum A(k,j)*x(k). 704* k<>j 705* 706 XJ = CABS1( X( J ) ) 707 USCAL = TSCAL 708 REC = ONE / MAX( XMAX, ONE ) 709 IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN 710* 711* If x(j) could overflow, scale x by 1/(2*XMAX). 712* 713 REC = REC*HALF 714 IF( NOUNIT ) THEN 715 TJJS = AB( MAIND, J )*TSCAL 716 ELSE 717 TJJS = TSCAL 718 END IF 719 TJJ = CABS1( TJJS ) 720 IF( TJJ.GT.ONE ) THEN 721* 722* Divide by A(j,j) when scaling x if A(j,j) > 1. 723* 724 REC = MIN( ONE, REC*TJJ ) 725 USCAL = CLADIV( USCAL, TJJS ) 726 END IF 727 IF( REC.LT.ONE ) THEN 728 CALL CSSCAL( N, REC, X, 1 ) 729 SCALE = SCALE*REC 730 XMAX = XMAX*REC 731 END IF 732 END IF 733* 734 CSUMJ = ZERO 735 IF( USCAL.EQ.CMPLX( ONE ) ) THEN 736* 737* If the scaling needed for A in the dot product is 1, 738* call CDOTU to perform the dot product. 739* 740 IF( UPPER ) THEN 741 JLEN = MIN( KD, J-1 ) 742 CSUMJ = CDOTU( JLEN, AB( KD+1-JLEN, J ), 1, 743 $ X( J-JLEN ), 1 ) 744 ELSE 745 JLEN = MIN( KD, N-J ) 746 IF( JLEN.GT.1 ) 747 $ CSUMJ = CDOTU( JLEN, AB( 2, J ), 1, X( J+1 ), 748 $ 1 ) 749 END IF 750 ELSE 751* 752* Otherwise, use in-line code for the dot product. 753* 754 IF( UPPER ) THEN 755 JLEN = MIN( KD, J-1 ) 756 DO 120 I = 1, JLEN 757 CSUMJ = CSUMJ + ( AB( KD+I-JLEN, J )*USCAL )* 758 $ X( J-JLEN-1+I ) 759 120 CONTINUE 760 ELSE 761 JLEN = MIN( KD, N-J ) 762 DO 130 I = 1, JLEN 763 CSUMJ = CSUMJ + ( AB( I+1, J )*USCAL )*X( J+I ) 764 130 CONTINUE 765 END IF 766 END IF 767* 768 IF( USCAL.EQ.CMPLX( TSCAL ) ) THEN 769* 770* Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j) 771* was not used to scale the dotproduct. 772* 773 X( J ) = X( J ) - CSUMJ 774 XJ = CABS1( X( J ) ) 775 IF( NOUNIT ) THEN 776* 777* Compute x(j) = x(j) / A(j,j), scaling if necessary. 778* 779 TJJS = AB( MAIND, J )*TSCAL 780 ELSE 781 TJJS = TSCAL 782 IF( TSCAL.EQ.ONE ) 783 $ GO TO 145 784 END IF 785 TJJ = CABS1( TJJS ) 786 IF( TJJ.GT.SMLNUM ) THEN 787* 788* abs(A(j,j)) > SMLNUM: 789* 790 IF( TJJ.LT.ONE ) THEN 791 IF( XJ.GT.TJJ*BIGNUM ) THEN 792* 793* Scale X by 1/abs(x(j)). 794* 795 REC = ONE / XJ 796 CALL CSSCAL( N, REC, X, 1 ) 797 SCALE = SCALE*REC 798 XMAX = XMAX*REC 799 END IF 800 END IF 801 X( J ) = CLADIV( X( J ), TJJS ) 802 ELSE IF( TJJ.GT.ZERO ) THEN 803* 804* 0 < abs(A(j,j)) <= SMLNUM: 805* 806 IF( XJ.GT.TJJ*BIGNUM ) THEN 807* 808* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM. 809* 810 REC = ( TJJ*BIGNUM ) / XJ 811 CALL CSSCAL( N, REC, X, 1 ) 812 SCALE = SCALE*REC 813 XMAX = XMAX*REC 814 END IF 815 X( J ) = CLADIV( X( J ), TJJS ) 816 ELSE 817* 818* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and 819* scale = 0 and compute a solution to A**T *x = 0. 820* 821 DO 140 I = 1, N 822 X( I ) = ZERO 823 140 CONTINUE 824 X( J ) = ONE 825 SCALE = ZERO 826 XMAX = ZERO 827 END IF 828 145 CONTINUE 829 ELSE 830* 831* Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot 832* product has already been divided by 1/A(j,j). 833* 834 X( J ) = CLADIV( X( J ), TJJS ) - CSUMJ 835 END IF 836 XMAX = MAX( XMAX, CABS1( X( J ) ) ) 837 150 CONTINUE 838* 839 ELSE 840* 841* Solve A**H * x = b 842* 843 DO 190 J = JFIRST, JLAST, JINC 844* 845* Compute x(j) = b(j) - sum A(k,j)*x(k). 846* k<>j 847* 848 XJ = CABS1( X( J ) ) 849 USCAL = TSCAL 850 REC = ONE / MAX( XMAX, ONE ) 851 IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN 852* 853* If x(j) could overflow, scale x by 1/(2*XMAX). 854* 855 REC = REC*HALF 856 IF( NOUNIT ) THEN 857 TJJS = CONJG( AB( MAIND, J ) )*TSCAL 858 ELSE 859 TJJS = TSCAL 860 END IF 861 TJJ = CABS1( TJJS ) 862 IF( TJJ.GT.ONE ) THEN 863* 864* Divide by A(j,j) when scaling x if A(j,j) > 1. 865* 866 REC = MIN( ONE, REC*TJJ ) 867 USCAL = CLADIV( USCAL, TJJS ) 868 END IF 869 IF( REC.LT.ONE ) THEN 870 CALL CSSCAL( N, REC, X, 1 ) 871 SCALE = SCALE*REC 872 XMAX = XMAX*REC 873 END IF 874 END IF 875* 876 CSUMJ = ZERO 877 IF( USCAL.EQ.CMPLX( ONE ) ) THEN 878* 879* If the scaling needed for A in the dot product is 1, 880* call CDOTC to perform the dot product. 881* 882 IF( UPPER ) THEN 883 JLEN = MIN( KD, J-1 ) 884 CSUMJ = CDOTC( JLEN, AB( KD+1-JLEN, J ), 1, 885 $ X( J-JLEN ), 1 ) 886 ELSE 887 JLEN = MIN( KD, N-J ) 888 IF( JLEN.GT.1 ) 889 $ CSUMJ = CDOTC( JLEN, AB( 2, J ), 1, X( J+1 ), 890 $ 1 ) 891 END IF 892 ELSE 893* 894* Otherwise, use in-line code for the dot product. 895* 896 IF( UPPER ) THEN 897 JLEN = MIN( KD, J-1 ) 898 DO 160 I = 1, JLEN 899 CSUMJ = CSUMJ + ( CONJG( AB( KD+I-JLEN, J ) )* 900 $ USCAL )*X( J-JLEN-1+I ) 901 160 CONTINUE 902 ELSE 903 JLEN = MIN( KD, N-J ) 904 DO 170 I = 1, JLEN 905 CSUMJ = CSUMJ + ( CONJG( AB( I+1, J ) )*USCAL )* 906 $ X( J+I ) 907 170 CONTINUE 908 END IF 909 END IF 910* 911 IF( USCAL.EQ.CMPLX( TSCAL ) ) THEN 912* 913* Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j) 914* was not used to scale the dotproduct. 915* 916 X( J ) = X( J ) - CSUMJ 917 XJ = CABS1( X( J ) ) 918 IF( NOUNIT ) THEN 919* 920* Compute x(j) = x(j) / A(j,j), scaling if necessary. 921* 922 TJJS = CONJG( AB( MAIND, J ) )*TSCAL 923 ELSE 924 TJJS = TSCAL 925 IF( TSCAL.EQ.ONE ) 926 $ GO TO 185 927 END IF 928 TJJ = CABS1( TJJS ) 929 IF( TJJ.GT.SMLNUM ) THEN 930* 931* abs(A(j,j)) > SMLNUM: 932* 933 IF( TJJ.LT.ONE ) THEN 934 IF( XJ.GT.TJJ*BIGNUM ) THEN 935* 936* Scale X by 1/abs(x(j)). 937* 938 REC = ONE / XJ 939 CALL CSSCAL( N, REC, X, 1 ) 940 SCALE = SCALE*REC 941 XMAX = XMAX*REC 942 END IF 943 END IF 944 X( J ) = CLADIV( X( J ), TJJS ) 945 ELSE IF( TJJ.GT.ZERO ) THEN 946* 947* 0 < abs(A(j,j)) <= SMLNUM: 948* 949 IF( XJ.GT.TJJ*BIGNUM ) THEN 950* 951* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM. 952* 953 REC = ( TJJ*BIGNUM ) / XJ 954 CALL CSSCAL( N, REC, X, 1 ) 955 SCALE = SCALE*REC 956 XMAX = XMAX*REC 957 END IF 958 X( J ) = CLADIV( X( J ), TJJS ) 959 ELSE 960* 961* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and 962* scale = 0 and compute a solution to A**H *x = 0. 963* 964 DO 180 I = 1, N 965 X( I ) = ZERO 966 180 CONTINUE 967 X( J ) = ONE 968 SCALE = ZERO 969 XMAX = ZERO 970 END IF 971 185 CONTINUE 972 ELSE 973* 974* Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot 975* product has already been divided by 1/A(j,j). 976* 977 X( J ) = CLADIV( X( J ), TJJS ) - CSUMJ 978 END IF 979 XMAX = MAX( XMAX, CABS1( X( J ) ) ) 980 190 CONTINUE 981 END IF 982 SCALE = SCALE / TSCAL 983 END IF 984* 985* Scale the column norms by 1/TSCAL for return. 986* 987 IF( TSCAL.NE.ONE ) THEN 988 CALL SSCAL( N, ONE / TSCAL, CNORM, 1 ) 989 END IF 990* 991 RETURN 992* 993* End of CLATBS 994* 995 END 996