1: *> \brief \b ZLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZLARFGP + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlarfgp.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlarfgp.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlarfgp.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZLARFGP( N, ALPHA, X, INCX, TAU )
22: *
23: * .. Scalar Arguments ..
24: * INTEGER INCX, N
25: * COMPLEX*16 ALPHA, TAU
26: * ..
27: * .. Array Arguments ..
28: * COMPLEX*16 X( * )
29: * ..
30: *
31: *
32: *> \par Purpose:
33: * =============
34: *>
35: *> \verbatim
36: *>
37: *> ZLARFGP generates a complex elementary reflector H of order n, such
38: *> that
39: *>
40: *> H**H * ( alpha ) = ( beta ), H**H * H = I.
41: *> ( x ) ( 0 )
42: *>
43: *> where alpha and beta are scalars, beta is real and non-negative, and
44: *> x is an (n-1)-element complex vector. H is represented in the form
45: *>
46: *> H = I - tau * ( 1 ) * ( 1 v**H ) ,
47: *> ( v )
48: *>
49: *> where tau is a complex scalar and v is a complex (n-1)-element
50: *> vector. Note that H is not hermitian.
51: *>
52: *> If the elements of x are all zero and alpha is real, then tau = 0
53: *> and H is taken to be the unit matrix.
54: *> \endverbatim
55: *
56: * Arguments:
57: * ==========
58: *
59: *> \param[in] N
60: *> \verbatim
61: *> N is INTEGER
62: *> The order of the elementary reflector.
63: *> \endverbatim
64: *>
65: *> \param[in,out] ALPHA
66: *> \verbatim
67: *> ALPHA is COMPLEX*16
68: *> On entry, the value alpha.
69: *> On exit, it is overwritten with the value beta.
70: *> \endverbatim
71: *>
72: *> \param[in,out] X
73: *> \verbatim
74: *> X is COMPLEX*16 array, dimension
75: *> (1+(N-2)*abs(INCX))
76: *> On entry, the vector x.
77: *> On exit, it is overwritten with the vector v.
78: *> \endverbatim
79: *>
80: *> \param[in] INCX
81: *> \verbatim
82: *> INCX is INTEGER
83: *> The increment between elements of X. INCX > 0.
84: *> \endverbatim
85: *>
86: *> \param[out] TAU
87: *> \verbatim
88: *> TAU is COMPLEX*16
89: *> The value tau.
90: *> \endverbatim
91: *
92: * Authors:
93: * ========
94: *
95: *> \author Univ. of Tennessee
96: *> \author Univ. of California Berkeley
97: *> \author Univ. of Colorado Denver
98: *> \author NAG Ltd.
99: *
100: *> \ingroup complex16OTHERauxiliary
101: *
102: * =====================================================================
103: SUBROUTINE ZLARFGP( N, ALPHA, X, INCX, TAU )
104: *
105: * -- LAPACK auxiliary routine --
106: * -- LAPACK is a software package provided by Univ. of Tennessee, --
107: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
108: *
109: * .. Scalar Arguments ..
110: INTEGER INCX, N
111: COMPLEX*16 ALPHA, TAU
112: * ..
113: * .. Array Arguments ..
114: COMPLEX*16 X( * )
115: * ..
116: *
117: * =====================================================================
118: *
119: * .. Parameters ..
120: DOUBLE PRECISION TWO, ONE, ZERO
121: PARAMETER ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 )
122: * ..
123: * .. Local Scalars ..
124: INTEGER J, KNT
125: DOUBLE PRECISION ALPHI, ALPHR, BETA, BIGNUM, SMLNUM, XNORM
126: COMPLEX*16 SAVEALPHA
127: * ..
128: * .. External Functions ..
129: DOUBLE PRECISION DLAMCH, DLAPY3, DLAPY2, DZNRM2
130: COMPLEX*16 ZLADIV
131: EXTERNAL DLAMCH, DLAPY3, DLAPY2, DZNRM2, ZLADIV
132: * ..
133: * .. Intrinsic Functions ..
134: INTRINSIC ABS, DBLE, DCMPLX, DIMAG, SIGN
135: * ..
136: * .. External Subroutines ..
137: EXTERNAL ZDSCAL, ZSCAL
138: * ..
139: * .. Executable Statements ..
140: *
141: IF( N.LE.0 ) THEN
142: TAU = ZERO
143: RETURN
144: END IF
145: *
146: XNORM = DZNRM2( N-1, X, INCX )
147: ALPHR = DBLE( ALPHA )
148: ALPHI = DIMAG( ALPHA )
149: *
150: IF( XNORM.EQ.ZERO ) THEN
151: *
152: * H = [1-alpha/abs(alpha) 0; 0 I], sign chosen so ALPHA >= 0.
153: *
154: IF( ALPHI.EQ.ZERO ) THEN
155: IF( ALPHR.GE.ZERO ) THEN
156: * When TAU.eq.ZERO, the vector is special-cased to be
157: * all zeros in the application routines. We do not need
158: * to clear it.
159: TAU = ZERO
160: ELSE
161: * However, the application routines rely on explicit
162: * zero checks when TAU.ne.ZERO, and we must clear X.
163: TAU = TWO
164: DO J = 1, N-1
165: X( 1 + (J-1)*INCX ) = ZERO
166: END DO
167: ALPHA = -ALPHA
168: END IF
169: ELSE
170: * Only "reflecting" the diagonal entry to be real and non-negative.
171: XNORM = DLAPY2( ALPHR, ALPHI )
172: TAU = DCMPLX( ONE - ALPHR / XNORM, -ALPHI / XNORM )
173: DO J = 1, N-1
174: X( 1 + (J-1)*INCX ) = ZERO
175: END DO
176: ALPHA = XNORM
177: END IF
178: ELSE
179: *
180: * general case
181: *
182: BETA = SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
183: SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'E' )
184: BIGNUM = ONE / SMLNUM
185: *
186: KNT = 0
187: IF( ABS( BETA ).LT.SMLNUM ) THEN
188: *
189: * XNORM, BETA may be inaccurate; scale X and recompute them
190: *
191: 10 CONTINUE
192: KNT = KNT + 1
193: CALL ZDSCAL( N-1, BIGNUM, X, INCX )
194: BETA = BETA*BIGNUM
195: ALPHI = ALPHI*BIGNUM
196: ALPHR = ALPHR*BIGNUM
197: IF( (ABS( BETA ).LT.SMLNUM) .AND. (KNT .LT. 20) )
198: $ GO TO 10
199: *
200: * New BETA is at most 1, at least SMLNUM
201: *
202: XNORM = DZNRM2( N-1, X, INCX )
203: ALPHA = DCMPLX( ALPHR, ALPHI )
204: BETA = SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
205: END IF
206: SAVEALPHA = ALPHA
207: ALPHA = ALPHA + BETA
208: IF( BETA.LT.ZERO ) THEN
209: BETA = -BETA
210: TAU = -ALPHA / BETA
211: ELSE
212: ALPHR = ALPHI * (ALPHI/DBLE( ALPHA ))
213: ALPHR = ALPHR + XNORM * (XNORM/DBLE( ALPHA ))
214: TAU = DCMPLX( ALPHR/BETA, -ALPHI/BETA )
215: ALPHA = DCMPLX( -ALPHR, ALPHI )
216: END IF
217: ALPHA = ZLADIV( DCMPLX( ONE ), ALPHA )
218: *
219: IF ( ABS(TAU).LE.SMLNUM ) THEN
220: *
221: * In the case where the computed TAU ends up being a denormalized number,
222: * it loses relative accuracy. This is a BIG problem. Solution: flush TAU
223: * to ZERO (or TWO or whatever makes a nonnegative real number for BETA).
224: *
225: * (Bug report provided by Pat Quillen from MathWorks on Jul 29, 2009.)
226: * (Thanks Pat. Thanks MathWorks.)
227: *
228: ALPHR = DBLE( SAVEALPHA )
229: ALPHI = DIMAG( SAVEALPHA )
230: IF( ALPHI.EQ.ZERO ) THEN
231: IF( ALPHR.GE.ZERO ) THEN
232: TAU = ZERO
233: ELSE
234: TAU = TWO
235: DO J = 1, N-1
236: X( 1 + (J-1)*INCX ) = ZERO
237: END DO
238: BETA = DBLE( -SAVEALPHA )
239: END IF
240: ELSE
241: XNORM = DLAPY2( ALPHR, ALPHI )
242: TAU = DCMPLX( ONE - ALPHR / XNORM, -ALPHI / XNORM )
243: DO J = 1, N-1
244: X( 1 + (J-1)*INCX ) = ZERO
245: END DO
246: BETA = XNORM
247: END IF
248: *
249: ELSE
250: *
251: * This is the general case.
252: *
253: CALL ZSCAL( N-1, ALPHA, X, INCX )
254: *
255: END IF
256: *
257: * If BETA is subnormal, it may lose relative accuracy
258: *
259: DO 20 J = 1, KNT
260: BETA = BETA*SMLNUM
261: 20 CONTINUE
262: ALPHA = BETA
263: END IF
264: *
265: RETURN
266: *
267: * End of ZLARFGP
268: *
269: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlarfgp.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack