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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dladiv.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack