1*> \brief \b DLARFGP generates an elementary reflector (Householder matrix) with non-negative beta. 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DLARFGP + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarfgp.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarfgp.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarfgp.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU ) 22* 23* .. Scalar Arguments .. 24* INTEGER INCX, N 25* DOUBLE PRECISION ALPHA, TAU 26* .. 27* .. Array Arguments .. 28* DOUBLE PRECISION X( * ) 29* .. 30* 31* 32*> \par Purpose: 33* ============= 34*> 35*> \verbatim 36*> 37*> DLARFGP generates a real elementary reflector H of order n, such 38*> that 39*> 40*> H * ( alpha ) = ( beta ), H**T * H = I. 41*> ( x ) ( 0 ) 42*> 43*> where alpha and beta are scalars, beta is non-negative, and x is 44*> an (n-1)-element real vector. H is represented in the form 45*> 46*> H = I - tau * ( 1 ) * ( 1 v**T ) , 47*> ( v ) 48*> 49*> where tau is a real scalar and v is a real (n-1)-element 50*> vector. 51*> 52*> If the elements of x are all zero, then tau = 0 and H is taken to be 53*> the unit matrix. 54*> \endverbatim 55* 56* Arguments: 57* ========== 58* 59*> \param[in] N 60*> \verbatim 61*> N is INTEGER 62*> The order of the elementary reflector. 63*> \endverbatim 64*> 65*> \param[in,out] ALPHA 66*> \verbatim 67*> ALPHA is DOUBLE PRECISION 68*> On entry, the value alpha. 69*> On exit, it is overwritten with the value beta. 70*> \endverbatim 71*> 72*> \param[in,out] X 73*> \verbatim 74*> X is DOUBLE PRECISION array, dimension 75*> (1+(N-2)*abs(INCX)) 76*> On entry, the vector x. 77*> On exit, it is overwritten with the vector v. 78*> \endverbatim 79*> 80*> \param[in] INCX 81*> \verbatim 82*> INCX is INTEGER 83*> The increment between elements of X. INCX > 0. 84*> \endverbatim 85*> 86*> \param[out] TAU 87*> \verbatim 88*> TAU is DOUBLE PRECISION 89*> The value tau. 90*> \endverbatim 91* 92* Authors: 93* ======== 94* 95*> \author Univ. of Tennessee 96*> \author Univ. of California Berkeley 97*> \author Univ. of Colorado Denver 98*> \author NAG Ltd. 99* 100*> \ingroup doubleOTHERauxiliary 101* 102* ===================================================================== 103 SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU ) 104* 105* -- LAPACK auxiliary routine -- 106* -- LAPACK is a software package provided by Univ. of Tennessee, -- 107* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 108* 109* .. Scalar Arguments .. 110 INTEGER INCX, N 111 DOUBLE PRECISION ALPHA, TAU 112* .. 113* .. Array Arguments .. 114 DOUBLE PRECISION X( * ) 115* .. 116* 117* ===================================================================== 118* 119* .. Parameters .. 120 DOUBLE PRECISION TWO, ONE, ZERO 121 PARAMETER ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 ) 122* .. 123* .. Local Scalars .. 124 INTEGER J, KNT 125 DOUBLE PRECISION BETA, BIGNUM, SAVEALPHA, SMLNUM, XNORM 126* .. 127* .. External Functions .. 128 DOUBLE PRECISION DLAMCH, DLAPY2, DNRM2 129 EXTERNAL DLAMCH, DLAPY2, DNRM2 130* .. 131* .. Intrinsic Functions .. 132 INTRINSIC ABS, SIGN 133* .. 134* .. External Subroutines .. 135 EXTERNAL DSCAL 136* .. 137* .. Executable Statements .. 138* 139 IF( N.LE.0 ) THEN 140 TAU = ZERO 141 RETURN 142 END IF 143* 144 XNORM = DNRM2( N-1, X, INCX ) 145* 146 IF( XNORM.EQ.ZERO ) THEN 147* 148* H = [+/-1, 0; I], sign chosen so ALPHA >= 0 149* 150 IF( ALPHA.GE.ZERO ) THEN 151* When TAU.eq.ZERO, the vector is special-cased to be 152* all zeros in the application routines. We do not need 153* to clear it. 154 TAU = ZERO 155 ELSE 156* However, the application routines rely on explicit 157* zero checks when TAU.ne.ZERO, and we must clear X. 158 TAU = TWO 159 DO J = 1, N-1 160 X( 1 + (J-1)*INCX ) = 0 161 END DO 162 ALPHA = -ALPHA 163 END IF 164 ELSE 165* 166* general case 167* 168 BETA = SIGN( DLAPY2( ALPHA, XNORM ), ALPHA ) 169 SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'E' ) 170 KNT = 0 171 IF( ABS( BETA ).LT.SMLNUM ) THEN 172* 173* XNORM, BETA may be inaccurate; scale X and recompute them 174* 175 BIGNUM = ONE / SMLNUM 176 10 CONTINUE 177 KNT = KNT + 1 178 CALL DSCAL( N-1, BIGNUM, X, INCX ) 179 BETA = BETA*BIGNUM 180 ALPHA = ALPHA*BIGNUM 181 IF( (ABS( BETA ).LT.SMLNUM) .AND. (KNT .LT. 20) ) 182 $ GO TO 10 183* 184* New BETA is at most 1, at least SMLNUM 185* 186 XNORM = DNRM2( N-1, X, INCX ) 187 BETA = SIGN( DLAPY2( ALPHA, XNORM ), ALPHA ) 188 END IF 189 SAVEALPHA = ALPHA 190 ALPHA = ALPHA + BETA 191 IF( BETA.LT.ZERO ) THEN 192 BETA = -BETA 193 TAU = -ALPHA / BETA 194 ELSE 195 ALPHA = XNORM * (XNORM/ALPHA) 196 TAU = ALPHA / BETA 197 ALPHA = -ALPHA 198 END IF 199* 200 IF ( ABS(TAU).LE.SMLNUM ) THEN 201* 202* In the case where the computed TAU ends up being a denormalized number, 203* it loses relative accuracy. This is a BIG problem. Solution: flush TAU 204* to ZERO. This explains the next IF statement. 205* 206* (Bug report provided by Pat Quillen from MathWorks on Jul 29, 2009.) 207* (Thanks Pat. Thanks MathWorks.) 208* 209 IF( SAVEALPHA.GE.ZERO ) THEN 210 TAU = ZERO 211 ELSE 212 TAU = TWO 213 DO J = 1, N-1 214 X( 1 + (J-1)*INCX ) = 0 215 END DO 216 BETA = -SAVEALPHA 217 END IF 218* 219 ELSE 220* 221* This is the general case. 222* 223 CALL DSCAL( N-1, ONE / ALPHA, X, INCX ) 224* 225 END IF 226* 227* If BETA is subnormal, it may lose relative accuracy 228* 229 DO 20 J = 1, KNT 230 BETA = BETA*SMLNUM 231 20 CONTINUE 232 ALPHA = BETA 233 END IF 234* 235 RETURN 236* 237* End of DLARFGP 238* 239 END 240