LAPACK  3.8.0
LAPACK: Linear Algebra PACKage
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 *> \date January 2013
88 *
89 *> \ingroup doubleOTHERauxiliary
90 *
91 * =====================================================================
92  SUBROUTINE dladiv( A, B, C, D, P, Q )
93 *
94 * -- LAPACK auxiliary routine (version 3.7.0) --
95 * -- LAPACK is a software package provided by Univ. of Tennessee, --
96 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
97 * January 2013
98 *
99 * .. Scalar Arguments ..
100  DOUBLE PRECISION A, B, C, D, P, Q
101 * ..
102 *
103 * =====================================================================
104 *
105 * .. Parameters ..
106  DOUBLE PRECISION BS
107  parameter( bs = 2.0d0 )
108  DOUBLE PRECISION HALF
109  parameter( half = 0.5d0 )
110  DOUBLE PRECISION TWO
111  parameter( two = 2.0d0 )
112 *
113 * .. Local Scalars ..
114  DOUBLE PRECISION AA, BB, CC, DD, AB, CD, S, OV, UN, BE, EPS
115 * ..
116 * .. External Functions ..
117  DOUBLE PRECISION DLAMCH
118  EXTERNAL dlamch
119 * ..
120 * .. External Subroutines ..
121  EXTERNAL dladiv1
122 * ..
123 * .. Intrinsic Functions ..
124  INTRINSIC abs, max
125 * ..
126 * .. Executable Statements ..
127 *
128  aa = a
129  bb = b
130  cc = c
131  dd = d
132  ab = max( abs(a), abs(b) )
133  cd = max( abs(c), abs(d) )
134  s = 1.0d0
135 
136  ov = dlamch( 'Overflow threshold' )
137  un = dlamch( 'Safe minimum' )
138  eps = dlamch( 'Epsilon' )
139  be = bs / (eps*eps)
140 
141  IF( ab >= half*ov ) THEN
142  aa = half * aa
143  bb = half * bb
144  s = two * s
145  END IF
146  IF( cd >= half*ov ) THEN
147  cc = half * cc
148  dd = half * dd
149  s = half * s
150  END IF
151  IF( ab <= un*bs/eps ) THEN
152  aa = aa * be
153  bb = bb * be
154  s = s / be
155  END IF
156  IF( cd <= un*bs/eps ) THEN
157  cc = cc * be
158  dd = dd * be
159  s = s * be
160  END IF
161  IF( abs( d ).LE.abs( c ) ) THEN
162  CALL dladiv1(aa, bb, cc, dd, p, q)
163  ELSE
164  CALL dladiv1(bb, aa, dd, cc, p, q)
165  q = -q
166  END IF
167  p = p * s
168  q = q * s
169 *
170  RETURN
171 *
172 * End of DLADIV
173 *
174  END
175 
176 *> \ingroup doubleOTHERauxiliary
177 
178 
179  SUBROUTINE dladiv1( A, B, C, D, P, Q )
180 *
181 * -- LAPACK auxiliary routine (version 3.7.0) --
182 * -- LAPACK is a software package provided by Univ. of Tennessee, --
183 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
184 * January 2013
185 *
186 * .. Scalar Arguments ..
187  DOUBLE PRECISION A, B, C, D, P, Q
188 * ..
189 *
190 * =====================================================================
191 *
192 * .. Parameters ..
193  DOUBLE PRECISION ONE
194  parameter( one = 1.0d0 )
195 *
196 * .. Local Scalars ..
197  DOUBLE PRECISION R, T
198 * ..
199 * .. External Functions ..
200  DOUBLE PRECISION DLADIV2
201  EXTERNAL dladiv2
202 * ..
203 * .. Executable Statements ..
204 *
205  r = d / c
206  t = one / (c + d * r)
207  p = dladiv2(a, b, c, d, r, t)
208  a = -a
209  q = dladiv2(b, a, c, d, r, t)
210 *
211  RETURN
212 *
213 * End of DLADIV1
214 *
215  END
216 
217 *> \ingroup doubleOTHERauxiliary
218 
219  DOUBLE PRECISION FUNCTION dladiv2( A, B, C, D, R, T )
220 *
221 * -- LAPACK auxiliary routine (version 3.7.0) --
222 * -- LAPACK is a software package provided by Univ. of Tennessee, --
223 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
224 * January 2013
225 *
226 * .. Scalar Arguments ..
227  DOUBLE PRECISION A, B, C, D, R, T
228 * ..
229 *
230 * =====================================================================
231 *
232 * .. Parameters ..
233  DOUBLE PRECISION ZERO
234  parameter( zero = 0.0d0 )
235 *
236 * .. Local Scalars ..
237  DOUBLE PRECISION BR
238 * ..
239 * .. Executable Statements ..
240 *
241  IF( r.NE.zero ) THEN
242  br = b * r
243  IF( br.NE.zero ) THEN
244  dladiv2 = (a + br) * t
245  ELSE
246  dladiv2 = a * t + (b * t) * r
247  END IF
248  ELSE
249  dladiv2 = (a + d * (b / c)) * t
250  END IF
251 *
252  RETURN
253 *
254 * End of DLADIV12
255 *
256  END
subroutine dladiv1(A, B, C, D, P, Q)
Definition: dladiv.f:180
subroutine dladiv(A, B, C, D, P, Q)
DLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
Definition: dladiv.f:93
double precision function dladiv2(A, B, C, D, R, T)
Definition: dladiv.f:220