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