LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zlarfgp.f
Go to the documentation of this file.
1*> \brief \b ZLARFGP 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 ZLARFGP + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlarfgp.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlarfgp.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlarfgp.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE ZLARFGP( N, ALPHA, X, INCX, TAU )
22*
23* .. Scalar Arguments ..
24* INTEGER INCX, N
25* COMPLEX*16 ALPHA, TAU
26* ..
27* .. Array Arguments ..
28* COMPLEX*16 X( * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> ZLARFGP generates a complex elementary reflector H of order n, such
38*> that
39*>
40*> H**H * ( alpha ) = ( beta ), H**H * H = I.
41*> ( x ) ( 0 )
42*>
43*> where alpha and beta are scalars, beta is real and non-negative, and
44*> x is an (n-1)-element complex vector. H is represented in the form
45*>
46*> H = I - tau * ( 1 ) * ( 1 v**H ) ,
47*> ( v )
48*>
49*> where tau is a complex scalar and v is a complex (n-1)-element
50*> vector. Note that H is not hermitian.
51*>
52*> If the elements of x are all zero and alpha is real, then tau = 0
53*> and H is taken to be 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 COMPLEX*16
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 COMPLEX*16 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 COMPLEX*16
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 larfgp
101*
102* =====================================================================
103 SUBROUTINE zlarfgp( 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 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, EPS, SMLNUM, XNORM
126 COMPLEX*16 SAVEALPHA
127* ..
128* .. External Functions ..
129 DOUBLE PRECISION DLAMCH, DLAPY3, DLAPY2, DZNRM2
130 COMPLEX*16 ZLADIV
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 eps = dlamch( 'Precision' )
147 xnorm = dznrm2( n-1, x, incx )
148 alphr = dble( alpha )
149 alphi = dimag( alpha )
150*
151 IF( xnorm.LE.eps*abs(alpha) ) THEN
152*
153* H = [1-alpha/abs(alpha) 0; 0 I], sign chosen so ALPHA >= 0.
154*
155 IF( alphi.EQ.zero ) THEN
156 IF( alphr.GE.zero ) THEN
157* When TAU.eq.ZERO, the vector is special-cased to be
158* all zeros in the application routines. We do not need
159* to clear it.
160 tau = zero
161 ELSE
162* However, the application routines rely on explicit
163* zero checks when TAU.ne.ZERO, and we must clear X.
164 tau = two
165 DO j = 1, n-1
166 x( 1 + (j-1)*incx ) = zero
167 END DO
168 alpha = -alpha
169 END IF
170 ELSE
171* Only "reflecting" the diagonal entry to be real and non-negative.
172 xnorm = dlapy2( alphr, alphi )
173 tau = dcmplx( one - alphr / xnorm, -alphi / xnorm )
174 DO j = 1, n-1
175 x( 1 + (j-1)*incx ) = zero
176 END DO
177 alpha = xnorm
178 END IF
179 ELSE
180*
181* general case
182*
183 beta = sign( dlapy3( alphr, alphi, xnorm ), alphr )
184 smlnum = dlamch( 'S' ) / dlamch( 'E' )
185 bignum = one / smlnum
186*
187 knt = 0
188 IF( abs( beta ).LT.smlnum ) THEN
189*
190* XNORM, BETA may be inaccurate; scale X and recompute them
191*
192 10 CONTINUE
193 knt = knt + 1
194 CALL zdscal( n-1, bignum, x, incx )
195 beta = beta*bignum
196 alphi = alphi*bignum
197 alphr = alphr*bignum
198 IF( (abs( beta ).LT.smlnum) .AND. (knt .LT. 20) )
199 $ GO TO 10
200*
201* New BETA is at most 1, at least SMLNUM
202*
203 xnorm = dznrm2( n-1, x, incx )
204 alpha = dcmplx( alphr, alphi )
205 beta = sign( dlapy3( alphr, alphi, xnorm ), alphr )
206 END IF
207 savealpha = alpha
208 alpha = alpha + beta
209 IF( beta.LT.zero ) THEN
210 beta = -beta
211 tau = -alpha / beta
212 ELSE
213 alphr = alphi * (alphi/dble( alpha ))
214 alphr = alphr + xnorm * (xnorm/dble( alpha ))
215 tau = dcmplx( alphr/beta, -alphi/beta )
216 alpha = dcmplx( -alphr, alphi )
217 END IF
218 alpha = zladiv( dcmplx( one ), alpha )
219*
220 IF ( abs(tau).LE.smlnum ) THEN
221*
222* In the case where the computed TAU ends up being a denormalized number,
223* it loses relative accuracy. This is a BIG problem. Solution: flush TAU
224* to ZERO (or TWO or whatever makes a nonnegative real number for BETA).
225*
226* (Bug report provided by Pat Quillen from MathWorks on Jul 29, 2009.)
227* (Thanks Pat. Thanks MathWorks.)
228*
229 alphr = dble( savealpha )
230 alphi = dimag( savealpha )
231 IF( alphi.EQ.zero ) THEN
232 IF( alphr.GE.zero ) THEN
233 tau = zero
234 ELSE
235 tau = two
236 DO j = 1, n-1
237 x( 1 + (j-1)*incx ) = zero
238 END DO
239 beta = dble( -savealpha )
240 END IF
241 ELSE
242 xnorm = dlapy2( alphr, alphi )
243 tau = dcmplx( one - alphr / xnorm, -alphi / xnorm )
244 DO j = 1, n-1
245 x( 1 + (j-1)*incx ) = zero
246 END DO
247 beta = xnorm
248 END IF
249*
250 ELSE
251*
252* This is the general case.
253*
254 CALL zscal( n-1, alpha, x, incx )
255*
256 END IF
257*
258* If BETA is subnormal, it may lose relative accuracy
259*
260 DO 20 j = 1, knt
261 beta = beta*smlnum
262 20 CONTINUE
263 alpha = beta
264 END IF
265*
266 RETURN
267*
268* End of ZLARFGP
269*
270 END
subroutine zlarfgp(n, alpha, x, incx, tau)
ZLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.
Definition zlarfgp.f:104
subroutine zdscal(n, da, zx, incx)
ZDSCAL
Definition zdscal.f:78
subroutine zscal(n, za, zx, incx)
ZSCAL
Definition zscal.f:78