113 COMPLEX*16 f, g, r, sn
119 DOUBLE PRECISION two, one, zero
120 parameter( two = 2.0d+0, one = 1.0d+0, zero = 0.0d+0 )
122 parameter( czero = ( 0.0d+0, 0.0d+0 ) )
127 DOUBLE PRECISION d, di, dr, eps, f2, f2s, g2, g2s, safmin,
128 $ safmn2, safmx2, scale
129 COMPLEX*16 ff, fs, gs
137 INTRINSIC abs, dble, dcmplx, dconjg, dimag, int, log,
141 DOUBLE PRECISION abs1, abssq
144 abs1( ff ) = max( abs( dble( ff ) ), abs( dimag( ff ) ) )
145 abssq( ff ) = dble( ff )**2 + dimag( ff )**2
151 safmn2 =
dlamch(
'B' )**int( log( safmin / eps ) /
152 $ log(
dlamch(
'B' ) ) / two )
153 safmx2 = one / safmn2
154 scale = max( abs1( f ), abs1( g ) )
158 IF( scale.GE.safmx2 )
THEN
164 IF( scale.GE.safmx2 )
166 ELSE IF( scale.LE.safmn2 )
THEN
167 IF( g.EQ.czero.OR.
disnan( abs( g ) ) )
THEN
178 IF( scale.LE.safmn2 )
183 IF( f2.LE.max( g2, one )*safmin )
THEN
187 IF( f.EQ.czero )
THEN
189 r =
dlapy2( dble( g ), dimag( g ) )
191 d =
dlapy2( dble( gs ), dimag( gs ) )
192 sn = dcmplx( dble( gs ) / d, -dimag( gs ) / d )
195 f2s =
dlapy2( dble( fs ), dimag( fs ) )
209 IF( abs1( f ).GT.one )
THEN
210 d =
dlapy2( dble( f ), dimag( f ) )
211 ff = dcmplx( dble( f ) / d, dimag( f ) / d )
213 dr = safmx2*dble( f )
214 di = safmx2*dimag( f )
216 ff = dcmplx( dr / d, di / d )
218 sn = ff*dcmplx( dble( gs ) / g2s, -dimag( gs ) / g2s )
226 f2s = sqrt( one+g2 / f2 )
228 r = dcmplx( f2s*dble( fs ), f2s*dimag( fs ) )
232 sn = dcmplx( dble( r ) / d, dimag( r ) / d )
234 IF( count.NE.0 )
THEN
235 IF( count.GT.0 )
THEN
subroutine zlartg(F, G, CS, SN, R)
ZLARTG generates a plane rotation with real cosine and complex sine.
double precision function dlapy2(X, Y)
DLAPY2 returns sqrt(x2+y2).
double precision function dlamch(CMACH)
DLAMCH
logical function disnan(DIN)
DISNAN tests input for NaN.