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.

ctpmlqt.f 9.8 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349
  1. * Definition:
  2. * ===========
  3. *
  4. * SUBROUTINE CTPMLQT( SIDE, TRANS, M, N, K, L, MB, V, LDV, T, LDT,
  5. * A, LDA, B, LDB, WORK, INFO )
  6. *
  7. * .. Scalar Arguments ..
  8. * CHARACTER SIDE, TRANS
  9. * INTEGER INFO, K, LDV, LDA, LDB, M, N, L, MB, LDT
  10. * ..
  11. * .. Array Arguments ..
  12. * COMPLEX V( LDV, * ), A( LDA, * ), B( LDB, * ),
  13. * $ T( LDT, * ), WORK( * )
  14. * ..
  15. *
  16. *
  17. *> \par Purpose:
  18. * =============
  19. *>
  20. *> \verbatim
  21. *>
  22. *> CTPMLQT applies a complex orthogonal matrix Q obtained from a
  23. *> "triangular-pentagonal" complex block reflector H to a general
  24. *> complex matrix C, which consists of two blocks A and B.
  25. *> \endverbatim
  26. *
  27. * Arguments:
  28. * ==========
  29. *
  30. *> \param[in] SIDE
  31. *> \verbatim
  32. *> SIDE is CHARACTER*1
  33. *> = 'L': apply Q or Q**H from the Left;
  34. *> = 'R': apply Q or Q**H from the Right.
  35. *> \endverbatim
  36. *>
  37. *> \param[in] TRANS
  38. *> \verbatim
  39. *> TRANS is CHARACTER*1
  40. *> = 'N': No transpose, apply Q;
  41. *> = 'C': Transpose, apply Q**H.
  42. *> \endverbatim
  43. *>
  44. *> \param[in] M
  45. *> \verbatim
  46. *> M is INTEGER
  47. *> The number of rows of the matrix B. M >= 0.
  48. *> \endverbatim
  49. *>
  50. *> \param[in] N
  51. *> \verbatim
  52. *> N is INTEGER
  53. *> The number of columns of the matrix B. N >= 0.
  54. *> \endverbatim
  55. *>
  56. *> \param[in] K
  57. *> \verbatim
  58. *> K is INTEGER
  59. *> The number of elementary reflectors whose product defines
  60. *> the matrix Q.
  61. *> \endverbatim
  62. *>
  63. *> \param[in] L
  64. *> \verbatim
  65. *> L is INTEGER
  66. *> The order of the trapezoidal part of V.
  67. *> K >= L >= 0. See Further Details.
  68. *> \endverbatim
  69. *>
  70. *> \param[in] MB
  71. *> \verbatim
  72. *> MB is INTEGER
  73. *> The block size used for the storage of T. K >= MB >= 1.
  74. *> This must be the same value of MB used to generate T
  75. *> in DTPLQT.
  76. *> \endverbatim
  77. *>
  78. *> \param[in] V
  79. *> \verbatim
  80. *> V is COMPLEX array, dimension (LDA,K)
  81. *> The i-th row must contain the vector which defines the
  82. *> elementary reflector H(i), for i = 1,2,...,k, as returned by
  83. *> DTPLQT in B. See Further Details.
  84. *> \endverbatim
  85. *>
  86. *> \param[in] LDV
  87. *> \verbatim
  88. *> LDV is INTEGER
  89. *> The leading dimension of the array V.
  90. *> If SIDE = 'L', LDV >= max(1,M);
  91. *> if SIDE = 'R', LDV >= max(1,N).
  92. *> \endverbatim
  93. *>
  94. *> \param[in] T
  95. *> \verbatim
  96. *> T is COMPLEX array, dimension (LDT,K)
  97. *> The upper triangular factors of the block reflectors
  98. *> as returned by DTPLQT, stored as a MB-by-K matrix.
  99. *> \endverbatim
  100. *>
  101. *> \param[in] LDT
  102. *> \verbatim
  103. *> LDT is INTEGER
  104. *> The leading dimension of the array T. LDT >= MB.
  105. *> \endverbatim
  106. *>
  107. *> \param[in,out] A
  108. *> \verbatim
  109. *> A is COMPLEX array, dimension
  110. *> (LDA,N) if SIDE = 'L' or
  111. *> (LDA,K) if SIDE = 'R'
  112. *> On entry, the K-by-N or M-by-K matrix A.
  113. *> On exit, A is overwritten by the corresponding block of
  114. *> Q*C or Q**H*C or C*Q or C*Q**H. See Further Details.
  115. *> \endverbatim
  116. *>
  117. *> \param[in] LDA
  118. *> \verbatim
  119. *> LDA is INTEGER
  120. *> The leading dimension of the array A.
  121. *> If SIDE = 'L', LDC >= max(1,K);
  122. *> If SIDE = 'R', LDC >= max(1,M).
  123. *> \endverbatim
  124. *>
  125. *> \param[in,out] B
  126. *> \verbatim
  127. *> B is COMPLEX array, dimension (LDB,N)
  128. *> On entry, the M-by-N matrix B.
  129. *> On exit, B is overwritten by the corresponding block of
  130. *> Q*C or Q**H*C or C*Q or C*Q**H. See Further Details.
  131. *> \endverbatim
  132. *>
  133. *> \param[in] LDB
  134. *> \verbatim
  135. *> LDB is INTEGER
  136. *> The leading dimension of the array B.
  137. *> LDB >= max(1,M).
  138. *> \endverbatim
  139. *>
  140. *> \param[out] WORK
  141. *> \verbatim
  142. *> WORK is COMPLEX array. The dimension of WORK is
  143. *> N*MB if SIDE = 'L', or M*MB if SIDE = 'R'.
  144. *> \endverbatim
  145. *>
  146. *> \param[out] INFO
  147. *> \verbatim
  148. *> INFO is INTEGER
  149. *> = 0: successful exit
  150. *> < 0: if INFO = -i, the i-th argument had an illegal value
  151. *> \endverbatim
  152. *
  153. * Authors:
  154. * ========
  155. *
  156. *> \author Univ. of Tennessee
  157. *> \author Univ. of California Berkeley
  158. *> \author Univ. of Colorado Denver
  159. *> \author NAG Ltd.
  160. *
  161. *> \date June 2017
  162. *
  163. *> \ingroup doubleOTHERcomputational
  164. *
  165. *> \par Further Details:
  166. * =====================
  167. *>
  168. *> \verbatim
  169. *>
  170. *> The columns of the pentagonal matrix V contain the elementary reflectors
  171. *> H(1), H(2), ..., H(K); V is composed of a rectangular block V1 and a
  172. *> trapezoidal block V2:
  173. *>
  174. *> V = [V1] [V2].
  175. *>
  176. *>
  177. *> The size of the trapezoidal block V2 is determined by the parameter L,
  178. *> where 0 <= L <= K; V2 is lower trapezoidal, consisting of the first L
  179. *> rows of a K-by-K upper triangular matrix. If L=K, V2 is lower triangular;
  180. *> if L=0, there is no trapezoidal block, hence V = V1 is rectangular.
  181. *>
  182. *> If SIDE = 'L': C = [A] where A is K-by-N, B is M-by-N and V is K-by-M.
  183. *> [B]
  184. *>
  185. *> If SIDE = 'R': C = [A B] where A is M-by-K, B is M-by-N and V is K-by-N.
  186. *>
  187. *> The real orthogonal matrix Q is formed from V and T.
  188. *>
  189. *> If TRANS='N' and SIDE='L', C is on exit replaced with Q * C.
  190. *>
  191. *> If TRANS='C' and SIDE='L', C is on exit replaced with Q**H * C.
  192. *>
  193. *> If TRANS='N' and SIDE='R', C is on exit replaced with C * Q.
  194. *>
  195. *> If TRANS='C' and SIDE='R', C is on exit replaced with C * Q**H.
  196. *> \endverbatim
  197. *>
  198. * =====================================================================
  199. SUBROUTINE CTPMLQT( SIDE, TRANS, M, N, K, L, MB, V, LDV, T, LDT,
  200. $ A, LDA, B, LDB, WORK, INFO )
  201. *
  202. * -- LAPACK computational routine (version 3.7.1) --
  203. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  204. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  205. * June 2017
  206. *
  207. * .. Scalar Arguments ..
  208. CHARACTER SIDE, TRANS
  209. INTEGER INFO, K, LDV, LDA, LDB, M, N, L, MB, LDT
  210. * ..
  211. * .. Array Arguments ..
  212. COMPLEX V( LDV, * ), A( LDA, * ), B( LDB, * ),
  213. $ T( LDT, * ), WORK( * )
  214. * ..
  215. *
  216. * =====================================================================
  217. *
  218. * ..
  219. * .. Local Scalars ..
  220. LOGICAL LEFT, RIGHT, TRAN, NOTRAN
  221. INTEGER I, IB, NB, LB, KF, LDAQ
  222. * ..
  223. * .. External Functions ..
  224. LOGICAL LSAME
  225. EXTERNAL LSAME
  226. * ..
  227. * .. External Subroutines ..
  228. EXTERNAL XERBLA, CTPRFB
  229. * ..
  230. * .. Intrinsic Functions ..
  231. INTRINSIC MAX, MIN
  232. * ..
  233. * .. Executable Statements ..
  234. *
  235. * .. Test the input arguments ..
  236. *
  237. INFO = 0
  238. LEFT = LSAME( SIDE, 'L' )
  239. RIGHT = LSAME( SIDE, 'R' )
  240. TRAN = LSAME( TRANS, 'C' )
  241. NOTRAN = LSAME( TRANS, 'N' )
  242. *
  243. IF ( LEFT ) THEN
  244. LDAQ = MAX( 1, K )
  245. ELSE IF ( RIGHT ) THEN
  246. LDAQ = MAX( 1, M )
  247. END IF
  248. IF( .NOT.LEFT .AND. .NOT.RIGHT ) THEN
  249. INFO = -1
  250. ELSE IF( .NOT.TRAN .AND. .NOT.NOTRAN ) THEN
  251. INFO = -2
  252. ELSE IF( M.LT.0 ) THEN
  253. INFO = -3
  254. ELSE IF( N.LT.0 ) THEN
  255. INFO = -4
  256. ELSE IF( K.LT.0 ) THEN
  257. INFO = -5
  258. ELSE IF( L.LT.0 .OR. L.GT.K ) THEN
  259. INFO = -6
  260. ELSE IF( MB.LT.1 .OR. (MB.GT.K .AND. K.GT.0) ) THEN
  261. INFO = -7
  262. ELSE IF( LDV.LT.K ) THEN
  263. INFO = -9
  264. ELSE IF( LDT.LT.MB ) THEN
  265. INFO = -11
  266. ELSE IF( LDA.LT.LDAQ ) THEN
  267. INFO = -13
  268. ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
  269. INFO = -15
  270. END IF
  271. *
  272. IF( INFO.NE.0 ) THEN
  273. CALL XERBLA( 'CTPMLQT', -INFO )
  274. RETURN
  275. END IF
  276. *
  277. * .. Quick return if possible ..
  278. *
  279. IF( M.EQ.0 .OR. N.EQ.0 .OR. K.EQ.0 ) RETURN
  280. *
  281. IF( LEFT .AND. NOTRAN ) THEN
  282. *
  283. DO I = 1, K, MB
  284. IB = MIN( MB, K-I+1 )
  285. NB = MIN( M-L+I+IB-1, M )
  286. IF( I.GE.L ) THEN
  287. LB = 0
  288. ELSE
  289. LB = 0
  290. END IF
  291. CALL CTPRFB( 'L', 'C', 'F', 'R', NB, N, IB, LB,
  292. $ V( I, 1 ), LDV, T( 1, I ), LDT,
  293. $ A( I, 1 ), LDA, B, LDB, WORK, IB )
  294. END DO
  295. *
  296. ELSE IF( RIGHT .AND. TRAN ) THEN
  297. *
  298. DO I = 1, K, MB
  299. IB = MIN( MB, K-I+1 )
  300. NB = MIN( N-L+I+IB-1, N )
  301. IF( I.GE.L ) THEN
  302. LB = 0
  303. ELSE
  304. LB = NB-N+L-I+1
  305. END IF
  306. CALL CTPRFB( 'R', 'N', 'F', 'R', M, NB, IB, LB,
  307. $ V( I, 1 ), LDV, T( 1, I ), LDT,
  308. $ A( 1, I ), LDA, B, LDB, WORK, M )
  309. END DO
  310. *
  311. ELSE IF( LEFT .AND. TRAN ) THEN
  312. *
  313. KF = ((K-1)/MB)*MB+1
  314. DO I = KF, 1, -MB
  315. IB = MIN( MB, K-I+1 )
  316. NB = MIN( M-L+I+IB-1, M )
  317. IF( I.GE.L ) THEN
  318. LB = 0
  319. ELSE
  320. LB = 0
  321. END IF
  322. CALL CTPRFB( 'L', 'N', 'F', 'R', NB, N, IB, LB,
  323. $ V( I, 1 ), LDV, T( 1, I ), LDT,
  324. $ A( I, 1 ), LDA, B, LDB, WORK, IB )
  325. END DO
  326. *
  327. ELSE IF( RIGHT .AND. NOTRAN ) THEN
  328. *
  329. KF = ((K-1)/MB)*MB+1
  330. DO I = KF, 1, -MB
  331. IB = MIN( MB, K-I+1 )
  332. NB = MIN( N-L+I+IB-1, N )
  333. IF( I.GE.L ) THEN
  334. LB = 0
  335. ELSE
  336. LB = NB-N+L-I+1
  337. END IF
  338. CALL CTPRFB( 'R', 'C', 'F', 'R', M, NB, IB, LB,
  339. $ V( I, 1 ), LDV, T( 1, I ), LDT,
  340. $ A( 1, I ), LDA, B, LDB, WORK, M )
  341. END DO
  342. *
  343. END IF
  344. *
  345. RETURN
  346. *
  347. * End of CTPMLQT
  348. *
  349. END