LAPACK  3.10.1 LAPACK: Linear Algebra PACKage
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
11 *> [TGZ]</a>
13 *> [ZIP]</a>
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 ..
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 *
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 ..
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 *
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 *