Télécharger dladiv.eso

Retour à la liste

Numérotation des lignes :

  1. C DLADIV SOURCE BP208322 15/10/13 21:15:25 8670
  2. *> \brief \b DLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
  3. *
  4. * =========== DOCUMENTATION ===========
  5. *
  6. * Online html documentation available at
  7. * http://www.netlib.org/lapack/explore-html/
  8. *
  9. *> \htmlonly
  10. *> Download DLADIV + dependencies
  11. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dladiv.f">
  12. *> [TGZ]</a>
  13. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dladiv.f">
  14. *> [ZIP]</a>
  15. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dladiv.f">
  16. *> [TXT]</a>
  17. *> \endhtmlonly
  18. *
  19. * Definition:
  20. * ===========
  21. *
  22. * SUBROUTINE DLADIV( A, B, C, D, P, Q )
  23. *
  24. * .. Scalar Arguments ..
  25. * REAL*8 A, B, C, D, P, Q
  26. * ..
  27. *
  28. *
  29. *> \par Purpose:
  30. * =============
  31. *>
  32. *> \verbatim
  33. *>
  34. *> DLADIV performs complex division in real arithmetic
  35. *>
  36. *> a + i*b
  37. *> p + i*q = ---------
  38. *> c + i*d
  39. *>
  40. *> The algorithm is due to Michael Baudin and Robert L. Smith
  41. *> and can be found in the paper
  42. *> "A Robust Complex Division in Scilab"
  43. *> \endverbatim
  44. *
  45. * Arguments:
  46. * ==========
  47. *
  48. *> \param[in] A
  49. *> \verbatim
  50. *> A is DOUBLE PRECISION
  51. *> \endverbatim
  52. *>
  53. *> \param[in] B
  54. *> \verbatim
  55. *> B is DOUBLE PRECISION
  56. *> \endverbatim
  57. *>
  58. *> \param[in] C
  59. *> \verbatim
  60. *> C is DOUBLE PRECISION
  61. *> \endverbatim
  62. *>
  63. *> \param[in] D
  64. *> \verbatim
  65. *> D is DOUBLE PRECISION
  66. *> The scalars a, b, c, and d in the above expression.
  67. *> \endverbatim
  68. *>
  69. *> \param[out] P
  70. *> \verbatim
  71. *> P is DOUBLE PRECISION
  72. *> \endverbatim
  73. *>
  74. *> \param[out] Q
  75. *> \verbatim
  76. *> Q is DOUBLE PRECISION
  77. *> The scalars p and q in the above expression.
  78. *> \endverbatim
  79. *
  80. * Authors:
  81. * ========
  82. *
  83. *> \author Univ. of Tennessee
  84. *> \author Univ. of California Berkeley
  85. *> \author Univ. of Colorado Denver
  86. *> \author NAG Ltd.
  87. *
  88. *> \date January 2013
  89. *
  90. *> \ingroup auxOTHERauxiliary
  91. *
  92. * =====================================================================
  93. SUBROUTINE DLADIV( A, B, C, D, P, Q )
  94. *
  95. * -- LAPACK auxiliary routine (version 3.5.0) --
  96. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  97. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  98. * January 2013
  99. *
  100. * .. Scalar Arguments ..
  101. REAL*8 A, B, C, D, P, Q
  102. * ..
  103. *
  104. * =====================================================================
  105. *
  106. * .. Parameters ..
  107. REAL*8 BS
  108. PARAMETER ( BS = 2.0D0 )
  109. REAL*8 HALF
  110. PARAMETER ( HALF = 0.5D0 )
  111. REAL*8 TWO
  112. PARAMETER ( TWO = 2.0D0 )
  113. *
  114. * .. Local Scalars ..
  115. REAL*8 AA, BB, CC, DD, AB, CD, S, OV, UN, BE, EPS
  116. * ..
  117. c * .. External Functions ..
  118. c REAL*8 DLAMCH
  119. c EXTERNAL DLAMCH
  120. c * ..
  121. c * .. External Subroutines ..
  122. c EXTERNAL DLADIV1
  123. * ..
  124. ** .. Intrinsic Functions ..
  125. * INTRINSIC ABS, MAX
  126. ** ..
  127. ** .. Executable Statements ..
  128. *
  129. AA = A
  130. BB = B
  131. CC = C
  132. DD = D
  133. AB = MAX( ABS(A), ABS(B) )
  134. CD = MAX( ABS(C), ABS(D) )
  135. S = 1.0D0
  136.  
  137. OV = DLAMCH( 'Overflow threshold' )
  138. UN = DLAMCH( 'Safe minimum' )
  139. EPS = DLAMCH( 'Epsilon' )
  140. BE = BS / (EPS*EPS)
  141.  
  142. IF( AB .GE. HALF*OV ) THEN
  143. AA = HALF * AA
  144. BB = HALF * BB
  145. S = TWO * S
  146. END IF
  147. IF( CD .GE. HALF*OV ) THEN
  148. CC = HALF * CC
  149. DD = HALF * DD
  150. S = HALF * S
  151. END IF
  152. IF( AB .LE. UN*BS/EPS ) THEN
  153. AA = AA * BE
  154. BB = BB * BE
  155. S = S / BE
  156. END IF
  157. IF( CD .LE. UN*BS/EPS ) THEN
  158. CC = CC * BE
  159. DD = DD * BE
  160. S = S * BE
  161. END IF
  162. IF( ABS( D ).LE.ABS( C ) ) THEN
  163. CALL DLADIV1(AA, BB, CC, DD, P, Q)
  164. ELSE
  165. CALL DLADIV1(BB, AA, DD, CC, P, Q)
  166. Q = -Q
  167. END IF
  168. P = P * S
  169. Q = Q * S
  170. *
  171. RETURN
  172. *
  173. * End of DLADIV
  174. *
  175. END
  176.  
  177.  
  178.  
  179. SUBROUTINE DLADIV1( A, B, C, D, P, Q )
  180. *
  181. * -- LAPACK auxiliary routine (version 3.5.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. REAL*8 A, B, C, D, P, Q
  188. * ..
  189. *
  190. * =====================================================================
  191. *
  192. * .. Parameters ..
  193. REAL*8 ONE
  194. PARAMETER ( ONE = 1.0D0 )
  195. *
  196. * .. Local Scalars ..
  197. REAL*8 R, T
  198. * ..
  199. * .. External Functions ..
  200. REAL*8 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.  
  218.  
  219.  
  220.  
  221.  
  222. FUNCTION DLADIV2( A, B, C, D, R, T )
  223. *
  224. * -- LAPACK auxiliary routine (version 3.5.0) --
  225. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  226. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  227. * January 2013
  228. *
  229. * .. Scalar Arguments ..
  230. REAL*8 A, B, C, D, R, T
  231. REAL*8 DLADIV2
  232. * ..
  233. *
  234. * =====================================================================
  235. *
  236. * .. Parameters ..
  237. REAL*8 ZERO
  238. PARAMETER ( ZERO = 0.0D0 )
  239. *
  240. * .. Local Scalars ..
  241. REAL*8 BR
  242. * ..
  243. * .. Executable Statements ..
  244. *
  245. IF ( R.NE.ZERO ) THEN
  246. BR = B * R
  247. if( BR.NE.ZERO ) THEN
  248. DLADIV2 = (A + BR) * T
  249. ELSE
  250. DLADIV2 = A * T + (B * T) * R
  251. END IF
  252. ELSE
  253. DLADIV2 = (A + D * (B / C)) * T
  254. END IF
  255. *
  256. RETURN
  257. *
  258. * End of DLADIV12
  259. *
  260. END
  261.  
  262.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales