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

◆ slarfgp()

subroutine slarfgp ( integer  N,
real  ALPHA,
real, dimension( * )  X,
integer  INCX,
real  TAU 
)

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

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

Purpose:
 SLARFGP generates a real elementary reflector H of order n, such
 that

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

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

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

 where tau is a real scalar and v is a real (n-1)-element
 vector.

 If the elements of x are all zero, 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 REAL
          On entry, the value alpha.
          On exit, it is overwritten with the value beta.
[in,out]X
          X is REAL 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 REAL
          The value tau.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 103 of file slarfgp.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 REAL ALPHA, TAU
112* ..
113* .. Array Arguments ..
114 REAL 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 BETA, BIGNUM, SAVEALPHA, SMLNUM, XNORM
126* ..
127* .. External Functions ..
128 REAL SLAMCH, SLAPY2, SNRM2
129 EXTERNAL slamch, slapy2, snrm2
130* ..
131* .. Intrinsic Functions ..
132 INTRINSIC abs, sign
133* ..
134* .. External Subroutines ..
135 EXTERNAL sscal
136* ..
137* .. Executable Statements ..
138*
139 IF( n.LE.0 ) THEN
140 tau = zero
141 RETURN
142 END IF
143*
144 xnorm = snrm2( n-1, x, incx )
145*
146 IF( xnorm.EQ.zero ) THEN
147*
148* H = [+/-1, 0; I], sign chosen so ALPHA >= 0.
149*
150 IF( alpha.GE.zero ) THEN
151* When TAU.eq.ZERO, the vector is special-cased to be
152* all zeros in the application routines. We do not need
153* to clear it.
154 tau = zero
155 ELSE
156* However, the application routines rely on explicit
157* zero checks when TAU.ne.ZERO, and we must clear X.
158 tau = two
159 DO j = 1, n-1
160 x( 1 + (j-1)*incx ) = 0
161 END DO
162 alpha = -alpha
163 END IF
164 ELSE
165*
166* general case
167*
168 beta = sign( slapy2( alpha, xnorm ), alpha )
169 smlnum = slamch( 'S' ) / slamch( 'E' )
170 knt = 0
171 IF( abs( beta ).LT.smlnum ) THEN
172*
173* XNORM, BETA may be inaccurate; scale X and recompute them
174*
175 bignum = one / smlnum
176 10 CONTINUE
177 knt = knt + 1
178 CALL sscal( n-1, bignum, x, incx )
179 beta = beta*bignum
180 alpha = alpha*bignum
181 IF( (abs( beta ).LT.smlnum) .AND. (knt .LT. 20) )
182 $ GO TO 10
183*
184* New BETA is at most 1, at least SMLNUM
185*
186 xnorm = snrm2( n-1, x, incx )
187 beta = sign( slapy2( alpha, xnorm ), alpha )
188 END IF
189 savealpha = alpha
190 alpha = alpha + beta
191 IF( beta.LT.zero ) THEN
192 beta = -beta
193 tau = -alpha / beta
194 ELSE
195 alpha = xnorm * (xnorm/alpha)
196 tau = alpha / beta
197 alpha = -alpha
198 END IF
199*
200 IF ( abs(tau).LE.smlnum ) THEN
201*
202* In the case where the computed TAU ends up being a denormalized number,
203* it loses relative accuracy. This is a BIG problem. Solution: flush TAU
204* to ZERO. This explains the next IF statement.
205*
206* (Bug report provided by Pat Quillen from MathWorks on Jul 29, 2009.)
207* (Thanks Pat. Thanks MathWorks.)
208*
209 IF( savealpha.GE.zero ) THEN
210 tau = zero
211 ELSE
212 tau = two
213 DO j = 1, n-1
214 x( 1 + (j-1)*incx ) = 0
215 END DO
216 beta = -savealpha
217 END IF
218*
219 ELSE
220*
221* This is the general case.
222*
223 CALL sscal( n-1, one / alpha, x, incx )
224*
225 END IF
226*
227* If BETA is subnormal, it may lose relative accuracy
228*
229 DO 20 j = 1, knt
230 beta = beta*smlnum
231 20 CONTINUE
232 alpha = beta
233 END IF
234*
235 RETURN
236*
237* End of SLARFGP
238*
real function slapy2(X, Y)
SLAPY2 returns sqrt(x2+y2).
Definition: slapy2.f:63
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
real(wp) function snrm2(n, x, incx)
SNRM2
Definition: snrm2.f90:89
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: