![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU ) 2: * 3: * -- LAPACK auxiliary routine (version 3.2.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: * June 2010 7: * 8: * .. Scalar Arguments .. 9: INTEGER INCX, N 10: DOUBLE PRECISION ALPHA, TAU 11: * .. 12: * .. Array Arguments .. 13: DOUBLE PRECISION X( * ) 14: * .. 15: * 16: * Purpose 17: * ======= 18: * 19: * DLARFGP generates a real 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, beta is non-negative, and x is 26: * an (n-1)-element real vector. H is represented in the form 27: * 28: * H = I - tau * ( 1 ) * ( 1 v' ) , 29: * ( v ) 30: * 31: * where tau is a real scalar and v is a real (n-1)-element 32: * vector. 33: * 34: * If the elements of x are all zero, then tau = 0 and H is taken to be 35: * the unit matrix. 36: * 37: * Arguments 38: * ========= 39: * 40: * N (input) INTEGER 41: * The order of the elementary reflector. 42: * 43: * ALPHA (input/output) DOUBLE PRECISION 44: * On entry, the value alpha. 45: * On exit, it is overwritten with the value beta. 46: * 47: * X (input/output) DOUBLE PRECISION array, dimension 48: * (1+(N-2)*abs(INCX)) 49: * On entry, the vector x. 50: * On exit, it is overwritten with the vector v. 51: * 52: * INCX (input) INTEGER 53: * The increment between elements of X. INCX > 0. 54: * 55: * TAU (output) DOUBLE PRECISION 56: * The value tau. 57: * 58: * ===================================================================== 59: * 60: * .. Parameters .. 61: DOUBLE PRECISION TWO, ONE, ZERO 62: PARAMETER ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 ) 63: * .. 64: * .. Local Scalars .. 65: INTEGER J, KNT 66: DOUBLE PRECISION BETA, BIGNUM, SAVEALPHA, SMLNUM, XNORM 67: * .. 68: * .. External Functions .. 69: DOUBLE PRECISION DLAMCH, DLAPY2, DNRM2 70: EXTERNAL DLAMCH, DLAPY2, DNRM2 71: * .. 72: * .. Intrinsic Functions .. 73: INTRINSIC ABS, SIGN 74: * .. 75: * .. External Subroutines .. 76: EXTERNAL DSCAL 77: * .. 78: * .. Executable Statements .. 79: * 80: IF( N.LE.0 ) THEN 81: TAU = ZERO 82: RETURN 83: END IF 84: * 85: XNORM = DNRM2( N-1, X, INCX ) 86: * 87: IF( XNORM.EQ.ZERO ) THEN 88: * 89: * H = [+/-1, 0; I], sign chosen so ALPHA >= 0 90: * 91: IF( ALPHA.GE.ZERO ) THEN 92: * When TAU.eq.ZERO, the vector is special-cased to be 93: * all zeros in the application routines. We do not need 94: * to clear it. 95: TAU = ZERO 96: ELSE 97: * However, the application routines rely on explicit 98: * zero checks when TAU.ne.ZERO, and we must clear X. 99: TAU = TWO 100: DO J = 1, N-1 101: X( 1 + (J-1)*INCX ) = 0 102: END DO 103: ALPHA = -ALPHA 104: END IF 105: ELSE 106: * 107: * general case 108: * 109: BETA = SIGN( DLAPY2( ALPHA, XNORM ), ALPHA ) 110: SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'E' ) 111: KNT = 0 112: IF( ABS( BETA ).LT.SMLNUM ) THEN 113: * 114: * XNORM, BETA may be inaccurate; scale X and recompute them 115: * 116: BIGNUM = ONE / SMLNUM 117: 10 CONTINUE 118: KNT = KNT + 1 119: CALL DSCAL( N-1, BIGNUM, X, INCX ) 120: BETA = BETA*BIGNUM 121: ALPHA = ALPHA*BIGNUM 122: IF( ABS( BETA ).LT.SMLNUM ) 123: $ GO TO 10 124: * 125: * New BETA is at most 1, at least SMLNUM 126: * 127: XNORM = DNRM2( N-1, X, INCX ) 128: BETA = SIGN( DLAPY2( ALPHA, XNORM ), ALPHA ) 129: END IF 130: SAVEALPHA = ALPHA 131: ALPHA = ALPHA + BETA 132: IF( BETA.LT.ZERO ) THEN 133: BETA = -BETA 134: TAU = -ALPHA / BETA 135: ELSE 136: ALPHA = XNORM * (XNORM/ALPHA) 137: TAU = ALPHA / BETA 138: ALPHA = -ALPHA 139: END IF 140: * 141: IF ( ABS(TAU).LE.SMLNUM ) THEN 142: * 143: * In the case where the computed TAU ends up being a denormalized number, 144: * it loses relative accuracy. This is a BIG problem. Solution: flush TAU 145: * to ZERO. This explains the next IF statement. 146: * 147: * (Bug report provided by Pat Quillen from MathWorks on Jul 29, 2009.) 148: * (Thanks Pat. Thanks MathWorks.) 149: * 150: IF( SAVEALPHA.GE.ZERO ) THEN 151: TAU = ZERO 152: ELSE 153: TAU = TWO 154: DO J = 1, N-1 155: X( 1 + (J-1)*INCX ) = 0 156: END DO 157: BETA = -SAVEALPHA 158: END IF 159: * 160: ELSE 161: * 162: * This is the general case. 163: * 164: CALL DSCAL( N-1, ONE / ALPHA, X, INCX ) 165: * 166: END IF 167: * 168: * If BETA is subnormal, it may lose relative accuracy 169: * 170: DO 20 J = 1, KNT 171: BETA = BETA*SMLNUM 172: 20 CONTINUE 173: ALPHA = BETA 174: END IF 175: * 176: RETURN 177: * 178: * End of DLARFGP 179: * 180: END