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*> \date November 2015
101*
102*> \ingroup realOTHERauxiliary
103*
104*  =====================================================================
105      SUBROUTINE SLARFGP( 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      REAL               ALPHA, TAU
115*     ..
116*     .. Array Arguments ..
117      REAL               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               BETA, BIGNUM, SAVEALPHA, SMLNUM, XNORM
129*     ..
130*     .. External Functions ..
131      REAL               SLAMCH, SLAPY2, SNRM2
132      EXTERNAL           SLAMCH, SLAPY2, SNRM2
133*     ..
134*     .. Intrinsic Functions ..
135      INTRINSIC          ABS, SIGN
136*     ..
137*     .. External Subroutines ..
138      EXTERNAL           SSCAL
139*     ..
140*     .. Executable Statements ..
141*
142      IF( N.LE.0 ) THEN
143         TAU = ZERO
144         RETURN
145      END IF
146*
147      XNORM = SNRM2( N-1, X, INCX )
148*
149      IF( XNORM.EQ.ZERO ) THEN
150*
151*        H  =  [+/-1, 0; I], sign chosen so ALPHA >= 0.
152*
153         IF( ALPHA.GE.ZERO ) THEN
154*           When TAU.eq.ZERO, the vector is special-cased to be
155*           all zeros in the application routines.  We do not need
156*           to clear it.
157            TAU = ZERO
158         ELSE
159*           However, the application routines rely on explicit
160*           zero checks when TAU.ne.ZERO, and we must clear X.
161            TAU = TWO
162            DO J = 1, N-1
163               X( 1 + (J-1)*INCX ) = 0
164            END DO
165            ALPHA = -ALPHA
166         END IF
167      ELSE
168*
169*        general case
170*
171         BETA = SIGN( SLAPY2( ALPHA, XNORM ), ALPHA )
172         SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'E' )
173         KNT = 0
174         IF( ABS( BETA ).LT.SMLNUM ) THEN
175*
176*           XNORM, BETA may be inaccurate; scale X and recompute them
177*
178            BIGNUM = ONE / SMLNUM
179   10       CONTINUE
180            KNT = KNT + 1
181            CALL SSCAL( N-1, BIGNUM, X, INCX )
182            BETA = BETA*BIGNUM
183            ALPHA = ALPHA*BIGNUM
184            IF( ABS( BETA ).LT.SMLNUM )
185     $         GO TO 10
186*
187*           New BETA is at most 1, at least SMLNUM
188*
189            XNORM = SNRM2( N-1, X, INCX )
190            BETA = SIGN( SLAPY2( ALPHA, XNORM ), ALPHA )
191         END IF
192         SAVEALPHA = ALPHA
193         ALPHA = ALPHA + BETA
194         IF( BETA.LT.ZERO ) THEN
195            BETA = -BETA
196            TAU = -ALPHA / BETA
197         ELSE
198            ALPHA = XNORM * (XNORM/ALPHA)
199            TAU = ALPHA / BETA
200            ALPHA = -ALPHA
201         END IF
202*
203         IF ( ABS(TAU).LE.SMLNUM ) THEN
204*
205*           In the case where the computed TAU ends up being a denormalized number,
206*           it loses relative accuracy. This is a BIG problem. Solution: flush TAU
207*           to ZERO. This explains the next IF statement.
208*
209*           (Bug report provided by Pat Quillen from MathWorks on Jul 29, 2009.)
210*           (Thanks Pat. Thanks MathWorks.)
211*
212            IF( SAVEALPHA.GE.ZERO ) THEN
213               TAU = ZERO
214            ELSE
215               TAU = TWO
216               DO J = 1, N-1
217                  X( 1 + (J-1)*INCX ) = 0
218               END DO
219               BETA = -SAVEALPHA
220            END IF
221*
222         ELSE
223*
224*           This is the general case.
225*
226            CALL SSCAL( N-1, ONE / ALPHA, X, INCX )
227*
228         END IF
229*
230*        If BETA is subnormal, it may lose relative accuracy
231*
232         DO 20 J = 1, KNT
233            BETA = BETA*SMLNUM
234 20      CONTINUE
235         ALPHA = BETA
236      END IF
237*
238      RETURN
239*
240*     End of SLARFGP
241*
242      END
243