|
2329
|
1 SUBROUTINE ZLARFG( 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 COMPLEX*16 ALPHA, TAU |
|
|
10 * .. |
|
|
11 * .. Array Arguments .. |
|
|
12 COMPLEX*16 X( * ) |
|
|
13 * .. |
|
|
14 * |
|
|
15 * Purpose |
|
|
16 * ======= |
|
|
17 * |
|
|
18 * ZLARFG generates a complex 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, with beta real, and x is an |
|
|
25 * (n-1)-element complex vector. H is represented in the form |
|
|
26 * |
|
|
27 * H = I - tau * ( 1 ) * ( 1 v' ) , |
|
|
28 * ( v ) |
|
|
29 * |
|
|
30 * where tau is a complex scalar and v is a complex (n-1)-element |
|
|
31 * vector. Note that H is not hermitian. |
|
|
32 * |
|
|
33 * If the elements of x are all zero and alpha is real, then tau = 0 |
|
|
34 * and H is taken to be the unit matrix. |
|
|
35 * |
|
|
36 * Otherwise 1 <= real(tau) <= 2 and abs(tau-1) <= 1 . |
|
|
37 * |
|
|
38 * Arguments |
|
|
39 * ========= |
|
|
40 * |
|
|
41 * N (input) INTEGER |
|
|
42 * The order of the elementary reflector. |
|
|
43 * |
|
|
44 * ALPHA (input/output) COMPLEX*16 |
|
|
45 * On entry, the value alpha. |
|
|
46 * On exit, it is overwritten with the value beta. |
|
|
47 * |
|
|
48 * X (input/output) COMPLEX*16 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) COMPLEX*16 |
|
|
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 ALPHI, ALPHR, BETA, RSAFMN, SAFMIN, XNORM |
|
|
68 * .. |
|
|
69 * .. External Functions .. |
|
|
70 DOUBLE PRECISION DLAMCH, DLAPY3, DZNRM2 |
|
|
71 COMPLEX*16 ZLADIV |
|
|
72 EXTERNAL DLAMCH, DLAPY3, 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 .AND. ALPHI.EQ.ZERO ) THEN |
|
|
92 * |
|
|
93 * H = I |
|
|
94 * |
|
|
95 TAU = ZERO |
|
|
96 ELSE |
|
|
97 * |
|
|
98 * general case |
|
|
99 * |
|
|
100 BETA = -SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) |
|
|
101 SAFMIN = DLAMCH( 'S' ) / DLAMCH( 'E' ) |
|
|
102 RSAFMN = ONE / SAFMIN |
|
|
103 * |
|
|
104 IF( ABS( BETA ).LT.SAFMIN ) THEN |
|
|
105 * |
|
|
106 * XNORM, BETA may be inaccurate; scale X and recompute them |
|
|
107 * |
|
|
108 KNT = 0 |
|
|
109 10 CONTINUE |
|
|
110 KNT = KNT + 1 |
|
|
111 CALL ZDSCAL( N-1, RSAFMN, X, INCX ) |
|
|
112 BETA = BETA*RSAFMN |
|
|
113 ALPHI = ALPHI*RSAFMN |
|
|
114 ALPHR = ALPHR*RSAFMN |
|
|
115 IF( ABS( BETA ).LT.SAFMIN ) |
|
|
116 $ GO TO 10 |
|
|
117 * |
|
|
118 * New BETA is at most 1, at least SAFMIN |
|
|
119 * |
|
|
120 XNORM = DZNRM2( N-1, X, INCX ) |
|
|
121 ALPHA = DCMPLX( ALPHR, ALPHI ) |
|
|
122 BETA = -SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) |
|
|
123 TAU = DCMPLX( ( BETA-ALPHR ) / BETA, -ALPHI / BETA ) |
|
|
124 ALPHA = ZLADIV( DCMPLX( ONE ), ALPHA-BETA ) |
|
|
125 CALL ZSCAL( N-1, ALPHA, X, INCX ) |
|
|
126 * |
|
|
127 * If ALPHA is subnormal, it may lose relative accuracy |
|
|
128 * |
|
|
129 ALPHA = BETA |
|
|
130 DO 20 J = 1, KNT |
|
|
131 ALPHA = ALPHA*SAFMIN |
|
|
132 20 CONTINUE |
|
|
133 ELSE |
|
|
134 TAU = DCMPLX( ( BETA-ALPHR ) / BETA, -ALPHI / BETA ) |
|
|
135 ALPHA = ZLADIV( DCMPLX( ONE ), ALPHA-BETA ) |
|
|
136 CALL ZSCAL( N-1, ALPHA, X, INCX ) |
|
|
137 ALPHA = BETA |
|
|
138 END IF |
|
|
139 END IF |
|
|
140 * |
|
|
141 RETURN |
|
|
142 * |
|
|
143 * End of ZLARFG |
|
|
144 * |
|
|
145 END |
