LAPACK  3.8.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.
Date
November 2017

Definition at line 106 of file dlarfgp.f.

106 *
107 * -- LAPACK auxiliary routine (version 3.8.0) --
108 * -- LAPACK is a software package provided by Univ. of Tennessee, --
109 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
110 * November 2017
111 *
112 * .. Scalar Arguments ..
113  INTEGER incx, n
114  DOUBLE PRECISION alpha, tau
115 * ..
116 * .. Array Arguments ..
117  DOUBLE PRECISION x( * )
118 * ..
119 *
120 * =====================================================================
121 *
122 * .. Parameters ..
123  DOUBLE PRECISION two, one, zero
124  parameter( two = 2.0d+0, one = 1.0d+0, zero = 0.0d+0 )
125 * ..
126 * .. Local Scalars ..
127  INTEGER j, knt
128  DOUBLE PRECISION beta, bignum, savealpha, smlnum, xnorm
129 * ..
130 * .. External Functions ..
131  DOUBLE PRECISION dlamch, dlapy2, dnrm2
132  EXTERNAL dlamch, dlapy2, dnrm2
133 * ..
134 * .. Intrinsic Functions ..
135  INTRINSIC abs, sign
136 * ..
137 * .. External Subroutines ..
138  EXTERNAL dscal
139 * ..
140 * .. Executable Statements ..
141 *
142  IF( n.LE.0 ) THEN
143  tau = zero
144  RETURN
145  END IF
146 *
147  xnorm = dnrm2( n-1, x, incx )
148 *
149  IF( xnorm.EQ.zero ) THEN
150 *
151 * H = [+/-1, 0; I], sign chosen so ALPHA >= 0
152 *
153  IF( alpha.GE.zero ) THEN
154 * When TAU.eq.ZERO, the vector is special-cased to be
155 * all zeros in the application routines. We do not need
156 * to clear it.
157  tau = zero
158  ELSE
159 * However, the application routines rely on explicit
160 * zero checks when TAU.ne.ZERO, and we must clear X.
161  tau = two
162  DO j = 1, n-1
163  x( 1 + (j-1)*incx ) = 0
164  END DO
165  alpha = -alpha
166  END IF
167  ELSE
168 *
169 * general case
170 *
171  beta = sign( dlapy2( alpha, xnorm ), alpha )
172  smlnum = dlamch( 'S' ) / dlamch( 'E' )
173  knt = 0
174  IF( abs( beta ).LT.smlnum ) THEN
175 *
176 * XNORM, BETA may be inaccurate; scale X and recompute them
177 *
178  bignum = one / smlnum
179  10 CONTINUE
180  knt = knt + 1
181  CALL dscal( n-1, bignum, x, incx )
182  beta = beta*bignum
183  alpha = alpha*bignum
184  IF( (abs( beta ).LT.smlnum) .AND. (knt .LT. 20) )
185  $ GO TO 10
186 *
187 * New BETA is at most 1, at least SMLNUM
188 *
189  xnorm = dnrm2( n-1, x, incx )
190  beta = sign( dlapy2( alpha, xnorm ), alpha )
191  END IF
192  savealpha = alpha
193  alpha = alpha + beta
194  IF( beta.LT.zero ) THEN
195  beta = -beta
196  tau = -alpha / beta
197  ELSE
198  alpha = xnorm * (xnorm/alpha)
199  tau = alpha / beta
200  alpha = -alpha
201  END IF
202 *
203  IF ( abs(tau).LE.smlnum ) THEN
204 *
205 * In the case where the computed TAU ends up being a denormalized number,
206 * it loses relative accuracy. This is a BIG problem. Solution: flush TAU
207 * to ZERO. This explains the next IF statement.
208 *
209 * (Bug report provided by Pat Quillen from MathWorks on Jul 29, 2009.)
210 * (Thanks Pat. Thanks MathWorks.)
211 *
212  IF( savealpha.GE.zero ) THEN
213  tau = zero
214  ELSE
215  tau = two
216  DO j = 1, n-1
217  x( 1 + (j-1)*incx ) = 0
218  END DO
219  beta = -savealpha
220  END IF
221 *
222  ELSE
223 *
224 * This is the general case.
225 *
226  CALL dscal( n-1, one / alpha, x, incx )
227 *
228  END IF
229 *
230 * If BETA is subnormal, it may lose relative accuracy
231 *
232  DO 20 j = 1, knt
233  beta = beta*smlnum
234  20 CONTINUE
235  alpha = beta
236  END IF
237 *
238  RETURN
239 *
240 * End of DLARFGP
241 *
double precision function dnrm2(N, X, INCX)
DNRM2
Definition: dnrm2.f:76
double precision function dlapy2(X, Y)
DLAPY2 returns sqrt(x2+y2).
Definition: dlapy2.f:65
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:81
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
Here is the call graph for this function:
Here is the caller graph for this function: