You can not select more than 25 topics Topics must start with a chinese character,a letter or number, can include dashes ('-') and can be up to 35 characters long.

clartg.f 7.5 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258
  1. *> \brief \b CLARTG generates a plane rotation with real cosine and complex sine.
  2. *
  3. * =========== DOCUMENTATION ===========
  4. *
  5. * Online html documentation available at
  6. * http://www.netlib.org/lapack/explore-html/
  7. *
  8. *> \htmlonly
  9. *> Download CLARTG + dependencies
  10. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clartg.f">
  11. *> [TGZ]</a>
  12. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clartg.f">
  13. *> [ZIP]</a>
  14. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clartg.f">
  15. *> [TXT]</a>
  16. *> \endhtmlonly
  17. *
  18. * Definition:
  19. * ===========
  20. *
  21. * SUBROUTINE CLARTG( F, G, CS, SN, R )
  22. *
  23. * .. Scalar Arguments ..
  24. * REAL CS
  25. * COMPLEX F, G, R, SN
  26. * ..
  27. *
  28. *
  29. *> \par Purpose:
  30. * =============
  31. *>
  32. *> \verbatim
  33. *>
  34. *> CLARTG generates a plane rotation so that
  35. *>
  36. *> [ CS SN ] [ F ] [ R ]
  37. *> [ __ ] . [ ] = [ ] where CS**2 + |SN|**2 = 1.
  38. *> [ -SN CS ] [ G ] [ 0 ]
  39. *>
  40. *> This is a faster version of the BLAS1 routine CROTG, except for
  41. *> the following differences:
  42. *> F and G are unchanged on return.
  43. *> If G=0, then CS=1 and SN=0.
  44. *> If F=0, then CS=0 and SN is chosen so that R is real.
  45. *> \endverbatim
  46. *
  47. * Arguments:
  48. * ==========
  49. *
  50. *> \param[in] F
  51. *> \verbatim
  52. *> F is COMPLEX
  53. *> The first component of vector to be rotated.
  54. *> \endverbatim
  55. *>
  56. *> \param[in] G
  57. *> \verbatim
  58. *> G is COMPLEX
  59. *> The second component of vector to be rotated.
  60. *> \endverbatim
  61. *>
  62. *> \param[out] CS
  63. *> \verbatim
  64. *> CS is REAL
  65. *> The cosine of the rotation.
  66. *> \endverbatim
  67. *>
  68. *> \param[out] SN
  69. *> \verbatim
  70. *> SN is COMPLEX
  71. *> The sine of the rotation.
  72. *> \endverbatim
  73. *>
  74. *> \param[out] R
  75. *> \verbatim
  76. *> R is COMPLEX
  77. *> The nonzero component of the rotated vector.
  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 September 2012
  89. *
  90. *> \ingroup complexOTHERauxiliary
  91. *
  92. *> \par Further Details:
  93. * =====================
  94. *>
  95. *> \verbatim
  96. *>
  97. *> 3-5-96 - Modified with a new algorithm by W. Kahan and J. Demmel
  98. *>
  99. *> This version has a few statements commented out for thread safety
  100. *> (machine parameters are computed on each entry). 10 feb 03, SJH.
  101. *> \endverbatim
  102. *>
  103. * =====================================================================
  104. SUBROUTINE CLARTG( F, G, CS, SN, R )
  105. *
  106. * -- LAPACK auxiliary routine (version 3.4.2) --
  107. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  108. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  109. * September 2012
  110. *
  111. * .. Scalar Arguments ..
  112. REAL CS
  113. COMPLEX F, G, R, SN
  114. * ..
  115. *
  116. * =====================================================================
  117. *
  118. * .. Parameters ..
  119. REAL TWO, ONE, ZERO
  120. PARAMETER ( TWO = 2.0E+0, ONE = 1.0E+0, ZERO = 0.0E+0 )
  121. COMPLEX CZERO
  122. PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ) )
  123. * ..
  124. * .. Local Scalars ..
  125. * LOGICAL FIRST
  126. INTEGER COUNT, I
  127. REAL D, DI, DR, EPS, F2, F2S, G2, G2S, SAFMIN,
  128. $ SAFMN2, SAFMX2, SCALE
  129. COMPLEX FF, FS, GS
  130. * ..
  131. * .. External Functions ..
  132. REAL SLAMCH, SLAPY2
  133. EXTERNAL SLAMCH, SLAPY2
  134. * ..
  135. * .. Intrinsic Functions ..
  136. INTRINSIC ABS, AIMAG, CMPLX, CONJG, INT, LOG, MAX, REAL,
  137. $ SQRT
  138. * ..
  139. * .. Statement Functions ..
  140. REAL ABS1, ABSSQ
  141. * ..
  142. * .. Save statement ..
  143. * SAVE FIRST, SAFMX2, SAFMIN, SAFMN2
  144. * ..
  145. * .. Data statements ..
  146. * DATA FIRST / .TRUE. /
  147. * ..
  148. * .. Statement Function definitions ..
  149. ABS1( FF ) = MAX( ABS( REAL( FF ) ), ABS( AIMAG( FF ) ) )
  150. ABSSQ( FF ) = REAL( FF )**2 + AIMAG( FF )**2
  151. * ..
  152. * .. Executable Statements ..
  153. *
  154. * IF( FIRST ) THEN
  155. SAFMIN = SLAMCH( 'S' )
  156. EPS = SLAMCH( 'E' )
  157. SAFMN2 = SLAMCH( 'B' )**INT( LOG( SAFMIN / EPS ) /
  158. $ LOG( SLAMCH( 'B' ) ) / TWO )
  159. SAFMX2 = ONE / SAFMN2
  160. * FIRST = .FALSE.
  161. * END IF
  162. SCALE = MAX( ABS1( F ), ABS1( G ) )
  163. FS = F
  164. GS = G
  165. COUNT = 0
  166. IF( SCALE.GE.SAFMX2 ) THEN
  167. 10 CONTINUE
  168. COUNT = COUNT + 1
  169. FS = FS*SAFMN2
  170. GS = GS*SAFMN2
  171. SCALE = SCALE*SAFMN2
  172. IF( SCALE.GE.SAFMX2 )
  173. $ GO TO 10
  174. ELSE IF( SCALE.LE.SAFMN2 ) THEN
  175. IF( G.EQ.CZERO ) THEN
  176. CS = ONE
  177. SN = CZERO
  178. R = F
  179. RETURN
  180. END IF
  181. 20 CONTINUE
  182. COUNT = COUNT - 1
  183. FS = FS*SAFMX2
  184. GS = GS*SAFMX2
  185. SCALE = SCALE*SAFMX2
  186. IF( SCALE.LE.SAFMN2 )
  187. $ GO TO 20
  188. END IF
  189. F2 = ABSSQ( FS )
  190. G2 = ABSSQ( GS )
  191. IF( F2.LE.MAX( G2, ONE )*SAFMIN ) THEN
  192. *
  193. * This is a rare case: F is very small.
  194. *
  195. IF( F.EQ.CZERO ) THEN
  196. CS = ZERO
  197. R = SLAPY2( REAL( G ), AIMAG( G ) )
  198. * Do complex/real division explicitly with two real divisions
  199. D = SLAPY2( REAL( GS ), AIMAG( GS ) )
  200. SN = CMPLX( REAL( GS ) / D, -AIMAG( GS ) / D )
  201. RETURN
  202. END IF
  203. F2S = SLAPY2( REAL( FS ), AIMAG( FS ) )
  204. * G2 and G2S are accurate
  205. * G2 is at least SAFMIN, and G2S is at least SAFMN2
  206. G2S = SQRT( G2 )
  207. * Error in CS from underflow in F2S is at most
  208. * UNFL / SAFMN2 .lt. sqrt(UNFL*EPS) .lt. EPS
  209. * If MAX(G2,ONE)=G2, then F2 .lt. G2*SAFMIN,
  210. * and so CS .lt. sqrt(SAFMIN)
  211. * If MAX(G2,ONE)=ONE, then F2 .lt. SAFMIN
  212. * and so CS .lt. sqrt(SAFMIN)/SAFMN2 = sqrt(EPS)
  213. * Therefore, CS = F2S/G2S / sqrt( 1 + (F2S/G2S)**2 ) = F2S/G2S
  214. CS = F2S / G2S
  215. * Make sure abs(FF) = 1
  216. * Do complex/real division explicitly with 2 real divisions
  217. IF( ABS1( F ).GT.ONE ) THEN
  218. D = SLAPY2( REAL( F ), AIMAG( F ) )
  219. FF = CMPLX( REAL( F ) / D, AIMAG( F ) / D )
  220. ELSE
  221. DR = SAFMX2*REAL( F )
  222. DI = SAFMX2*AIMAG( F )
  223. D = SLAPY2( DR, DI )
  224. FF = CMPLX( DR / D, DI / D )
  225. END IF
  226. SN = FF*CMPLX( REAL( GS ) / G2S, -AIMAG( GS ) / G2S )
  227. R = CS*F + SN*G
  228. ELSE
  229. *
  230. * This is the most common case.
  231. * Neither F2 nor F2/G2 are less than SAFMIN
  232. * F2S cannot overflow, and it is accurate
  233. *
  234. F2S = SQRT( ONE+G2 / F2 )
  235. * Do the F2S(real)*FS(complex) multiply with two real multiplies
  236. R = CMPLX( F2S*REAL( FS ), F2S*AIMAG( FS ) )
  237. CS = ONE / F2S
  238. D = F2 + G2
  239. * Do complex/real division explicitly with two real divisions
  240. SN = CMPLX( REAL( R ) / D, AIMAG( R ) / D )
  241. SN = SN*CONJG( GS )
  242. IF( COUNT.NE.0 ) THEN
  243. IF( COUNT.GT.0 ) THEN
  244. DO 30 I = 1, COUNT
  245. R = R*SAFMX2
  246. 30 CONTINUE
  247. ELSE
  248. DO 40 I = 1, -COUNT
  249. R = R*SAFMN2
  250. 40 CONTINUE
  251. END IF
  252. END IF
  253. END IF
  254. RETURN
  255. *
  256. * End of CLARTG
  257. *
  258. END