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
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