LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

## ◆ clarfgp()

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

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

Purpose:
``` CLARFGP 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 On entry, the value alpha. On exit, it is overwritten with the value beta.``` [in,out] X ``` X is COMPLEX 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 The value tau.```

Definition at line 103 of file clarfgp.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 ALPHA, TAU
112 * ..
113 * .. Array Arguments ..
114  COMPLEX X( * )
115 * ..
116 *
117 * =====================================================================
118 *
119 * .. Parameters ..
120  REAL TWO, ONE, ZERO
121  parameter( two = 2.0e+0, one = 1.0e+0, zero = 0.0e+0 )
122 * ..
123 * .. Local Scalars ..
124  INTEGER J, KNT
125  REAL ALPHI, ALPHR, BETA, BIGNUM, SMLNUM, XNORM
126  COMPLEX SAVEALPHA
127 * ..
128 * .. External Functions ..
129  REAL SCNRM2, SLAMCH, SLAPY3, SLAPY2
131  EXTERNAL scnrm2, slamch, slapy3, slapy2, cladiv
132 * ..
133 * .. Intrinsic Functions ..
134  INTRINSIC abs, aimag, cmplx, real, sign
135 * ..
136 * .. External Subroutines ..
137  EXTERNAL cscal, csscal
138 * ..
139 * .. Executable Statements ..
140 *
141  IF( n.LE.0 ) THEN
142  tau = zero
143  RETURN
144  END IF
145 *
146  xnorm = scnrm2( n-1, x, incx )
147  alphr = real( alpha )
148  alphi = aimag( 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 = slapy2( alphr, alphi )
172  tau = cmplx( 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( slapy3( alphr, alphi, xnorm ), alphr )
183  smlnum = slamch( 'S' ) / slamch( '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 csscal( 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 = scnrm2( n-1, x, incx )
203  alpha = cmplx( alphr, alphi )
204  beta = sign( slapy3( 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/real( alpha ))
213  alphr = alphr + xnorm * (xnorm/real( alpha ))
214  tau = cmplx( alphr/beta, -alphi/beta )
215  alpha = cmplx( -alphr, alphi )
216  END IF
217  alpha = cladiv( cmplx( 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 = real( savealpha )
229  alphi = aimag( 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 = real( -savealpha )
239  END IF
240  ELSE
241  xnorm = slapy2( alphr, alphi )
242  tau = cmplx( 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 cscal( 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 CLARFGP
268 *
real function slapy2(X, Y)
SLAPY2 returns sqrt(x2+y2).
Definition: slapy2.f:63
real function slapy3(X, Y, Z)
SLAPY3 returns sqrt(x2+y2+z2).
Definition: slapy3.f:68
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:78
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:78