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.

ctpmqrt.f 11 kB

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