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