1
SUBROUTINE ZLARTG( F, G, CS, SN, R )
3
* -- LAPACK auxiliary routine (version 3.0) --
4
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
5
* Courant Institute, Argonne National Lab, and Rice University
8
* .. Scalar Arguments ..
10
COMPLEX*16 F, G, R, SN
16
* ZLARTG generates a plane rotation so that
18
* [ CS SN ] [ F ] [ R ]
19
* [ __ ] . [ ] = [ ] where CS**2 + |SN|**2 = 1.
20
* [ -SN CS ] [ G ] [ 0 ]
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.
31
* F (input) COMPLEX*16
32
* The first component of vector to be rotated.
34
* G (input) COMPLEX*16
35
* The second component of vector to be rotated.
37
* CS (output) DOUBLE PRECISION
38
* The cosine of the rotation.
40
* SN (output) COMPLEX*16
41
* The sine of the rotation.
43
* R (output) COMPLEX*16
44
* The nonzero component of the rotated vector.
49
* 3-5-96 - Modified with a new algorithm by W. Kahan and J. Demmel
51
* =====================================================================
54
DOUBLE PRECISION TWO, ONE, ZERO
55
PARAMETER ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 )
57
PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) )
62
DOUBLE PRECISION D, DI, DR, EPS, F2, F2S, G2, G2S, SAFMIN,
63
$ SAFMN2, SAFMX2, SCALE
66
* .. External Functions ..
67
DOUBLE PRECISION DLAMCH, DLAPY2
68
EXTERNAL DLAMCH, DLAPY2
70
* .. Intrinsic Functions ..
71
INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, INT, LOG,
74
* .. Statement Functions ..
75
DOUBLE PRECISION ABS1, ABSSQ
77
* .. Save statement ..
78
SAVE FIRST, SAFMX2, SAFMIN, SAFMN2
80
* .. Data statements ..
83
* .. Statement Function definitions ..
84
ABS1( FF ) = MAX( ABS( DBLE( FF ) ), ABS( DIMAG( FF ) ) )
85
ABSSQ( FF ) = DBLE( FF )**2 + DIMAG( FF )**2
87
* .. Executable Statements ..
91
SAFMIN = DLAMCH( 'S' )
93
SAFMN2 = DLAMCH( 'B' )**INT( LOG( SAFMIN / EPS ) /
94
$ LOG( DLAMCH( 'B' ) ) / TWO )
97
SCALE = MAX( ABS1( F ), ABS1( G ) )
101
IF( SCALE.GE.SAFMX2 ) THEN
107
IF( SCALE.GE.SAFMX2 )
109
ELSE IF( SCALE.LE.SAFMN2 ) THEN
110
IF( G.EQ.CZERO ) THEN
121
IF( SCALE.LE.SAFMN2 )
126
IF( F2.LE.MAX( G2, ONE )*SAFMIN ) THEN
128
* This is a rare case: F is very small.
130
IF( F.EQ.CZERO ) THEN
132
R = DLAPY2( DBLE( G ), DIMAG( G ) )
133
* Do complex/real division explicitly with two real divisions
134
D = DLAPY2( DBLE( GS ), DIMAG( GS ) )
135
SN = DCMPLX( DBLE( GS ) / D, -DIMAG( GS ) / D )
138
F2S = DLAPY2( DBLE( FS ), DIMAG( FS ) )
139
* G2 and G2S are accurate
140
* G2 is at least SAFMIN, and G2S is at least SAFMN2
142
* Error in CS from underflow in F2S is at most
143
* UNFL / SAFMN2 .lt. sqrt(UNFL*EPS) .lt. EPS
144
* If MAX(G2,ONE)=G2, then F2 .lt. G2*SAFMIN,
145
* and so CS .lt. sqrt(SAFMIN)
146
* If MAX(G2,ONE)=ONE, then F2 .lt. SAFMIN
147
* and so CS .lt. sqrt(SAFMIN)/SAFMN2 = sqrt(EPS)
148
* Therefore, CS = F2S/G2S / sqrt( 1 + (F2S/G2S)**2 ) = F2S/G2S
150
* Make sure abs(FF) = 1
151
* Do complex/real division explicitly with 2 real divisions
152
IF( ABS1( F ).GT.ONE ) THEN
153
D = DLAPY2( DBLE( F ), DIMAG( F ) )
154
FF = DCMPLX( DBLE( F ) / D, DIMAG( F ) / D )
156
DR = SAFMX2*DBLE( F )
157
DI = SAFMX2*DIMAG( F )
159
FF = DCMPLX( DR / D, DI / D )
161
SN = FF*DCMPLX( DBLE( GS ) / G2S, -DIMAG( GS ) / G2S )
165
* This is the most common case.
166
* Neither F2 nor F2/G2 are less than SAFMIN
167
* F2S cannot overflow, and it is accurate
169
F2S = SQRT( ONE+G2 / F2 )
170
* Do the F2S(real)*FS(complex) multiply with two real multiplies
171
R = DCMPLX( F2S*DBLE( FS ), F2S*DIMAG( FS ) )
174
* Do complex/real division explicitly with two real divisions
175
SN = DCMPLX( DBLE( R ) / D, DIMAG( R ) / D )
177
IF( COUNT.NE.0 ) THEN
178
IF( COUNT.GT.0 ) THEN