![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZLARFG( N, ALPHA, X, INCX, TAU ) 2: * 3: * -- LAPACK auxiliary routine (version 3.2) -- 4: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 6: * November 2006 7: * 8: * .. Scalar Arguments .. 9: INTEGER INCX, N 10: COMPLEX*16 ALPHA, TAU 11: * .. 12: * .. Array Arguments .. 13: COMPLEX*16 X( * ) 14: * .. 15: * 16: * Purpose 17: * ======= 18: * 19: * ZLARFG generates a complex elementary reflector H of order n, such 20: * that 21: * 22: * H' * ( alpha ) = ( beta ), H' * H = I. 23: * ( x ) ( 0 ) 24: * 25: * where alpha and beta are scalars, with beta real, and x is an 26: * (n-1)-element complex vector. H is represented in the form 27: * 28: * H = I - tau * ( 1 ) * ( 1 v' ) , 29: * ( v ) 30: * 31: * where tau is a complex scalar and v is a complex (n-1)-element 32: * vector. Note that H is not hermitian. 33: * 34: * If the elements of x are all zero and alpha is real, then tau = 0 35: * and H is taken to be the unit matrix. 36: * 37: * Otherwise 1 <= real(tau) <= 2 and abs(tau-1) <= 1 . 38: * 39: * Arguments 40: * ========= 41: * 42: * N (input) INTEGER 43: * The order of the elementary reflector. 44: * 45: * ALPHA (input/output) COMPLEX*16 46: * On entry, the value alpha. 47: * On exit, it is overwritten with the value beta. 48: * 49: * X (input/output) COMPLEX*16 array, dimension 50: * (1+(N-2)*abs(INCX)) 51: * On entry, the vector x. 52: * On exit, it is overwritten with the vector v. 53: * 54: * INCX (input) INTEGER 55: * The increment between elements of X. INCX > 0. 56: * 57: * TAU (output) COMPLEX*16 58: * The value tau. 59: * 60: * ===================================================================== 61: * 62: * .. Parameters .. 63: DOUBLE PRECISION ONE, ZERO 64: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 65: * .. 66: * .. Local Scalars .. 67: INTEGER J, KNT 68: DOUBLE PRECISION ALPHI, ALPHR, BETA, RSAFMN, SAFMIN, XNORM 69: * .. 70: * .. External Functions .. 71: DOUBLE PRECISION DLAMCH, DLAPY3, DZNRM2 72: COMPLEX*16 ZLADIV 73: EXTERNAL DLAMCH, DLAPY3, DZNRM2, ZLADIV 74: * .. 75: * .. Intrinsic Functions .. 76: INTRINSIC ABS, DBLE, DCMPLX, DIMAG, SIGN 77: * .. 78: * .. External Subroutines .. 79: EXTERNAL ZDSCAL, ZSCAL 80: * .. 81: * .. Executable Statements .. 82: * 83: IF( N.LE.0 ) THEN 84: TAU = ZERO 85: RETURN 86: END IF 87: * 88: XNORM = DZNRM2( N-1, X, INCX ) 89: ALPHR = DBLE( ALPHA ) 90: ALPHI = DIMAG( ALPHA ) 91: * 92: IF( XNORM.EQ.ZERO .AND. ALPHI.EQ.ZERO ) THEN 93: * 94: * H = I 95: * 96: TAU = ZERO 97: ELSE 98: * 99: * general case 100: * 101: BETA = -SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) 102: SAFMIN = DLAMCH( 'S' ) / DLAMCH( 'E' ) 103: RSAFMN = ONE / SAFMIN 104: * 105: KNT = 0 106: IF( ABS( BETA ).LT.SAFMIN ) THEN 107: * 108: * XNORM, BETA may be inaccurate; scale X and recompute them 109: * 110: 10 CONTINUE 111: KNT = KNT + 1 112: CALL ZDSCAL( N-1, RSAFMN, X, INCX ) 113: BETA = BETA*RSAFMN 114: ALPHI = ALPHI*RSAFMN 115: ALPHR = ALPHR*RSAFMN 116: IF( ABS( BETA ).LT.SAFMIN ) 117: $ GO TO 10 118: * 119: * New BETA is at most 1, at least SAFMIN 120: * 121: XNORM = DZNRM2( N-1, X, INCX ) 122: ALPHA = DCMPLX( ALPHR, ALPHI ) 123: BETA = -SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) 124: END IF 125: TAU = DCMPLX( ( BETA-ALPHR ) / BETA, -ALPHI / BETA ) 126: ALPHA = ZLADIV( DCMPLX( ONE ), ALPHA-BETA ) 127: CALL ZSCAL( N-1, ALPHA, X, INCX ) 128: * 129: * If ALPHA is subnormal, it may lose relative accuracy 130: * 131: DO 20 J = 1, KNT 132: BETA = BETA*SAFMIN 133: 20 CONTINUE 134: ALPHA = BETA 135: END IF 136: * 137: RETURN 138: * 139: * End of ZLARFG 140: * 141: END