LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ zlarfgp()

 subroutine zlarfgp ( integer N, complex*16 ALPHA, complex*16, dimension( * ) X, integer INCX, complex*16 TAU )

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

Purpose:
``` ZLARFGP generates a complex elementary reflector H of order n, such
that

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

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

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

where tau is a complex scalar and v is a complex (n-1)-element
vector. Note that H is not hermitian.

If the elements of x are all zero and alpha is real, 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 COMPLEX*16 On entry, the value alpha. On exit, it is overwritten with the value beta.``` [in,out] X ``` X is COMPLEX*16 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 COMPLEX*16 The value tau.```

Definition at line 103 of file zlarfgp.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 COMPLEX*16 ALPHA, TAU
112* ..
113* .. Array Arguments ..
114 COMPLEX*16 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 ALPHI, ALPHR, BETA, BIGNUM, SMLNUM, XNORM
126 COMPLEX*16 SAVEALPHA
127* ..
128* .. External Functions ..
129 DOUBLE PRECISION DLAMCH, DLAPY3, DLAPY2, DZNRM2
131 EXTERNAL dlamch, dlapy3, dlapy2, dznrm2, zladiv
132* ..
133* .. Intrinsic Functions ..
134 INTRINSIC abs, dble, dcmplx, dimag, sign
135* ..
136* .. External Subroutines ..
137 EXTERNAL zdscal, zscal
138* ..
139* .. Executable Statements ..
140*
141 IF( n.LE.0 ) THEN
142 tau = zero
143 RETURN
144 END IF
145*
146 xnorm = dznrm2( n-1, x, incx )
147 alphr = dble( alpha )
148 alphi = dimag( alpha )
149*
150 IF( xnorm.EQ.zero ) THEN
151*
152* H = [1-alpha/abs(alpha) 0; 0 I], sign chosen so ALPHA >= 0.
153*
154 IF( alphi.EQ.zero ) THEN
155 IF( alphr.GE.zero ) THEN
156* When TAU.eq.ZERO, the vector is special-cased to be
157* all zeros in the application routines. We do not need
158* to clear it.
159 tau = zero
160 ELSE
161* However, the application routines rely on explicit
162* zero checks when TAU.ne.ZERO, and we must clear X.
163 tau = two
164 DO j = 1, n-1
165 x( 1 + (j-1)*incx ) = zero
166 END DO
167 alpha = -alpha
168 END IF
169 ELSE
170* Only "reflecting" the diagonal entry to be real and non-negative.
171 xnorm = dlapy2( alphr, alphi )
172 tau = dcmplx( one - alphr / xnorm, -alphi / xnorm )
173 DO j = 1, n-1
174 x( 1 + (j-1)*incx ) = zero
175 END DO
176 alpha = xnorm
177 END IF
178 ELSE
179*
180* general case
181*
182 beta = sign( dlapy3( alphr, alphi, xnorm ), alphr )
183 smlnum = dlamch( 'S' ) / dlamch( 'E' )
184 bignum = one / smlnum
185*
186 knt = 0
187 IF( abs( beta ).LT.smlnum ) THEN
188*
189* XNORM, BETA may be inaccurate; scale X and recompute them
190*
191 10 CONTINUE
192 knt = knt + 1
193 CALL zdscal( n-1, bignum, x, incx )
194 beta = beta*bignum
195 alphi = alphi*bignum
196 alphr = alphr*bignum
197 IF( (abs( beta ).LT.smlnum) .AND. (knt .LT. 20) )
198 \$ GO TO 10
199*
200* New BETA is at most 1, at least SMLNUM
201*
202 xnorm = dznrm2( n-1, x, incx )
203 alpha = dcmplx( alphr, alphi )
204 beta = sign( dlapy3( alphr, alphi, xnorm ), alphr )
205 END IF
206 savealpha = alpha
207 alpha = alpha + beta
208 IF( beta.LT.zero ) THEN
209 beta = -beta
210 tau = -alpha / beta
211 ELSE
212 alphr = alphi * (alphi/dble( alpha ))
213 alphr = alphr + xnorm * (xnorm/dble( alpha ))
214 tau = dcmplx( alphr/beta, -alphi/beta )
215 alpha = dcmplx( -alphr, alphi )
216 END IF
217 alpha = zladiv( dcmplx( one ), alpha )
218*
219 IF ( abs(tau).LE.smlnum ) THEN
220*
221* In the case where the computed TAU ends up being a denormalized number,
222* it loses relative accuracy. This is a BIG problem. Solution: flush TAU
223* to ZERO (or TWO or whatever makes a nonnegative real number for BETA).
224*
225* (Bug report provided by Pat Quillen from MathWorks on Jul 29, 2009.)
226* (Thanks Pat. Thanks MathWorks.)
227*
228 alphr = dble( savealpha )
229 alphi = dimag( savealpha )
230 IF( alphi.EQ.zero ) THEN
231 IF( alphr.GE.zero ) THEN
232 tau = zero
233 ELSE
234 tau = two
235 DO j = 1, n-1
236 x( 1 + (j-1)*incx ) = zero
237 END DO
238 beta = dble( -savealpha )
239 END IF
240 ELSE
241 xnorm = dlapy2( alphr, alphi )
242 tau = dcmplx( one - alphr / xnorm, -alphi / xnorm )
243 DO j = 1, n-1
244 x( 1 + (j-1)*incx ) = zero
245 END DO
246 beta = xnorm
247 END IF
248*
249 ELSE
250*
251* This is the general case.
252*
253 CALL zscal( n-1, alpha, x, incx )
254*
255 END IF
256*
257* If BETA is subnormal, it may lose relative accuracy
258*
259 DO 20 j = 1, knt
260 beta = beta*smlnum
261 20 CONTINUE
262 alpha = beta
263 END IF
264*
265 RETURN
266*
267* End of ZLARFGP
268*
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
double precision function dlapy3(X, Y, Z)
DLAPY3 returns sqrt(x2+y2+z2).
Definition: dlapy3.f:68
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:78
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:78