1 SUBROUTINE DLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, 2 $ CNORM, INFO ) 3* 4* -- LAPACK auxiliary routine (version 3.0) -- 5* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6* Courant Institute, Argonne National Lab, and Rice University 7* June 30, 1992 8* 9* .. Scalar Arguments .. 10 CHARACTER DIAG, NORMIN, TRANS, UPLO 11 INTEGER INFO, LDA, N 12 DOUBLE PRECISION SCALE 13* .. 14* .. Array Arguments .. 15 DOUBLE PRECISION A( LDA, * ), CNORM( * ), X( * ) 16* .. 17* 18* Purpose 19* ======= 20* 21* DLATRS solves one of the triangular systems 22* 23* A *x = s*b or A'*x = s*b 24* 25* with scaling to prevent overflow. Here A is an upper or lower 26* triangular matrix, A' denotes the transpose of A, x and b are 27* n-element vectors, and s is a scaling factor, usually less than 28* or equal to 1, chosen so that the components of x will be less than 29* the overflow threshold. If the unscaled problem will not cause 30* overflow, the Level 2 BLAS routine DTRSV is called. If the matrix A 31* is singular (A(j,j) = 0 for some j), then s is set to 0 and a 32* non-trivial solution to A*x = 0 is returned. 33* 34* Arguments 35* ========= 36* 37* UPLO (input) CHARACTER*1 38* Specifies whether the matrix A is upper or lower triangular. 39* = 'U': Upper triangular 40* = 'L': Lower triangular 41* 42* TRANS (input) CHARACTER*1 43* Specifies the operation applied to A. 44* = 'N': Solve A * x = s*b (No transpose) 45* = 'T': Solve A'* x = s*b (Transpose) 46* = 'C': Solve A'* x = s*b (Conjugate transpose = Transpose) 47* 48* DIAG (input) CHARACTER*1 49* Specifies whether or not the matrix A is unit triangular. 50* = 'N': Non-unit triangular 51* = 'U': Unit triangular 52* 53* NORMIN (input) CHARACTER*1 54* Specifies whether CNORM has been set or not. 55* = 'Y': CNORM contains the column norms on entry 56* = 'N': CNORM is not set on entry. On exit, the norms will 57* be computed and stored in CNORM. 58* 59* N (input) INTEGER 60* The order of the matrix A. N >= 0. 61* 62* A (input) DOUBLE PRECISION array, dimension (LDA,N) 63* The triangular matrix A. If UPLO = 'U', the leading n by n 64* upper triangular part of the array A contains the upper 65* triangular matrix, and the strictly lower triangular part of 66* A is not referenced. If UPLO = 'L', the leading n by n lower 67* triangular part of the array A contains the lower triangular 68* matrix, and the strictly upper triangular part of A is not 69* referenced. If DIAG = 'U', the diagonal elements of A are 70* also not referenced and are assumed to be 1. 71* 72* LDA (input) INTEGER 73* The leading dimension of the array A. LDA >= max (1,N). 74* 75* X (input/output) DOUBLE PRECISION array, dimension (N) 76* On entry, the right hand side b of the triangular system. 77* On exit, X is overwritten by the solution vector x. 78* 79* SCALE (output) DOUBLE PRECISION 80* The scaling factor s for the triangular system 81* A * x = s*b or A'* x = s*b. 82* If SCALE = 0, the matrix A is singular or badly scaled, and 83* the vector x is an exact or approximate solution to A*x = 0. 84* 85* CNORM (input or output) DOUBLE PRECISION array, dimension (N) 86* 87* If NORMIN = 'Y', CNORM is an input argument and CNORM(j) 88* contains the norm of the off-diagonal part of the j-th column 89* of A. If TRANS = 'N', CNORM(j) must be greater than or equal 90* to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) 91* must be greater than or equal to the 1-norm. 92* 93* If NORMIN = 'N', CNORM is an output argument and CNORM(j) 94* returns the 1-norm of the offdiagonal part of the j-th column 95* of A. 96* 97* INFO (output) INTEGER 98* = 0: successful exit 99* < 0: if INFO = -k, the k-th argument had an illegal value 100* 101* Further Details 102* ======= ======= 103* 104* A rough bound on x is computed; if that is less than overflow, DTRSV 105* is called, otherwise, specific code is used which checks for possible 106* overflow or divide-by-zero at every operation. 107* 108* A columnwise scheme is used for solving A*x = b. The basic algorithm 109* if A is lower triangular is 110* 111* x[1:n] := b[1:n] 112* for j = 1, ..., n 113* x(j) := x(j) / A(j,j) 114* x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j] 115* end 116* 117* Define bounds on the components of x after j iterations of the loop: 118* M(j) = bound on x[1:j] 119* G(j) = bound on x[j+1:n] 120* Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}. 121* 122* Then for iteration j+1 we have 123* M(j+1) <= G(j) / | A(j+1,j+1) | 124* G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] | 125* <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | ) 126* 127* where CNORM(j+1) is greater than or equal to the infinity-norm of 128* column j+1 of A, not counting the diagonal. Hence 129* 130* G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | ) 131* 1<=i<=j 132* and 133* 134* |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| ) 135* 1<=i< j 136* 137* Since |x(j)| <= M(j), we use the Level 2 BLAS routine DTRSV if the 138* reciprocal of the largest M(j), j=1,..,n, is larger than 139* max(underflow, 1/overflow). 140* 141* The bound on x(j) is also used to determine when a step in the 142* columnwise method can be performed without fear of overflow. If 143* the computed bound is greater than a large constant, x is scaled to 144* prevent overflow, but if the bound overflows, x is set to 0, x(j) to 145* 1, and scale to 0, and a non-trivial solution to A*x = 0 is found. 146* 147* Similarly, a row-wise scheme is used to solve A'*x = b. The basic 148* algorithm for A upper triangular is 149* 150* for j = 1, ..., n 151* x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j) 152* end 153* 154* We simultaneously compute two bounds 155* G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j 156* M(j) = bound on x(i), 1<=i<=j 157* 158* The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we 159* add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1. 160* Then the bound on x(j) is 161* 162* M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) | 163* 164* <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| ) 165* 1<=i<=j 166* 167* and we can safely call DTRSV if 1/M(n) and 1/G(n) are both greater 168* than max(underflow, 1/overflow). 169* 170* ===================================================================== 171* 172* .. Parameters .. 173 DOUBLE PRECISION ZERO, HALF, ONE 174 PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0 ) 175* .. 176* .. Local Scalars .. 177 LOGICAL NOTRAN, NOUNIT, UPPER 178 INTEGER I, IMAX, J, JFIRST, JINC, JLAST 179 DOUBLE PRECISION BIGNUM, GROW, REC, SMLNUM, SUMJ, TJJ, TJJS, 180 $ TMAX, TSCAL, USCAL, XBND, XJ, XMAX 181* .. 182* .. External Functions .. 183 LOGICAL LSAME 184 INTEGER IDAMAX 185 DOUBLE PRECISION DASUM, DDOT, DLAMCH 186 EXTERNAL LSAME, IDAMAX, DASUM, DDOT, DLAMCH 187* .. 188* .. External Subroutines .. 189 EXTERNAL DAXPY, DSCAL, DTRSV, XERBLA 190* .. 191* .. Intrinsic Functions .. 192 INTRINSIC ABS, MAX, MIN 193* .. 194* .. Executable Statements .. 195* 196 INFO = 0 197 UPPER = LSAME( UPLO, 'U' ) 198 NOTRAN = LSAME( TRANS, 'N' ) 199 NOUNIT = LSAME( DIAG, 'N' ) 200* 201* Test the input parameters. 202* 203 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 204 INFO = -1 205 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. 206 $ LSAME( TRANS, 'C' ) ) THEN 207 INFO = -2 208 ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN 209 INFO = -3 210 ELSE IF( .NOT.LSAME( NORMIN, 'Y' ) .AND. .NOT. 211 $ LSAME( NORMIN, 'N' ) ) THEN 212 INFO = -4 213 ELSE IF( N.LT.0 ) THEN 214 INFO = -5 215 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 216 INFO = -7 217 END IF 218 IF( INFO.NE.0 ) THEN 219 CALL XERBLA( 'DLATRS', -INFO ) 220 RETURN 221 END IF 222* 223* Quick return if possible 224* 225 IF( N.EQ.0 ) 226 $ RETURN 227* 228* Determine machine dependent parameters to control overflow. 229* 230 SMLNUM = DLAMCH( 'Safe minimum' ) / DLAMCH( 'Precision' ) 231 BIGNUM = ONE / SMLNUM 232 SCALE = ONE 233* 234 IF( LSAME( NORMIN, 'N' ) ) THEN 235* 236* Compute the 1-norm of each column, not including the diagonal. 237* 238 IF( UPPER ) THEN 239* 240* A is upper triangular. 241* 242 DO 10 J = 1, N 243 CNORM( J ) = DASUM( J-1, A( 1, J ), 1 ) 244 10 CONTINUE 245 ELSE 246* 247* A is lower triangular. 248* 249 DO 20 J = 1, N - 1 250 CNORM( J ) = DASUM( N-J, A( J+1, J ), 1 ) 251 20 CONTINUE 252 CNORM( N ) = ZERO 253 END IF 254 END IF 255* 256* Scale the column norms by TSCAL if the maximum element in CNORM is 257* greater than BIGNUM. 258* 259 IMAX = IDAMAX( N, CNORM, 1 ) 260 TMAX = CNORM( IMAX ) 261 IF( TMAX.LE.BIGNUM ) THEN 262 TSCAL = ONE 263 ELSE 264 TSCAL = ONE / ( SMLNUM*TMAX ) 265 CALL DSCAL( N, TSCAL, CNORM, 1 ) 266 END IF 267* 268* Compute a bound on the computed solution vector to see if the 269* Level 2 BLAS routine DTRSV can be used. 270* 271 J = IDAMAX( N, X, 1 ) 272 XMAX = ABS( X( J ) ) 273 XBND = XMAX 274 IF( NOTRAN ) THEN 275* 276* Compute the growth in A * x = b. 277* 278 IF( UPPER ) THEN 279 JFIRST = N 280 JLAST = 1 281 JINC = -1 282 ELSE 283 JFIRST = 1 284 JLAST = N 285 JINC = 1 286 END IF 287* 288 IF( TSCAL.NE.ONE ) THEN 289 GROW = ZERO 290 GO TO 50 291 END IF 292* 293 IF( NOUNIT ) THEN 294* 295* A is non-unit triangular. 296* 297* Compute GROW = 1/G(j) and XBND = 1/M(j). 298* Initially, G(0) = max{x(i), i=1,...,n}. 299* 300 GROW = ONE / MAX( XBND, SMLNUM ) 301 XBND = GROW 302 DO 30 J = JFIRST, JLAST, JINC 303* 304* Exit the loop if the growth factor is too small. 305* 306 IF( GROW.LE.SMLNUM ) 307 $ GO TO 50 308* 309* M(j) = G(j-1) / abs(A(j,j)) 310* 311 TJJ = ABS( A( J, J ) ) 312 XBND = MIN( XBND, MIN( ONE, TJJ )*GROW ) 313 IF( TJJ+CNORM( J ).GE.SMLNUM ) THEN 314* 315* G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) ) 316* 317 GROW = GROW*( TJJ / ( TJJ+CNORM( J ) ) ) 318 ELSE 319* 320* G(j) could overflow, set GROW to 0. 321* 322 GROW = ZERO 323 END IF 324 30 CONTINUE 325 GROW = XBND 326 ELSE 327* 328* A is unit triangular. 329* 330* Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}. 331* 332 GROW = MIN( ONE, ONE / MAX( XBND, SMLNUM ) ) 333 DO 40 J = JFIRST, JLAST, JINC 334* 335* Exit the loop if the growth factor is too small. 336* 337 IF( GROW.LE.SMLNUM ) 338 $ GO TO 50 339* 340* G(j) = G(j-1)*( 1 + CNORM(j) ) 341* 342 GROW = GROW*( ONE / ( ONE+CNORM( J ) ) ) 343 40 CONTINUE 344 END IF 345 50 CONTINUE 346* 347 ELSE 348* 349* Compute the growth in A' * x = b. 350* 351 IF( UPPER ) THEN 352 JFIRST = 1 353 JLAST = N 354 JINC = 1 355 ELSE 356 JFIRST = N 357 JLAST = 1 358 JINC = -1 359 END IF 360* 361 IF( TSCAL.NE.ONE ) THEN 362 GROW = ZERO 363 GO TO 80 364 END IF 365* 366 IF( NOUNIT ) THEN 367* 368* A is non-unit triangular. 369* 370* Compute GROW = 1/G(j) and XBND = 1/M(j). 371* Initially, M(0) = max{x(i), i=1,...,n}. 372* 373 GROW = ONE / MAX( XBND, SMLNUM ) 374 XBND = GROW 375 DO 60 J = JFIRST, JLAST, JINC 376* 377* Exit the loop if the growth factor is too small. 378* 379 IF( GROW.LE.SMLNUM ) 380 $ GO TO 80 381* 382* G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) ) 383* 384 XJ = ONE + CNORM( J ) 385 GROW = MIN( GROW, XBND / XJ ) 386* 387* M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j)) 388* 389 TJJ = ABS( A( J, J ) ) 390 IF( XJ.GT.TJJ ) 391 $ XBND = XBND*( TJJ / XJ ) 392 60 CONTINUE 393 GROW = MIN( GROW, XBND ) 394 ELSE 395* 396* A is unit triangular. 397* 398* Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}. 399* 400 GROW = MIN( ONE, ONE / MAX( XBND, SMLNUM ) ) 401 DO 70 J = JFIRST, JLAST, JINC 402* 403* Exit the loop if the growth factor is too small. 404* 405 IF( GROW.LE.SMLNUM ) 406 $ GO TO 80 407* 408* G(j) = ( 1 + CNORM(j) )*G(j-1) 409* 410 XJ = ONE + CNORM( J ) 411 GROW = GROW / XJ 412 70 CONTINUE 413 END IF 414 80 CONTINUE 415 END IF 416* 417 IF( ( GROW*TSCAL ).GT.SMLNUM ) THEN 418* 419* Use the Level 2 BLAS solve if the reciprocal of the bound on 420* elements of X is not too small. 421* 422 CALL DTRSV( UPLO, TRANS, DIAG, N, A, LDA, X, 1 ) 423 ELSE 424* 425* Use a Level 1 BLAS solve, scaling intermediate results. 426* 427 IF( XMAX.GT.BIGNUM ) THEN 428* 429* Scale X so that its components are less than or equal to 430* BIGNUM in absolute value. 431* 432 SCALE = BIGNUM / XMAX 433 CALL DSCAL( N, SCALE, X, 1 ) 434 XMAX = BIGNUM 435 END IF 436* 437 IF( NOTRAN ) THEN 438* 439* Solve A * x = b 440* 441 DO 110 J = JFIRST, JLAST, JINC 442* 443* Compute x(j) = b(j) / A(j,j), scaling x if necessary. 444* 445 XJ = ABS( X( J ) ) 446 IF( NOUNIT ) THEN 447 TJJS = A( J, J )*TSCAL 448 ELSE 449 TJJS = TSCAL 450 IF( TSCAL.EQ.ONE ) 451 $ GO TO 100 452 END IF 453 TJJ = ABS( TJJS ) 454 IF( TJJ.GT.SMLNUM ) THEN 455* 456* abs(A(j,j)) > SMLNUM: 457* 458 IF( TJJ.LT.ONE ) THEN 459 IF( XJ.GT.TJJ*BIGNUM ) THEN 460* 461* Scale x by 1/b(j). 462* 463 REC = ONE / XJ 464 CALL DSCAL( N, REC, X, 1 ) 465 SCALE = SCALE*REC 466 XMAX = XMAX*REC 467 END IF 468 END IF 469 X( J ) = X( J ) / TJJS 470 XJ = ABS( X( J ) ) 471 ELSE IF( TJJ.GT.ZERO ) THEN 472* 473* 0 < abs(A(j,j)) <= SMLNUM: 474* 475 IF( XJ.GT.TJJ*BIGNUM ) THEN 476* 477* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM 478* to avoid overflow when dividing by A(j,j). 479* 480 REC = ( TJJ*BIGNUM ) / XJ 481 IF( CNORM( J ).GT.ONE ) THEN 482* 483* Scale by 1/CNORM(j) to avoid overflow when 484* multiplying x(j) times column j. 485* 486 REC = REC / CNORM( J ) 487 END IF 488 CALL DSCAL( N, REC, X, 1 ) 489 SCALE = SCALE*REC 490 XMAX = XMAX*REC 491 END IF 492 X( J ) = X( J ) / TJJS 493 XJ = ABS( X( J ) ) 494 ELSE 495* 496* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and 497* scale = 0, and compute a solution to A*x = 0. 498* 499 DO 90 I = 1, N 500 X( I ) = ZERO 501 90 CONTINUE 502 X( J ) = ONE 503 XJ = ONE 504 SCALE = ZERO 505 XMAX = ZERO 506 END IF 507 100 CONTINUE 508* 509* Scale x if necessary to avoid overflow when adding a 510* multiple of column j of A. 511* 512 IF( XJ.GT.ONE ) THEN 513 REC = ONE / XJ 514 IF( CNORM( J ).GT.( BIGNUM-XMAX )*REC ) THEN 515* 516* Scale x by 1/(2*abs(x(j))). 517* 518 REC = REC*HALF 519 CALL DSCAL( N, REC, X, 1 ) 520 SCALE = SCALE*REC 521 END IF 522 ELSE IF( XJ*CNORM( J ).GT.( BIGNUM-XMAX ) ) THEN 523* 524* Scale x by 1/2. 525* 526 CALL DSCAL( N, HALF, X, 1 ) 527 SCALE = SCALE*HALF 528 END IF 529* 530 IF( UPPER ) THEN 531 IF( J.GT.1 ) THEN 532* 533* Compute the update 534* x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j) 535* 536 CALL DAXPY( J-1, -X( J )*TSCAL, A( 1, J ), 1, X, 537 $ 1 ) 538 I = IDAMAX( J-1, X, 1 ) 539 XMAX = ABS( X( I ) ) 540 END IF 541 ELSE 542 IF( J.LT.N ) THEN 543* 544* Compute the update 545* x(j+1:n) := x(j+1:n) - x(j) * A(j+1:n,j) 546* 547 CALL DAXPY( N-J, -X( J )*TSCAL, A( J+1, J ), 1, 548 $ X( J+1 ), 1 ) 549 I = J + IDAMAX( N-J, X( J+1 ), 1 ) 550 XMAX = ABS( X( I ) ) 551 END IF 552 END IF 553 110 CONTINUE 554* 555 ELSE 556* 557* Solve A' * x = b 558* 559 DO 160 J = JFIRST, JLAST, JINC 560* 561* Compute x(j) = b(j) - sum A(k,j)*x(k). 562* k<>j 563* 564 XJ = ABS( X( J ) ) 565 USCAL = TSCAL 566 REC = ONE / MAX( XMAX, ONE ) 567 IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN 568* 569* If x(j) could overflow, scale x by 1/(2*XMAX). 570* 571 REC = REC*HALF 572 IF( NOUNIT ) THEN 573 TJJS = A( J, J )*TSCAL 574 ELSE 575 TJJS = TSCAL 576 END IF 577 TJJ = ABS( TJJS ) 578 IF( TJJ.GT.ONE ) THEN 579* 580* Divide by A(j,j) when scaling x if A(j,j) > 1. 581* 582 REC = MIN( ONE, REC*TJJ ) 583 USCAL = USCAL / TJJS 584 END IF 585 IF( REC.LT.ONE ) THEN 586 CALL DSCAL( N, REC, X, 1 ) 587 SCALE = SCALE*REC 588 XMAX = XMAX*REC 589 END IF 590 END IF 591* 592 SUMJ = ZERO 593 IF( USCAL.EQ.ONE ) THEN 594* 595* If the scaling needed for A in the dot product is 1, 596* call DDOT to perform the dot product. 597* 598 IF( UPPER ) THEN 599 SUMJ = DDOT( J-1, A( 1, J ), 1, X, 1 ) 600 ELSE IF( J.LT.N ) THEN 601 SUMJ = DDOT( N-J, A( J+1, J ), 1, X( J+1 ), 1 ) 602 END IF 603 ELSE 604* 605* Otherwise, use in-line code for the dot product. 606* 607 IF( UPPER ) THEN 608 DO 120 I = 1, J - 1 609 SUMJ = SUMJ + ( A( I, J )*USCAL )*X( I ) 610 120 CONTINUE 611 ELSE IF( J.LT.N ) THEN 612 DO 130 I = J + 1, N 613 SUMJ = SUMJ + ( A( I, J )*USCAL )*X( I ) 614 130 CONTINUE 615 END IF 616 END IF 617* 618 IF( USCAL.EQ.TSCAL ) THEN 619* 620* Compute x(j) := ( x(j) - sumj ) / A(j,j) if 1/A(j,j) 621* was not used to scale the dotproduct. 622* 623 X( J ) = X( J ) - SUMJ 624 XJ = ABS( X( J ) ) 625 IF( NOUNIT ) THEN 626 TJJS = A( J, J )*TSCAL 627 ELSE 628 TJJS = TSCAL 629 IF( TSCAL.EQ.ONE ) 630 $ GO TO 150 631 END IF 632* 633* Compute x(j) = x(j) / A(j,j), scaling if necessary. 634* 635 TJJ = ABS( TJJS ) 636 IF( TJJ.GT.SMLNUM ) THEN 637* 638* abs(A(j,j)) > SMLNUM: 639* 640 IF( TJJ.LT.ONE ) THEN 641 IF( XJ.GT.TJJ*BIGNUM ) THEN 642* 643* Scale X by 1/abs(x(j)). 644* 645 REC = ONE / XJ 646 CALL DSCAL( N, REC, X, 1 ) 647 SCALE = SCALE*REC 648 XMAX = XMAX*REC 649 END IF 650 END IF 651 X( J ) = X( J ) / TJJS 652 ELSE IF( TJJ.GT.ZERO ) THEN 653* 654* 0 < abs(A(j,j)) <= SMLNUM: 655* 656 IF( XJ.GT.TJJ*BIGNUM ) THEN 657* 658* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM. 659* 660 REC = ( TJJ*BIGNUM ) / XJ 661 CALL DSCAL( N, REC, X, 1 ) 662 SCALE = SCALE*REC 663 XMAX = XMAX*REC 664 END IF 665 X( J ) = X( J ) / TJJS 666 ELSE 667* 668* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and 669* scale = 0, and compute a solution to A'*x = 0. 670* 671 DO 140 I = 1, N 672 X( I ) = ZERO 673 140 CONTINUE 674 X( J ) = ONE 675 SCALE = ZERO 676 XMAX = ZERO 677 END IF 678 150 CONTINUE 679 ELSE 680* 681* Compute x(j) := x(j) / A(j,j) - sumj if the dot 682* product has already been divided by 1/A(j,j). 683* 684 X( J ) = X( J ) / TJJS - SUMJ 685 END IF 686 XMAX = MAX( XMAX, ABS( X( J ) ) ) 687 160 CONTINUE 688 END IF 689 SCALE = SCALE / TSCAL 690 END IF 691* 692* Scale the column norms by 1/TSCAL for return. 693* 694 IF( TSCAL.NE.ONE ) THEN 695 CALL DSCAL( N, ONE / TSCAL, CNORM, 1 ) 696 END IF 697* 698 RETURN 699* 700* End of DLATRS 701* 702 END 703