Annotation of rpl/lapack/lapack/zlatbs.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE ZLATBS( 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 CNORM( * )
! 16: COMPLEX*16 AB( LDAB, * ), X( * )
! 17: * ..
! 18: *
! 19: * Purpose
! 20: * =======
! 21: *
! 22: * ZLATBS solves one of the triangular systems
! 23: *
! 24: * A * x = s*b, A**T * x = s*b, or A**H * x = s*b,
! 25: *
! 26: * with scaling to prevent overflow, where A is an upper or lower
! 27: * triangular band matrix. Here A' denotes the transpose of A, x and b
! 28: * are n-element vectors, and s is a scaling factor, usually less than
! 29: * or equal to 1, chosen so that the components of x will be less than
! 30: * the overflow threshold. If the unscaled problem will not cause
! 31: * overflow, the Level 2 BLAS routine ZTBSV is called. If the matrix A
! 32: * is singular (A(j,j) = 0 for some j), then s is set to 0 and a
! 33: * non-trivial solution to A*x = 0 is returned.
! 34: *
! 35: * Arguments
! 36: * =========
! 37: *
! 38: * UPLO (input) CHARACTER*1
! 39: * Specifies whether the matrix A is upper or lower triangular.
! 40: * = 'U': Upper triangular
! 41: * = 'L': Lower triangular
! 42: *
! 43: * TRANS (input) CHARACTER*1
! 44: * Specifies the operation applied to A.
! 45: * = 'N': Solve A * x = s*b (No transpose)
! 46: * = 'T': Solve A**T * x = s*b (Transpose)
! 47: * = 'C': Solve A**H * x = s*b (Conjugate transpose)
! 48: *
! 49: * DIAG (input) CHARACTER*1
! 50: * Specifies whether or not the matrix A is unit triangular.
! 51: * = 'N': Non-unit triangular
! 52: * = 'U': Unit triangular
! 53: *
! 54: * NORMIN (input) CHARACTER*1
! 55: * Specifies whether CNORM has been set or not.
! 56: * = 'Y': CNORM contains the column norms on entry
! 57: * = 'N': CNORM is not set on entry. On exit, the norms will
! 58: * be computed and stored in CNORM.
! 59: *
! 60: * N (input) INTEGER
! 61: * The order of the matrix A. N >= 0.
! 62: *
! 63: * KD (input) INTEGER
! 64: * The number of subdiagonals or superdiagonals in the
! 65: * triangular matrix A. KD >= 0.
! 66: *
! 67: * AB (input) COMPLEX*16 array, dimension (LDAB,N)
! 68: * The upper or lower triangular band matrix A, stored in the
! 69: * first KD+1 rows of the array. The j-th column of A is stored
! 70: * in the j-th column of the array AB as follows:
! 71: * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
! 72: * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
! 73: *
! 74: * LDAB (input) INTEGER
! 75: * The leading dimension of the array AB. LDAB >= KD+1.
! 76: *
! 77: * X (input/output) COMPLEX*16 array, dimension (N)
! 78: * On entry, the right hand side b of the triangular system.
! 79: * On exit, X is overwritten by the solution vector x.
! 80: *
! 81: * SCALE (output) DOUBLE PRECISION
! 82: * The scaling factor s for the triangular system
! 83: * A * x = s*b, A**T * x = s*b, or A**H * x = s*b.
! 84: * If SCALE = 0, the matrix A is singular or badly scaled, and
! 85: * the vector x is an exact or approximate solution to A*x = 0.
! 86: *
! 87: * CNORM (input or output) DOUBLE PRECISION array, dimension (N)
! 88: *
! 89: * If NORMIN = 'Y', CNORM is an input argument and CNORM(j)
! 90: * contains the norm of the off-diagonal part of the j-th column
! 91: * of A. If TRANS = 'N', CNORM(j) must be greater than or equal
! 92: * to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)
! 93: * must be greater than or equal to the 1-norm.
! 94: *
! 95: * If NORMIN = 'N', CNORM is an output argument and CNORM(j)
! 96: * returns the 1-norm of the offdiagonal part of the j-th column
! 97: * of A.
! 98: *
! 99: * INFO (output) INTEGER
! 100: * = 0: successful exit
! 101: * < 0: if INFO = -k, the k-th argument had an illegal value
! 102: *
! 103: * Further Details
! 104: * ======= =======
! 105: *
! 106: * A rough bound on x is computed; if that is less than overflow, ZTBSV
! 107: * is called, otherwise, specific code is used which checks for possible
! 108: * overflow or divide-by-zero at every operation.
! 109: *
! 110: * A columnwise scheme is used for solving A*x = b. The basic algorithm
! 111: * if A is lower triangular is
! 112: *
! 113: * x[1:n] := b[1:n]
! 114: * for j = 1, ..., n
! 115: * x(j) := x(j) / A(j,j)
! 116: * x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
! 117: * end
! 118: *
! 119: * Define bounds on the components of x after j iterations of the loop:
! 120: * M(j) = bound on x[1:j]
! 121: * G(j) = bound on x[j+1:n]
! 122: * Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
! 123: *
! 124: * Then for iteration j+1 we have
! 125: * M(j+1) <= G(j) / | A(j+1,j+1) |
! 126: * G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
! 127: * <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
! 128: *
! 129: * where CNORM(j+1) is greater than or equal to the infinity-norm of
! 130: * column j+1 of A, not counting the diagonal. Hence
! 131: *
! 132: * G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
! 133: * 1<=i<=j
! 134: * and
! 135: *
! 136: * |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
! 137: * 1<=i< j
! 138: *
! 139: * Since |x(j)| <= M(j), we use the Level 2 BLAS routine ZTBSV if the
! 140: * reciprocal of the largest M(j), j=1,..,n, is larger than
! 141: * max(underflow, 1/overflow).
! 142: *
! 143: * The bound on x(j) is also used to determine when a step in the
! 144: * columnwise method can be performed without fear of overflow. If
! 145: * the computed bound is greater than a large constant, x is scaled to
! 146: * prevent overflow, but if the bound overflows, x is set to 0, x(j) to
! 147: * 1, and scale to 0, and a non-trivial solution to A*x = 0 is found.
! 148: *
! 149: * Similarly, a row-wise scheme is used to solve A**T *x = b or
! 150: * A**H *x = b. The basic algorithm for A upper triangular is
! 151: *
! 152: * for j = 1, ..., n
! 153: * x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)
! 154: * end
! 155: *
! 156: * We simultaneously compute two bounds
! 157: * G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j
! 158: * M(j) = bound on x(i), 1<=i<=j
! 159: *
! 160: * The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we
! 161: * add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.
! 162: * Then the bound on x(j) is
! 163: *
! 164: * M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
! 165: *
! 166: * <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
! 167: * 1<=i<=j
! 168: *
! 169: * and we can safely call ZTBSV if 1/M(n) and 1/G(n) are both greater
! 170: * than max(underflow, 1/overflow).
! 171: *
! 172: * =====================================================================
! 173: *
! 174: * .. Parameters ..
! 175: DOUBLE PRECISION ZERO, HALF, ONE, TWO
! 176: PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0,
! 177: $ TWO = 2.0D+0 )
! 178: * ..
! 179: * .. Local Scalars ..
! 180: LOGICAL NOTRAN, NOUNIT, UPPER
! 181: INTEGER I, IMAX, J, JFIRST, JINC, JLAST, JLEN, MAIND
! 182: DOUBLE PRECISION BIGNUM, GROW, REC, SMLNUM, TJJ, TMAX, TSCAL,
! 183: $ XBND, XJ, XMAX
! 184: COMPLEX*16 CSUMJ, TJJS, USCAL, ZDUM
! 185: * ..
! 186: * .. External Functions ..
! 187: LOGICAL LSAME
! 188: INTEGER IDAMAX, IZAMAX
! 189: DOUBLE PRECISION DLAMCH, DZASUM
! 190: COMPLEX*16 ZDOTC, ZDOTU, ZLADIV
! 191: EXTERNAL LSAME, IDAMAX, IZAMAX, DLAMCH, DZASUM, ZDOTC,
! 192: $ ZDOTU, ZLADIV
! 193: * ..
! 194: * .. External Subroutines ..
! 195: EXTERNAL DSCAL, XERBLA, ZAXPY, ZDSCAL, ZTBSV
! 196: * ..
! 197: * .. Intrinsic Functions ..
! 198: INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN
! 199: * ..
! 200: * .. Statement Functions ..
! 201: DOUBLE PRECISION CABS1, CABS2
! 202: * ..
! 203: * .. Statement Function definitions ..
! 204: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
! 205: CABS2( ZDUM ) = ABS( DBLE( ZDUM ) / 2.D0 ) +
! 206: $ ABS( DIMAG( ZDUM ) / 2.D0 )
! 207: * ..
! 208: * .. Executable Statements ..
! 209: *
! 210: INFO = 0
! 211: UPPER = LSAME( UPLO, 'U' )
! 212: NOTRAN = LSAME( TRANS, 'N' )
! 213: NOUNIT = LSAME( DIAG, 'N' )
! 214: *
! 215: * Test the input parameters.
! 216: *
! 217: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
! 218: INFO = -1
! 219: ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
! 220: $ LSAME( TRANS, 'C' ) ) THEN
! 221: INFO = -2
! 222: ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
! 223: INFO = -3
! 224: ELSE IF( .NOT.LSAME( NORMIN, 'Y' ) .AND. .NOT.
! 225: $ LSAME( NORMIN, 'N' ) ) THEN
! 226: INFO = -4
! 227: ELSE IF( N.LT.0 ) THEN
! 228: INFO = -5
! 229: ELSE IF( KD.LT.0 ) THEN
! 230: INFO = -6
! 231: ELSE IF( LDAB.LT.KD+1 ) THEN
! 232: INFO = -8
! 233: END IF
! 234: IF( INFO.NE.0 ) THEN
! 235: CALL XERBLA( 'ZLATBS', -INFO )
! 236: RETURN
! 237: END IF
! 238: *
! 239: * Quick return if possible
! 240: *
! 241: IF( N.EQ.0 )
! 242: $ RETURN
! 243: *
! 244: * Determine machine dependent parameters to control overflow.
! 245: *
! 246: SMLNUM = DLAMCH( 'Safe minimum' )
! 247: BIGNUM = ONE / SMLNUM
! 248: CALL DLABAD( SMLNUM, BIGNUM )
! 249: SMLNUM = SMLNUM / DLAMCH( 'Precision' )
! 250: BIGNUM = ONE / SMLNUM
! 251: SCALE = ONE
! 252: *
! 253: IF( LSAME( NORMIN, 'N' ) ) THEN
! 254: *
! 255: * Compute the 1-norm of each column, not including the diagonal.
! 256: *
! 257: IF( UPPER ) THEN
! 258: *
! 259: * A is upper triangular.
! 260: *
! 261: DO 10 J = 1, N
! 262: JLEN = MIN( KD, J-1 )
! 263: CNORM( J ) = DZASUM( JLEN, AB( KD+1-JLEN, J ), 1 )
! 264: 10 CONTINUE
! 265: ELSE
! 266: *
! 267: * A is lower triangular.
! 268: *
! 269: DO 20 J = 1, N
! 270: JLEN = MIN( KD, N-J )
! 271: IF( JLEN.GT.0 ) THEN
! 272: CNORM( J ) = DZASUM( JLEN, AB( 2, J ), 1 )
! 273: ELSE
! 274: CNORM( J ) = ZERO
! 275: END IF
! 276: 20 CONTINUE
! 277: END IF
! 278: END IF
! 279: *
! 280: * Scale the column norms by TSCAL if the maximum element in CNORM is
! 281: * greater than BIGNUM/2.
! 282: *
! 283: IMAX = IDAMAX( N, CNORM, 1 )
! 284: TMAX = CNORM( IMAX )
! 285: IF( TMAX.LE.BIGNUM*HALF ) THEN
! 286: TSCAL = ONE
! 287: ELSE
! 288: TSCAL = HALF / ( SMLNUM*TMAX )
! 289: CALL DSCAL( N, TSCAL, CNORM, 1 )
! 290: END IF
! 291: *
! 292: * Compute a bound on the computed solution vector to see if the
! 293: * Level 2 BLAS routine ZTBSV can be used.
! 294: *
! 295: XMAX = ZERO
! 296: DO 30 J = 1, N
! 297: XMAX = MAX( XMAX, CABS2( X( J ) ) )
! 298: 30 CONTINUE
! 299: XBND = XMAX
! 300: IF( NOTRAN ) THEN
! 301: *
! 302: * Compute the growth in A * x = b.
! 303: *
! 304: IF( UPPER ) THEN
! 305: JFIRST = N
! 306: JLAST = 1
! 307: JINC = -1
! 308: MAIND = KD + 1
! 309: ELSE
! 310: JFIRST = 1
! 311: JLAST = N
! 312: JINC = 1
! 313: MAIND = 1
! 314: END IF
! 315: *
! 316: IF( TSCAL.NE.ONE ) THEN
! 317: GROW = ZERO
! 318: GO TO 60
! 319: END IF
! 320: *
! 321: IF( NOUNIT ) THEN
! 322: *
! 323: * A is non-unit triangular.
! 324: *
! 325: * Compute GROW = 1/G(j) and XBND = 1/M(j).
! 326: * Initially, G(0) = max{x(i), i=1,...,n}.
! 327: *
! 328: GROW = HALF / MAX( XBND, SMLNUM )
! 329: XBND = GROW
! 330: DO 40 J = JFIRST, JLAST, JINC
! 331: *
! 332: * Exit the loop if the growth factor is too small.
! 333: *
! 334: IF( GROW.LE.SMLNUM )
! 335: $ GO TO 60
! 336: *
! 337: TJJS = AB( MAIND, J )
! 338: TJJ = CABS1( TJJS )
! 339: *
! 340: IF( TJJ.GE.SMLNUM ) THEN
! 341: *
! 342: * M(j) = G(j-1) / abs(A(j,j))
! 343: *
! 344: XBND = MIN( XBND, MIN( ONE, TJJ )*GROW )
! 345: ELSE
! 346: *
! 347: * M(j) could overflow, set XBND to 0.
! 348: *
! 349: XBND = ZERO
! 350: END IF
! 351: *
! 352: IF( TJJ+CNORM( J ).GE.SMLNUM ) THEN
! 353: *
! 354: * G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
! 355: *
! 356: GROW = GROW*( TJJ / ( TJJ+CNORM( J ) ) )
! 357: ELSE
! 358: *
! 359: * G(j) could overflow, set GROW to 0.
! 360: *
! 361: GROW = ZERO
! 362: END IF
! 363: 40 CONTINUE
! 364: GROW = XBND
! 365: ELSE
! 366: *
! 367: * A is unit triangular.
! 368: *
! 369: * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
! 370: *
! 371: GROW = MIN( ONE, HALF / MAX( XBND, SMLNUM ) )
! 372: DO 50 J = JFIRST, JLAST, JINC
! 373: *
! 374: * Exit the loop if the growth factor is too small.
! 375: *
! 376: IF( GROW.LE.SMLNUM )
! 377: $ GO TO 60
! 378: *
! 379: * G(j) = G(j-1)*( 1 + CNORM(j) )
! 380: *
! 381: GROW = GROW*( ONE / ( ONE+CNORM( J ) ) )
! 382: 50 CONTINUE
! 383: END IF
! 384: 60 CONTINUE
! 385: *
! 386: ELSE
! 387: *
! 388: * Compute the growth in A**T * x = b or A**H * x = b.
! 389: *
! 390: IF( UPPER ) THEN
! 391: JFIRST = 1
! 392: JLAST = N
! 393: JINC = 1
! 394: MAIND = KD + 1
! 395: ELSE
! 396: JFIRST = N
! 397: JLAST = 1
! 398: JINC = -1
! 399: MAIND = 1
! 400: END IF
! 401: *
! 402: IF( TSCAL.NE.ONE ) THEN
! 403: GROW = ZERO
! 404: GO TO 90
! 405: END IF
! 406: *
! 407: IF( NOUNIT ) THEN
! 408: *
! 409: * A is non-unit triangular.
! 410: *
! 411: * Compute GROW = 1/G(j) and XBND = 1/M(j).
! 412: * Initially, M(0) = max{x(i), i=1,...,n}.
! 413: *
! 414: GROW = HALF / MAX( XBND, SMLNUM )
! 415: XBND = GROW
! 416: DO 70 J = JFIRST, JLAST, JINC
! 417: *
! 418: * Exit the loop if the growth factor is too small.
! 419: *
! 420: IF( GROW.LE.SMLNUM )
! 421: $ GO TO 90
! 422: *
! 423: * G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
! 424: *
! 425: XJ = ONE + CNORM( J )
! 426: GROW = MIN( GROW, XBND / XJ )
! 427: *
! 428: TJJS = AB( MAIND, J )
! 429: TJJ = CABS1( TJJS )
! 430: *
! 431: IF( TJJ.GE.SMLNUM ) THEN
! 432: *
! 433: * M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
! 434: *
! 435: IF( XJ.GT.TJJ )
! 436: $ XBND = XBND*( TJJ / XJ )
! 437: ELSE
! 438: *
! 439: * M(j) could overflow, set XBND to 0.
! 440: *
! 441: XBND = ZERO
! 442: END IF
! 443: 70 CONTINUE
! 444: GROW = MIN( GROW, XBND )
! 445: ELSE
! 446: *
! 447: * A is unit triangular.
! 448: *
! 449: * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
! 450: *
! 451: GROW = MIN( ONE, HALF / MAX( XBND, SMLNUM ) )
! 452: DO 80 J = JFIRST, JLAST, JINC
! 453: *
! 454: * Exit the loop if the growth factor is too small.
! 455: *
! 456: IF( GROW.LE.SMLNUM )
! 457: $ GO TO 90
! 458: *
! 459: * G(j) = ( 1 + CNORM(j) )*G(j-1)
! 460: *
! 461: XJ = ONE + CNORM( J )
! 462: GROW = GROW / XJ
! 463: 80 CONTINUE
! 464: END IF
! 465: 90 CONTINUE
! 466: END IF
! 467: *
! 468: IF( ( GROW*TSCAL ).GT.SMLNUM ) THEN
! 469: *
! 470: * Use the Level 2 BLAS solve if the reciprocal of the bound on
! 471: * elements of X is not too small.
! 472: *
! 473: CALL ZTBSV( UPLO, TRANS, DIAG, N, KD, AB, LDAB, X, 1 )
! 474: ELSE
! 475: *
! 476: * Use a Level 1 BLAS solve, scaling intermediate results.
! 477: *
! 478: IF( XMAX.GT.BIGNUM*HALF ) THEN
! 479: *
! 480: * Scale X so that its components are less than or equal to
! 481: * BIGNUM in absolute value.
! 482: *
! 483: SCALE = ( BIGNUM*HALF ) / XMAX
! 484: CALL ZDSCAL( N, SCALE, X, 1 )
! 485: XMAX = BIGNUM
! 486: ELSE
! 487: XMAX = XMAX*TWO
! 488: END IF
! 489: *
! 490: IF( NOTRAN ) THEN
! 491: *
! 492: * Solve A * x = b
! 493: *
! 494: DO 120 J = JFIRST, JLAST, JINC
! 495: *
! 496: * Compute x(j) = b(j) / A(j,j), scaling x if necessary.
! 497: *
! 498: XJ = CABS1( X( J ) )
! 499: IF( NOUNIT ) THEN
! 500: TJJS = AB( MAIND, J )*TSCAL
! 501: ELSE
! 502: TJJS = TSCAL
! 503: IF( TSCAL.EQ.ONE )
! 504: $ GO TO 110
! 505: END IF
! 506: TJJ = CABS1( TJJS )
! 507: IF( TJJ.GT.SMLNUM ) THEN
! 508: *
! 509: * abs(A(j,j)) > SMLNUM:
! 510: *
! 511: IF( TJJ.LT.ONE ) THEN
! 512: IF( XJ.GT.TJJ*BIGNUM ) THEN
! 513: *
! 514: * Scale x by 1/b(j).
! 515: *
! 516: REC = ONE / XJ
! 517: CALL ZDSCAL( N, REC, X, 1 )
! 518: SCALE = SCALE*REC
! 519: XMAX = XMAX*REC
! 520: END IF
! 521: END IF
! 522: X( J ) = ZLADIV( X( J ), TJJS )
! 523: XJ = CABS1( X( J ) )
! 524: ELSE IF( TJJ.GT.ZERO ) THEN
! 525: *
! 526: * 0 < abs(A(j,j)) <= SMLNUM:
! 527: *
! 528: IF( XJ.GT.TJJ*BIGNUM ) THEN
! 529: *
! 530: * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
! 531: * to avoid overflow when dividing by A(j,j).
! 532: *
! 533: REC = ( TJJ*BIGNUM ) / XJ
! 534: IF( CNORM( J ).GT.ONE ) THEN
! 535: *
! 536: * Scale by 1/CNORM(j) to avoid overflow when
! 537: * multiplying x(j) times column j.
! 538: *
! 539: REC = REC / CNORM( J )
! 540: END IF
! 541: CALL ZDSCAL( N, REC, X, 1 )
! 542: SCALE = SCALE*REC
! 543: XMAX = XMAX*REC
! 544: END IF
! 545: X( J ) = ZLADIV( X( J ), TJJS )
! 546: XJ = CABS1( X( J ) )
! 547: ELSE
! 548: *
! 549: * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
! 550: * scale = 0, and compute a solution to A*x = 0.
! 551: *
! 552: DO 100 I = 1, N
! 553: X( I ) = ZERO
! 554: 100 CONTINUE
! 555: X( J ) = ONE
! 556: XJ = ONE
! 557: SCALE = ZERO
! 558: XMAX = ZERO
! 559: END IF
! 560: 110 CONTINUE
! 561: *
! 562: * Scale x if necessary to avoid overflow when adding a
! 563: * multiple of column j of A.
! 564: *
! 565: IF( XJ.GT.ONE ) THEN
! 566: REC = ONE / XJ
! 567: IF( CNORM( J ).GT.( BIGNUM-XMAX )*REC ) THEN
! 568: *
! 569: * Scale x by 1/(2*abs(x(j))).
! 570: *
! 571: REC = REC*HALF
! 572: CALL ZDSCAL( N, REC, X, 1 )
! 573: SCALE = SCALE*REC
! 574: END IF
! 575: ELSE IF( XJ*CNORM( J ).GT.( BIGNUM-XMAX ) ) THEN
! 576: *
! 577: * Scale x by 1/2.
! 578: *
! 579: CALL ZDSCAL( N, HALF, X, 1 )
! 580: SCALE = SCALE*HALF
! 581: END IF
! 582: *
! 583: IF( UPPER ) THEN
! 584: IF( J.GT.1 ) THEN
! 585: *
! 586: * Compute the update
! 587: * x(max(1,j-kd):j-1) := x(max(1,j-kd):j-1) -
! 588: * x(j)* A(max(1,j-kd):j-1,j)
! 589: *
! 590: JLEN = MIN( KD, J-1 )
! 591: CALL ZAXPY( JLEN, -X( J )*TSCAL,
! 592: $ AB( KD+1-JLEN, J ), 1, X( J-JLEN ), 1 )
! 593: I = IZAMAX( J-1, X, 1 )
! 594: XMAX = CABS1( X( I ) )
! 595: END IF
! 596: ELSE IF( J.LT.N ) THEN
! 597: *
! 598: * Compute the update
! 599: * x(j+1:min(j+kd,n)) := x(j+1:min(j+kd,n)) -
! 600: * x(j) * A(j+1:min(j+kd,n),j)
! 601: *
! 602: JLEN = MIN( KD, N-J )
! 603: IF( JLEN.GT.0 )
! 604: $ CALL ZAXPY( JLEN, -X( J )*TSCAL, AB( 2, J ), 1,
! 605: $ X( J+1 ), 1 )
! 606: I = J + IZAMAX( N-J, X( J+1 ), 1 )
! 607: XMAX = CABS1( X( I ) )
! 608: END IF
! 609: 120 CONTINUE
! 610: *
! 611: ELSE IF( LSAME( TRANS, 'T' ) ) THEN
! 612: *
! 613: * Solve A**T * x = b
! 614: *
! 615: DO 170 J = JFIRST, JLAST, JINC
! 616: *
! 617: * Compute x(j) = b(j) - sum A(k,j)*x(k).
! 618: * k<>j
! 619: *
! 620: XJ = CABS1( X( J ) )
! 621: USCAL = TSCAL
! 622: REC = ONE / MAX( XMAX, ONE )
! 623: IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN
! 624: *
! 625: * If x(j) could overflow, scale x by 1/(2*XMAX).
! 626: *
! 627: REC = REC*HALF
! 628: IF( NOUNIT ) THEN
! 629: TJJS = AB( MAIND, J )*TSCAL
! 630: ELSE
! 631: TJJS = TSCAL
! 632: END IF
! 633: TJJ = CABS1( TJJS )
! 634: IF( TJJ.GT.ONE ) THEN
! 635: *
! 636: * Divide by A(j,j) when scaling x if A(j,j) > 1.
! 637: *
! 638: REC = MIN( ONE, REC*TJJ )
! 639: USCAL = ZLADIV( USCAL, TJJS )
! 640: END IF
! 641: IF( REC.LT.ONE ) THEN
! 642: CALL ZDSCAL( N, REC, X, 1 )
! 643: SCALE = SCALE*REC
! 644: XMAX = XMAX*REC
! 645: END IF
! 646: END IF
! 647: *
! 648: CSUMJ = ZERO
! 649: IF( USCAL.EQ.DCMPLX( ONE ) ) THEN
! 650: *
! 651: * If the scaling needed for A in the dot product is 1,
! 652: * call ZDOTU to perform the dot product.
! 653: *
! 654: IF( UPPER ) THEN
! 655: JLEN = MIN( KD, J-1 )
! 656: CSUMJ = ZDOTU( JLEN, AB( KD+1-JLEN, J ), 1,
! 657: $ X( J-JLEN ), 1 )
! 658: ELSE
! 659: JLEN = MIN( KD, N-J )
! 660: IF( JLEN.GT.1 )
! 661: $ CSUMJ = ZDOTU( JLEN, AB( 2, J ), 1, X( J+1 ),
! 662: $ 1 )
! 663: END IF
! 664: ELSE
! 665: *
! 666: * Otherwise, use in-line code for the dot product.
! 667: *
! 668: IF( UPPER ) THEN
! 669: JLEN = MIN( KD, J-1 )
! 670: DO 130 I = 1, JLEN
! 671: CSUMJ = CSUMJ + ( AB( KD+I-JLEN, J )*USCAL )*
! 672: $ X( J-JLEN-1+I )
! 673: 130 CONTINUE
! 674: ELSE
! 675: JLEN = MIN( KD, N-J )
! 676: DO 140 I = 1, JLEN
! 677: CSUMJ = CSUMJ + ( AB( I+1, J )*USCAL )*X( J+I )
! 678: 140 CONTINUE
! 679: END IF
! 680: END IF
! 681: *
! 682: IF( USCAL.EQ.DCMPLX( TSCAL ) ) THEN
! 683: *
! 684: * Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
! 685: * was not used to scale the dotproduct.
! 686: *
! 687: X( J ) = X( J ) - CSUMJ
! 688: XJ = CABS1( X( J ) )
! 689: IF( NOUNIT ) THEN
! 690: *
! 691: * Compute x(j) = x(j) / A(j,j), scaling if necessary.
! 692: *
! 693: TJJS = AB( MAIND, J )*TSCAL
! 694: ELSE
! 695: TJJS = TSCAL
! 696: IF( TSCAL.EQ.ONE )
! 697: $ GO TO 160
! 698: END IF
! 699: TJJ = CABS1( TJJS )
! 700: IF( TJJ.GT.SMLNUM ) THEN
! 701: *
! 702: * abs(A(j,j)) > SMLNUM:
! 703: *
! 704: IF( TJJ.LT.ONE ) THEN
! 705: IF( XJ.GT.TJJ*BIGNUM ) THEN
! 706: *
! 707: * Scale X by 1/abs(x(j)).
! 708: *
! 709: REC = ONE / XJ
! 710: CALL ZDSCAL( N, REC, X, 1 )
! 711: SCALE = SCALE*REC
! 712: XMAX = XMAX*REC
! 713: END IF
! 714: END IF
! 715: X( J ) = ZLADIV( X( J ), TJJS )
! 716: ELSE IF( TJJ.GT.ZERO ) THEN
! 717: *
! 718: * 0 < abs(A(j,j)) <= SMLNUM:
! 719: *
! 720: IF( XJ.GT.TJJ*BIGNUM ) THEN
! 721: *
! 722: * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
! 723: *
! 724: REC = ( TJJ*BIGNUM ) / XJ
! 725: CALL ZDSCAL( N, REC, X, 1 )
! 726: SCALE = SCALE*REC
! 727: XMAX = XMAX*REC
! 728: END IF
! 729: X( J ) = ZLADIV( X( J ), TJJS )
! 730: ELSE
! 731: *
! 732: * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
! 733: * scale = 0 and compute a solution to A**T *x = 0.
! 734: *
! 735: DO 150 I = 1, N
! 736: X( I ) = ZERO
! 737: 150 CONTINUE
! 738: X( J ) = ONE
! 739: SCALE = ZERO
! 740: XMAX = ZERO
! 741: END IF
! 742: 160 CONTINUE
! 743: ELSE
! 744: *
! 745: * Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
! 746: * product has already been divided by 1/A(j,j).
! 747: *
! 748: X( J ) = ZLADIV( X( J ), TJJS ) - CSUMJ
! 749: END IF
! 750: XMAX = MAX( XMAX, CABS1( X( J ) ) )
! 751: 170 CONTINUE
! 752: *
! 753: ELSE
! 754: *
! 755: * Solve A**H * x = b
! 756: *
! 757: DO 220 J = JFIRST, JLAST, JINC
! 758: *
! 759: * Compute x(j) = b(j) - sum A(k,j)*x(k).
! 760: * k<>j
! 761: *
! 762: XJ = CABS1( X( J ) )
! 763: USCAL = TSCAL
! 764: REC = ONE / MAX( XMAX, ONE )
! 765: IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN
! 766: *
! 767: * If x(j) could overflow, scale x by 1/(2*XMAX).
! 768: *
! 769: REC = REC*HALF
! 770: IF( NOUNIT ) THEN
! 771: TJJS = DCONJG( AB( MAIND, J ) )*TSCAL
! 772: ELSE
! 773: TJJS = TSCAL
! 774: END IF
! 775: TJJ = CABS1( TJJS )
! 776: IF( TJJ.GT.ONE ) THEN
! 777: *
! 778: * Divide by A(j,j) when scaling x if A(j,j) > 1.
! 779: *
! 780: REC = MIN( ONE, REC*TJJ )
! 781: USCAL = ZLADIV( USCAL, TJJS )
! 782: END IF
! 783: IF( REC.LT.ONE ) THEN
! 784: CALL ZDSCAL( N, REC, X, 1 )
! 785: SCALE = SCALE*REC
! 786: XMAX = XMAX*REC
! 787: END IF
! 788: END IF
! 789: *
! 790: CSUMJ = ZERO
! 791: IF( USCAL.EQ.DCMPLX( ONE ) ) THEN
! 792: *
! 793: * If the scaling needed for A in the dot product is 1,
! 794: * call ZDOTC to perform the dot product.
! 795: *
! 796: IF( UPPER ) THEN
! 797: JLEN = MIN( KD, J-1 )
! 798: CSUMJ = ZDOTC( JLEN, AB( KD+1-JLEN, J ), 1,
! 799: $ X( J-JLEN ), 1 )
! 800: ELSE
! 801: JLEN = MIN( KD, N-J )
! 802: IF( JLEN.GT.1 )
! 803: $ CSUMJ = ZDOTC( JLEN, AB( 2, J ), 1, X( J+1 ),
! 804: $ 1 )
! 805: END IF
! 806: ELSE
! 807: *
! 808: * Otherwise, use in-line code for the dot product.
! 809: *
! 810: IF( UPPER ) THEN
! 811: JLEN = MIN( KD, J-1 )
! 812: DO 180 I = 1, JLEN
! 813: CSUMJ = CSUMJ + ( DCONJG( AB( KD+I-JLEN, J ) )*
! 814: $ USCAL )*X( J-JLEN-1+I )
! 815: 180 CONTINUE
! 816: ELSE
! 817: JLEN = MIN( KD, N-J )
! 818: DO 190 I = 1, JLEN
! 819: CSUMJ = CSUMJ + ( DCONJG( AB( I+1, J ) )*USCAL )
! 820: $ *X( J+I )
! 821: 190 CONTINUE
! 822: END IF
! 823: END IF
! 824: *
! 825: IF( USCAL.EQ.DCMPLX( TSCAL ) ) THEN
! 826: *
! 827: * Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
! 828: * was not used to scale the dotproduct.
! 829: *
! 830: X( J ) = X( J ) - CSUMJ
! 831: XJ = CABS1( X( J ) )
! 832: IF( NOUNIT ) THEN
! 833: *
! 834: * Compute x(j) = x(j) / A(j,j), scaling if necessary.
! 835: *
! 836: TJJS = DCONJG( AB( MAIND, J ) )*TSCAL
! 837: ELSE
! 838: TJJS = TSCAL
! 839: IF( TSCAL.EQ.ONE )
! 840: $ GO TO 210
! 841: END IF
! 842: TJJ = CABS1( TJJS )
! 843: IF( TJJ.GT.SMLNUM ) THEN
! 844: *
! 845: * abs(A(j,j)) > SMLNUM:
! 846: *
! 847: IF( TJJ.LT.ONE ) THEN
! 848: IF( XJ.GT.TJJ*BIGNUM ) THEN
! 849: *
! 850: * Scale X by 1/abs(x(j)).
! 851: *
! 852: REC = ONE / XJ
! 853: CALL ZDSCAL( N, REC, X, 1 )
! 854: SCALE = SCALE*REC
! 855: XMAX = XMAX*REC
! 856: END IF
! 857: END IF
! 858: X( J ) = ZLADIV( X( J ), TJJS )
! 859: ELSE IF( TJJ.GT.ZERO ) THEN
! 860: *
! 861: * 0 < abs(A(j,j)) <= SMLNUM:
! 862: *
! 863: IF( XJ.GT.TJJ*BIGNUM ) THEN
! 864: *
! 865: * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
! 866: *
! 867: REC = ( TJJ*BIGNUM ) / XJ
! 868: CALL ZDSCAL( N, REC, X, 1 )
! 869: SCALE = SCALE*REC
! 870: XMAX = XMAX*REC
! 871: END IF
! 872: X( J ) = ZLADIV( X( J ), TJJS )
! 873: ELSE
! 874: *
! 875: * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
! 876: * scale = 0 and compute a solution to A**H *x = 0.
! 877: *
! 878: DO 200 I = 1, N
! 879: X( I ) = ZERO
! 880: 200 CONTINUE
! 881: X( J ) = ONE
! 882: SCALE = ZERO
! 883: XMAX = ZERO
! 884: END IF
! 885: 210 CONTINUE
! 886: ELSE
! 887: *
! 888: * Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
! 889: * product has already been divided by 1/A(j,j).
! 890: *
! 891: X( J ) = ZLADIV( X( J ), TJJS ) - CSUMJ
! 892: END IF
! 893: XMAX = MAX( XMAX, CABS1( X( J ) ) )
! 894: 220 CONTINUE
! 895: END IF
! 896: SCALE = SCALE / TSCAL
! 897: END IF
! 898: *
! 899: * Scale the column norms by 1/TSCAL for return.
! 900: *
! 901: IF( TSCAL.NE.ONE ) THEN
! 902: CALL DSCAL( N, ONE / TSCAL, CNORM, 1 )
! 903: END IF
! 904: *
! 905: RETURN
! 906: *
! 907: * End of ZLATBS
! 908: *
! 909: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlatbs.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack