LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine clarfgp ( integer  N,
complex  ALPHA,
complex, dimension( * )  X,
integer  INCX,
complex  TAU 
)

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

Download CLARFGP + dependencies [TGZ] [ZIP] [TXT]

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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2015

Definition at line 106 of file clarfgp.f.

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 alpha, tau
115 * ..
116 * .. Array Arguments ..
117  COMPLEX x( * )
118 * ..
119 *
120 * =====================================================================
121 *
122 * .. Parameters ..
123  REAL two, one, zero
124  parameter ( two = 2.0e+0, one = 1.0e+0, zero = 0.0e+0 )
125 * ..
126 * .. Local Scalars ..
127  INTEGER j, knt
128  REAL alphi, alphr, beta, bignum, smlnum, xnorm
129  COMPLEX savealpha
130 * ..
131 * .. External Functions ..
132  REAL scnrm2, slamch, slapy3, slapy2
133  COMPLEX cladiv
134  EXTERNAL scnrm2, slamch, slapy3, slapy2, cladiv
135 * ..
136 * .. Intrinsic Functions ..
137  INTRINSIC abs, aimag, cmplx, REAL, sign
138 * ..
139 * .. External Subroutines ..
140  EXTERNAL cscal, csscal
141 * ..
142 * .. Executable Statements ..
143 *
144  IF( n.LE.0 ) THEN
145  tau = zero
146  RETURN
147  END IF
148 *
149  xnorm = scnrm2( n-1, x, incx )
150  alphr = REAL( alpha )
151  alphi = aimag( 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 = slapy2( alphr, alphi )
175  tau = cmplx( 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( slapy3( alphr, alphi, xnorm ), alphr )
186  smlnum = slamch( 'S' ) / slamch( '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 csscal( 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 = scnrm2( n-1, x, incx )
206  alpha = cmplx( alphr, alphi )
207  beta = sign( slapy3( 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/REAL( alpha ))
216  alphr = alphr + xnorm * (xnorm/REAL( alpha ))
217  tau = cmplx( alphr/beta, -alphi/beta )
218  alpha = cmplx( -alphr, alphi )
219  END IF
220  alpha = cladiv( cmplx( 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 = REAL( savealpha )
232  alphi = aimag( 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 = slapy2( alphr, alphi )
245  tau = cmplx( 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 cscal( 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 CLARFGP
271 *
scnrm2
real function scnrm2(N, X, INCX)
SCNRM2
Definition: scnrm2.f:56
slapy3
real function slapy3(X, Y, Z)
SLAPY3 returns sqrt(x2+y2+z2).
Definition: slapy3.f:70
slapy2
real function slapy2(X, Y)
SLAPY2 returns sqrt(x2+y2).
Definition: slapy2.f:65
cscal
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:54
cladiv
complex function cladiv(X, Y)
CLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
Definition: cladiv.f:66
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
csscal
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:54

Here is the call graph for this function:

Here is the caller graph for this function: