1 SUBROUTINE DLARFG( N, ALPHA, X, INCX, TAU ) 2* 3* -- LAPACK auxiliary routine (version 3.0) -- 4* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 5* Courant Institute, Argonne National Lab, and Rice University 6* September 30, 1994 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* DLARFG 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, and x is an (n-1)-element real 26* 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* Otherwise 1 <= tau <= 2. 38* 39* Arguments 40* ========= 41* 42* N (input) INTEGER 43* The order of the elementary reflector. 44* 45* ALPHA (input/output) DOUBLE PRECISION 46* On entry, the value alpha. 47* On exit, it is overwritten with the value beta. 48* 49* X (input/output) DOUBLE PRECISION 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) DOUBLE PRECISION 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 BETA, RSAFMN, SAFMIN, XNORM 69* .. 70* .. External Functions .. 71 DOUBLE PRECISION DLAMCH, DLAPY2, DNRM2 72 EXTERNAL DLAMCH, DLAPY2, DNRM2 73* .. 74* .. Intrinsic Functions .. 75 INTRINSIC ABS, SIGN 76* .. 77* .. External Subroutines .. 78 EXTERNAL DSCAL 79* .. 80* .. Executable Statements .. 81* 82 IF( N.LE.1 ) THEN 83 TAU = ZERO 84 RETURN 85 END IF 86* 87 XNORM = DNRM2( N-1, X, INCX ) 88* 89 IF( XNORM.EQ.ZERO ) THEN 90* 91* H = I 92* 93 TAU = ZERO 94 ELSE 95* 96* general case 97* 98 BETA = -SIGN( DLAPY2( ALPHA, XNORM ), ALPHA ) 99 SAFMIN = DLAMCH( 'S' ) / DLAMCH( 'E' ) 100 IF( ABS( BETA ).LT.SAFMIN ) THEN 101* 102* XNORM, BETA may be inaccurate; scale X and recompute them 103* 104 RSAFMN = ONE / SAFMIN 105 KNT = 0 106 10 CONTINUE 107 KNT = KNT + 1 108 CALL DSCAL( N-1, RSAFMN, X, INCX ) 109 BETA = BETA*RSAFMN 110 ALPHA = ALPHA*RSAFMN 111 IF( ABS( BETA ).LT.SAFMIN ) 112 $ GO TO 10 113* 114* New BETA is at most 1, at least SAFMIN 115* 116 XNORM = DNRM2( N-1, X, INCX ) 117 BETA = -SIGN( DLAPY2( ALPHA, XNORM ), ALPHA ) 118 TAU = ( BETA-ALPHA ) / BETA 119 CALL DSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX ) 120* 121* If ALPHA is subnormal, it may lose relative accuracy 122* 123 ALPHA = BETA 124 DO 20 J = 1, KNT 125 ALPHA = ALPHA*SAFMIN 126 20 CONTINUE 127 ELSE 128 TAU = ( BETA-ALPHA ) / BETA 129 CALL DSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX ) 130 ALPHA = BETA 131 END IF 132 END IF 133* 134 RETURN 135* 136* End of DLARFG 137* 138 END 139