Télécharger dlarft.eso

Retour à la liste

Numérotation des lignes :

  1. C DLARFT SOURCE BP208322 20/09/18 21:16:06 10718
  2. *> \brief \b DLARFT forms the triangular factor T of a block reflector H = I - vtvH
  3. *
  4. * =========== DOCUMENTATION ===========
  5. *
  6. * Online html documentation available at
  7. * http://www.netlib.org/lapack/explore-html/
  8. *
  9. *> \htmlonly
  10. *> Download DLARFT + dependencies
  11. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarft.f">
  12. *> [TGZ]</a>
  13. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarft.f">
  14. *> [ZIP]</a>
  15. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarft.f">
  16. *> [TXT]</a>
  17. *> \endhtmlonly
  18. *
  19. * Definition:
  20. * ===========
  21. *
  22. * SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
  23. *
  24. * .. Scalar Arguments ..
  25. * CHARACTER DIRECT, STOREV
  26. * INTEGER K, LDT, LDV, N
  27. * ..
  28. * .. Array Arguments ..
  29. * REAL*8 T( LDT, * ), TAU( * ), V( LDV, * )
  30. * ..
  31. *
  32. *
  33. *> \par Purpose:
  34. * =============
  35. *>
  36. *> \verbatim
  37. *>
  38. *> DLARFT forms the triangular factor T of a real block reflector H
  39. *> of order n, which is defined as a product of k elementary reflectors.
  40. *>
  41. *> If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular;
  42. *>
  43. *> If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular.
  44. *>
  45. *> If STOREV = 'C', the vector which defines the elementary reflector
  46. *> H(i) is stored in the i-th column of the array V, and
  47. *>
  48. *> H = I - V * T * V**T
  49. *>
  50. *> If STOREV = 'R', the vector which defines the elementary reflector
  51. *> H(i) is stored in the i-th row of the array V, and
  52. *>
  53. *> H = I - V**T * T * V
  54. *> \endverbatim
  55. *
  56. * Arguments:
  57. * ==========
  58. *
  59. *> \param[in] DIRECT
  60. *> \verbatim
  61. *> DIRECT is CHARACTER*1
  62. *> Specifies the order in which the elementary reflectors are
  63. *> multiplied to form the block reflector:
  64. *> = 'F': H = H(1) H(2) . . . H(k) (Forward)
  65. *> = 'B': H = H(k) . . . H(2) H(1) (Backward)
  66. *> \endverbatim
  67. *>
  68. *> \param[in] STOREV
  69. *> \verbatim
  70. *> STOREV is CHARACTER*1
  71. *> Specifies how the vectors which define the elementary
  72. *> reflectors are stored (see also Further Details):
  73. *> = 'C': columnwise
  74. *> = 'R': rowwise
  75. *> \endverbatim
  76. *>
  77. *> \param[in] N
  78. *> \verbatim
  79. *> N is INTEGER
  80. *> The order of the block reflector H. N >= 0.
  81. *> \endverbatim
  82. *>
  83. *> \param[in] K
  84. *> \verbatim
  85. *> K is INTEGER
  86. *> The order of the triangular factor T (= the number of
  87. *> elementary reflectors). K >= 1.
  88. *> \endverbatim
  89. *>
  90. *> \param[in] V
  91. *> \verbatim
  92. *> V is REAL*8 array, dimension
  93. *> (LDV,K) if STOREV = 'C'
  94. *> (LDV,N) if STOREV = 'R'
  95. *> The matrix V. See further details.
  96. *> \endverbatim
  97. *>
  98. *> \param[in] LDV
  99. *> \verbatim
  100. *> LDV is INTEGER
  101. *> The leading dimension of the array V.
  102. *> If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K.
  103. *> \endverbatim
  104. *>
  105. *> \param[in] TAU
  106. *> \verbatim
  107. *> TAU is REAL*8 array, dimension (K)
  108. *> TAU(i) must contain the scalar factor of the elementary
  109. *> reflector H(i).
  110. *> \endverbatim
  111. *>
  112. *> \param[out] T
  113. *> \verbatim
  114. *> T is REAL*8 array, dimension (LDT,K)
  115. *> The k by k triangular factor T of the block reflector.
  116. *> If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is
  117. *> lower triangular. The rest of the array is not used.
  118. *> \endverbatim
  119. *>
  120. *> \param[in] LDT
  121. *> \verbatim
  122. *> LDT is INTEGER
  123. *> The leading dimension of the array T. LDT >= K.
  124. *> \endverbatim
  125. *
  126. * Authors:
  127. * ========
  128. *
  129. *> \author Univ. of Tennessee
  130. *> \author Univ. of California Berkeley
  131. *> \author Univ. of Colorado Denver
  132. *> \author NAG Ltd.
  133. *
  134. *> \date December 2016
  135. *
  136. *> \ingroup doubleOTHERauxiliary
  137. *
  138. *> \par Further Details:
  139. * =====================
  140. *>
  141. *> \verbatim
  142. *>
  143. *> The shape of the matrix V and the storage of the vectors which define
  144. *> the H(i) is best illustrated by the following example with n = 5 and
  145. *> k = 3. The elements equal to 1 are not stored.
  146. *>
  147. *> DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R':
  148. *>
  149. *> V = ( 1 ) V = ( 1 v1 v1 v1 v1 )
  150. *> ( v1 1 ) ( 1 v2 v2 v2 )
  151. *> ( v1 v2 1 ) ( 1 v3 v3 )
  152. *> ( v1 v2 v3 )
  153. *> ( v1 v2 v3 )
  154. *>
  155. *> DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R':
  156. *>
  157. *> V = ( v1 v2 v3 ) V = ( v1 v1 1 )
  158. *> ( v1 v2 v3 ) ( v2 v2 v2 1 )
  159. *> ( 1 v2 v3 ) ( v3 v3 v3 v3 1 )
  160. *> ( 1 v3 )
  161. *> ( 1 )
  162. *> \endverbatim
  163. *>
  164. * =====================================================================
  165. SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
  166. *
  167. * -- LAPACK auxiliary routine (version 3.7.0) --
  168. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  169. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  170. * December 2016
  171.  
  172. IMPLICIT INTEGER(I-N)
  173. IMPLICIT REAL*8(A-H,O-Z)
  174. *
  175. * .. Scalar Arguments ..
  176. CHARACTER DIRECT, STOREV
  177. INTEGER K, LDT, LDV, N
  178. * ..
  179. * .. Array Arguments ..
  180. REAL*8 T( LDT, * ), TAU( * ), V( LDV, * )
  181. * ..
  182. *
  183. * =====================================================================
  184. *
  185. * .. Parameters ..
  186. REAL*8 ONE, ZERO
  187. PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  188. * ..
  189. * .. Local Scalars ..
  190. INTEGER I, J, PREVLASTV, LASTV
  191. * ..
  192. * .. External Subroutines ..
  193. * EXTERNAL DGEMV, DTRMV
  194. * ..
  195. * .. External Functions ..
  196. LOGICAL LSAME
  197. * EXTERNAL LSAME
  198. * ..
  199. * .. Executable Statements ..
  200. *
  201. * Quick return if possible
  202. *
  203. IF( N.EQ.0 )
  204. $ RETURN
  205. *
  206. IF( LSAME( DIRECT, 'F' ) ) THEN
  207. PREVLASTV = N
  208. DO I = 1, K
  209. PREVLASTV = MAX( I, PREVLASTV )
  210. IF( TAU( I ).EQ.ZERO ) THEN
  211. *
  212. * H(i) = I
  213. *
  214. DO J = 1, I
  215. T( J, I ) = ZERO
  216. END DO
  217. ELSE
  218. *
  219. * general case
  220. *
  221. IF( LSAME( STOREV, 'C' ) ) THEN
  222. * Skip any trailing zeros.
  223. DO LASTV = N, I+1, -1
  224. IF( V( LASTV, I ).NE.ZERO ) RETURN
  225. END DO
  226. DO J = 1, I-1
  227. T( J, I ) = -TAU( I ) * V( I , J )
  228. END DO
  229. J = MIN( LASTV, PREVLASTV )
  230. *
  231. * T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**T * V(i:j,i)
  232. *
  233. CALL DGEMV( 'Transpose', J-I, I-1, -TAU( I ),
  234. $ V( I+1, 1 ), LDV, V( I+1, I ), 1, ONE,
  235. $ T( 1, I ), 1 )
  236. ELSE
  237. * Skip any trailing zeros.
  238. DO LASTV = N, I+1, -1
  239. IF( V( I, LASTV ).NE.ZERO ) RETURN
  240. END DO
  241. DO J = 1, I-1
  242. T( J, I ) = -TAU( I ) * V( J , I )
  243. END DO
  244. J = MIN( LASTV, PREVLASTV )
  245. *
  246. * T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**T
  247. *
  248. CALL DGEMV( 'No transpose', I-1, J-I, -TAU( I ),
  249. $ V( 1, I+1 ), LDV, V( I, I+1 ), LDV, ONE,
  250. $ T( 1, I ), 1 )
  251. END IF
  252. *
  253. * T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i)
  254. *
  255. CALL DTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T,
  256. $ LDT, T( 1, I ), 1 )
  257. T( I, I ) = TAU( I )
  258. IF( I.GT.1 ) THEN
  259. PREVLASTV = MAX( PREVLASTV, LASTV )
  260. ELSE
  261. PREVLASTV = LASTV
  262. END IF
  263. END IF
  264. END DO
  265. ELSE
  266. PREVLASTV = 1
  267. DO I = K, 1, -1
  268. IF( TAU( I ).EQ.ZERO ) THEN
  269. *
  270. * H(i) = I
  271. *
  272. DO J = I, K
  273. T( J, I ) = ZERO
  274. END DO
  275. ELSE
  276. *
  277. * general case
  278. *
  279. IF( I.LT.K ) THEN
  280. IF( LSAME( STOREV, 'C' ) ) THEN
  281. * Skip any leading zeros.
  282. DO LASTV = 1, I-1
  283. IF( V( LASTV, I ).NE.ZERO ) RETURN
  284. END DO
  285. DO J = I+1, K
  286. T( J, I ) = -TAU( I ) * V( N-K+I , J )
  287. END DO
  288. J = MAX( LASTV, PREVLASTV )
  289. *
  290. * T(i+1:k,i) = -tau(i) * V(j:n-k+i,i+1:k)**T * V(j:n-k+i,i)
  291. *
  292. CALL DGEMV( 'Transpose', N-K+I-J, K-I, -TAU( I ),
  293. $ V( J, I+1 ), LDV, V( J, I ), 1, ONE,
  294. $ T( I+1, I ), 1 )
  295. ELSE
  296. * Skip any leading zeros.
  297. DO LASTV = 1, I-1
  298. IF( V( I, LASTV ).NE.ZERO ) RETURN
  299. END DO
  300. DO J = I+1, K
  301. T( J, I ) = -TAU( I ) * V( J, N-K+I )
  302. END DO
  303. J = MAX( LASTV, PREVLASTV )
  304. *
  305. * T(i+1:k,i) = -tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**T
  306. *
  307. CALL DGEMV( 'No transpose', K-I, N-K+I-J,
  308. $ -TAU( I ), V( I+1, J ), LDV, V( I, J ), LDV,
  309. $ ONE, T( I+1, I ), 1 )
  310. END IF
  311. *
  312. * T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i)
  313. *
  314. CALL DTRMV( 'Lower', 'No transpose', 'Non-unit', K-I,
  315. $ T( I+1, I+1 ), LDT, T( I+1, I ), 1 )
  316. IF( I.GT.1 ) THEN
  317. PREVLASTV = MIN( PREVLASTV, LASTV )
  318. ELSE
  319. PREVLASTV = LASTV
  320. END IF
  321. END IF
  322. T( I, I ) = TAU( I )
  323. END IF
  324. END DO
  325. END IF
  326. RETURN
  327. *
  328. * End of DLARFT
  329. *
  330. END
  331.  
  332.  
  333.  

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