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 auxOTHERauxiliary
90: *
91: * =====================================================================
92: SUBROUTINE DLADIV( A, B, C, D, P, Q )
93: *
94: * -- LAPACK auxiliary routine (version 3.5.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:
177:
178: SUBROUTINE DLADIV1( A, B, C, D, P, Q )
179: *
180: * -- LAPACK auxiliary routine (version 3.5.0) --
181: * -- LAPACK is a software package provided by Univ. of Tennessee, --
182: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
183: * January 2013
184: *
185: * .. Scalar Arguments ..
186: DOUBLE PRECISION A, B, C, D, P, Q
187: * ..
188: *
189: * =====================================================================
190: *
191: * .. Parameters ..
192: DOUBLE PRECISION ONE
193: PARAMETER ( ONE = 1.0D0 )
194: *
195: * .. Local Scalars ..
196: DOUBLE PRECISION R, T
197: * ..
198: * .. External Functions ..
199: DOUBLE PRECISION DLADIV2
200: EXTERNAL DLADIV2
201: * ..
202: * .. Executable Statements ..
203: *
204: R = D / C
205: T = ONE / (C + D * R)
206: P = DLADIV2(A, B, C, D, R, T)
207: A = -A
208: Q = DLADIV2(B, A, C, D, R, T)
209: *
210: RETURN
211: *
212: * End of DLADIV1
213: *
214: END
215:
216: DOUBLE PRECISION FUNCTION DLADIV2( A, B, C, D, R, T )
217: *
218: * -- LAPACK auxiliary routine (version 3.5.0) --
219: * -- LAPACK is a software package provided by Univ. of Tennessee, --
220: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
221: * January 2013
222: *
223: * .. Scalar Arguments ..
224: DOUBLE PRECISION A, B, C, D, R, T
225: * ..
226: *
227: * =====================================================================
228: *
229: * .. Parameters ..
230: DOUBLE PRECISION ZERO
231: PARAMETER ( ZERO = 0.0D0 )
232: *
233: * .. Local Scalars ..
234: DOUBLE PRECISION BR
235: * ..
236: * .. Executable Statements ..
237: *
238: IF( R.NE.ZERO ) THEN
239: BR = B * R
240: if( BR.NE.ZERO ) THEN
241: DLADIV2 = (A + BR) * T
242: ELSE
243: DLADIV2 = A * T + (B * T) * R
244: END IF
245: ELSE
246: DLADIV2 = (A + D * (B / C)) * T
247: END IF
248: *
249: RETURN
250: *
251: * End of DLADIV12
252: *
253: 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