![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlarfg.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 dlarfg.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: SUBROUTINE DLARFG( N, ALPHA, X, INCX, TAU )
2: *
3: * -- LAPACK auxiliary routine (version 3.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: * November 2006
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: * DLARFG 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, and x is an (n-1)-element real
26: * 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: * Otherwise 1 <= tau <= 2.
38: *
39: * Arguments
40: * =========
41: *
42: * N (input) INTEGER
43: * The order of the elementary reflector.
44: *
45: * ALPHA (input/output) DOUBLE PRECISION
46: * On entry, the value alpha.
47: * On exit, it is overwritten with the value beta.
48: *
49: * X (input/output) DOUBLE PRECISION 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) DOUBLE PRECISION
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 BETA, RSAFMN, SAFMIN, XNORM
69: * ..
70: * .. External Functions ..
71: DOUBLE PRECISION DLAMCH, DLAPY2, DNRM2
72: EXTERNAL DLAMCH, DLAPY2, DNRM2
73: * ..
74: * .. Intrinsic Functions ..
75: INTRINSIC ABS, SIGN
76: * ..
77: * .. External Subroutines ..
78: EXTERNAL DSCAL
79: * ..
80: * .. Executable Statements ..
81: *
82: IF( N.LE.1 ) THEN
83: TAU = ZERO
84: RETURN
85: END IF
86: *
87: XNORM = DNRM2( N-1, X, INCX )
88: *
89: IF( XNORM.EQ.ZERO ) THEN
90: *
91: * H = I
92: *
93: TAU = ZERO
94: ELSE
95: *
96: * general case
97: *
98: BETA = -SIGN( DLAPY2( ALPHA, XNORM ), ALPHA )
99: SAFMIN = DLAMCH( 'S' ) / DLAMCH( 'E' )
100: KNT = 0
101: IF( ABS( BETA ).LT.SAFMIN ) THEN
102: *
103: * XNORM, BETA may be inaccurate; scale X and recompute them
104: *
105: RSAFMN = ONE / SAFMIN
106: 10 CONTINUE
107: KNT = KNT + 1
108: CALL DSCAL( N-1, RSAFMN, X, INCX )
109: BETA = BETA*RSAFMN
110: ALPHA = ALPHA*RSAFMN
111: IF( ABS( BETA ).LT.SAFMIN )
112: $ GO TO 10
113: *
114: * New BETA is at most 1, at least SAFMIN
115: *
116: XNORM = DNRM2( N-1, X, INCX )
117: BETA = -SIGN( DLAPY2( ALPHA, XNORM ), ALPHA )
118: END IF
119: TAU = ( BETA-ALPHA ) / BETA
120: CALL DSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX )
121: *
122: * If ALPHA is subnormal, it may lose relative accuracy
123: *
124: DO 20 J = 1, KNT
125: BETA = BETA*SAFMIN
126: 20 CONTINUE
127: ALPHA = BETA
128: END IF
129: *
130: RETURN
131: *
132: * End of DLARFG
133: *
134: END