1*> \brief \b SLATPS 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 SLATPS + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slatps.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slatps.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slatps.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE SLATPS( 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 AP( * ), CNORM( * ), X( * ) 31* .. 32* 33* 34*> \par Purpose: 35* ============= 36*> 37*> \verbatim 38*> 39*> SLATPS solves one of the triangular systems 40*> 41*> A *x = s*b or A**T*x = s*b 42*> 43*> with scaling to prevent overflow, where A is an upper or lower 44*> triangular matrix stored in packed form. Here A**T denotes the 45*> transpose of A, x and b are n-element vectors, and s is a scaling 46*> factor, usually less than or equal to 1, chosen so that the 47*> components of x will be less than the overflow threshold. If the 48*> unscaled problem will not cause overflow, the Level 2 BLAS routine 49*> STPSV is called. If the matrix A is singular (A(j,j) = 0 for some j), 50*> then s is set to 0 and a non-trivial solution to A*x = 0 is returned. 51*> \endverbatim 52* 53* Arguments: 54* ========== 55* 56*> \param[in] UPLO 57*> \verbatim 58*> UPLO is CHARACTER*1 59*> Specifies whether the matrix A is upper or lower triangular. 60*> = 'U': Upper triangular 61*> = 'L': Lower triangular 62*> \endverbatim 63*> 64*> \param[in] TRANS 65*> \verbatim 66*> TRANS is CHARACTER*1 67*> Specifies the operation applied to A. 68*> = 'N': Solve A * x = s*b (No transpose) 69*> = 'T': Solve A**T* x = s*b (Transpose) 70*> = 'C': Solve A**T* x = s*b (Conjugate transpose = Transpose) 71*> \endverbatim 72*> 73*> \param[in] DIAG 74*> \verbatim 75*> DIAG is CHARACTER*1 76*> Specifies whether or not the matrix A is unit triangular. 77*> = 'N': Non-unit triangular 78*> = 'U': Unit triangular 79*> \endverbatim 80*> 81*> \param[in] NORMIN 82*> \verbatim 83*> NORMIN is CHARACTER*1 84*> Specifies whether CNORM has been set or not. 85*> = 'Y': CNORM contains the column norms on entry 86*> = 'N': CNORM is not set on entry. On exit, the norms will 87*> be computed and stored in CNORM. 88*> \endverbatim 89*> 90*> \param[in] N 91*> \verbatim 92*> N is INTEGER 93*> The order of the matrix A. N >= 0. 94*> \endverbatim 95*> 96*> \param[in] AP 97*> \verbatim 98*> AP is REAL array, dimension (N*(N+1)/2) 99*> The upper or lower triangular matrix A, packed columnwise in 100*> a linear array. The j-th column of A is stored in the array 101*> AP as follows: 102*> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; 103*> if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. 104*> \endverbatim 105*> 106*> \param[in,out] X 107*> \verbatim 108*> X is REAL array, dimension (N) 109*> On entry, the right hand side b of the triangular system. 110*> On exit, X is overwritten by the solution vector x. 111*> \endverbatim 112*> 113*> \param[out] SCALE 114*> \verbatim 115*> SCALE is REAL 116*> The scaling factor s for the triangular system 117*> A * x = s*b or A**T* x = s*b. 118*> If SCALE = 0, the matrix A is singular or badly scaled, and 119*> the vector x is an exact or approximate solution to A*x = 0. 120*> \endverbatim 121*> 122*> \param[in,out] CNORM 123*> \verbatim 124*> CNORM is REAL array, dimension (N) 125*> 126*> If NORMIN = 'Y', CNORM is an input argument and CNORM(j) 127*> contains the norm of the off-diagonal part of the j-th column 128*> of A. If TRANS = 'N', CNORM(j) must be greater than or equal 129*> to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) 130*> must be greater than or equal to the 1-norm. 131*> 132*> If NORMIN = 'N', CNORM is an output argument and CNORM(j) 133*> returns the 1-norm of the offdiagonal part of the j-th column 134*> of A. 135*> \endverbatim 136*> 137*> \param[out] INFO 138*> \verbatim 139*> INFO is INTEGER 140*> = 0: successful exit 141*> < 0: if INFO = -k, the k-th argument had an illegal value 142*> \endverbatim 143* 144* Authors: 145* ======== 146* 147*> \author Univ. of Tennessee 148*> \author Univ. of California Berkeley 149*> \author Univ. of Colorado Denver 150*> \author NAG Ltd. 151* 152*> \date December 2016 153* 154*> \ingroup realOTHERauxiliary 155* 156*> \par Further Details: 157* ===================== 158*> 159*> \verbatim 160*> 161*> A rough bound on x is computed; if that is less than overflow, STPSV 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 STPSV 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. The basic 205*> algorithm for A upper triangular is 206*> 207*> for j = 1, ..., n 208*> x(j) := ( b(j) - A[1:j-1,j]**T * 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]**T * 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 STPSV if 1/M(n) and 1/G(n) are both greater 225*> than max(underflow, 1/overflow). 226*> \endverbatim 227*> 228* ===================================================================== 229 SUBROUTINE SLATPS( UPLO, TRANS, DIAG, NORMIN, N, AP, X, SCALE, 230 $ CNORM, INFO ) 231* 232* -- LAPACK auxiliary routine (version 3.7.0) -- 233* -- LAPACK is a software package provided by Univ. of Tennessee, -- 234* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 235* December 2016 236* 237* .. Scalar Arguments .. 238 CHARACTER DIAG, NORMIN, TRANS, UPLO 239 INTEGER INFO, N 240 REAL SCALE 241* .. 242* .. Array Arguments .. 243 REAL AP( * ), CNORM( * ), X( * ) 244* .. 245* 246* ===================================================================== 247* 248* .. Parameters .. 249 REAL ZERO, HALF, ONE 250 PARAMETER ( ZERO = 0.0E+0, HALF = 0.5E+0, ONE = 1.0E+0 ) 251* .. 252* .. Local Scalars .. 253 LOGICAL NOTRAN, NOUNIT, UPPER 254 INTEGER I, IMAX, IP, J, JFIRST, JINC, JLAST, JLEN 255 REAL BIGNUM, GROW, REC, SMLNUM, SUMJ, TJJ, TJJS, 256 $ TMAX, TSCAL, USCAL, XBND, XJ, XMAX 257* .. 258* .. External Functions .. 259 LOGICAL LSAME 260 INTEGER ISAMAX 261 REAL SASUM, SDOT, SLAMCH 262 EXTERNAL LSAME, ISAMAX, SASUM, SDOT, SLAMCH 263* .. 264* .. External Subroutines .. 265 EXTERNAL SAXPY, SSCAL, STPSV, XERBLA 266* .. 267* .. Intrinsic Functions .. 268 INTRINSIC ABS, MAX, MIN 269* .. 270* .. Executable Statements .. 271* 272 INFO = 0 273 UPPER = LSAME( UPLO, 'U' ) 274 NOTRAN = LSAME( TRANS, 'N' ) 275 NOUNIT = LSAME( DIAG, 'N' ) 276* 277* Test the input parameters. 278* 279 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 280 INFO = -1 281 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. 282 $ LSAME( TRANS, 'C' ) ) THEN 283 INFO = -2 284 ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN 285 INFO = -3 286 ELSE IF( .NOT.LSAME( NORMIN, 'Y' ) .AND. .NOT. 287 $ LSAME( NORMIN, 'N' ) ) THEN 288 INFO = -4 289 ELSE IF( N.LT.0 ) THEN 290 INFO = -5 291 END IF 292 IF( INFO.NE.0 ) THEN 293 CALL XERBLA( 'SLATPS', -INFO ) 294 RETURN 295 END IF 296* 297* Quick return if possible 298* 299 IF( N.EQ.0 ) 300 $ RETURN 301* 302* Determine machine dependent parameters to control overflow. 303* 304 SMLNUM = SLAMCH( 'Safe minimum' ) / SLAMCH( 'Precision' ) 305 BIGNUM = ONE / SMLNUM 306 SCALE = ONE 307* 308 IF( LSAME( NORMIN, 'N' ) ) THEN 309* 310* Compute the 1-norm of each column, not including the diagonal. 311* 312 IF( UPPER ) THEN 313* 314* A is upper triangular. 315* 316 IP = 1 317 DO 10 J = 1, N 318 CNORM( J ) = SASUM( J-1, AP( IP ), 1 ) 319 IP = IP + J 320 10 CONTINUE 321 ELSE 322* 323* A is lower triangular. 324* 325 IP = 1 326 DO 20 J = 1, N - 1 327 CNORM( J ) = SASUM( N-J, AP( IP+1 ), 1 ) 328 IP = IP + N - J + 1 329 20 CONTINUE 330 CNORM( N ) = ZERO 331 END IF 332 END IF 333* 334* Scale the column norms by TSCAL if the maximum element in CNORM is 335* greater than BIGNUM. 336* 337 IMAX = ISAMAX( N, CNORM, 1 ) 338 TMAX = CNORM( IMAX ) 339 IF( TMAX.LE.BIGNUM ) THEN 340 TSCAL = ONE 341 ELSE 342 TSCAL = ONE / ( SMLNUM*TMAX ) 343 CALL SSCAL( N, TSCAL, CNORM, 1 ) 344 END IF 345* 346* Compute a bound on the computed solution vector to see if the 347* Level 2 BLAS routine STPSV can be used. 348* 349 J = ISAMAX( N, X, 1 ) 350 XMAX = ABS( X( J ) ) 351 XBND = XMAX 352 IF( NOTRAN ) THEN 353* 354* Compute the growth in A * x = b. 355* 356 IF( UPPER ) THEN 357 JFIRST = N 358 JLAST = 1 359 JINC = -1 360 ELSE 361 JFIRST = 1 362 JLAST = N 363 JINC = 1 364 END IF 365* 366 IF( TSCAL.NE.ONE ) THEN 367 GROW = ZERO 368 GO TO 50 369 END IF 370* 371 IF( NOUNIT ) THEN 372* 373* A is non-unit triangular. 374* 375* Compute GROW = 1/G(j) and XBND = 1/M(j). 376* Initially, G(0) = max{x(i), i=1,...,n}. 377* 378 GROW = ONE / MAX( XBND, SMLNUM ) 379 XBND = GROW 380 IP = JFIRST*( JFIRST+1 ) / 2 381 JLEN = N 382 DO 30 J = JFIRST, JLAST, JINC 383* 384* Exit the loop if the growth factor is too small. 385* 386 IF( GROW.LE.SMLNUM ) 387 $ GO TO 50 388* 389* M(j) = G(j-1) / abs(A(j,j)) 390* 391 TJJ = ABS( AP( IP ) ) 392 XBND = MIN( XBND, MIN( ONE, TJJ )*GROW ) 393 IF( TJJ+CNORM( J ).GE.SMLNUM ) THEN 394* 395* G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) ) 396* 397 GROW = GROW*( TJJ / ( TJJ+CNORM( J ) ) ) 398 ELSE 399* 400* G(j) could overflow, set GROW to 0. 401* 402 GROW = ZERO 403 END IF 404 IP = IP + JINC*JLEN 405 JLEN = JLEN - 1 406 30 CONTINUE 407 GROW = XBND 408 ELSE 409* 410* A is unit triangular. 411* 412* Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}. 413* 414 GROW = MIN( ONE, ONE / MAX( XBND, SMLNUM ) ) 415 DO 40 J = JFIRST, JLAST, JINC 416* 417* Exit the loop if the growth factor is too small. 418* 419 IF( GROW.LE.SMLNUM ) 420 $ GO TO 50 421* 422* G(j) = G(j-1)*( 1 + CNORM(j) ) 423* 424 GROW = GROW*( ONE / ( ONE+CNORM( J ) ) ) 425 40 CONTINUE 426 END IF 427 50 CONTINUE 428* 429 ELSE 430* 431* Compute the growth in A**T * x = b. 432* 433 IF( UPPER ) THEN 434 JFIRST = 1 435 JLAST = N 436 JINC = 1 437 ELSE 438 JFIRST = N 439 JLAST = 1 440 JINC = -1 441 END IF 442* 443 IF( TSCAL.NE.ONE ) THEN 444 GROW = ZERO 445 GO TO 80 446 END IF 447* 448 IF( NOUNIT ) THEN 449* 450* A is non-unit triangular. 451* 452* Compute GROW = 1/G(j) and XBND = 1/M(j). 453* Initially, M(0) = max{x(i), i=1,...,n}. 454* 455 GROW = ONE / MAX( XBND, SMLNUM ) 456 XBND = GROW 457 IP = JFIRST*( JFIRST+1 ) / 2 458 JLEN = 1 459 DO 60 J = JFIRST, JLAST, JINC 460* 461* Exit the loop if the growth factor is too small. 462* 463 IF( GROW.LE.SMLNUM ) 464 $ GO TO 80 465* 466* G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) ) 467* 468 XJ = ONE + CNORM( J ) 469 GROW = MIN( GROW, XBND / XJ ) 470* 471* M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j)) 472* 473 TJJ = ABS( AP( IP ) ) 474 IF( XJ.GT.TJJ ) 475 $ XBND = XBND*( TJJ / XJ ) 476 JLEN = JLEN + 1 477 IP = IP + JINC*JLEN 478 60 CONTINUE 479 GROW = MIN( GROW, XBND ) 480 ELSE 481* 482* A is unit triangular. 483* 484* Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}. 485* 486 GROW = MIN( ONE, ONE / MAX( XBND, SMLNUM ) ) 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 80 493* 494* G(j) = ( 1 + CNORM(j) )*G(j-1) 495* 496 XJ = ONE + CNORM( J ) 497 GROW = GROW / XJ 498 70 CONTINUE 499 END IF 500 80 CONTINUE 501 END IF 502* 503 IF( ( GROW*TSCAL ).GT.SMLNUM ) THEN 504* 505* Use the Level 2 BLAS solve if the reciprocal of the bound on 506* elements of X is not too small. 507* 508 CALL STPSV( UPLO, TRANS, DIAG, N, AP, X, 1 ) 509 ELSE 510* 511* Use a Level 1 BLAS solve, scaling intermediate results. 512* 513 IF( XMAX.GT.BIGNUM ) THEN 514* 515* Scale X so that its components are less than or equal to 516* BIGNUM in absolute value. 517* 518 SCALE = BIGNUM / XMAX 519 CALL SSCAL( N, SCALE, X, 1 ) 520 XMAX = BIGNUM 521 END IF 522* 523 IF( NOTRAN ) THEN 524* 525* Solve A * x = b 526* 527 IP = JFIRST*( JFIRST+1 ) / 2 528 DO 100 J = JFIRST, JLAST, JINC 529* 530* Compute x(j) = b(j) / A(j,j), scaling x if necessary. 531* 532 XJ = ABS( X( J ) ) 533 IF( NOUNIT ) THEN 534 TJJS = AP( IP )*TSCAL 535 ELSE 536 TJJS = TSCAL 537 IF( TSCAL.EQ.ONE ) 538 $ GO TO 95 539 END IF 540 TJJ = ABS( TJJS ) 541 IF( TJJ.GT.SMLNUM ) THEN 542* 543* abs(A(j,j)) > SMLNUM: 544* 545 IF( TJJ.LT.ONE ) THEN 546 IF( XJ.GT.TJJ*BIGNUM ) THEN 547* 548* Scale x by 1/b(j). 549* 550 REC = ONE / XJ 551 CALL SSCAL( N, REC, X, 1 ) 552 SCALE = SCALE*REC 553 XMAX = XMAX*REC 554 END IF 555 END IF 556 X( J ) = X( J ) / TJJS 557 XJ = ABS( X( J ) ) 558 ELSE IF( TJJ.GT.ZERO ) THEN 559* 560* 0 < abs(A(j,j)) <= SMLNUM: 561* 562 IF( XJ.GT.TJJ*BIGNUM ) THEN 563* 564* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM 565* to avoid overflow when dividing by A(j,j). 566* 567 REC = ( TJJ*BIGNUM ) / XJ 568 IF( CNORM( J ).GT.ONE ) THEN 569* 570* Scale by 1/CNORM(j) to avoid overflow when 571* multiplying x(j) times column j. 572* 573 REC = REC / CNORM( J ) 574 END IF 575 CALL SSCAL( N, REC, X, 1 ) 576 SCALE = SCALE*REC 577 XMAX = XMAX*REC 578 END IF 579 X( J ) = X( J ) / TJJS 580 XJ = ABS( X( J ) ) 581 ELSE 582* 583* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and 584* scale = 0, and compute a solution to A*x = 0. 585* 586 DO 90 I = 1, N 587 X( I ) = ZERO 588 90 CONTINUE 589 X( J ) = ONE 590 XJ = ONE 591 SCALE = ZERO 592 XMAX = ZERO 593 END IF 594 95 CONTINUE 595* 596* Scale x if necessary to avoid overflow when adding a 597* multiple of column j of A. 598* 599 IF( XJ.GT.ONE ) THEN 600 REC = ONE / XJ 601 IF( CNORM( J ).GT.( BIGNUM-XMAX )*REC ) THEN 602* 603* Scale x by 1/(2*abs(x(j))). 604* 605 REC = REC*HALF 606 CALL SSCAL( N, REC, X, 1 ) 607 SCALE = SCALE*REC 608 END IF 609 ELSE IF( XJ*CNORM( J ).GT.( BIGNUM-XMAX ) ) THEN 610* 611* Scale x by 1/2. 612* 613 CALL SSCAL( N, HALF, X, 1 ) 614 SCALE = SCALE*HALF 615 END IF 616* 617 IF( UPPER ) THEN 618 IF( J.GT.1 ) THEN 619* 620* Compute the update 621* x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j) 622* 623 CALL SAXPY( J-1, -X( J )*TSCAL, AP( IP-J+1 ), 1, X, 624 $ 1 ) 625 I = ISAMAX( J-1, X, 1 ) 626 XMAX = ABS( X( I ) ) 627 END IF 628 IP = IP - J 629 ELSE 630 IF( J.LT.N ) THEN 631* 632* Compute the update 633* x(j+1:n) := x(j+1:n) - x(j) * A(j+1:n,j) 634* 635 CALL SAXPY( N-J, -X( J )*TSCAL, AP( IP+1 ), 1, 636 $ X( J+1 ), 1 ) 637 I = J + ISAMAX( N-J, X( J+1 ), 1 ) 638 XMAX = ABS( X( I ) ) 639 END IF 640 IP = IP + N - J + 1 641 END IF 642 100 CONTINUE 643* 644 ELSE 645* 646* Solve A**T * x = b 647* 648 IP = JFIRST*( JFIRST+1 ) / 2 649 JLEN = 1 650 DO 140 J = JFIRST, JLAST, JINC 651* 652* Compute x(j) = b(j) - sum A(k,j)*x(k). 653* k<>j 654* 655 XJ = ABS( X( J ) ) 656 USCAL = TSCAL 657 REC = ONE / MAX( XMAX, ONE ) 658 IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN 659* 660* If x(j) could overflow, scale x by 1/(2*XMAX). 661* 662 REC = REC*HALF 663 IF( NOUNIT ) THEN 664 TJJS = AP( IP )*TSCAL 665 ELSE 666 TJJS = TSCAL 667 END IF 668 TJJ = ABS( TJJS ) 669 IF( TJJ.GT.ONE ) THEN 670* 671* Divide by A(j,j) when scaling x if A(j,j) > 1. 672* 673 REC = MIN( ONE, REC*TJJ ) 674 USCAL = USCAL / TJJS 675 END IF 676 IF( REC.LT.ONE ) THEN 677 CALL SSCAL( N, REC, X, 1 ) 678 SCALE = SCALE*REC 679 XMAX = XMAX*REC 680 END IF 681 END IF 682* 683 SUMJ = ZERO 684 IF( USCAL.EQ.ONE ) THEN 685* 686* If the scaling needed for A in the dot product is 1, 687* call SDOT to perform the dot product. 688* 689 IF( UPPER ) THEN 690 SUMJ = SDOT( J-1, AP( IP-J+1 ), 1, X, 1 ) 691 ELSE IF( J.LT.N ) THEN 692 SUMJ = SDOT( N-J, AP( IP+1 ), 1, X( J+1 ), 1 ) 693 END IF 694 ELSE 695* 696* Otherwise, use in-line code for the dot product. 697* 698 IF( UPPER ) THEN 699 DO 110 I = 1, J - 1 700 SUMJ = SUMJ + ( AP( IP-J+I )*USCAL )*X( I ) 701 110 CONTINUE 702 ELSE IF( J.LT.N ) THEN 703 DO 120 I = 1, N - J 704 SUMJ = SUMJ + ( AP( IP+I )*USCAL )*X( J+I ) 705 120 CONTINUE 706 END IF 707 END IF 708* 709 IF( USCAL.EQ.TSCAL ) THEN 710* 711* Compute x(j) := ( x(j) - sumj ) / A(j,j) if 1/A(j,j) 712* was not used to scale the dotproduct. 713* 714 X( J ) = X( J ) - SUMJ 715 XJ = ABS( X( J ) ) 716 IF( NOUNIT ) THEN 717* 718* Compute x(j) = x(j) / A(j,j), scaling if necessary. 719* 720 TJJS = AP( IP )*TSCAL 721 ELSE 722 TJJS = TSCAL 723 IF( TSCAL.EQ.ONE ) 724 $ GO TO 135 725 END IF 726 TJJ = ABS( TJJS ) 727 IF( TJJ.GT.SMLNUM ) THEN 728* 729* abs(A(j,j)) > SMLNUM: 730* 731 IF( TJJ.LT.ONE ) THEN 732 IF( XJ.GT.TJJ*BIGNUM ) THEN 733* 734* Scale X by 1/abs(x(j)). 735* 736 REC = ONE / XJ 737 CALL SSCAL( N, REC, X, 1 ) 738 SCALE = SCALE*REC 739 XMAX = XMAX*REC 740 END IF 741 END IF 742 X( J ) = X( J ) / TJJS 743 ELSE IF( TJJ.GT.ZERO ) THEN 744* 745* 0 < abs(A(j,j)) <= SMLNUM: 746* 747 IF( XJ.GT.TJJ*BIGNUM ) THEN 748* 749* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM. 750* 751 REC = ( TJJ*BIGNUM ) / XJ 752 CALL SSCAL( N, REC, X, 1 ) 753 SCALE = SCALE*REC 754 XMAX = XMAX*REC 755 END IF 756 X( J ) = X( J ) / TJJS 757 ELSE 758* 759* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and 760* scale = 0, and compute a solution to A**T*x = 0. 761* 762 DO 130 I = 1, N 763 X( I ) = ZERO 764 130 CONTINUE 765 X( J ) = ONE 766 SCALE = ZERO 767 XMAX = ZERO 768 END IF 769 135 CONTINUE 770 ELSE 771* 772* Compute x(j) := x(j) / A(j,j) - sumj if the dot 773* product has already been divided by 1/A(j,j). 774* 775 X( J ) = X( J ) / TJJS - SUMJ 776 END IF 777 XMAX = MAX( XMAX, ABS( X( J ) ) ) 778 JLEN = JLEN + 1 779 IP = IP + JINC*JLEN 780 140 CONTINUE 781 END IF 782 SCALE = SCALE / TSCAL 783 END IF 784* 785* Scale the column norms by 1/TSCAL for return. 786* 787 IF( TSCAL.NE.ONE ) THEN 788 CALL SSCAL( N, ONE / TSCAL, CNORM, 1 ) 789 END IF 790* 791 RETURN 792* 793* End of SLATPS 794* 795 END 796