1 SUBROUTINE ZLARFGP( 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 * ZLARFGP 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, beta is real and non-negative, and
26 * x is an (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 * Arguments
38 * =========
39 *
40 * N (input) INTEGER
41 * The order of the elementary reflector.
42 *
43 * ALPHA (input/output) COMPLEX*16
44 * On entry, the value alpha.
45 * On exit, it is overwritten with the value beta.
46 *
47 * X (input/output) COMPLEX*16 array, dimension
48 * (1+(N-2)*abs(INCX))
49 * On entry, the vector x.
50 * On exit, it is overwritten with the vector v.
51 *
52 * INCX (input) INTEGER
53 * The increment between elements of X. INCX > 0.
54 *
55 * TAU (output) COMPLEX*16
56 * The value tau.
57 *
58 * =====================================================================
59 *
60 * .. Parameters ..
61 DOUBLE PRECISION TWO, ONE, ZERO
62 PARAMETER ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 )
63 * ..
64 * .. Local Scalars ..
65 INTEGER J, KNT
66 DOUBLE PRECISION ALPHI, ALPHR, BETA, BIGNUM, SMLNUM, XNORM
67 COMPLEX*16 SAVEALPHA
68 * ..
69 * .. External Functions ..
70 DOUBLE PRECISION DLAMCH, DLAPY3, DLAPY2, DZNRM2
71 COMPLEX*16 ZLADIV
72 EXTERNAL DLAMCH, DLAPY3, DLAPY2, DZNRM2, ZLADIV
73 * ..
74 * .. Intrinsic Functions ..
75 INTRINSIC ABS, DBLE, DCMPLX, DIMAG, SIGN
76 * ..
77 * .. External Subroutines ..
78 EXTERNAL ZDSCAL, ZSCAL
79 * ..
80 * .. Executable Statements ..
81 *
82 IF( N.LE.0 ) THEN
83 TAU = ZERO
84 RETURN
85 END IF
86 *
87 XNORM = DZNRM2( N-1, X, INCX )
88 ALPHR = DBLE( ALPHA )
89 ALPHI = DIMAG( ALPHA )
90 *
91 IF( XNORM.EQ.ZERO ) THEN
92 *
93 * H = [1-alpha/abs(alpha) 0; 0 I], sign chosen so ALPHA >= 0.
94 *
95 IF( ALPHI.EQ.ZERO ) THEN
96 IF( ALPHR.GE.ZERO ) THEN
97 * When TAU.eq.ZERO, the vector is special-cased to be
98 * all zeros in the application routines. We do not need
99 * 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 ) = ZERO
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 = DCMPLX( ONE - ALPHR / XNORM, -ALPHI / XNORM )
114 DO J = 1, N-1
115 X( 1 + (J-1)*INCX ) = ZERO
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 SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'E' )
125 BIGNUM = ONE / SMLNUM
126 *
127 KNT = 0
128 IF( ABS( BETA ).LT.SMLNUM ) 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, BIGNUM, X, INCX )
135 BETA = BETA*BIGNUM
136 ALPHI = ALPHI*BIGNUM
137 ALPHR = ALPHR*BIGNUM
138 IF( ABS( BETA ).LT.SMLNUM )
139 $ GO TO 10
140 *
141 * New BETA is at most 1, at least SMLNUM
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 SAVEALPHA = ALPHA
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 *
160 IF ( ABS(TAU).LE.SMLNUM ) THEN
161 *
162 * In the case where the computed TAU ends up being a denormalized number,
163 * it loses relative accuracy. This is a BIG problem. Solution: flush TAU
164 * to ZERO (or TWO or whatever makes a nonnegative real number for BETA).
165 *
166 * (Bug report provided by Pat Quillen from MathWorks on Jul 29, 2009.)
167 * (Thanks Pat. Thanks MathWorks.)
168 *
169 ALPHR = DBLE( SAVEALPHA )
170 ALPHI = DIMAG( SAVEALPHA )
171 IF( ALPHI.EQ.ZERO ) THEN
172 IF( ALPHR.GE.ZERO ) THEN
173 TAU = ZERO
174 ELSE
175 TAU = TWO
176 DO J = 1, N-1
177 X( 1 + (J-1)*INCX ) = ZERO
178 END DO
179 BETA = -SAVEALPHA
180 END IF
181 ELSE
182 XNORM = DLAPY2( ALPHR, ALPHI )
183 TAU = DCMPLX( ONE - ALPHR / XNORM, -ALPHI / XNORM )
184 DO J = 1, N-1
185 X( 1 + (J-1)*INCX ) = ZERO
186 END DO
187 BETA = XNORM
188 END IF
189 *
190 ELSE
191 *
192 * This is the general case.
193 *
194 CALL ZSCAL( N-1, ALPHA, X, INCX )
195 *
196 END IF
197 *
198 * If BETA is subnormal, it may lose relative accuracy
199 *
200 DO 20 J = 1, KNT
201 BETA = BETA*SMLNUM
202 20 CONTINUE
203 ALPHA = BETA
204 END IF
205 *
206 RETURN
207 *
208 * End of ZLARFGP
209 *
210 END
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 * ZLARFGP 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, beta is real and non-negative, and
26 * x is an (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 * Arguments
38 * =========
39 *
40 * N (input) INTEGER
41 * The order of the elementary reflector.
42 *
43 * ALPHA (input/output) COMPLEX*16
44 * On entry, the value alpha.
45 * On exit, it is overwritten with the value beta.
46 *
47 * X (input/output) COMPLEX*16 array, dimension
48 * (1+(N-2)*abs(INCX))
49 * On entry, the vector x.
50 * On exit, it is overwritten with the vector v.
51 *
52 * INCX (input) INTEGER
53 * The increment between elements of X. INCX > 0.
54 *
55 * TAU (output) COMPLEX*16
56 * The value tau.
57 *
58 * =====================================================================
59 *
60 * .. Parameters ..
61 DOUBLE PRECISION TWO, ONE, ZERO
62 PARAMETER ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 )
63 * ..
64 * .. Local Scalars ..
65 INTEGER J, KNT
66 DOUBLE PRECISION ALPHI, ALPHR, BETA, BIGNUM, SMLNUM, XNORM
67 COMPLEX*16 SAVEALPHA
68 * ..
69 * .. External Functions ..
70 DOUBLE PRECISION DLAMCH, DLAPY3, DLAPY2, DZNRM2
71 COMPLEX*16 ZLADIV
72 EXTERNAL DLAMCH, DLAPY3, DLAPY2, DZNRM2, ZLADIV
73 * ..
74 * .. Intrinsic Functions ..
75 INTRINSIC ABS, DBLE, DCMPLX, DIMAG, SIGN
76 * ..
77 * .. External Subroutines ..
78 EXTERNAL ZDSCAL, ZSCAL
79 * ..
80 * .. Executable Statements ..
81 *
82 IF( N.LE.0 ) THEN
83 TAU = ZERO
84 RETURN
85 END IF
86 *
87 XNORM = DZNRM2( N-1, X, INCX )
88 ALPHR = DBLE( ALPHA )
89 ALPHI = DIMAG( ALPHA )
90 *
91 IF( XNORM.EQ.ZERO ) THEN
92 *
93 * H = [1-alpha/abs(alpha) 0; 0 I], sign chosen so ALPHA >= 0.
94 *
95 IF( ALPHI.EQ.ZERO ) THEN
96 IF( ALPHR.GE.ZERO ) THEN
97 * When TAU.eq.ZERO, the vector is special-cased to be
98 * all zeros in the application routines. We do not need
99 * 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 ) = ZERO
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 = DCMPLX( ONE - ALPHR / XNORM, -ALPHI / XNORM )
114 DO J = 1, N-1
115 X( 1 + (J-1)*INCX ) = ZERO
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 SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'E' )
125 BIGNUM = ONE / SMLNUM
126 *
127 KNT = 0
128 IF( ABS( BETA ).LT.SMLNUM ) 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, BIGNUM, X, INCX )
135 BETA = BETA*BIGNUM
136 ALPHI = ALPHI*BIGNUM
137 ALPHR = ALPHR*BIGNUM
138 IF( ABS( BETA ).LT.SMLNUM )
139 $ GO TO 10
140 *
141 * New BETA is at most 1, at least SMLNUM
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 SAVEALPHA = ALPHA
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 *
160 IF ( ABS(TAU).LE.SMLNUM ) THEN
161 *
162 * In the case where the computed TAU ends up being a denormalized number,
163 * it loses relative accuracy. This is a BIG problem. Solution: flush TAU
164 * to ZERO (or TWO or whatever makes a nonnegative real number for BETA).
165 *
166 * (Bug report provided by Pat Quillen from MathWorks on Jul 29, 2009.)
167 * (Thanks Pat. Thanks MathWorks.)
168 *
169 ALPHR = DBLE( SAVEALPHA )
170 ALPHI = DIMAG( SAVEALPHA )
171 IF( ALPHI.EQ.ZERO ) THEN
172 IF( ALPHR.GE.ZERO ) THEN
173 TAU = ZERO
174 ELSE
175 TAU = TWO
176 DO J = 1, N-1
177 X( 1 + (J-1)*INCX ) = ZERO
178 END DO
179 BETA = -SAVEALPHA
180 END IF
181 ELSE
182 XNORM = DLAPY2( ALPHR, ALPHI )
183 TAU = DCMPLX( ONE - ALPHR / XNORM, -ALPHI / XNORM )
184 DO J = 1, N-1
185 X( 1 + (J-1)*INCX ) = ZERO
186 END DO
187 BETA = XNORM
188 END IF
189 *
190 ELSE
191 *
192 * This is the general case.
193 *
194 CALL ZSCAL( N-1, ALPHA, X, INCX )
195 *
196 END IF
197 *
198 * If BETA is subnormal, it may lose relative accuracy
199 *
200 DO 20 J = 1, KNT
201 BETA = BETA*SMLNUM
202 20 CONTINUE
203 ALPHA = BETA
204 END IF
205 *
206 RETURN
207 *
208 * End of ZLARFGP
209 *
210 END