1       SUBROUTINE ZLARFG( N, ALPHA, X, INCX, TAU )
  2 *
  3 *  -- LAPACK auxiliary routine (version 3.3.1) --
  4 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  5 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  6 *  -- April 2011                                                      --
  7 *
  8 *     .. Scalar Arguments ..
  9       INTEGER            INCX, N
 10       COMPLEX*16         ALPHA, TAU
 11 *     ..
 12 *     .. Array Arguments ..
 13       COMPLEX*16         X( * )
 14 *     ..
 15 *
 16 *  Purpose
 17 *  =======
 18 *
 19 *  ZLARFG generates a complex elementary reflector H of order n, such
 20 *  that
 21 *
 22 *        H**H * ( alpha ) = ( beta ),   H**H * H = I.
 23 *               (   x   )   (   0  )
 24 *
 25 *  where alpha and beta are scalars, with beta real, and x is an
 26 *  (n-1)-element complex vector. H is represented in the form
 27 *
 28 *        H = I - tau * ( 1 ) * ( 1 v**H ) ,
 29 *                      ( v )
 30 *
 31 *  where tau is a complex scalar and v is a complex (n-1)-element
 32 *  vector. Note that H is not hermitian.
 33 *
 34 *  If the elements of x are all zero and alpha is real, then tau = 0
 35 *  and H is taken to be the unit matrix.
 36 *
 37 *  Otherwise  1 <= real(tau) <= 2  and  abs(tau-1) <= 1 .
 38 *
 39 *  Arguments
 40 *  =========
 41 *
 42 *  N       (input) INTEGER
 43 *          The order of the elementary reflector.
 44 *
 45 *  ALPHA   (input/output) COMPLEX*16
 46 *          On entry, the value alpha.
 47 *          On exit, it is overwritten with the value beta.
 48 *
 49 *  X       (input/output) COMPLEX*16 array, dimension
 50 *                         (1+(N-2)*abs(INCX))
 51 *          On entry, the vector x.
 52 *          On exit, it is overwritten with the vector v.
 53 *
 54 *  INCX    (input) INTEGER
 55 *          The increment between elements of X. INCX > 0.
 56 *
 57 *  TAU     (output) COMPLEX*16
 58 *          The value tau.
 59 *
 60 *  =====================================================================
 61 *
 62 *     .. Parameters ..
 63       DOUBLE PRECISION   ONE, ZERO
 64       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 65 *     ..
 66 *     .. Local Scalars ..
 67       INTEGER            J, KNT
 68       DOUBLE PRECISION   ALPHI, ALPHR, BETA, RSAFMN, SAFMIN, XNORM
 69 *     ..
 70 *     .. External Functions ..
 71       DOUBLE PRECISION   DLAMCH, DLAPY3, DZNRM2
 72       COMPLEX*16         ZLADIV
 73       EXTERNAL           DLAMCH, DLAPY3, DZNRM2, ZLADIV
 74 *     ..
 75 *     .. Intrinsic Functions ..
 76       INTRINSIC          ABSDBLEDCMPLXDIMAGSIGN
 77 *     ..
 78 *     .. External Subroutines ..
 79       EXTERNAL           ZDSCAL, ZSCAL
 80 *     ..
 81 *     .. Executable Statements ..
 82 *
 83       IF( N.LE.0 ) THEN
 84          TAU = ZERO
 85          RETURN
 86       END IF
 87 *
 88       XNORM = DZNRM2( N-1, X, INCX )
 89       ALPHR = DBLE( ALPHA )
 90       ALPHI = DIMAG( ALPHA )
 91 *
 92       IF( XNORM.EQ.ZERO .AND. ALPHI.EQ.ZERO ) THEN
 93 *
 94 *        H  =  I
 95 *
 96          TAU = ZERO
 97       ELSE
 98 *
 99 *        general case
100 *
101          BETA = -SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
102          SAFMIN = DLAMCH( 'S' ) / DLAMCH( 'E' )
103          RSAFMN = ONE / SAFMIN
104 *
105          KNT = 0
106          IFABS( BETA ).LT.SAFMIN ) THEN
107 *
108 *           XNORM, BETA may be inaccurate; scale X and recompute them
109 *
110    10       CONTINUE
111             KNT = KNT + 1
112             CALL ZDSCAL( N-1, RSAFMN, X, INCX )
113             BETA = BETA*RSAFMN
114             ALPHI = ALPHI*RSAFMN
115             ALPHR = ALPHR*RSAFMN
116             IFABS( BETA ).LT.SAFMIN )
117      $         GO TO 10
118 *
119 *           New BETA is at most 1, at least SAFMIN
120 *
121             XNORM = DZNRM2( N-1, X, INCX )
122             ALPHA = DCMPLX( ALPHR, ALPHI )
123             BETA = -SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
124          END IF
125          TAU = DCMPLX( ( BETA-ALPHR ) / BETA, -ALPHI / BETA )
126          ALPHA = ZLADIV( DCMPLX( ONE ), ALPHA-BETA )
127          CALL ZSCAL( N-1, ALPHA, X, INCX )
128 *
129 *        If ALPHA is subnormal, it may lose relative accuracy
130 *
131          DO 20 J = 1, KNT
132             BETA = BETA*SAFMIN
133  20      CONTINUE
134          ALPHA = BETA
135       END IF
136 *
137       RETURN
138 *
139 *     End of ZLARFG
140 *
141       END