Annotation of rpl/lapack/lapack/dlatbs.f, revision 1.1
1.1 ! bertrand 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlatbs.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack