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.

zlaswlq.f 7.8 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270
  1. *> \brief \b ZLASWLQ
  2. *
  3. * Definition:
  4. * ===========
  5. *
  6. * SUBROUTINE ZLASWLQ( M, N, MB, NB, A, LDA, T, LDT, WORK,
  7. * LWORK, INFO)
  8. *
  9. * .. Scalar Arguments ..
  10. * INTEGER INFO, LDA, M, N, MB, NB, LDT, LWORK
  11. * ..
  12. * .. Array Arguments ..
  13. * COMPLEX*16 A( LDA, * ), T( LDT, * ), WORK( * )
  14. * ..
  15. *
  16. *
  17. *> \par Purpose:
  18. * =============
  19. *>
  20. *> \verbatim
  21. *>
  22. *> ZLASWLQ computes a blocked Tall-Skinny LQ factorization of
  23. *> a complexx M-by-N matrix A for M <= N:
  24. *>
  25. *> A = ( L 0 ) * Q,
  26. *>
  27. *> where:
  28. *>
  29. *> Q is a n-by-N orthogonal matrix, stored on exit in an implicit
  30. *> form in the elements above the digonal of the array A and in
  31. *> the elemenst of the array T;
  32. *> L is an lower-triangular M-by-M matrix stored on exit in
  33. *> the elements on and below the diagonal of the array A.
  34. *> 0 is a M-by-(N-M) zero matrix, if M < N, and is not stored.
  35. *>
  36. *> \endverbatim
  37. *
  38. * Arguments:
  39. * ==========
  40. *
  41. *> \param[in] M
  42. *> \verbatim
  43. *> M is INTEGER
  44. *> The number of rows of the matrix A. M >= 0.
  45. *> \endverbatim
  46. *>
  47. *> \param[in] N
  48. *> \verbatim
  49. *> N is INTEGER
  50. *> The number of columns of the matrix A. N >= M >= 0.
  51. *> \endverbatim
  52. *>
  53. *> \param[in] MB
  54. *> \verbatim
  55. *> MB is INTEGER
  56. *> The row block size to be used in the blocked QR.
  57. *> M >= MB >= 1
  58. *> \endverbatim
  59. *> \param[in] NB
  60. *> \verbatim
  61. *> NB is INTEGER
  62. *> The column block size to be used in the blocked QR.
  63. *> NB > M.
  64. *> \endverbatim
  65. *>
  66. *> \param[in,out] A
  67. *> \verbatim
  68. *> A is COMPLEX*16 array, dimension (LDA,N)
  69. *> On entry, the M-by-N matrix A.
  70. *> On exit, the elements on and below the diagonal
  71. *> of the array contain the N-by-N lower triangular matrix L;
  72. *> the elements above the diagonal represent Q by the rows
  73. *> of blocked V (see Further Details).
  74. *>
  75. *> \endverbatim
  76. *>
  77. *> \param[in] LDA
  78. *> \verbatim
  79. *> LDA is INTEGER
  80. *> The leading dimension of the array A. LDA >= max(1,M).
  81. *> \endverbatim
  82. *>
  83. *> \param[out] T
  84. *> \verbatim
  85. *> T is COMPLEX*16 array,
  86. *> dimension (LDT, N * Number_of_row_blocks)
  87. *> where Number_of_row_blocks = CEIL((N-M)/(NB-M))
  88. *> The blocked upper triangular block reflectors stored in compact form
  89. *> as a sequence of upper triangular blocks.
  90. *> See Further Details below.
  91. *> \endverbatim
  92. *>
  93. *> \param[in] LDT
  94. *> \verbatim
  95. *> LDT is INTEGER
  96. *> The leading dimension of the array T. LDT >= MB.
  97. *> \endverbatim
  98. *>
  99. *>
  100. *> \param[out] WORK
  101. *> \verbatim
  102. *> (workspace) COMPLEX*16 array, dimension (MAX(1,LWORK))
  103. *>
  104. *> \endverbatim
  105. *> \param[in] LWORK
  106. *> \verbatim
  107. *> The dimension of the array WORK. LWORK >= MB*M.
  108. *> If LWORK = -1, then a workspace query is assumed; the routine
  109. *> only calculates the optimal size of the WORK array, returns
  110. *> this value as the first entry of the WORK array, and no error
  111. *> message related to LWORK is issued by XERBLA.
  112. *>
  113. *> \endverbatim
  114. *> \param[out] INFO
  115. *> \verbatim
  116. *> INFO is INTEGER
  117. *> = 0: successful exit
  118. *> < 0: if INFO = -i, the i-th argument had an illegal value
  119. *> \endverbatim
  120. *
  121. * Authors:
  122. * ========
  123. *
  124. *> \author Univ. of Tennessee
  125. *> \author Univ. of California Berkeley
  126. *> \author Univ. of Colorado Denver
  127. *> \author NAG Ltd.
  128. *
  129. *> \par Further Details:
  130. * =====================
  131. *>
  132. *> \verbatim
  133. *> Short-Wide LQ (SWLQ) performs LQ by a sequence of orthogonal transformations,
  134. *> representing Q as a product of other orthogonal matrices
  135. *> Q = Q(1) * Q(2) * . . . * Q(k)
  136. *> where each Q(i) zeros out upper diagonal entries of a block of NB rows of A:
  137. *> Q(1) zeros out the upper diagonal entries of rows 1:NB of A
  138. *> Q(2) zeros out the bottom MB-N rows of rows [1:M,NB+1:2*NB-M] of A
  139. *> Q(3) zeros out the bottom MB-N rows of rows [1:M,2*NB-M+1:3*NB-2*M] of A
  140. *> . . .
  141. *>
  142. *> Q(1) is computed by GELQT, which represents Q(1) by Householder vectors
  143. *> stored under the diagonal of rows 1:MB of A, and by upper triangular
  144. *> block reflectors, stored in array T(1:LDT,1:N).
  145. *> For more information see Further Details in GELQT.
  146. *>
  147. *> Q(i) for i>1 is computed by TPLQT, which represents Q(i) by Householder vectors
  148. *> stored in columns [(i-1)*(NB-M)+M+1:i*(NB-M)+M] of A, and by upper triangular
  149. *> block reflectors, stored in array T(1:LDT,(i-1)*M+1:i*M).
  150. *> The last Q(k) may use fewer rows.
  151. *> For more information see Further Details in TPQRT.
  152. *>
  153. *> For more details of the overall algorithm, see the description of
  154. *> Sequential TSQR in Section 2.2 of [1].
  155. *>
  156. *> [1] “Communication-Optimal Parallel and Sequential QR and LU Factorizations,”
  157. *> J. Demmel, L. Grigori, M. Hoemmen, J. Langou,
  158. *> SIAM J. Sci. Comput, vol. 34, no. 1, 2012
  159. *> \endverbatim
  160. *>
  161. * =====================================================================
  162. SUBROUTINE ZLASWLQ( M, N, MB, NB, A, LDA, T, LDT, WORK, LWORK,
  163. $ INFO)
  164. *
  165. * -- LAPACK computational routine (version 3.9.0) --
  166. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  167. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd. --
  168. * June 2017
  169. *
  170. * .. Scalar Arguments ..
  171. INTEGER INFO, LDA, M, N, MB, NB, LWORK, LDT
  172. * ..
  173. * .. Array Arguments ..
  174. COMPLEX*16 A( LDA, * ), WORK( * ), T( LDT, *)
  175. * ..
  176. *
  177. * =====================================================================
  178. *
  179. * ..
  180. * .. Local Scalars ..
  181. LOGICAL LQUERY
  182. INTEGER I, II, KK, CTR
  183. * ..
  184. * .. EXTERNAL FUNCTIONS ..
  185. LOGICAL LSAME
  186. EXTERNAL LSAME
  187. * .. EXTERNAL SUBROUTINES ..
  188. EXTERNAL ZGELQT, ZTPLQT, XERBLA
  189. * .. INTRINSIC FUNCTIONS ..
  190. INTRINSIC MAX, MIN, MOD
  191. * ..
  192. * .. EXECUTABLE STATEMENTS ..
  193. *
  194. * TEST THE INPUT ARGUMENTS
  195. *
  196. INFO = 0
  197. *
  198. LQUERY = ( LWORK.EQ.-1 )
  199. *
  200. IF( M.LT.0 ) THEN
  201. INFO = -1
  202. ELSE IF( N.LT.0 .OR. N.LT.M ) THEN
  203. INFO = -2
  204. ELSE IF( MB.LT.1 .OR. ( MB.GT.M .AND. M.GT.0 )) THEN
  205. INFO = -3
  206. ELSE IF( NB.LE.M ) THEN
  207. INFO = -4
  208. ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  209. INFO = -5
  210. ELSE IF( LDT.LT.MB ) THEN
  211. INFO = -8
  212. ELSE IF( ( LWORK.LT.M*MB) .AND. (.NOT.LQUERY) ) THEN
  213. INFO = -10
  214. END IF
  215. IF( INFO.EQ.0) THEN
  216. WORK(1) = MB*M
  217. END IF
  218. *
  219. IF( INFO.NE.0 ) THEN
  220. CALL XERBLA( 'ZLASWLQ', -INFO )
  221. RETURN
  222. ELSE IF (LQUERY) THEN
  223. RETURN
  224. END IF
  225. *
  226. * Quick return if possible
  227. *
  228. IF( MIN(M,N).EQ.0 ) THEN
  229. RETURN
  230. END IF
  231. *
  232. * The LQ Decomposition
  233. *
  234. IF((M.GE.N).OR.(NB.LE.M).OR.(NB.GE.N)) THEN
  235. CALL ZGELQT( M, N, MB, A, LDA, T, LDT, WORK, INFO)
  236. RETURN
  237. END IF
  238. *
  239. KK = MOD((N-M),(NB-M))
  240. II=N-KK+1
  241. *
  242. * Compute the LQ factorization of the first block A(1:M,1:NB)
  243. *
  244. CALL ZGELQT( M, NB, MB, A(1,1), LDA, T, LDT, WORK, INFO)
  245. CTR = 1
  246. *
  247. DO I = NB+1, II-NB+M , (NB-M)
  248. *
  249. * Compute the QR factorization of the current block A(1:M,I:I+NB-M)
  250. *
  251. CALL ZTPLQT( M, NB-M, 0, MB, A(1,1), LDA, A( 1, I ),
  252. $ LDA, T(1, CTR * M + 1),
  253. $ LDT, WORK, INFO )
  254. CTR = CTR + 1
  255. END DO
  256. *
  257. * Compute the QR factorization of the last block A(1:M,II:N)
  258. *
  259. IF (II.LE.N) THEN
  260. CALL ZTPLQT( M, KK, 0, MB, A(1,1), LDA, A( 1, II ),
  261. $ LDA, T(1, CTR * M + 1), LDT,
  262. $ WORK, INFO )
  263. END IF
  264. *
  265. WORK( 1 ) = M * MB
  266. RETURN
  267. *
  268. * End of ZLASWLQ
  269. *
  270. END