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