Annotation of rpl/lapack/lapack/zlatrs.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE ZLATRS( 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 CNORM( * )
! 16: COMPLEX*16 A( LDA, * ), X( * )
! 17: * ..
! 18: *
! 19: * Purpose
! 20: * =======
! 21: *
! 22: * ZLATRS 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. Here A is an upper or lower
! 27: * triangular matrix, A**T denotes the transpose of A, A**H denotes the
! 28: * conjugate transpose of A, x and b are n-element vectors, and s is a
! 29: * scaling factor, usually less than or equal to 1, chosen so that the
! 30: * components of x will be less than the overflow threshold. If the
! 31: * unscaled problem will not cause overflow, the Level 2 BLAS routine
! 32: * ZTRSV is called. If the matrix A is singular (A(j,j) = 0 for some j),
! 33: * then s is set to 0 and a 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: * A (input) COMPLEX*16 array, dimension (LDA,N)
! 64: * The triangular matrix A. If UPLO = 'U', the leading n by n
! 65: * upper triangular part of the array A contains the upper
! 66: * triangular matrix, and the strictly lower triangular part of
! 67: * A is not referenced. If UPLO = 'L', the leading n by n lower
! 68: * triangular part of the array A contains the lower triangular
! 69: * matrix, and the strictly upper triangular part of A is not
! 70: * referenced. If DIAG = 'U', the diagonal elements of A are
! 71: * also not referenced and are assumed to be 1.
! 72: *
! 73: * LDA (input) INTEGER
! 74: * The leading dimension of the array A. LDA >= max (1,N).
! 75: *
! 76: * X (input/output) COMPLEX*16 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, A**T * x = s*b, or A**H * 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, ZTRSV
! 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 ZTRSV 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**T *x = b or
! 149: * A**H *x = b. The basic 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 ZTRSV 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, TWO
! 175: PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0,
! 176: $ TWO = 2.0D+0 )
! 177: * ..
! 178: * .. Local Scalars ..
! 179: LOGICAL NOTRAN, NOUNIT, UPPER
! 180: INTEGER I, IMAX, J, JFIRST, JINC, JLAST
! 181: DOUBLE PRECISION BIGNUM, GROW, REC, SMLNUM, TJJ, TMAX, TSCAL,
! 182: $ XBND, XJ, XMAX
! 183: COMPLEX*16 CSUMJ, TJJS, USCAL, ZDUM
! 184: * ..
! 185: * .. External Functions ..
! 186: LOGICAL LSAME
! 187: INTEGER IDAMAX, IZAMAX
! 188: DOUBLE PRECISION DLAMCH, DZASUM
! 189: COMPLEX*16 ZDOTC, ZDOTU, ZLADIV
! 190: EXTERNAL LSAME, IDAMAX, IZAMAX, DLAMCH, DZASUM, ZDOTC,
! 191: $ ZDOTU, ZLADIV
! 192: * ..
! 193: * .. External Subroutines ..
! 194: EXTERNAL DSCAL, XERBLA, ZAXPY, ZDSCAL, ZTRSV
! 195: * ..
! 196: * .. Intrinsic Functions ..
! 197: INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN
! 198: * ..
! 199: * .. Statement Functions ..
! 200: DOUBLE PRECISION CABS1, CABS2
! 201: * ..
! 202: * .. Statement Function definitions ..
! 203: CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
! 204: CABS2( ZDUM ) = ABS( DBLE( ZDUM ) / 2.D0 ) +
! 205: $ ABS( DIMAG( ZDUM ) / 2.D0 )
! 206: * ..
! 207: * .. Executable Statements ..
! 208: *
! 209: INFO = 0
! 210: UPPER = LSAME( UPLO, 'U' )
! 211: NOTRAN = LSAME( TRANS, 'N' )
! 212: NOUNIT = LSAME( DIAG, 'N' )
! 213: *
! 214: * Test the input parameters.
! 215: *
! 216: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
! 217: INFO = -1
! 218: ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
! 219: $ LSAME( TRANS, 'C' ) ) THEN
! 220: INFO = -2
! 221: ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
! 222: INFO = -3
! 223: ELSE IF( .NOT.LSAME( NORMIN, 'Y' ) .AND. .NOT.
! 224: $ LSAME( NORMIN, 'N' ) ) THEN
! 225: INFO = -4
! 226: ELSE IF( N.LT.0 ) THEN
! 227: INFO = -5
! 228: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
! 229: INFO = -7
! 230: END IF
! 231: IF( INFO.NE.0 ) THEN
! 232: CALL XERBLA( 'ZLATRS', -INFO )
! 233: RETURN
! 234: END IF
! 235: *
! 236: * Quick return if possible
! 237: *
! 238: IF( N.EQ.0 )
! 239: $ RETURN
! 240: *
! 241: * Determine machine dependent parameters to control overflow.
! 242: *
! 243: SMLNUM = DLAMCH( 'Safe minimum' )
! 244: BIGNUM = ONE / SMLNUM
! 245: CALL DLABAD( SMLNUM, BIGNUM )
! 246: SMLNUM = SMLNUM / DLAMCH( 'Precision' )
! 247: BIGNUM = ONE / SMLNUM
! 248: SCALE = ONE
! 249: *
! 250: IF( LSAME( NORMIN, 'N' ) ) THEN
! 251: *
! 252: * Compute the 1-norm of each column, not including the diagonal.
! 253: *
! 254: IF( UPPER ) THEN
! 255: *
! 256: * A is upper triangular.
! 257: *
! 258: DO 10 J = 1, N
! 259: CNORM( J ) = DZASUM( J-1, A( 1, J ), 1 )
! 260: 10 CONTINUE
! 261: ELSE
! 262: *
! 263: * A is lower triangular.
! 264: *
! 265: DO 20 J = 1, N - 1
! 266: CNORM( J ) = DZASUM( N-J, A( J+1, J ), 1 )
! 267: 20 CONTINUE
! 268: CNORM( N ) = ZERO
! 269: END IF
! 270: END IF
! 271: *
! 272: * Scale the column norms by TSCAL if the maximum element in CNORM is
! 273: * greater than BIGNUM/2.
! 274: *
! 275: IMAX = IDAMAX( N, CNORM, 1 )
! 276: TMAX = CNORM( IMAX )
! 277: IF( TMAX.LE.BIGNUM*HALF ) THEN
! 278: TSCAL = ONE
! 279: ELSE
! 280: TSCAL = HALF / ( SMLNUM*TMAX )
! 281: CALL DSCAL( N, TSCAL, CNORM, 1 )
! 282: END IF
! 283: *
! 284: * Compute a bound on the computed solution vector to see if the
! 285: * Level 2 BLAS routine ZTRSV can be used.
! 286: *
! 287: XMAX = ZERO
! 288: DO 30 J = 1, N
! 289: XMAX = MAX( XMAX, CABS2( X( J ) ) )
! 290: 30 CONTINUE
! 291: XBND = XMAX
! 292: *
! 293: IF( NOTRAN ) THEN
! 294: *
! 295: * Compute the growth in A * x = b.
! 296: *
! 297: IF( UPPER ) THEN
! 298: JFIRST = N
! 299: JLAST = 1
! 300: JINC = -1
! 301: ELSE
! 302: JFIRST = 1
! 303: JLAST = N
! 304: JINC = 1
! 305: END IF
! 306: *
! 307: IF( TSCAL.NE.ONE ) THEN
! 308: GROW = ZERO
! 309: GO TO 60
! 310: END IF
! 311: *
! 312: IF( NOUNIT ) THEN
! 313: *
! 314: * A is non-unit triangular.
! 315: *
! 316: * Compute GROW = 1/G(j) and XBND = 1/M(j).
! 317: * Initially, G(0) = max{x(i), i=1,...,n}.
! 318: *
! 319: GROW = HALF / MAX( XBND, SMLNUM )
! 320: XBND = GROW
! 321: DO 40 J = JFIRST, JLAST, JINC
! 322: *
! 323: * Exit the loop if the growth factor is too small.
! 324: *
! 325: IF( GROW.LE.SMLNUM )
! 326: $ GO TO 60
! 327: *
! 328: TJJS = A( J, J )
! 329: TJJ = CABS1( TJJS )
! 330: *
! 331: IF( TJJ.GE.SMLNUM ) THEN
! 332: *
! 333: * M(j) = G(j-1) / abs(A(j,j))
! 334: *
! 335: XBND = MIN( XBND, MIN( ONE, TJJ )*GROW )
! 336: ELSE
! 337: *
! 338: * M(j) could overflow, set XBND to 0.
! 339: *
! 340: XBND = ZERO
! 341: END IF
! 342: *
! 343: IF( TJJ+CNORM( J ).GE.SMLNUM ) THEN
! 344: *
! 345: * G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
! 346: *
! 347: GROW = GROW*( TJJ / ( TJJ+CNORM( J ) ) )
! 348: ELSE
! 349: *
! 350: * G(j) could overflow, set GROW to 0.
! 351: *
! 352: GROW = ZERO
! 353: END IF
! 354: 40 CONTINUE
! 355: GROW = XBND
! 356: ELSE
! 357: *
! 358: * A is unit triangular.
! 359: *
! 360: * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
! 361: *
! 362: GROW = MIN( ONE, HALF / MAX( XBND, SMLNUM ) )
! 363: DO 50 J = JFIRST, JLAST, JINC
! 364: *
! 365: * Exit the loop if the growth factor is too small.
! 366: *
! 367: IF( GROW.LE.SMLNUM )
! 368: $ GO TO 60
! 369: *
! 370: * G(j) = G(j-1)*( 1 + CNORM(j) )
! 371: *
! 372: GROW = GROW*( ONE / ( ONE+CNORM( J ) ) )
! 373: 50 CONTINUE
! 374: END IF
! 375: 60 CONTINUE
! 376: *
! 377: ELSE
! 378: *
! 379: * Compute the growth in A**T * x = b or A**H * x = b.
! 380: *
! 381: IF( UPPER ) THEN
! 382: JFIRST = 1
! 383: JLAST = N
! 384: JINC = 1
! 385: ELSE
! 386: JFIRST = N
! 387: JLAST = 1
! 388: JINC = -1
! 389: END IF
! 390: *
! 391: IF( TSCAL.NE.ONE ) THEN
! 392: GROW = ZERO
! 393: GO TO 90
! 394: END IF
! 395: *
! 396: IF( NOUNIT ) THEN
! 397: *
! 398: * A is non-unit triangular.
! 399: *
! 400: * Compute GROW = 1/G(j) and XBND = 1/M(j).
! 401: * Initially, M(0) = max{x(i), i=1,...,n}.
! 402: *
! 403: GROW = HALF / MAX( XBND, SMLNUM )
! 404: XBND = GROW
! 405: DO 70 J = JFIRST, JLAST, JINC
! 406: *
! 407: * Exit the loop if the growth factor is too small.
! 408: *
! 409: IF( GROW.LE.SMLNUM )
! 410: $ GO TO 90
! 411: *
! 412: * G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
! 413: *
! 414: XJ = ONE + CNORM( J )
! 415: GROW = MIN( GROW, XBND / XJ )
! 416: *
! 417: TJJS = A( J, J )
! 418: TJJ = CABS1( TJJS )
! 419: *
! 420: IF( TJJ.GE.SMLNUM ) THEN
! 421: *
! 422: * M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
! 423: *
! 424: IF( XJ.GT.TJJ )
! 425: $ XBND = XBND*( TJJ / XJ )
! 426: ELSE
! 427: *
! 428: * M(j) could overflow, set XBND to 0.
! 429: *
! 430: XBND = ZERO
! 431: END IF
! 432: 70 CONTINUE
! 433: GROW = MIN( GROW, XBND )
! 434: ELSE
! 435: *
! 436: * A is unit triangular.
! 437: *
! 438: * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
! 439: *
! 440: GROW = MIN( ONE, HALF / MAX( XBND, SMLNUM ) )
! 441: DO 80 J = JFIRST, JLAST, JINC
! 442: *
! 443: * Exit the loop if the growth factor is too small.
! 444: *
! 445: IF( GROW.LE.SMLNUM )
! 446: $ GO TO 90
! 447: *
! 448: * G(j) = ( 1 + CNORM(j) )*G(j-1)
! 449: *
! 450: XJ = ONE + CNORM( J )
! 451: GROW = GROW / XJ
! 452: 80 CONTINUE
! 453: END IF
! 454: 90 CONTINUE
! 455: END IF
! 456: *
! 457: IF( ( GROW*TSCAL ).GT.SMLNUM ) THEN
! 458: *
! 459: * Use the Level 2 BLAS solve if the reciprocal of the bound on
! 460: * elements of X is not too small.
! 461: *
! 462: CALL ZTRSV( UPLO, TRANS, DIAG, N, A, LDA, X, 1 )
! 463: ELSE
! 464: *
! 465: * Use a Level 1 BLAS solve, scaling intermediate results.
! 466: *
! 467: IF( XMAX.GT.BIGNUM*HALF ) THEN
! 468: *
! 469: * Scale X so that its components are less than or equal to
! 470: * BIGNUM in absolute value.
! 471: *
! 472: SCALE = ( BIGNUM*HALF ) / XMAX
! 473: CALL ZDSCAL( N, SCALE, X, 1 )
! 474: XMAX = BIGNUM
! 475: ELSE
! 476: XMAX = XMAX*TWO
! 477: END IF
! 478: *
! 479: IF( NOTRAN ) THEN
! 480: *
! 481: * Solve A * x = b
! 482: *
! 483: DO 120 J = JFIRST, JLAST, JINC
! 484: *
! 485: * Compute x(j) = b(j) / A(j,j), scaling x if necessary.
! 486: *
! 487: XJ = CABS1( X( J ) )
! 488: IF( NOUNIT ) THEN
! 489: TJJS = A( J, J )*TSCAL
! 490: ELSE
! 491: TJJS = TSCAL
! 492: IF( TSCAL.EQ.ONE )
! 493: $ GO TO 110
! 494: END IF
! 495: TJJ = CABS1( TJJS )
! 496: IF( TJJ.GT.SMLNUM ) THEN
! 497: *
! 498: * abs(A(j,j)) > SMLNUM:
! 499: *
! 500: IF( TJJ.LT.ONE ) THEN
! 501: IF( XJ.GT.TJJ*BIGNUM ) THEN
! 502: *
! 503: * Scale x by 1/b(j).
! 504: *
! 505: REC = ONE / XJ
! 506: CALL ZDSCAL( N, REC, X, 1 )
! 507: SCALE = SCALE*REC
! 508: XMAX = XMAX*REC
! 509: END IF
! 510: END IF
! 511: X( J ) = ZLADIV( X( J ), TJJS )
! 512: XJ = CABS1( X( J ) )
! 513: ELSE IF( TJJ.GT.ZERO ) THEN
! 514: *
! 515: * 0 < abs(A(j,j)) <= SMLNUM:
! 516: *
! 517: IF( XJ.GT.TJJ*BIGNUM ) THEN
! 518: *
! 519: * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
! 520: * to avoid overflow when dividing by A(j,j).
! 521: *
! 522: REC = ( TJJ*BIGNUM ) / XJ
! 523: IF( CNORM( J ).GT.ONE ) THEN
! 524: *
! 525: * Scale by 1/CNORM(j) to avoid overflow when
! 526: * multiplying x(j) times column j.
! 527: *
! 528: REC = REC / CNORM( J )
! 529: END IF
! 530: CALL ZDSCAL( N, REC, X, 1 )
! 531: SCALE = SCALE*REC
! 532: XMAX = XMAX*REC
! 533: END IF
! 534: X( J ) = ZLADIV( X( J ), TJJS )
! 535: XJ = CABS1( X( J ) )
! 536: ELSE
! 537: *
! 538: * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
! 539: * scale = 0, and compute a solution to A*x = 0.
! 540: *
! 541: DO 100 I = 1, N
! 542: X( I ) = ZERO
! 543: 100 CONTINUE
! 544: X( J ) = ONE
! 545: XJ = ONE
! 546: SCALE = ZERO
! 547: XMAX = ZERO
! 548: END IF
! 549: 110 CONTINUE
! 550: *
! 551: * Scale x if necessary to avoid overflow when adding a
! 552: * multiple of column j of A.
! 553: *
! 554: IF( XJ.GT.ONE ) THEN
! 555: REC = ONE / XJ
! 556: IF( CNORM( J ).GT.( BIGNUM-XMAX )*REC ) THEN
! 557: *
! 558: * Scale x by 1/(2*abs(x(j))).
! 559: *
! 560: REC = REC*HALF
! 561: CALL ZDSCAL( N, REC, X, 1 )
! 562: SCALE = SCALE*REC
! 563: END IF
! 564: ELSE IF( XJ*CNORM( J ).GT.( BIGNUM-XMAX ) ) THEN
! 565: *
! 566: * Scale x by 1/2.
! 567: *
! 568: CALL ZDSCAL( N, HALF, X, 1 )
! 569: SCALE = SCALE*HALF
! 570: END IF
! 571: *
! 572: IF( UPPER ) THEN
! 573: IF( J.GT.1 ) THEN
! 574: *
! 575: * Compute the update
! 576: * x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j)
! 577: *
! 578: CALL ZAXPY( J-1, -X( J )*TSCAL, A( 1, J ), 1, X,
! 579: $ 1 )
! 580: I = IZAMAX( J-1, X, 1 )
! 581: XMAX = CABS1( X( I ) )
! 582: END IF
! 583: ELSE
! 584: IF( J.LT.N ) THEN
! 585: *
! 586: * Compute the update
! 587: * x(j+1:n) := x(j+1:n) - x(j) * A(j+1:n,j)
! 588: *
! 589: CALL ZAXPY( N-J, -X( J )*TSCAL, A( J+1, J ), 1,
! 590: $ X( J+1 ), 1 )
! 591: I = J + IZAMAX( N-J, X( J+1 ), 1 )
! 592: XMAX = CABS1( X( I ) )
! 593: END IF
! 594: END IF
! 595: 120 CONTINUE
! 596: *
! 597: ELSE IF( LSAME( TRANS, 'T' ) ) THEN
! 598: *
! 599: * Solve A**T * x = b
! 600: *
! 601: DO 170 J = JFIRST, JLAST, JINC
! 602: *
! 603: * Compute x(j) = b(j) - sum A(k,j)*x(k).
! 604: * k<>j
! 605: *
! 606: XJ = CABS1( X( J ) )
! 607: USCAL = TSCAL
! 608: REC = ONE / MAX( XMAX, ONE )
! 609: IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN
! 610: *
! 611: * If x(j) could overflow, scale x by 1/(2*XMAX).
! 612: *
! 613: REC = REC*HALF
! 614: IF( NOUNIT ) THEN
! 615: TJJS = A( J, J )*TSCAL
! 616: ELSE
! 617: TJJS = TSCAL
! 618: END IF
! 619: TJJ = CABS1( TJJS )
! 620: IF( TJJ.GT.ONE ) THEN
! 621: *
! 622: * Divide by A(j,j) when scaling x if A(j,j) > 1.
! 623: *
! 624: REC = MIN( ONE, REC*TJJ )
! 625: USCAL = ZLADIV( USCAL, TJJS )
! 626: END IF
! 627: IF( REC.LT.ONE ) THEN
! 628: CALL ZDSCAL( N, REC, X, 1 )
! 629: SCALE = SCALE*REC
! 630: XMAX = XMAX*REC
! 631: END IF
! 632: END IF
! 633: *
! 634: CSUMJ = ZERO
! 635: IF( USCAL.EQ.DCMPLX( ONE ) ) THEN
! 636: *
! 637: * If the scaling needed for A in the dot product is 1,
! 638: * call ZDOTU to perform the dot product.
! 639: *
! 640: IF( UPPER ) THEN
! 641: CSUMJ = ZDOTU( J-1, A( 1, J ), 1, X, 1 )
! 642: ELSE IF( J.LT.N ) THEN
! 643: CSUMJ = ZDOTU( N-J, A( J+1, J ), 1, X( J+1 ), 1 )
! 644: END IF
! 645: ELSE
! 646: *
! 647: * Otherwise, use in-line code for the dot product.
! 648: *
! 649: IF( UPPER ) THEN
! 650: DO 130 I = 1, J - 1
! 651: CSUMJ = CSUMJ + ( A( I, J )*USCAL )*X( I )
! 652: 130 CONTINUE
! 653: ELSE IF( J.LT.N ) THEN
! 654: DO 140 I = J + 1, N
! 655: CSUMJ = CSUMJ + ( A( I, J )*USCAL )*X( I )
! 656: 140 CONTINUE
! 657: END IF
! 658: END IF
! 659: *
! 660: IF( USCAL.EQ.DCMPLX( TSCAL ) ) THEN
! 661: *
! 662: * Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
! 663: * was not used to scale the dotproduct.
! 664: *
! 665: X( J ) = X( J ) - CSUMJ
! 666: XJ = CABS1( X( J ) )
! 667: IF( NOUNIT ) THEN
! 668: TJJS = A( J, J )*TSCAL
! 669: ELSE
! 670: TJJS = TSCAL
! 671: IF( TSCAL.EQ.ONE )
! 672: $ GO TO 160
! 673: END IF
! 674: *
! 675: * Compute x(j) = x(j) / A(j,j), scaling if necessary.
! 676: *
! 677: TJJ = CABS1( TJJS )
! 678: IF( TJJ.GT.SMLNUM ) THEN
! 679: *
! 680: * abs(A(j,j)) > SMLNUM:
! 681: *
! 682: IF( TJJ.LT.ONE ) THEN
! 683: IF( XJ.GT.TJJ*BIGNUM ) THEN
! 684: *
! 685: * Scale X by 1/abs(x(j)).
! 686: *
! 687: REC = ONE / XJ
! 688: CALL ZDSCAL( N, REC, X, 1 )
! 689: SCALE = SCALE*REC
! 690: XMAX = XMAX*REC
! 691: END IF
! 692: END IF
! 693: X( J ) = ZLADIV( X( J ), TJJS )
! 694: ELSE IF( TJJ.GT.ZERO ) THEN
! 695: *
! 696: * 0 < abs(A(j,j)) <= SMLNUM:
! 697: *
! 698: IF( XJ.GT.TJJ*BIGNUM ) THEN
! 699: *
! 700: * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
! 701: *
! 702: REC = ( TJJ*BIGNUM ) / XJ
! 703: CALL ZDSCAL( N, REC, X, 1 )
! 704: SCALE = SCALE*REC
! 705: XMAX = XMAX*REC
! 706: END IF
! 707: X( J ) = ZLADIV( X( J ), TJJS )
! 708: ELSE
! 709: *
! 710: * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
! 711: * scale = 0 and compute a solution to A**T *x = 0.
! 712: *
! 713: DO 150 I = 1, N
! 714: X( I ) = ZERO
! 715: 150 CONTINUE
! 716: X( J ) = ONE
! 717: SCALE = ZERO
! 718: XMAX = ZERO
! 719: END IF
! 720: 160 CONTINUE
! 721: ELSE
! 722: *
! 723: * Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
! 724: * product has already been divided by 1/A(j,j).
! 725: *
! 726: X( J ) = ZLADIV( X( J ), TJJS ) - CSUMJ
! 727: END IF
! 728: XMAX = MAX( XMAX, CABS1( X( J ) ) )
! 729: 170 CONTINUE
! 730: *
! 731: ELSE
! 732: *
! 733: * Solve A**H * x = b
! 734: *
! 735: DO 220 J = JFIRST, JLAST, JINC
! 736: *
! 737: * Compute x(j) = b(j) - sum A(k,j)*x(k).
! 738: * k<>j
! 739: *
! 740: XJ = CABS1( X( J ) )
! 741: USCAL = TSCAL
! 742: REC = ONE / MAX( XMAX, ONE )
! 743: IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN
! 744: *
! 745: * If x(j) could overflow, scale x by 1/(2*XMAX).
! 746: *
! 747: REC = REC*HALF
! 748: IF( NOUNIT ) THEN
! 749: TJJS = DCONJG( A( J, J ) )*TSCAL
! 750: ELSE
! 751: TJJS = TSCAL
! 752: END IF
! 753: TJJ = CABS1( TJJS )
! 754: IF( TJJ.GT.ONE ) THEN
! 755: *
! 756: * Divide by A(j,j) when scaling x if A(j,j) > 1.
! 757: *
! 758: REC = MIN( ONE, REC*TJJ )
! 759: USCAL = ZLADIV( USCAL, TJJS )
! 760: END IF
! 761: IF( REC.LT.ONE ) THEN
! 762: CALL ZDSCAL( N, REC, X, 1 )
! 763: SCALE = SCALE*REC
! 764: XMAX = XMAX*REC
! 765: END IF
! 766: END IF
! 767: *
! 768: CSUMJ = ZERO
! 769: IF( USCAL.EQ.DCMPLX( ONE ) ) THEN
! 770: *
! 771: * If the scaling needed for A in the dot product is 1,
! 772: * call ZDOTC to perform the dot product.
! 773: *
! 774: IF( UPPER ) THEN
! 775: CSUMJ = ZDOTC( J-1, A( 1, J ), 1, X, 1 )
! 776: ELSE IF( J.LT.N ) THEN
! 777: CSUMJ = ZDOTC( N-J, A( J+1, J ), 1, X( J+1 ), 1 )
! 778: END IF
! 779: ELSE
! 780: *
! 781: * Otherwise, use in-line code for the dot product.
! 782: *
! 783: IF( UPPER ) THEN
! 784: DO 180 I = 1, J - 1
! 785: CSUMJ = CSUMJ + ( DCONJG( A( I, J ) )*USCAL )*
! 786: $ X( I )
! 787: 180 CONTINUE
! 788: ELSE IF( J.LT.N ) THEN
! 789: DO 190 I = J + 1, N
! 790: CSUMJ = CSUMJ + ( DCONJG( A( I, J ) )*USCAL )*
! 791: $ X( I )
! 792: 190 CONTINUE
! 793: END IF
! 794: END IF
! 795: *
! 796: IF( USCAL.EQ.DCMPLX( TSCAL ) ) THEN
! 797: *
! 798: * Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
! 799: * was not used to scale the dotproduct.
! 800: *
! 801: X( J ) = X( J ) - CSUMJ
! 802: XJ = CABS1( X( J ) )
! 803: IF( NOUNIT ) THEN
! 804: TJJS = DCONJG( A( J, J ) )*TSCAL
! 805: ELSE
! 806: TJJS = TSCAL
! 807: IF( TSCAL.EQ.ONE )
! 808: $ GO TO 210
! 809: END IF
! 810: *
! 811: * Compute x(j) = x(j) / A(j,j), scaling if necessary.
! 812: *
! 813: TJJ = CABS1( TJJS )
! 814: IF( TJJ.GT.SMLNUM ) THEN
! 815: *
! 816: * abs(A(j,j)) > SMLNUM:
! 817: *
! 818: IF( TJJ.LT.ONE ) THEN
! 819: IF( XJ.GT.TJJ*BIGNUM ) THEN
! 820: *
! 821: * Scale X by 1/abs(x(j)).
! 822: *
! 823: REC = ONE / XJ
! 824: CALL ZDSCAL( N, REC, X, 1 )
! 825: SCALE = SCALE*REC
! 826: XMAX = XMAX*REC
! 827: END IF
! 828: END IF
! 829: X( J ) = ZLADIV( X( J ), TJJS )
! 830: ELSE IF( TJJ.GT.ZERO ) THEN
! 831: *
! 832: * 0 < abs(A(j,j)) <= SMLNUM:
! 833: *
! 834: IF( XJ.GT.TJJ*BIGNUM ) THEN
! 835: *
! 836: * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
! 837: *
! 838: REC = ( TJJ*BIGNUM ) / XJ
! 839: CALL ZDSCAL( N, REC, X, 1 )
! 840: SCALE = SCALE*REC
! 841: XMAX = XMAX*REC
! 842: END IF
! 843: X( J ) = ZLADIV( X( J ), TJJS )
! 844: ELSE
! 845: *
! 846: * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
! 847: * scale = 0 and compute a solution to A**H *x = 0.
! 848: *
! 849: DO 200 I = 1, N
! 850: X( I ) = ZERO
! 851: 200 CONTINUE
! 852: X( J ) = ONE
! 853: SCALE = ZERO
! 854: XMAX = ZERO
! 855: END IF
! 856: 210 CONTINUE
! 857: ELSE
! 858: *
! 859: * Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
! 860: * product has already been divided by 1/A(j,j).
! 861: *
! 862: X( J ) = ZLADIV( X( J ), TJJS ) - CSUMJ
! 863: END IF
! 864: XMAX = MAX( XMAX, CABS1( X( J ) ) )
! 865: 220 CONTINUE
! 866: END IF
! 867: SCALE = SCALE / TSCAL
! 868: END IF
! 869: *
! 870: * Scale the column norms by 1/TSCAL for return.
! 871: *
! 872: IF( TSCAL.NE.ONE ) THEN
! 873: CALL DSCAL( N, ONE / TSCAL, CNORM, 1 )
! 874: END IF
! 875: *
! 876: RETURN
! 877: *
! 878: * End of ZLATRS
! 879: *
! 880: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlatrs.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack