![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlarfgp.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 dlarfgp.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: SUBROUTINE DLARFGP( 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: DOUBLE PRECISION ALPHA, TAU
11: * ..
12: * .. Array Arguments ..
13: DOUBLE PRECISION X( * )
14: * ..
15: *
16: * Purpose
17: * =======
18: *
19: * DLARFGP generates a real 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 non-negative, and x is
26: * an (n-1)-element real vector. H is represented in the form
27: *
28: * H = I - tau * ( 1 ) * ( 1 v' ) ,
29: * ( v )
30: *
31: * where tau is a real scalar and v is a real (n-1)-element
32: * vector.
33: *
34: * If the elements of x are all zero, then tau = 0 and H is taken to be
35: * the unit matrix.
36: *
37: * Arguments
38: * =========
39: *
40: * N (input) INTEGER
41: * The order of the elementary reflector.
42: *
43: * ALPHA (input/output) DOUBLE PRECISION
44: * On entry, the value alpha.
45: * On exit, it is overwritten with the value beta.
46: *
47: * X (input/output) DOUBLE PRECISION 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) DOUBLE PRECISION
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 BETA, BIGNUM, SAVEALPHA, SMLNUM, XNORM
67: * ..
68: * .. External Functions ..
69: DOUBLE PRECISION DLAMCH, DLAPY2, DNRM2
70: EXTERNAL DLAMCH, DLAPY2, DNRM2
71: * ..
72: * .. Intrinsic Functions ..
73: INTRINSIC ABS, SIGN
74: * ..
75: * .. External Subroutines ..
76: EXTERNAL DSCAL
77: * ..
78: * .. Executable Statements ..
79: *
80: IF( N.LE.0 ) THEN
81: TAU = ZERO
82: RETURN
83: END IF
84: *
85: XNORM = DNRM2( N-1, X, INCX )
86: *
87: IF( XNORM.EQ.ZERO ) THEN
88: *
89: * H = [+/-1, 0; I], sign chosen so ALPHA >= 0
90: *
91: IF( ALPHA.GE.ZERO ) THEN
92: * When TAU.eq.ZERO, the vector is special-cased to be
93: * all zeros in the application routines. We do not need
94: * to clear it.
95: TAU = ZERO
96: ELSE
97: * However, the application routines rely on explicit
98: * zero checks when TAU.ne.ZERO, and we must clear X.
99: TAU = TWO
100: DO J = 1, N-1
101: X( 1 + (J-1)*INCX ) = 0
102: END DO
103: ALPHA = -ALPHA
104: END IF
105: ELSE
106: *
107: * general case
108: *
109: BETA = SIGN( DLAPY2( ALPHA, XNORM ), ALPHA )
110: SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'E' )
111: KNT = 0
112: IF( ABS( BETA ).LT.SMLNUM ) THEN
113: *
114: * XNORM, BETA may be inaccurate; scale X and recompute them
115: *
116: BIGNUM = ONE / SMLNUM
117: 10 CONTINUE
118: KNT = KNT + 1
119: CALL DSCAL( N-1, BIGNUM, X, INCX )
120: BETA = BETA*BIGNUM
121: ALPHA = ALPHA*BIGNUM
122: IF( ABS( BETA ).LT.SMLNUM )
123: $ GO TO 10
124: *
125: * New BETA is at most 1, at least SMLNUM
126: *
127: XNORM = DNRM2( N-1, X, INCX )
128: BETA = SIGN( DLAPY2( ALPHA, XNORM ), ALPHA )
129: END IF
130: SAVEALPHA = ALPHA
131: ALPHA = ALPHA + BETA
132: IF( BETA.LT.ZERO ) THEN
133: BETA = -BETA
134: TAU = -ALPHA / BETA
135: ELSE
136: ALPHA = XNORM * (XNORM/ALPHA)
137: TAU = ALPHA / BETA
138: ALPHA = -ALPHA
139: END IF
140: *
141: IF ( ABS(TAU).LE.SMLNUM ) THEN
142: *
143: * In the case where the computed TAU ends up being a denormalized number,
144: * it loses relative accuracy. This is a BIG problem. Solution: flush TAU
145: * to ZERO. This explains the next IF statement.
146: *
147: * (Bug report provided by Pat Quillen from MathWorks on Jul 29, 2009.)
148: * (Thanks Pat. Thanks MathWorks.)
149: *
150: IF( SAVEALPHA.GE.ZERO ) THEN
151: TAU = ZERO
152: ELSE
153: TAU = TWO
154: DO J = 1, N-1
155: X( 1 + (J-1)*INCX ) = 0
156: END DO
157: BETA = -SAVEALPHA
158: END IF
159: *
160: ELSE
161: *
162: * This is the general case.
163: *
164: CALL DSCAL( N-1, ONE / ALPHA, X, INCX )
165: *
166: END IF
167: *
168: * If BETA is subnormal, it may lose relative accuracy
169: *
170: DO 20 J = 1, KNT
171: BETA = BETA*SMLNUM
172: 20 CONTINUE
173: ALPHA = BETA
174: END IF
175: *
176: RETURN
177: *
178: * End of DLARFGP
179: *
180: END