92 SUBROUTINE dladiv( A, B, C, D, P, Q )
100 DOUBLE PRECISION A, B, C, D, P, Q
107 parameter ( bs = 2.0d0 )
108 DOUBLE PRECISION HALF
109 parameter ( half = 0.5d0 )
111 parameter ( two = 2.0d0 )
114 DOUBLE PRECISION AA, BB, CC, DD, AB, CD, S, OV, UN, BE, EPS
117 DOUBLE PRECISION DLAMCH
132 ab = max( abs(a), abs(b) )
133 cd = max( abs(c), abs(d) )
136 ov = dlamch(
'Overflow threshold' )
137 un = dlamch(
'Safe minimum' )
138 eps = dlamch(
'Epsilon' )
141 IF( ab >= half*ov )
THEN
146 IF( cd >= half*ov )
THEN
151 IF( ab <= un*bs/eps )
THEN
156 IF( cd <= un*bs/eps )
THEN
161 IF( abs( d ).LE.abs( c ) )
THEN
162 CALL dladiv1(aa, bb, cc, dd, p, q)
164 CALL dladiv1(bb, aa, dd, cc, p, q)
179 SUBROUTINE dladiv1( A, B, C, D, P, Q )
187 DOUBLE PRECISION A, B, C, D, P, Q
194 parameter ( one = 1.0d0 )
197 DOUBLE PRECISION R, T
200 DOUBLE PRECISION DLADIV2
206 t = one / (c + d * r)
207 p = dladiv2(a, b, c, d, r, t)
209 q = dladiv2(b, a, c, d, r, t)
219 DOUBLE PRECISION FUNCTION dladiv2( A, B, C, D, R, T )
227 DOUBLE PRECISION A, B, C, D, R, T
233 DOUBLE PRECISION ZERO
234 parameter ( zero = 0.0d0 )
243 IF( br.NE.zero )
THEN
249 dladiv2 = (a + d * (b / c)) * t
subroutine dladiv1(A, B, C, D, P, Q)
subroutine dladiv(A, B, C, D, P, Q)
DLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
double precision function dladiv2(A, B, C, D, R, T)