Télécharger dlaev2.eso

Retour à la liste

Numérotation des lignes :

  1. C DLAEV2 SOURCE BP208322 15/10/13 21:15:26 8670
  2. *> \brief \b DLAEV2 computes the eigenvalues and eigenvectors of a 2-by-2 symmetric/Hermitian matrix.
  3. *
  4. * =========== DOCUMENTATION ===========
  5. *
  6. * Online html documentation available at
  7. * http://www.netlib.org/lapack/explore-html/
  8. *
  9. *> \htmlonly
  10. *> Download DLAEV2 + dependencies
  11. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaev2.f">
  12. *> [TGZ]</a>
  13. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaev2.f">
  14. *> [ZIP]</a>
  15. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaev2.f">
  16. *> [TXT]</a>
  17. *> \endhtmlonly
  18. *
  19. * Definition:
  20. * ===========
  21. *
  22. * SUBROUTINE DLAEV2( A, B, C, RT1, RT2, CS1, SN1 )
  23. *
  24. * .. Scalar Arguments ..
  25. * REAL*8 A, B, C, CS1, RT1, RT2, SN1
  26. * ..
  27. *
  28. *
  29. *> \par Purpose:
  30. * =============
  31. *>
  32. *> \verbatim
  33. *>
  34. *> DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix
  35. *> [ A B ]
  36. *> [ B C ].
  37. *> On return, RT1 is the eigenvalue of larger absolute value, RT2 is the
  38. *> eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right
  39. *> eigenvector for RT1, giving the decomposition
  40. *>
  41. *> [ CS1 SN1 ] [ A B ] [ CS1 -SN1 ] = [ RT1 0 ]
  42. *> [-SN1 CS1 ] [ B C ] [ SN1 CS1 ] [ 0 RT2 ].
  43. *> \endverbatim
  44. *
  45. * Arguments:
  46. * ==========
  47. *
  48. *> \param[in] A
  49. *> \verbatim
  50. *> A is DOUBLE PRECISION
  51. *> The (1,1) element of the 2-by-2 matrix.
  52. *> \endverbatim
  53. *>
  54. *> \param[in] B
  55. *> \verbatim
  56. *> B is DOUBLE PRECISION
  57. *> The (1,2) element and the conjugate of the (2,1) element of
  58. *> the 2-by-2 matrix.
  59. *> \endverbatim
  60. *>
  61. *> \param[in] C
  62. *> \verbatim
  63. *> C is DOUBLE PRECISION
  64. *> The (2,2) element of the 2-by-2 matrix.
  65. *> \endverbatim
  66. *>
  67. *> \param[out] RT1
  68. *> \verbatim
  69. *> RT1 is DOUBLE PRECISION
  70. *> The eigenvalue of larger absolute value.
  71. *> \endverbatim
  72. *>
  73. *> \param[out] RT2
  74. *> \verbatim
  75. *> RT2 is DOUBLE PRECISION
  76. *> The eigenvalue of smaller absolute value.
  77. *> \endverbatim
  78. *>
  79. *> \param[out] CS1
  80. *> \verbatim
  81. *> CS1 is DOUBLE PRECISION
  82. *> \endverbatim
  83. *>
  84. *> \param[out] SN1
  85. *> \verbatim
  86. *> SN1 is DOUBLE PRECISION
  87. *> The vector (CS1, SN1) is a unit right eigenvector for RT1.
  88. *> \endverbatim
  89. *
  90. * Authors:
  91. * ========
  92. *
  93. *> \author Univ. of Tennessee
  94. *> \author Univ. of California Berkeley
  95. *> \author Univ. of Colorado Denver
  96. *> \author NAG Ltd.
  97. *
  98. *> \date September 2012
  99. *
  100. *> \ingroup auxOTHERauxiliary
  101. *
  102. *> \par Further Details:
  103. * =====================
  104. *>
  105. *> \verbatim
  106. *>
  107. *> RT1 is accurate to a few ulps barring over/underflow.
  108. *>
  109. *> RT2 may be inaccurate if there is massive cancellation in the
  110. *> determinant A*C-B*B; higher precision or correctly rounded or
  111. *> correctly truncated arithmetic would be needed to compute RT2
  112. *> accurately in all cases.
  113. *>
  114. *> CS1 and SN1 are accurate to a few ulps barring over/underflow.
  115. *>
  116. *> Overflow is possible only if RT1 is within a factor of 5 of overflow.
  117. *> Underflow is harmless if the input data is 0 or exceeds
  118. *> underflow_threshold / macheps.
  119. *> \endverbatim
  120. *>
  121. * =====================================================================
  122. SUBROUTINE DLAEV2( A, B, C, RT1, RT2, CS1, SN1 )
  123. *
  124. * -- LAPACK auxiliary routine (version 3.4.2) --
  125. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  126. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  127. * September 2012
  128. *
  129. * .. Scalar Arguments ..
  130. REAL*8 A, B, C, CS1, RT1, RT2, SN1
  131. * ..
  132. *
  133. * =====================================================================
  134. *
  135. * .. Parameters ..
  136. REAL*8 ONE
  137. PARAMETER ( ONE = 1.0D0 )
  138. REAL*8 TWO
  139. PARAMETER ( TWO = 2.0D0 )
  140. REAL*8 ZERO
  141. PARAMETER ( ZERO = 0.0D0 )
  142. REAL*8 HALF
  143. PARAMETER ( HALF = 0.5D0 )
  144. * ..
  145. * .. Local Scalars ..
  146. INTEGER SGN1, SGN2
  147. REAL*8 AB, ACMN, ACMX, ACS, ADF, CS, CT, DF, RT, SM,
  148. $ TB, TN
  149. * ..
  150. ** .. Intrinsic Functions ..
  151. * INTRINSIC ABS, SQRT
  152. ** ..
  153. ** .. Executable Statements ..
  154. *
  155. * Compute the eigenvalues
  156. *
  157. SM = A + C
  158. DF = A - C
  159. ADF = ABS( DF )
  160. TB = B + B
  161. AB = ABS( TB )
  162. IF( ABS( A ).GT.ABS( C ) ) THEN
  163. ACMX = A
  164. ACMN = C
  165. ELSE
  166. ACMX = C
  167. ACMN = A
  168. END IF
  169. IF( ADF.GT.AB ) THEN
  170. RT = ADF*SQRT( ONE+( AB / ADF )**2 )
  171. ELSE IF( ADF.LT.AB ) THEN
  172. RT = AB*SQRT( ONE+( ADF / AB )**2 )
  173. ELSE
  174. *
  175. * Includes case AB=ADF=0
  176. *
  177. RT = AB*SQRT( TWO )
  178. END IF
  179. IF( SM.LT.ZERO ) THEN
  180. RT1 = HALF*( SM-RT )
  181. SGN1 = -1
  182. *
  183. * Order of execution important.
  184. * To get fully accurate smaller eigenvalue,
  185. * next line needs to be executed in higher precision.
  186. *
  187. RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
  188. ELSE IF( SM.GT.ZERO ) THEN
  189. RT1 = HALF*( SM+RT )
  190. SGN1 = 1
  191. *
  192. * Order of execution important.
  193. * To get fully accurate smaller eigenvalue,
  194. * next line needs to be executed in higher precision.
  195. *
  196. RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
  197. ELSE
  198. *
  199. * Includes case RT1 = RT2 = 0
  200. *
  201. RT1 = HALF*RT
  202. RT2 = -HALF*RT
  203. SGN1 = 1
  204. END IF
  205. *
  206. * Compute the eigenvector
  207. *
  208. IF( DF.GE.ZERO ) THEN
  209. CS = DF + RT
  210. SGN2 = 1
  211. ELSE
  212. CS = DF - RT
  213. SGN2 = -1
  214. END IF
  215. ACS = ABS( CS )
  216. IF( ACS.GT.AB ) THEN
  217. CT = -TB / CS
  218. SN1 = ONE / SQRT( ONE+CT*CT )
  219. CS1 = CT*SN1
  220. ELSE
  221. IF( AB.EQ.ZERO ) THEN
  222. CS1 = ONE
  223. SN1 = ZERO
  224. ELSE
  225. TN = -CS / TB
  226. CS1 = ONE / SQRT( ONE+TN*TN )
  227. SN1 = TN*CS1
  228. END IF
  229. END IF
  230. IF( SGN1.EQ.SGN2 ) THEN
  231. TN = CS1
  232. CS1 = -SN1
  233. SN1 = TN
  234. END IF
  235. RETURN
  236. *
  237. * End of DLAEV2
  238. *
  239. END
  240.  
  241.  

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