|
2329
|
1 SUBROUTINE DLARFG( N, ALPHA, X, INCX, TAU ) |
|
|
2 * |
|
7034
|
3 * -- LAPACK auxiliary routine (version 3.1) -- |
|
|
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. |
|
|
5 * November 2006 |
|
2329
|
6 * |
|
|
7 * .. Scalar Arguments .. |
|
|
8 INTEGER INCX, N |
|
|
9 DOUBLE PRECISION ALPHA, TAU |
|
|
10 * .. |
|
|
11 * .. Array Arguments .. |
|
|
12 DOUBLE PRECISION X( * ) |
|
|
13 * .. |
|
|
14 * |
|
|
15 * Purpose |
|
|
16 * ======= |
|
|
17 * |
|
|
18 * DLARFG generates a real elementary reflector H of order n, such |
|
|
19 * that |
|
|
20 * |
|
|
21 * H * ( alpha ) = ( beta ), H' * H = I. |
|
|
22 * ( x ) ( 0 ) |
|
|
23 * |
|
|
24 * where alpha and beta are scalars, and x is an (n-1)-element real |
|
|
25 * vector. H is represented in the form |
|
|
26 * |
|
|
27 * H = I - tau * ( 1 ) * ( 1 v' ) , |
|
|
28 * ( v ) |
|
|
29 * |
|
|
30 * where tau is a real scalar and v is a real (n-1)-element |
|
|
31 * vector. |
|
|
32 * |
|
|
33 * If the elements of x are all zero, then tau = 0 and H is taken to be |
|
|
34 * the unit matrix. |
|
|
35 * |
|
|
36 * Otherwise 1 <= tau <= 2. |
|
|
37 * |
|
|
38 * Arguments |
|
|
39 * ========= |
|
|
40 * |
|
|
41 * N (input) INTEGER |
|
|
42 * The order of the elementary reflector. |
|
|
43 * |
|
|
44 * ALPHA (input/output) DOUBLE PRECISION |
|
|
45 * On entry, the value alpha. |
|
|
46 * On exit, it is overwritten with the value beta. |
|
|
47 * |
|
|
48 * X (input/output) DOUBLE PRECISION array, dimension |
|
|
49 * (1+(N-2)*abs(INCX)) |
|
|
50 * On entry, the vector x. |
|
|
51 * On exit, it is overwritten with the vector v. |
|
|
52 * |
|
|
53 * INCX (input) INTEGER |
|
|
54 * The increment between elements of X. INCX > 0. |
|
|
55 * |
|
|
56 * TAU (output) DOUBLE PRECISION |
|
|
57 * The value tau. |
|
|
58 * |
|
|
59 * ===================================================================== |
|
|
60 * |
|
|
61 * .. Parameters .. |
|
|
62 DOUBLE PRECISION ONE, ZERO |
|
|
63 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) |
|
|
64 * .. |
|
|
65 * .. Local Scalars .. |
|
|
66 INTEGER J, KNT |
|
|
67 DOUBLE PRECISION BETA, RSAFMN, SAFMIN, XNORM |
|
|
68 * .. |
|
|
69 * .. External Functions .. |
|
|
70 DOUBLE PRECISION DLAMCH, DLAPY2, DNRM2 |
|
|
71 EXTERNAL DLAMCH, DLAPY2, DNRM2 |
|
|
72 * .. |
|
|
73 * .. Intrinsic Functions .. |
|
|
74 INTRINSIC ABS, SIGN |
|
|
75 * .. |
|
|
76 * .. External Subroutines .. |
|
|
77 EXTERNAL DSCAL |
|
|
78 * .. |
|
|
79 * .. Executable Statements .. |
|
|
80 * |
|
|
81 IF( N.LE.1 ) THEN |
|
|
82 TAU = ZERO |
|
|
83 RETURN |
|
|
84 END IF |
|
|
85 * |
|
|
86 XNORM = DNRM2( N-1, X, INCX ) |
|
|
87 * |
|
|
88 IF( XNORM.EQ.ZERO ) THEN |
|
|
89 * |
|
|
90 * H = I |
|
|
91 * |
|
|
92 TAU = ZERO |
|
|
93 ELSE |
|
|
94 * |
|
|
95 * general case |
|
|
96 * |
|
|
97 BETA = -SIGN( DLAPY2( ALPHA, XNORM ), ALPHA ) |
|
|
98 SAFMIN = DLAMCH( 'S' ) / DLAMCH( 'E' ) |
|
|
99 IF( ABS( BETA ).LT.SAFMIN ) THEN |
|
|
100 * |
|
|
101 * XNORM, BETA may be inaccurate; scale X and recompute them |
|
|
102 * |
|
|
103 RSAFMN = ONE / SAFMIN |
|
|
104 KNT = 0 |
|
|
105 10 CONTINUE |
|
|
106 KNT = KNT + 1 |
|
|
107 CALL DSCAL( N-1, RSAFMN, X, INCX ) |
|
|
108 BETA = BETA*RSAFMN |
|
|
109 ALPHA = ALPHA*RSAFMN |
|
|
110 IF( ABS( BETA ).LT.SAFMIN ) |
|
|
111 $ GO TO 10 |
|
|
112 * |
|
|
113 * New BETA is at most 1, at least SAFMIN |
|
|
114 * |
|
|
115 XNORM = DNRM2( N-1, X, INCX ) |
|
|
116 BETA = -SIGN( DLAPY2( ALPHA, XNORM ), ALPHA ) |
|
|
117 TAU = ( BETA-ALPHA ) / BETA |
|
|
118 CALL DSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX ) |
|
|
119 * |
|
|
120 * If ALPHA is subnormal, it may lose relative accuracy |
|
|
121 * |
|
|
122 ALPHA = BETA |
|
|
123 DO 20 J = 1, KNT |
|
|
124 ALPHA = ALPHA*SAFMIN |
|
|
125 20 CONTINUE |
|
|
126 ELSE |
|
|
127 TAU = ( BETA-ALPHA ) / BETA |
|
|
128 CALL DSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX ) |
|
|
129 ALPHA = BETA |
|
|
130 END IF |
|
|
131 END IF |
|
|
132 * |
|
|
133 RETURN |
|
|
134 * |
|
|
135 * End of DLARFG |
|
|
136 * |
|
|
137 END |
