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