Télécharger dlahr2.eso

Retour à la liste

Numérotation des lignes :

  1. C DLAHR2 SOURCE BP208322 20/09/18 21:15:56 10718
  2. *> \brief \b DLAHR2 reduces the specified number of first columns of a general rectangular matrix A so that elements below the specified subdiagonal are zero, and returns auxiliary matrices which are needed to apply the transformation to the unreduced part of A.
  3. *
  4. * =========== DOCUMENTATION ===========
  5. *
  6. * Online html documentation available at
  7. * http://www.netlib.org/lapack/explore-html/
  8. *
  9. *> \htmlonly
  10. *> Download DLAHR2 + dependencies
  11. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlahr2.f">
  12. *> [TGZ]</a>
  13. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlahr2.f">
  14. *> [ZIP]</a>
  15. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlahr2.f">
  16. *> [TXT]</a>
  17. *> \endhtmlonly
  18. *
  19. * Definition:
  20. * ===========
  21. *
  22. * SUBROUTINE DLAHR2( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY )
  23. *
  24. * .. Scalar Arguments ..
  25. * INTEGER K, LDA, LDT, LDY, N, NB
  26. * ..
  27. * .. Array Arguments ..
  28. * REAL*8 A( LDA, * ), T( LDT, NB ), TAU( NB ),
  29. * $ Y( LDY, NB )
  30. * ..
  31. *
  32. *
  33. *> \par Purpose:
  34. * =============
  35. *>
  36. *> \verbatim
  37. *>
  38. *> DLAHR2 reduces the first NB columns of A real general n-BY-(n-k+1)
  39. *> matrix A so that elements below the k-th subdiagonal are zero. The
  40. *> reduction is performed by an orthogonal similarity transformation
  41. *> Q**T * A * Q. The routine returns the matrices V and T which determine
  42. *> Q as a block reflector I - V*T*V**T, and also the matrix Y = A * V * T.
  43. *>
  44. *> This is an auxiliary routine called by DGEHRD.
  45. *> \endverbatim
  46. *
  47. * Arguments:
  48. * ==========
  49. *
  50. *> \param[in] N
  51. *> \verbatim
  52. *> N is INTEGER
  53. *> The order of the matrix A.
  54. *> \endverbatim
  55. *>
  56. *> \param[in] K
  57. *> \verbatim
  58. *> K is INTEGER
  59. *> The offset for the reduction. Elements below the k-th
  60. *> subdiagonal in the first NB columns are reduced to zero.
  61. *> K < N.
  62. *> \endverbatim
  63. *>
  64. *> \param[in] NB
  65. *> \verbatim
  66. *> NB is INTEGER
  67. *> The number of columns to be reduced.
  68. *> \endverbatim
  69. *>
  70. *> \param[in,out] A
  71. *> \verbatim
  72. *> A is REAL*8 array, dimension (LDA,N-K+1)
  73. *> On entry, the n-by-(n-k+1) general matrix A.
  74. *> On exit, the elements on and above the k-th subdiagonal in
  75. *> the first NB columns are overwritten with the corresponding
  76. *> elements of the reduced matrix; the elements below the k-th
  77. *> subdiagonal, with the array TAU, represent the matrix Q as a
  78. *> product of elementary reflectors. The other columns of A are
  79. *> unchanged. See Further Details.
  80. *> \endverbatim
  81. *>
  82. *> \param[in] LDA
  83. *> \verbatim
  84. *> LDA is INTEGER
  85. *> The leading dimension of the array A. LDA >= max(1,N).
  86. *> \endverbatim
  87. *>
  88. *> \param[out] TAU
  89. *> \verbatim
  90. *> TAU is REAL*8 array, dimension (NB)
  91. *> The scalar factors of the elementary reflectors. See Further
  92. *> Details.
  93. *> \endverbatim
  94. *>
  95. *> \param[out] T
  96. *> \verbatim
  97. *> T is REAL*8 array, dimension (LDT,NB)
  98. *> The upper triangular matrix T.
  99. *> \endverbatim
  100. *>
  101. *> \param[in] LDT
  102. *> \verbatim
  103. *> LDT is INTEGER
  104. *> The leading dimension of the array T. LDT >= NB.
  105. *> \endverbatim
  106. *>
  107. *> \param[out] Y
  108. *> \verbatim
  109. *> Y is REAL*8 array, dimension (LDY,NB)
  110. *> The n-by-nb matrix Y.
  111. *> \endverbatim
  112. *>
  113. *> \param[in] LDY
  114. *> \verbatim
  115. *> LDY is INTEGER
  116. *> The leading dimension of the array Y. LDY >= N.
  117. *> \endverbatim
  118. *
  119. * Authors:
  120. * ========
  121. *
  122. *> \author Univ. of Tennessee
  123. *> \author Univ. of California Berkeley
  124. *> \author Univ. of Colorado Denver
  125. *> \author NAG Ltd.
  126. *
  127. *> \date December 2016
  128. *
  129. *> \ingroup doubleOTHERauxiliary
  130. *
  131. *> \par Further Details:
  132. * =====================
  133. *>
  134. *> \verbatim
  135. *>
  136. *> The matrix Q is represented as a product of nb elementary reflectors
  137. *>
  138. *> Q = H(1) H(2) . . . H(nb).
  139. *>
  140. *> Each H(i) has the form
  141. *>
  142. *> H(i) = I - tau * v * v**T
  143. *>
  144. *> where tau is a real scalar, and v is a real vector with
  145. *> v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in
  146. *> A(i+k+1:n,i), and tau in TAU(i).
  147. *>
  148. *> The elements of the vectors v together form the (n-k+1)-by-nb matrix
  149. *> V which is needed, with T and Y, to apply the transformation to the
  150. *> unreduced part of the matrix, using an update of the form:
  151. *> A := (I - V*T*V**T) * (A - Y*V**T).
  152. *>
  153. *> The contents of A on exit are illustrated by the following example
  154. *> with n = 7, k = 3 and nb = 2:
  155. *>
  156. *> ( a a a a a )
  157. *> ( a a a a a )
  158. *> ( a a a a a )
  159. *> ( h h a a a )
  160. *> ( v1 h a a a )
  161. *> ( v1 v2 a a a )
  162. *> ( v1 v2 a a a )
  163. *>
  164. *> where a denotes an element of the original matrix A, h denotes a
  165. *> modified element of the upper Hessenberg matrix H, and vi denotes an
  166. *> element of the vector defining H(i).
  167. *>
  168. *> This subroutine is a slight modification of LAPACK-3.0's DLAHRD
  169. *> incorporating improvements proposed by Quintana-Orti and Van de
  170. *> Gejin. Note that the entries of A(1:K,2:NB) differ from those
  171. *> returned by the original LAPACK-3.0's DLAHRD routine. (This
  172. *> subroutine is not backward compatible with LAPACK-3.0's DLAHRD.)
  173. *> \endverbatim
  174. *
  175. *> \par References:
  176. * ================
  177. *>
  178. *> Gregorio Quintana-Orti and Robert van de Geijn, "Improving the
  179. *> performance of reduction to Hessenberg form," ACM Transactions on
  180. *> Mathematical Software, 32(2):180-194, June 2006.
  181. *>
  182. * =====================================================================
  183. SUBROUTINE DLAHR2( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY )
  184. *
  185. * -- LAPACK auxiliary routine (version 3.7.0) --
  186. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  187. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  188. * December 2016
  189.  
  190. IMPLICIT INTEGER(I-N)
  191. IMPLICIT REAL*8(A-H,O-Z)
  192. *
  193. * .. Scalar Arguments ..
  194. INTEGER K, LDA, LDT, LDY, N, NB
  195. * ..
  196. * .. Array Arguments ..
  197. REAL*8 A( LDA, * ), T( LDT, NB ), TAU( NB ),
  198. $ Y( LDY, NB )
  199. * ..
  200. *
  201. * =====================================================================
  202. *
  203. * .. Parameters ..
  204. REAL*8 ZERO, ONE
  205. PARAMETER ( ZERO = 0.0D+0,
  206. $ ONE = 1.0D+0 )
  207. * ..
  208. * .. Local Scalars ..
  209. INTEGER I
  210. REAL*8 EI
  211. * ..
  212. * .. External Subroutines ..
  213. * EXTERNAL DAXPY, DCOPY, DGEMM, DGEMV, DLACPY,
  214. * $ DLARFG, DSCAL, DTRMM, DTRMV
  215. * ..
  216. * .. Intrinsic Functions ..
  217. * INTRINSIC MIN
  218. * ..
  219. * .. Executable Statements ..
  220. *
  221. * Quick return if possible
  222. *
  223. IF( N.LE.1 )
  224. $ RETURN
  225. *
  226. DO 10 I = 1, NB
  227. IF( I.GT.1 ) THEN
  228. *
  229. * Update A(K+1:N,I)
  230. *
  231. * Update I-th column of A - Y * V**T
  232. *
  233. CALL DGEMV( 'NO TRANSPOSE', N-K, I-1, -ONE, Y(K+1,1), LDY,
  234. $ A( K+I-1, 1 ), LDA, ONE, A( K+1, I ), 1 )
  235. *
  236. * Apply I - V * T**T * V**T to this column (call it b) from the
  237. * left, using the last column of T as workspace
  238. *
  239. * Let V = ( V1 ) and b = ( b1 ) (first I-1 rows)
  240. * ( V2 ) ( b2 )
  241. *
  242. * where V1 is unit lower triangular
  243. *
  244. * w := V1**T * b1
  245. *
  246. CALL DCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 )
  247. CALL DTRMV( 'Lower', 'Transpose', 'UNIT',
  248. $ I-1, A( K+1, 1 ),
  249. $ LDA, T( 1, NB ), 1 )
  250. *
  251. * w := w + V2**T * b2
  252. *
  253. CALL DGEMV( 'Transpose', N-K-I+1, I-1,
  254. $ ONE, A( K+I, 1 ),
  255. $ LDA, A( K+I, I ), 1, ONE, T( 1, NB ), 1 )
  256. *
  257. * w := T**T * w
  258. *
  259. CALL DTRMV( 'Upper', 'Transpose', 'NON-UNIT',
  260. $ I-1, T, LDT,
  261. $ T( 1, NB ), 1 )
  262. *
  263. * b2 := b2 - V2*w
  264. *
  265. CALL DGEMV( 'NO TRANSPOSE', N-K-I+1, I-1, -ONE,
  266. $ A( K+I, 1 ),
  267. $ LDA, T( 1, NB ), 1, ONE, A( K+I, I ), 1 )
  268. *
  269. * b1 := b1 - V1*w
  270. *
  271. CALL DTRMV( 'Lower', 'NO TRANSPOSE',
  272. $ 'UNIT', I-1,
  273. $ A( K+1, 1 ), LDA, T( 1, NB ), 1 )
  274. CALL DAXPY( I-1, -ONE, T( 1, NB ), 1, A( K+1, I ), 1 )
  275. *
  276. A( K+I-1, I-1 ) = EI
  277. END IF
  278. *
  279. * Generate the elementary reflector H(I) to annihilate
  280. * A(K+I+1:N,I)
  281. *
  282. CALL DLARFG( N-K-I+1, A( K+I, I ), A( MIN( K+I+1, N ), I ), 1,
  283. $ TAU( I ) )
  284. EI = A( K+I, I )
  285. A( K+I, I ) = ONE
  286. *
  287. * Compute Y(K+1:N,I)
  288. *
  289. CALL DGEMV( 'NO TRANSPOSE', N-K, N-K-I+1,
  290. $ ONE, A( K+1, I+1 ),
  291. $ LDA, A( K+I, I ), 1, ZERO, Y( K+1, I ), 1 )
  292. CALL DGEMV( 'Transpose', N-K-I+1, I-1,
  293. $ ONE, A( K+I, 1 ), LDA,
  294. $ A( K+I, I ), 1, ZERO, T( 1, I ), 1 )
  295. CALL DGEMV( 'NO TRANSPOSE', N-K, I-1, -ONE,
  296. $ Y( K+1, 1 ), LDY,
  297. $ T( 1, I ), 1, ONE, Y( K+1, I ), 1 )
  298. CALL DSCAL( N-K, TAU( I ), Y( K+1, I ), 1 )
  299. *
  300. * Compute T(1:I,I)
  301. *
  302. CALL DSCAL( I-1, -TAU( I ), T( 1, I ), 1 )
  303. CALL DTRMV( 'Upper', 'No Transpose', 'NON-UNIT',
  304. $ I-1, T, LDT,
  305. $ T( 1, I ), 1 )
  306. T( I, I ) = TAU( I )
  307. *
  308. 10 CONTINUE
  309. A( K+NB, NB ) = EI
  310. *
  311. * Compute Y(1:K,1:NB)
  312. *
  313. CALL DLACPY( 'ALL', K, NB, A( 1, 2 ), LDA, Y, LDY )
  314. CALL DTRMM( 'RIGHT', 'Lower', 'NO TRANSPOSE',
  315. $ 'UNIT', K, NB,
  316. $ ONE, A( K+1, 1 ), LDA, Y, LDY )
  317. IF( N.GT.K+NB )
  318. $ CALL DGEMM( 'NO TRANSPOSE', 'NO TRANSPOSE', K,
  319. $ NB, N-K-NB, ONE,
  320. $ A( 1, 2+NB ), LDA, A( K+1+NB, 1 ), LDA, ONE, Y,
  321. $ LDY )
  322. CALL DTRMM( 'RIGHT', 'Upper', 'NO TRANSPOSE',
  323. $ 'NON-UNIT', K, NB,
  324. $ ONE, T, LDT, Y, LDY )
  325. *
  326. RETURN
  327. *
  328. * End of DLAHR2
  329. *
  330. END
  331.  
  332.  
  333.  

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