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