LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
slarfgp.f
Go to the documentation of this file.
1*> \brief \b SLARFGP 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 SLARFGP + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarfgp.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarfgp.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarfgp.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SLARFGP( N, ALPHA, X, INCX, TAU )
22*
23* .. Scalar Arguments ..
24* INTEGER INCX, N
25* REAL ALPHA, TAU
26* ..
27* .. Array Arguments ..
28* REAL X( * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> SLARFGP generates a real elementary reflector H of order n, such
38*> that
39*>
40*> H * ( alpha ) = ( beta ), H**T * H = I.
41*> ( x ) ( 0 )
42*>
43*> where alpha and beta are scalars, beta is non-negative, and x is
44*> an (n-1)-element real vector. H is represented in the form
45*>
46*> H = I - tau * ( 1 ) * ( 1 v**T ) ,
47*> ( v )
48*>
49*> where tau is a real scalar and v is a real (n-1)-element
50*> vector.
51*>
52*> If the elements of x are all zero, then tau = 0 and H is taken to be
53*> 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 REAL
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 REAL 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 REAL
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 slarfgp( 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 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, EPS, 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 eps = slamch( 'Precision' )
145 xnorm = snrm2( n-1, x, incx )
146*
147 IF( xnorm.LE.eps*abs(alpha) ) THEN
148*
149* H = [+/-1, 0; I], sign chosen so ALPHA >= 0.
150*
151 IF( alpha.GE.zero ) THEN
152* When TAU.eq.ZERO, the vector is special-cased to be
153* all zeros in the application routines. We do not need
154* to clear it.
155 tau = zero
156 ELSE
157* However, the application routines rely on explicit
158* zero checks when TAU.ne.ZERO, and we must clear X.
159 tau = two
160 DO j = 1, n-1
161 x( 1 + (j-1)*incx ) = 0
162 END DO
163 alpha = -alpha
164 END IF
165 ELSE
166*
167* general case
168*
169 beta = sign( slapy2( alpha, xnorm ), alpha )
170 smlnum = slamch( 'S' ) / slamch( 'E' )
171 knt = 0
172 IF( abs( beta ).LT.smlnum ) THEN
173*
174* XNORM, BETA may be inaccurate; scale X and recompute them
175*
176 bignum = one / smlnum
177 10 CONTINUE
178 knt = knt + 1
179 CALL sscal( n-1, bignum, x, incx )
180 beta = beta*bignum
181 alpha = alpha*bignum
182 IF( (abs( beta ).LT.smlnum) .AND. (knt .LT. 20) )
183 $ GO TO 10
184*
185* New BETA is at most 1, at least SMLNUM
186*
187 xnorm = snrm2( n-1, x, incx )
188 beta = sign( slapy2( alpha, xnorm ), alpha )
189 END IF
190 savealpha = alpha
191 alpha = alpha + beta
192 IF( beta.LT.zero ) THEN
193 beta = -beta
194 tau = -alpha / beta
195 ELSE
196 alpha = xnorm * (xnorm/alpha)
197 tau = alpha / beta
198 alpha = -alpha
199 END IF
200*
201 IF ( abs(tau).LE.smlnum ) THEN
202*
203* In the case where the computed TAU ends up being a denormalized number,
204* it loses relative accuracy. This is a BIG problem. Solution: flush TAU
205* to ZERO. This explains the next IF statement.
206*
207* (Bug report provided by Pat Quillen from MathWorks on Jul 29, 2009.)
208* (Thanks Pat. Thanks MathWorks.)
209*
210 IF( savealpha.GE.zero ) THEN
211 tau = zero
212 ELSE
213 tau = two
214 DO j = 1, n-1
215 x( 1 + (j-1)*incx ) = 0
216 END DO
217 beta = -savealpha
218 END IF
219*
220 ELSE
221*
222* This is the general case.
223*
224 CALL sscal( n-1, one / alpha, x, incx )
225*
226 END IF
227*
228* If BETA is subnormal, it may lose relative accuracy
229*
230 DO 20 j = 1, knt
231 beta = beta*smlnum
232 20 CONTINUE
233 alpha = beta
234 END IF
235*
236 RETURN
237*
238* End of SLARFGP
239*
240 END
subroutine slarfgp(n, alpha, x, incx, tau)
SLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.
Definition slarfgp.f:104
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79