LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dladiv.f
Go to the documentation of this file.
1*> \brief \b DLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DLADIV + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dladiv.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dladiv.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dladiv.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DLADIV( A, B, C, D, P, Q )
22*
23* .. Scalar Arguments ..
24* DOUBLE PRECISION A, B, C, D, P, Q
25* ..
26*
27*
28*> \par Purpose:
29* =============
30*>
31*> \verbatim
32*>
33*> DLADIV performs complex division in real arithmetic
34*>
35*> a + i*b
36*> p + i*q = ---------
37*> c + i*d
38*>
39*> The algorithm is due to Michael Baudin and Robert L. Smith
40*> and can be found in the paper
41*> "A Robust Complex Division in Scilab"
42*> \endverbatim
43*
44* Arguments:
45* ==========
46*
47*> \param[in] A
48*> \verbatim
49*> A is DOUBLE PRECISION
50*> \endverbatim
51*>
52*> \param[in] B
53*> \verbatim
54*> B is DOUBLE PRECISION
55*> \endverbatim
56*>
57*> \param[in] C
58*> \verbatim
59*> C is DOUBLE PRECISION
60*> \endverbatim
61*>
62*> \param[in] D
63*> \verbatim
64*> D is DOUBLE PRECISION
65*> The scalars a, b, c, and d in the above expression.
66*> \endverbatim
67*>
68*> \param[out] P
69*> \verbatim
70*> P is DOUBLE PRECISION
71*> \endverbatim
72*>
73*> \param[out] Q
74*> \verbatim
75*> Q is DOUBLE PRECISION
76*> The scalars p and q in the above expression.
77*> \endverbatim
78*
79* Authors:
80* ========
81*
82*> \author Univ. of Tennessee
83*> \author Univ. of California Berkeley
84*> \author Univ. of Colorado Denver
85*> \author NAG Ltd.
86*
87*> \ingroup doubleOTHERauxiliary
88*
89* =====================================================================
90 SUBROUTINE dladiv( A, B, C, D, P, Q )
91*
92* -- LAPACK auxiliary routine --
93* -- LAPACK is a software package provided by Univ. of Tennessee, --
94* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
95*
96* .. Scalar Arguments ..
97 DOUBLE PRECISION A, B, C, D, P, Q
98* ..
99*
100* =====================================================================
101*
102* .. Parameters ..
103 DOUBLE PRECISION BS
104 parameter( bs = 2.0d0 )
105 DOUBLE PRECISION HALF
106 parameter( half = 0.5d0 )
107 DOUBLE PRECISION TWO
108 parameter( two = 2.0d0 )
109*
110* .. Local Scalars ..
111 DOUBLE PRECISION AA, BB, CC, DD, AB, CD, S, OV, UN, BE, EPS
112* ..
113* .. External Functions ..
114 DOUBLE PRECISION DLAMCH
115 EXTERNAL dlamch
116* ..
117* .. External Subroutines ..
118 EXTERNAL dladiv1
119* ..
120* .. Intrinsic Functions ..
121 INTRINSIC abs, max
122* ..
123* .. Executable Statements ..
124*
125 aa = a
126 bb = b
127 cc = c
128 dd = d
129 ab = max( abs(a), abs(b) )
130 cd = max( abs(c), abs(d) )
131 s = 1.0d0
132
133 ov = dlamch( 'Overflow threshold' )
134 un = dlamch( 'Safe minimum' )
135 eps = dlamch( 'Epsilon' )
136 be = bs / (eps*eps)
137
138 IF( ab >= half*ov ) THEN
139 aa = half * aa
140 bb = half * bb
141 s = two * s
142 END IF
143 IF( cd >= half*ov ) THEN
144 cc = half * cc
145 dd = half * dd
146 s = half * s
147 END IF
148 IF( ab <= un*bs/eps ) THEN
149 aa = aa * be
150 bb = bb * be
151 s = s / be
152 END IF
153 IF( cd <= un*bs/eps ) THEN
154 cc = cc * be
155 dd = dd * be
156 s = s * be
157 END IF
158 IF( abs( d ).LE.abs( c ) ) THEN
159 CALL dladiv1(aa, bb, cc, dd, p, q)
160 ELSE
161 CALL dladiv1(bb, aa, dd, cc, p, q)
162 q = -q
163 END IF
164 p = p * s
165 q = q * s
166*
167 RETURN
168*
169* End of DLADIV
170*
171 END
172
173*> \ingroup doubleOTHERauxiliary
174
175
176 SUBROUTINE dladiv1( A, B, C, D, P, Q )
177*
178* -- LAPACK auxiliary routine --
179* -- LAPACK is a software package provided by Univ. of Tennessee, --
180* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
181*
182* .. Scalar Arguments ..
183 DOUBLE PRECISION A, B, C, D, P, Q
184* ..
185*
186* =====================================================================
187*
188* .. Parameters ..
189 DOUBLE PRECISION ONE
190 parameter( one = 1.0d0 )
191*
192* .. Local Scalars ..
193 DOUBLE PRECISION R, T
194* ..
195* .. External Functions ..
196 DOUBLE PRECISION DLADIV2
197 EXTERNAL dladiv2
198* ..
199* .. Executable Statements ..
200*
201 r = d / c
202 t = one / (c + d * r)
203 p = dladiv2(a, b, c, d, r, t)
204 a = -a
205 q = dladiv2(b, a, c, d, r, t)
206*
207 RETURN
208*
209* End of DLADIV1
210*
211 END
212
213*> \ingroup doubleOTHERauxiliary
214
215 DOUBLE PRECISION FUNCTION dladiv2( A, B, C, D, R, T )
216*
217* -- LAPACK auxiliary routine --
218* -- LAPACK is a software package provided by Univ. of Tennessee, --
219* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
220*
221* .. Scalar Arguments ..
222 DOUBLE PRECISION a, b, c, d, r, t
223* ..
224*
225* =====================================================================
226*
227* .. Parameters ..
228 DOUBLE PRECISION zero
229 parameter( zero = 0.0d0 )
230*
231* .. Local Scalars ..
232 DOUBLE PRECISION br
233* ..
234* .. Executable Statements ..
235*
236 IF( r.NE.zero ) THEN
237 br = b * r
238 IF( br.NE.zero ) THEN
239 dladiv2 = (a + br) * t
240 ELSE
241 dladiv2 = a * t + (b * t) * r
242 END IF
243 ELSE
244 dladiv2 = (a + d * (b / c)) * t
245 END IF
246*
247 RETURN
248*
249* End of DLADIV2
250*
251 END
double precision function dladiv2(A, B, C, D, R, T)
Definition: dladiv.f:216
subroutine dladiv1(A, B, C, D, P, Q)
Definition: dladiv.f:177
subroutine dladiv(A, B, C, D, P, Q)
DLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
Definition: dladiv.f:91