001:       SUBROUTINE SLARFP( N, ALPHA, X, INCX, TAU )
002: *
003: *  -- LAPACK auxiliary routine (version 3.2) --
004: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
005: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       INTEGER            INCX, N
010:       REAL               ALPHA, TAU
011: *     ..
012: *     .. Array Arguments ..
013:       REAL               X( * )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  SLARFP generates a real elementary reflector H of order n, such
020: *  that
021: *
022: *        H * ( alpha ) = ( beta ),   H' * H = I.
023: *            (   x   )   (   0  )
024: *
025: *  where alpha and beta are scalars, beta is non-negative, and x is
026: *  an (n-1)-element real vector.  H is represented in the form
027: *
028: *        H = I - tau * ( 1 ) * ( 1 v' ) ,
029: *                      ( v )
030: *
031: *  where tau is a real scalar and v is a real (n-1)-element
032: *  vector.
033: *
034: *  If the elements of x are all zero, then tau = 0 and H is taken to be
035: *  the unit matrix.
036: *
037: *  Otherwise  1 <= tau <= 2.
038: *
039: *  Arguments
040: *  =========
041: *
042: *  N       (input) INTEGER
043: *          The order of the elementary reflector.
044: *
045: *  ALPHA   (input/output) REAL
046: *          On entry, the value alpha.
047: *          On exit, it is overwritten with the value beta.
048: *
049: *  X       (input/output) REAL array, dimension
050: *                         (1+(N-2)*abs(INCX))
051: *          On entry, the vector x.
052: *          On exit, it is overwritten with the vector v.
053: *
054: *  INCX    (input) INTEGER
055: *          The increment between elements of X. INCX > 0.
056: *
057: *  TAU     (output) REAL
058: *          The value tau.
059: *
060: *  =====================================================================
061: *
062: *     .. Parameters ..
063:       REAL               TWO, ONE, ZERO
064:       PARAMETER          ( TWO = 2.0E+0, ONE = 1.0E+0, ZERO = 0.0E+0 )
065: *     ..
066: *     .. Local Scalars ..
067:       INTEGER            J, KNT
068:       REAL               BETA, RSAFMN, SAFMIN, XNORM
069: *     ..
070: *     .. External Functions ..
071:       REAL               SLAMCH, SLAPY2, SNRM2
072:       EXTERNAL           SLAMCH, SLAPY2, SNRM2
073: *     ..
074: *     .. Intrinsic Functions ..
075:       INTRINSIC          ABS, SIGN
076: *     ..
077: *     .. External Subroutines ..
078:       EXTERNAL           SSCAL
079: *     ..
080: *     .. Executable Statements ..
081: *
082:       IF( N.LE.0 ) THEN
083:          TAU = ZERO
084:          RETURN
085:       END IF
086: *
087:       XNORM = SNRM2( N-1, X, INCX )
088: *
089:       IF( XNORM.EQ.ZERO ) THEN
090: *
091: *        H  =  [+/-1, 0; I], sign chosen so ALPHA >= 0.
092: *
093:          IF( ALPHA.GE.ZERO ) THEN
094: !           When TAU.eq.ZERO, the vector is special-cased to be
095: !           all zeros in the application routines.  We do not need
096: !           to clear it.
097:             TAU = ZERO
098:          ELSE
099: !           However, the application routines rely on explicit
100: !           zero checks when TAU.ne.ZERO, and we must clear X.
101:             TAU = TWO
102:             DO J = 1, N-1
103:                X( 1 + (J-1)*INCX ) = 0
104:             END DO
105:             ALPHA = -ALPHA
106:          END IF
107:       ELSE
108: *
109: *        general case
110: *
111:          BETA = SIGN( SLAPY2( ALPHA, XNORM ), ALPHA )
112:          SAFMIN = SLAMCH( 'S' ) / SLAMCH( 'E' )
113:          KNT = 0
114:          IF( ABS( BETA ).LT.SAFMIN ) THEN
115: *
116: *           XNORM, BETA may be inaccurate; scale X and recompute them
117: *
118:             RSAFMN = ONE / SAFMIN
119:    10       CONTINUE
120:             KNT = KNT + 1
121:             CALL SSCAL( N-1, RSAFMN, X, INCX )
122:             BETA = BETA*RSAFMN
123:             ALPHA = ALPHA*RSAFMN
124:             IF( ABS( BETA ).LT.SAFMIN )
125:      $         GO TO 10
126: *
127: *           New BETA is at most 1, at least SAFMIN
128: *
129:             XNORM = SNRM2( N-1, X, INCX )
130:             BETA = SIGN( SLAPY2( ALPHA, XNORM ), ALPHA )
131:          END IF
132:          ALPHA = ALPHA + BETA
133:          IF( BETA.LT.ZERO ) THEN
134:             BETA = -BETA
135:             TAU = -ALPHA / BETA
136:          ELSE
137:             ALPHA = XNORM * (XNORM/ALPHA)
138:             TAU = ALPHA / BETA
139:             ALPHA = -ALPHA
140:          END IF
141:          CALL SSCAL( N-1, ONE / ALPHA, X, INCX )
142: *
143: *        If BETA is subnormal, it may lose relative accuracy
144: *
145:          DO 20 J = 1, KNT
146:             BETA = BETA*SAFMIN
147:  20      CONTINUE
148:          ALPHA = BETA
149:       END IF
150: *
151:       RETURN
152: *
153: *     End of SLARFP
154: *
155:       END
156: