![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
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