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