![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, 2: $ CNORM, INFO ) 3: * 4: * -- LAPACK auxiliary routine (version 3.2) -- 5: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 7: * November 2006 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