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)
186 DOUBLE PRECISION A, B, C, D, P, Q
193 parameter( one = 1.0d0 )
196 DOUBLE PRECISION R, T
199 DOUBLE PRECISION DLADIV2
205 t = one / (c + d * r)
206 p = dladiv2(a, b, c, d, r, t)
208 q = dladiv2(b, a, c, d, r, t)
216 DOUBLE PRECISION FUNCTION dladiv2( A, B, C, D, R, T )
224 DOUBLE PRECISION A, B, C, D, R, T
230 DOUBLE PRECISION ZERO
231 parameter( zero = 0.0d0 )
240 if( br.NE.zero )
THEN
241 dladiv2 = (a + br) * t
243 dladiv2 = a * t + (b * t) * r
246 dladiv2 = (a + d * (b / c)) * t
subroutine dladiv(A, B, C, D, P, Q)
DLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
subroutine dladiv1(A, B, C, D, P, Q)