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