```001:       SUBROUTINE ZLARFG( 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: *  ZLARFG 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, with beta real, and x is an
025: *  (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   ONE, ZERO
063:       PARAMETER          ( 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, DZNRM2
072:       EXTERNAL           DLAMCH, DLAPY3, 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  =  I
094: *
095:          TAU = ZERO
096:       ELSE
097: *
098: *        general case
099: *
100:          BETA = -SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
101:          SAFMIN = DLAMCH( 'S' ) / DLAMCH( 'E' )
102:          RSAFMN = ONE / SAFMIN
103: *
104:          KNT = 0
105:          IF( ABS( BETA ).LT.SAFMIN ) THEN
106: *
107: *           XNORM, BETA may be inaccurate; scale X and recompute them
108: *
109:    10       CONTINUE
110:             KNT = KNT + 1
111:             CALL ZDSCAL( N-1, RSAFMN, X, INCX )
112:             BETA = BETA*RSAFMN
113:             ALPHI = ALPHI*RSAFMN
114:             ALPHR = ALPHR*RSAFMN
115:             IF( ABS( BETA ).LT.SAFMIN )
116:      \$         GO TO 10
117: *
118: *           New BETA is at most 1, at least SAFMIN
119: *
120:             XNORM = DZNRM2( N-1, X, INCX )
121:             ALPHA = DCMPLX( ALPHR, ALPHI )
122:             BETA = -SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
123:          END IF
124:          TAU = DCMPLX( ( BETA-ALPHR ) / BETA, -ALPHI / BETA )
125:          ALPHA = ZLADIV( DCMPLX( ONE ), ALPHA-BETA )
126:          CALL ZSCAL( N-1, ALPHA, X, INCX )
127: *
128: *        If ALPHA is subnormal, it may lose relative accuracy
129: *
130:          DO 20 J = 1, KNT
131:             BETA = BETA*SAFMIN
132:  20      CONTINUE
133:          ALPHA = BETA
134:       END IF
135: *
136:       RETURN
137: *
138: *     End of ZLARFG
139: *
140:       END
141: ```