001:       SUBROUTINE ZLARFP( 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:       COMPLEX*16         ALPHA, TAU
010: *     ..
011: *     .. Array Arguments ..
012:       COMPLEX*16         X( * )
013: *     ..
014: *
015: *  Purpose
016: *  =======
017: *
018: *  ZLARFP generates a complex 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 real and non-negative, and
025: *  x is an (n-1)-element complex vector.  H is represented in the form
026: *
027: *        H = I - tau * ( 1 ) * ( 1 v' ) ,
028: *                      ( v )
029: *
030: *  where tau is a complex scalar and v is a complex (n-1)-element
031: *  vector. Note that H is not hermitian.
032: *
033: *  If the elements of x are all zero and alpha is real, then tau = 0
034: *  and H is taken to be the unit matrix.
035: *
036: *  Otherwise  1 <= real(tau) <= 2  and  abs(tau-1) <= 1 .
037: *
038: *  Arguments
039: *  =========
040: *
041: *  N       (input) INTEGER
042: *          The order of the elementary reflector.
043: *
044: *  ALPHA   (input/output) COMPLEX*16
045: *          On entry, the value alpha.
046: *          On exit, it is overwritten with the value beta.
047: *
048: *  X       (input/output) COMPLEX*16 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) COMPLEX*16
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   ALPHI, ALPHR, BETA, RSAFMN, SAFMIN, XNORM
068: *     ..
069: *     .. External Functions ..
070:       DOUBLE PRECISION   DLAMCH, DLAPY3, DLAPY2, DZNRM2
071:       COMPLEX*16         ZLADIV
072:       EXTERNAL           DLAMCH, DLAPY3, DLAPY2, DZNRM2, ZLADIV
073: *     ..
074: *     .. Intrinsic Functions ..
075:       INTRINSIC          ABS, DBLE, DCMPLX, DIMAG, SIGN
076: *     ..
077: *     .. External Subroutines ..
078:       EXTERNAL           ZDSCAL, ZSCAL
079: *     ..
080: *     .. Executable Statements ..
081: *
082:       IF( N.LE.0 ) THEN
083:          TAU = ZERO
084:          RETURN
085:       END IF
086: *
087:       XNORM = DZNRM2( N-1, X, INCX )
088:       ALPHR = DBLE( ALPHA )
089:       ALPHI = DIMAG( ALPHA )
090: *
091:       IF( XNORM.EQ.ZERO .AND. ALPHI.EQ.ZERO ) THEN
092: *
093: *        H  =  [1-alpha/abs(alpha) 0; 0 I], sign chosen so ALPHA >= 0.
094: *
095:          IF( ALPHI.EQ.ZERO ) THEN
096:             IF( ALPHR.GE.ZERO ) THEN
097: !              When TAU.eq.ZERO, the vector is special-cased to be
098: !              all zeros in the application routines.  We do not need
099: !              to clear it.
100:                TAU = ZERO
101:             ELSE
102: !              However, the application routines rely on explicit
103: !              zero checks when TAU.ne.ZERO, and we must clear X.
104:                TAU = TWO
105:                DO J = 1, N-1
106:                   X( 1 + (J-1)*INCX ) = 0
107:                END DO
108:                ALPHA = -ALPHA
109:             END IF
110:          ELSE
111: !           Only "reflecting" the diagonal entry to be real and non-negative.
112:             XNORM = DLAPY2( ALPHR, ALPHI )
113:             TAU = CMPLX( ONE - ALPHR / XNORM, -ALPHI / XNORM )
114:             DO J = 1, N-1
115:                X( 1 + (J-1)*INCX ) = 0
116:             END DO
117:             ALPHA = XNORM
118:          END IF
119:       ELSE
120: *
121: *        general case
122: *
123:          BETA = SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
124:          SAFMIN = DLAMCH( 'S' ) / DLAMCH( 'E' )
125:          RSAFMN = ONE / SAFMIN
126: *
127:          KNT = 0
128:          IF( ABS( BETA ).LT.SAFMIN ) THEN
129: *
130: *           XNORM, BETA may be inaccurate; scale X and recompute them
131: *
132:    10       CONTINUE
133:             KNT = KNT + 1
134:             CALL ZDSCAL( N-1, RSAFMN, X, INCX )
135:             BETA = BETA*RSAFMN
136:             ALPHI = ALPHI*RSAFMN
137:             ALPHR = ALPHR*RSAFMN
138:             IF( ABS( BETA ).LT.SAFMIN )
139:      $         GO TO 10
140: *
141: *           New BETA is at most 1, at least SAFMIN
142: *
143:             XNORM = DZNRM2( N-1, X, INCX )
144:             ALPHA = DCMPLX( ALPHR, ALPHI )
145:             BETA = SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
146:          END IF
147:          ALPHA = ALPHA + BETA
148:          IF( BETA.LT.ZERO ) THEN
149:             BETA = -BETA
150:             TAU = -ALPHA / BETA
151:          ELSE
152:             ALPHR = ALPHI * (ALPHI/DBLE( ALPHA ))
153:             ALPHR = ALPHR + XNORM * (XNORM/DBLE( ALPHA ))
154:             TAU = DCMPLX( ALPHR/BETA, -ALPHI/BETA )
155:             ALPHA = DCMPLX( -ALPHR, ALPHI )
156:          END IF
157:          ALPHA = ZLADIV( DCMPLX( ONE ), ALPHA )
158:          CALL ZSCAL( N-1, ALPHA, X, INCX )
159: *
160: *        If BETA is subnormal, it may lose relative accuracy
161: *
162:          DO 20 J = 1, KNT
163:             BETA = BETA*SAFMIN
164:  20      CONTINUE
165:          ALPHA = BETA
166:       END IF
167: *
168:       RETURN
169: *
170: *     End of ZLARFP
171: *
172:       END
173: