1 SUBROUTINE ZLARTG( F, G, CS, SN, R )
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 DOUBLE PRECISION CS
10 COMPLEX*16 F, G, R, SN
11 * ..
12 *
13 * Purpose
14 * =======
15 *
16 * ZLARTG generates a plane rotation so that
17 *
18 * [ CS SN ] [ F ] [ R ]
19 * [ __ ] . [ ] = [ ] where CS**2 + |SN|**2 = 1.
20 * [ -SN CS ] [ G ] [ 0 ]
21 *
22 * This is a faster version of the BLAS1 routine ZROTG, except for
23 * the following differences:
24 * F and G are unchanged on return.
25 * If G=0, then CS=1 and SN=0.
26 * If F=0, then CS=0 and SN is chosen so that R is real.
27 *
28 * Arguments
29 * =========
30 *
31 * F (input) COMPLEX*16
32 * The first component of vector to be rotated.
33 *
34 * G (input) COMPLEX*16
35 * The second component of vector to be rotated.
36 *
37 * CS (output) DOUBLE PRECISION
38 * The cosine of the rotation.
39 *
40 * SN (output) COMPLEX*16
41 * The sine of the rotation.
42 *
43 * R (output) COMPLEX*16
44 * The nonzero component of the rotated vector.
45 *
46 * Further Details
47 * ======= =======
48 *
49 * 3-5-96 - Modified with a new algorithm by W. Kahan and J. Demmel
50 *
51 * This version has a few statements commented out for thread safety
52 * (machine parameters are computed on each entry). 10 feb 03, SJH.
53 *
54 * =====================================================================
55 *
56 * .. Parameters ..
57 DOUBLE PRECISION TWO, ONE, ZERO
58 PARAMETER ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 )
59 COMPLEX*16 CZERO
60 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) )
61 * ..
62 * .. Local Scalars ..
63 * LOGICAL FIRST
64 INTEGER COUNT, I
65 DOUBLE PRECISION D, DI, DR, EPS, F2, F2S, G2, G2S, SAFMIN,
66 $ SAFMN2, SAFMX2, SCALE
67 COMPLEX*16 FF, FS, GS
68 * ..
69 * .. External Functions ..
70 DOUBLE PRECISION DLAMCH, DLAPY2
71 EXTERNAL DLAMCH, DLAPY2
72 * ..
73 * .. Intrinsic Functions ..
74 INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, INT, LOG,
75 $ MAX, SQRT
76 * ..
77 * .. Statement Functions ..
78 DOUBLE PRECISION ABS1, ABSSQ
79 * ..
80 * .. Save statement ..
81 * SAVE FIRST, SAFMX2, SAFMIN, SAFMN2
82 * ..
83 * .. Data statements ..
84 * DATA FIRST / .TRUE. /
85 * ..
86 * .. Statement Function definitions ..
87 ABS1( FF ) = MAX( ABS( DBLE( FF ) ), ABS( DIMAG( FF ) ) )
88 ABSSQ( FF ) = DBLE( FF )**2 + DIMAG( FF )**2
89 * ..
90 * .. Executable Statements ..
91 *
92 * IF( FIRST ) THEN
93 SAFMIN = DLAMCH( 'S' )
94 EPS = DLAMCH( 'E' )
95 SAFMN2 = DLAMCH( 'B' )**INT( LOG( SAFMIN / EPS ) /
96 $ LOG( DLAMCH( 'B' ) ) / TWO )
97 SAFMX2 = ONE / SAFMN2
98 * FIRST = .FALSE.
99 * END IF
100 SCALE = MAX( ABS1( F ), ABS1( G ) )
101 FS = F
102 GS = G
103 COUNT = 0
104 IF( SCALE.GE.SAFMX2 ) THEN
105 10 CONTINUE
106 COUNT = COUNT + 1
107 FS = FS*SAFMN2
108 GS = GS*SAFMN2
109 SCALE = SCALE*SAFMN2
110 IF( SCALE.GE.SAFMX2 )
111 $ GO TO 10
112 ELSE IF( SCALE.LE.SAFMN2 ) THEN
113 IF( G.EQ.CZERO ) THEN
114 CS = ONE
115 SN = CZERO
116 R = F
117 RETURN
118 END IF
119 20 CONTINUE
120 COUNT = COUNT - 1
121 FS = FS*SAFMX2
122 GS = GS*SAFMX2
123 SCALE = SCALE*SAFMX2
124 IF( SCALE.LE.SAFMN2 )
125 $ GO TO 20
126 END IF
127 F2 = ABSSQ( FS )
128 G2 = ABSSQ( GS )
129 IF( F2.LE.MAX( G2, ONE )*SAFMIN ) THEN
130 *
131 * This is a rare case: F is very small.
132 *
133 IF( F.EQ.CZERO ) THEN
134 CS = ZERO
135 R = DLAPY2( DBLE( G ), DIMAG( G ) )
136 * Do complex/real division explicitly with two real divisions
137 D = DLAPY2( DBLE( GS ), DIMAG( GS ) )
138 SN = DCMPLX( DBLE( GS ) / D, -DIMAG( GS ) / D )
139 RETURN
140 END IF
141 F2S = DLAPY2( DBLE( FS ), DIMAG( FS ) )
142 * G2 and G2S are accurate
143 * G2 is at least SAFMIN, and G2S is at least SAFMN2
144 G2S = SQRT( G2 )
145 * Error in CS from underflow in F2S is at most
146 * UNFL / SAFMN2 .lt. sqrt(UNFL*EPS) .lt. EPS
147 * If MAX(G2,ONE)=G2, then F2 .lt. G2*SAFMIN,
148 * and so CS .lt. sqrt(SAFMIN)
149 * If MAX(G2,ONE)=ONE, then F2 .lt. SAFMIN
150 * and so CS .lt. sqrt(SAFMIN)/SAFMN2 = sqrt(EPS)
151 * Therefore, CS = F2S/G2S / sqrt( 1 + (F2S/G2S)**2 ) = F2S/G2S
152 CS = F2S / G2S
153 * Make sure abs(FF) = 1
154 * Do complex/real division explicitly with 2 real divisions
155 IF( ABS1( F ).GT.ONE ) THEN
156 D = DLAPY2( DBLE( F ), DIMAG( F ) )
157 FF = DCMPLX( DBLE( F ) / D, DIMAG( F ) / D )
158 ELSE
159 DR = SAFMX2*DBLE( F )
160 DI = SAFMX2*DIMAG( F )
161 D = DLAPY2( DR, DI )
162 FF = DCMPLX( DR / D, DI / D )
163 END IF
164 SN = FF*DCMPLX( DBLE( GS ) / G2S, -DIMAG( GS ) / G2S )
165 R = CS*F + SN*G
166 ELSE
167 *
168 * This is the most common case.
169 * Neither F2 nor F2/G2 are less than SAFMIN
170 * F2S cannot overflow, and it is accurate
171 *
172 F2S = SQRT( ONE+G2 / F2 )
173 * Do the F2S(real)*FS(complex) multiply with two real multiplies
174 R = DCMPLX( F2S*DBLE( FS ), F2S*DIMAG( FS ) )
175 CS = ONE / F2S
176 D = F2 + G2
177 * Do complex/real division explicitly with two real divisions
178 SN = DCMPLX( DBLE( R ) / D, DIMAG( R ) / D )
179 SN = SN*DCONJG( GS )
180 IF( COUNT.NE.0 ) THEN
181 IF( COUNT.GT.0 ) THEN
182 DO 30 I = 1, COUNT
183 R = R*SAFMX2
184 30 CONTINUE
185 ELSE
186 DO 40 I = 1, -COUNT
187 R = R*SAFMN2
188 40 CONTINUE
189 END IF
190 END IF
191 END IF
192 RETURN
193 *
194 * End of ZLARTG
195 *
196 END
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 DOUBLE PRECISION CS
10 COMPLEX*16 F, G, R, SN
11 * ..
12 *
13 * Purpose
14 * =======
15 *
16 * ZLARTG generates a plane rotation so that
17 *
18 * [ CS SN ] [ F ] [ R ]
19 * [ __ ] . [ ] = [ ] where CS**2 + |SN|**2 = 1.
20 * [ -SN CS ] [ G ] [ 0 ]
21 *
22 * This is a faster version of the BLAS1 routine ZROTG, except for
23 * the following differences:
24 * F and G are unchanged on return.
25 * If G=0, then CS=1 and SN=0.
26 * If F=0, then CS=0 and SN is chosen so that R is real.
27 *
28 * Arguments
29 * =========
30 *
31 * F (input) COMPLEX*16
32 * The first component of vector to be rotated.
33 *
34 * G (input) COMPLEX*16
35 * The second component of vector to be rotated.
36 *
37 * CS (output) DOUBLE PRECISION
38 * The cosine of the rotation.
39 *
40 * SN (output) COMPLEX*16
41 * The sine of the rotation.
42 *
43 * R (output) COMPLEX*16
44 * The nonzero component of the rotated vector.
45 *
46 * Further Details
47 * ======= =======
48 *
49 * 3-5-96 - Modified with a new algorithm by W. Kahan and J. Demmel
50 *
51 * This version has a few statements commented out for thread safety
52 * (machine parameters are computed on each entry). 10 feb 03, SJH.
53 *
54 * =====================================================================
55 *
56 * .. Parameters ..
57 DOUBLE PRECISION TWO, ONE, ZERO
58 PARAMETER ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 )
59 COMPLEX*16 CZERO
60 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) )
61 * ..
62 * .. Local Scalars ..
63 * LOGICAL FIRST
64 INTEGER COUNT, I
65 DOUBLE PRECISION D, DI, DR, EPS, F2, F2S, G2, G2S, SAFMIN,
66 $ SAFMN2, SAFMX2, SCALE
67 COMPLEX*16 FF, FS, GS
68 * ..
69 * .. External Functions ..
70 DOUBLE PRECISION DLAMCH, DLAPY2
71 EXTERNAL DLAMCH, DLAPY2
72 * ..
73 * .. Intrinsic Functions ..
74 INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, INT, LOG,
75 $ MAX, SQRT
76 * ..
77 * .. Statement Functions ..
78 DOUBLE PRECISION ABS1, ABSSQ
79 * ..
80 * .. Save statement ..
81 * SAVE FIRST, SAFMX2, SAFMIN, SAFMN2
82 * ..
83 * .. Data statements ..
84 * DATA FIRST / .TRUE. /
85 * ..
86 * .. Statement Function definitions ..
87 ABS1( FF ) = MAX( ABS( DBLE( FF ) ), ABS( DIMAG( FF ) ) )
88 ABSSQ( FF ) = DBLE( FF )**2 + DIMAG( FF )**2
89 * ..
90 * .. Executable Statements ..
91 *
92 * IF( FIRST ) THEN
93 SAFMIN = DLAMCH( 'S' )
94 EPS = DLAMCH( 'E' )
95 SAFMN2 = DLAMCH( 'B' )**INT( LOG( SAFMIN / EPS ) /
96 $ LOG( DLAMCH( 'B' ) ) / TWO )
97 SAFMX2 = ONE / SAFMN2
98 * FIRST = .FALSE.
99 * END IF
100 SCALE = MAX( ABS1( F ), ABS1( G ) )
101 FS = F
102 GS = G
103 COUNT = 0
104 IF( SCALE.GE.SAFMX2 ) THEN
105 10 CONTINUE
106 COUNT = COUNT + 1
107 FS = FS*SAFMN2
108 GS = GS*SAFMN2
109 SCALE = SCALE*SAFMN2
110 IF( SCALE.GE.SAFMX2 )
111 $ GO TO 10
112 ELSE IF( SCALE.LE.SAFMN2 ) THEN
113 IF( G.EQ.CZERO ) THEN
114 CS = ONE
115 SN = CZERO
116 R = F
117 RETURN
118 END IF
119 20 CONTINUE
120 COUNT = COUNT - 1
121 FS = FS*SAFMX2
122 GS = GS*SAFMX2
123 SCALE = SCALE*SAFMX2
124 IF( SCALE.LE.SAFMN2 )
125 $ GO TO 20
126 END IF
127 F2 = ABSSQ( FS )
128 G2 = ABSSQ( GS )
129 IF( F2.LE.MAX( G2, ONE )*SAFMIN ) THEN
130 *
131 * This is a rare case: F is very small.
132 *
133 IF( F.EQ.CZERO ) THEN
134 CS = ZERO
135 R = DLAPY2( DBLE( G ), DIMAG( G ) )
136 * Do complex/real division explicitly with two real divisions
137 D = DLAPY2( DBLE( GS ), DIMAG( GS ) )
138 SN = DCMPLX( DBLE( GS ) / D, -DIMAG( GS ) / D )
139 RETURN
140 END IF
141 F2S = DLAPY2( DBLE( FS ), DIMAG( FS ) )
142 * G2 and G2S are accurate
143 * G2 is at least SAFMIN, and G2S is at least SAFMN2
144 G2S = SQRT( G2 )
145 * Error in CS from underflow in F2S is at most
146 * UNFL / SAFMN2 .lt. sqrt(UNFL*EPS) .lt. EPS
147 * If MAX(G2,ONE)=G2, then F2 .lt. G2*SAFMIN,
148 * and so CS .lt. sqrt(SAFMIN)
149 * If MAX(G2,ONE)=ONE, then F2 .lt. SAFMIN
150 * and so CS .lt. sqrt(SAFMIN)/SAFMN2 = sqrt(EPS)
151 * Therefore, CS = F2S/G2S / sqrt( 1 + (F2S/G2S)**2 ) = F2S/G2S
152 CS = F2S / G2S
153 * Make sure abs(FF) = 1
154 * Do complex/real division explicitly with 2 real divisions
155 IF( ABS1( F ).GT.ONE ) THEN
156 D = DLAPY2( DBLE( F ), DIMAG( F ) )
157 FF = DCMPLX( DBLE( F ) / D, DIMAG( F ) / D )
158 ELSE
159 DR = SAFMX2*DBLE( F )
160 DI = SAFMX2*DIMAG( F )
161 D = DLAPY2( DR, DI )
162 FF = DCMPLX( DR / D, DI / D )
163 END IF
164 SN = FF*DCMPLX( DBLE( GS ) / G2S, -DIMAG( GS ) / G2S )
165 R = CS*F + SN*G
166 ELSE
167 *
168 * This is the most common case.
169 * Neither F2 nor F2/G2 are less than SAFMIN
170 * F2S cannot overflow, and it is accurate
171 *
172 F2S = SQRT( ONE+G2 / F2 )
173 * Do the F2S(real)*FS(complex) multiply with two real multiplies
174 R = DCMPLX( F2S*DBLE( FS ), F2S*DIMAG( FS ) )
175 CS = ONE / F2S
176 D = F2 + G2
177 * Do complex/real division explicitly with two real divisions
178 SN = DCMPLX( DBLE( R ) / D, DIMAG( R ) / D )
179 SN = SN*DCONJG( GS )
180 IF( COUNT.NE.0 ) THEN
181 IF( COUNT.GT.0 ) THEN
182 DO 30 I = 1, COUNT
183 R = R*SAFMX2
184 30 CONTINUE
185 ELSE
186 DO 40 I = 1, -COUNT
187 R = R*SAFMN2
188 40 CONTINUE
189 END IF
190 END IF
191 END IF
192 RETURN
193 *
194 * End of ZLARTG
195 *
196 END