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