clarfg.f 5.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203
  1. *> \brief \b CLARFG
  2. *
  3. * =========== DOCUMENTATION ===========
  4. *
  5. * Online html documentation available at
  6. * http://www.netlib.org/lapack/explore-html/
  7. *
  8. *> \htmlonly
  9. *> Download CLARFG + dependencies
  10. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clarfg.f">
  11. *> [TGZ]</a>
  12. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clarfg.f">
  13. *> [ZIP]</a>
  14. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clarfg.f">
  15. *> [TXT]</a>
  16. *> \endhtmlonly
  17. *
  18. * Definition:
  19. * ===========
  20. *
  21. * SUBROUTINE CLARFG( N, ALPHA, X, INCX, TAU )
  22. *
  23. * .. Scalar Arguments ..
  24. * INTEGER INCX, N
  25. * COMPLEX ALPHA, TAU
  26. * ..
  27. * .. Array Arguments ..
  28. * COMPLEX X( * )
  29. * ..
  30. *
  31. *
  32. *> \par Purpose:
  33. * =============
  34. *>
  35. *> \verbatim
  36. *>
  37. *> CLARFG generates a complex elementary reflector H of order n, such
  38. *> that
  39. *>
  40. *> H**H * ( alpha ) = ( beta ), H**H * H = I.
  41. *> ( x ) ( 0 )
  42. *>
  43. *> where alpha and beta are scalars, with beta real, and x is an
  44. *> (n-1)-element complex vector. H is represented in the form
  45. *>
  46. *> H = I - tau * ( 1 ) * ( 1 v**H ) ,
  47. *> ( v )
  48. *>
  49. *> where tau is a complex scalar and v is a complex (n-1)-element
  50. *> vector. Note that H is not hermitian.
  51. *>
  52. *> If the elements of x are all zero and alpha is real, then tau = 0
  53. *> and H is taken to be the unit matrix.
  54. *>
  55. *> Otherwise 1 <= real(tau) <= 2 and abs(tau-1) <= 1 .
  56. *> \endverbatim
  57. *
  58. * Arguments:
  59. * ==========
  60. *
  61. *> \param[in] N
  62. *> \verbatim
  63. *> N is INTEGER
  64. *> The order of the elementary reflector.
  65. *> \endverbatim
  66. *>
  67. *> \param[in,out] ALPHA
  68. *> \verbatim
  69. *> ALPHA is COMPLEX
  70. *> On entry, the value alpha.
  71. *> On exit, it is overwritten with the value beta.
  72. *> \endverbatim
  73. *>
  74. *> \param[in,out] X
  75. *> \verbatim
  76. *> X is COMPLEX array, dimension
  77. *> (1+(N-2)*abs(INCX))
  78. *> On entry, the vector x.
  79. *> On exit, it is overwritten with the vector v.
  80. *> \endverbatim
  81. *>
  82. *> \param[in] INCX
  83. *> \verbatim
  84. *> INCX is INTEGER
  85. *> The increment between elements of X. INCX > 0.
  86. *> \endverbatim
  87. *>
  88. *> \param[out] TAU
  89. *> \verbatim
  90. *> TAU is COMPLEX
  91. *> The value tau.
  92. *> \endverbatim
  93. *
  94. * Authors:
  95. * ========
  96. *
  97. *> \author Univ. of Tennessee
  98. *> \author Univ. of California Berkeley
  99. *> \author Univ. of Colorado Denver
  100. *> \author NAG Ltd.
  101. *
  102. *> \date November 2011
  103. *
  104. *> \ingroup complexOTHERauxiliary
  105. *
  106. * =====================================================================
  107. SUBROUTINE CLARFG( N, ALPHA, X, INCX, TAU )
  108. *
  109. * -- LAPACK auxiliary routine (version 3.4.0) --
  110. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  111. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  112. * November 2011
  113. *
  114. * .. Scalar Arguments ..
  115. INTEGER INCX, N
  116. COMPLEX ALPHA, TAU
  117. * ..
  118. * .. Array Arguments ..
  119. COMPLEX X( * )
  120. * ..
  121. *
  122. * =====================================================================
  123. *
  124. * .. Parameters ..
  125. REAL ONE, ZERO
  126. PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
  127. * ..
  128. * .. Local Scalars ..
  129. INTEGER J, KNT
  130. REAL ALPHI, ALPHR, BETA, RSAFMN, SAFMIN, XNORM
  131. * ..
  132. * .. External Functions ..
  133. REAL SCNRM2, SLAMCH, SLAPY3
  134. COMPLEX CLADIV
  135. EXTERNAL SCNRM2, SLAMCH, SLAPY3, CLADIV
  136. * ..
  137. * .. Intrinsic Functions ..
  138. INTRINSIC ABS, AIMAG, CMPLX, REAL, SIGN
  139. * ..
  140. * .. External Subroutines ..
  141. EXTERNAL CSCAL, CSSCAL
  142. * ..
  143. * .. Executable Statements ..
  144. *
  145. IF( N.LE.0 ) THEN
  146. TAU = ZERO
  147. RETURN
  148. END IF
  149. *
  150. XNORM = SCNRM2( N-1, X, INCX )
  151. ALPHR = REAL( ALPHA )
  152. ALPHI = AIMAG( ALPHA )
  153. *
  154. IF( XNORM.EQ.ZERO .AND. ALPHI.EQ.ZERO ) THEN
  155. *
  156. * H = I
  157. *
  158. TAU = ZERO
  159. ELSE
  160. *
  161. * general case
  162. *
  163. BETA = -SIGN( SLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
  164. SAFMIN = SLAMCH( 'S' ) / SLAMCH( 'E' )
  165. RSAFMN = ONE / SAFMIN
  166. *
  167. KNT = 0
  168. IF( ABS( BETA ).LT.SAFMIN ) THEN
  169. *
  170. * XNORM, BETA may be inaccurate; scale X and recompute them
  171. *
  172. 10 CONTINUE
  173. KNT = KNT + 1
  174. CALL CSSCAL( N-1, RSAFMN, X, INCX )
  175. BETA = BETA*RSAFMN
  176. ALPHI = ALPHI*RSAFMN
  177. ALPHR = ALPHR*RSAFMN
  178. IF( ABS( BETA ).LT.SAFMIN )
  179. $ GO TO 10
  180. *
  181. * New BETA is at most 1, at least SAFMIN
  182. *
  183. XNORM = SCNRM2( N-1, X, INCX )
  184. ALPHA = CMPLX( ALPHR, ALPHI )
  185. BETA = -SIGN( SLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
  186. END IF
  187. TAU = CMPLX( ( BETA-ALPHR ) / BETA, -ALPHI / BETA )
  188. ALPHA = CLADIV( CMPLX( ONE ), ALPHA-BETA )
  189. CALL CSCAL( N-1, ALPHA, X, INCX )
  190. *
  191. * If ALPHA is subnormal, it may lose relative accuracy
  192. *
  193. DO 20 J = 1, KNT
  194. BETA = BETA*SAFMIN
  195. 20 CONTINUE
  196. ALPHA = BETA
  197. END IF
  198. *
  199. RETURN
  200. *
  201. * End of CLARFG
  202. *
  203. END