LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ dlarfgp()

subroutine dlarfgp ( integer  N,
double precision  ALPHA,
double precision, dimension( * )  X,
integer  INCX,
double precision  TAU 
)

DLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.

Download DLARFGP + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 DLARFGP generates a real elementary reflector H of order n, such
 that

       H * ( alpha ) = ( beta ),   H**T * H = I.
           (   x   )   (   0  )

 where alpha and beta are scalars, beta is non-negative, and x is
 an (n-1)-element real vector.  H is represented in the form

       H = I - tau * ( 1 ) * ( 1 v**T ) ,
                     ( v )

 where tau is a real scalar and v is a real (n-1)-element
 vector.

 If the elements of x are all zero, then tau = 0 and H is taken to be
 the unit matrix.
Parameters
[in]N
          N is INTEGER
          The order of the elementary reflector.
[in,out]ALPHA
          ALPHA is DOUBLE PRECISION
          On entry, the value alpha.
          On exit, it is overwritten with the value beta.
[in,out]X
          X is DOUBLE PRECISION array, dimension
                         (1+(N-2)*abs(INCX))
          On entry, the vector x.
          On exit, it is overwritten with the vector v.
[in]INCX
          INCX is INTEGER
          The increment between elements of X. INCX > 0.
[out]TAU
          TAU is DOUBLE PRECISION
          The value tau.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 103 of file dlarfgp.f.

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 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
double precision function dlapy2(X, Y)
DLAPY2 returns sqrt(x2+y2).
Definition: dlapy2.f:63
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
real(wp) function dnrm2(n, x, incx)
DNRM2
Definition: dnrm2.f90:89
Here is the call graph for this function:
Here is the caller graph for this function: