![[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|
Return to zlarfgp.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: SUBROUTINE ZLARFGP( N, ALPHA, X, INCX, TAU )
2: *
3: * -- LAPACK auxiliary routine (version 3.2.2) --
4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6: * June 2010
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' * ( alpha ) = ( beta ), 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' ) ,
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