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