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