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

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