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.

cgelqt.f 5.2 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194
  1. *
  2. * Definition:
  3. * ===========
  4. *
  5. * SUBROUTINE CGELQT( M, N, MB, A, LDA, T, LDT, WORK, INFO )
  6. *
  7. * .. Scalar Arguments ..
  8. * INTEGER INFO, LDA, LDT, M, N, MB
  9. * ..
  10. * .. Array Arguments ..
  11. * COMPLEX A( LDA, * ), T( LDT, * ), WORK( * )
  12. * ..
  13. *
  14. *
  15. *> \par Purpose:
  16. * =============
  17. *>
  18. *> \verbatim
  19. *>
  20. *> CGELQT computes a blocked LQ factorization of a complex M-by-N matrix A
  21. *> using the compact WY representation of Q.
  22. *> \endverbatim
  23. *
  24. * Arguments:
  25. * ==========
  26. *
  27. *> \param[in] M
  28. *> \verbatim
  29. *> M is INTEGER
  30. *> The number of rows of the matrix A. M >= 0.
  31. *> \endverbatim
  32. *>
  33. *> \param[in] N
  34. *> \verbatim
  35. *> N is INTEGER
  36. *> The number of columns of the matrix A. N >= 0.
  37. *> \endverbatim
  38. *>
  39. *> \param[in] MB
  40. *> \verbatim
  41. *> MB is INTEGER
  42. *> The block size to be used in the blocked QR. MIN(M,N) >= MB >= 1.
  43. *> \endverbatim
  44. *>
  45. *> \param[in,out] A
  46. *> \verbatim
  47. *> A is COMPLEX array, dimension (LDA,N)
  48. *> On entry, the M-by-N matrix A.
  49. *> On exit, the elements on and below the diagonal of the array
  50. *> contain the M-by-MIN(M,N) lower trapezoidal matrix L (L is
  51. *> lower triangular if M <= N); the elements above the diagonal
  52. *> are the rows of V.
  53. *> \endverbatim
  54. *>
  55. *> \param[in] LDA
  56. *> \verbatim
  57. *> LDA is INTEGER
  58. *> The leading dimension of the array A. LDA >= max(1,M).
  59. *> \endverbatim
  60. *>
  61. *> \param[out] T
  62. *> \verbatim
  63. *> T is COMPLEX array, dimension (LDT,MIN(M,N))
  64. *> The upper triangular block reflectors stored in compact form
  65. *> as a sequence of upper triangular blocks. See below
  66. *> for further details.
  67. *> \endverbatim
  68. *>
  69. *> \param[in] LDT
  70. *> \verbatim
  71. *> LDT is INTEGER
  72. *> The leading dimension of the array T. LDT >= MB.
  73. *> \endverbatim
  74. *>
  75. *> \param[out] WORK
  76. *> \verbatim
  77. *> WORK is COMPLEX array, dimension (MB*N)
  78. *> \endverbatim
  79. *>
  80. *> \param[out] INFO
  81. *> \verbatim
  82. *> INFO is INTEGER
  83. *> = 0: successful exit
  84. *> < 0: if INFO = -i, the i-th argument had an illegal value
  85. *> \endverbatim
  86. *
  87. * Authors:
  88. * ========
  89. *
  90. *> \author Univ. of Tennessee
  91. *> \author Univ. of California Berkeley
  92. *> \author Univ. of Colorado Denver
  93. *> \author NAG Ltd.
  94. *
  95. *> \date December 2016
  96. *
  97. *> \ingroup doubleGEcomputational
  98. *
  99. *> \par Further Details:
  100. * =====================
  101. *>
  102. *> \verbatim
  103. *>
  104. *> The matrix V stores the elementary reflectors H(i) in the i-th column
  105. *> below the diagonal. For example, if M=5 and N=3, the matrix V is
  106. *>
  107. *> V = ( 1 v1 v1 v1 v1 )
  108. *> ( 1 v2 v2 v2 )
  109. *> ( 1 v3 v3 )
  110. *>
  111. *>
  112. *> where the vi's represent the vectors which define H(i), which are returned
  113. *> in the matrix A. The 1's along the diagonal of V are not stored in A.
  114. *> Let K=MIN(M,N). The number of blocks is B = ceiling(K/NB), where each
  115. *> block is of order NB except for the last block, which is of order
  116. *> IB = K - (B-1)*NB. For each of the B blocks, a upper triangular block
  117. *> reflector factor is computed: T1, T2, ..., TB. The NB-by-NB (and IB-by-IB
  118. *> for the last block) T's are stored in the NB-by-N matrix T as
  119. *>
  120. *> T = (T1 T2 ... TB).
  121. *> \endverbatim
  122. *>
  123. * =====================================================================
  124. SUBROUTINE CGELQT( M, N, MB, A, LDA, T, LDT, WORK, INFO )
  125. *
  126. * -- LAPACK computational routine (version 3.7.0) --
  127. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  128. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  129. * December 2016
  130. *
  131. * .. Scalar Arguments ..
  132. INTEGER INFO, LDA, LDT, M, N, MB
  133. * ..
  134. * .. Array Arguments ..
  135. COMPLEX A( LDA, * ), T( LDT, * ), WORK( * )
  136. * ..
  137. *
  138. * =====================================================================
  139. *
  140. * ..
  141. * .. Local Scalars ..
  142. INTEGER I, IB, IINFO, K
  143. * ..
  144. * .. External Subroutines ..
  145. EXTERNAL CGELQT3, CLARFB, XERBLA
  146. * ..
  147. * .. Executable Statements ..
  148. *
  149. * Test the input arguments
  150. *
  151. INFO = 0
  152. IF( M.LT.0 ) THEN
  153. INFO = -1
  154. ELSE IF( N.LT.0 ) THEN
  155. INFO = -2
  156. ELSE IF( MB.LT.1 .OR. (MB.GT.MIN(M,N) .AND. MIN(M,N).GT.0 ))THEN
  157. INFO = -3
  158. ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  159. INFO = -5
  160. ELSE IF( LDT.LT.MB ) THEN
  161. INFO = -7
  162. END IF
  163. IF( INFO.NE.0 ) THEN
  164. CALL XERBLA( 'CGELQT', -INFO )
  165. RETURN
  166. END IF
  167. *
  168. * Quick return if possible
  169. *
  170. K = MIN( M, N )
  171. IF( K.EQ.0 ) RETURN
  172. *
  173. * Blocked loop of length K
  174. *
  175. DO I = 1, K, MB
  176. IB = MIN( K-I+1, MB )
  177. *
  178. * Compute the LQ factorization of the current block A(I:M,I:I+IB-1)
  179. *
  180. CALL CGELQT3( IB, N-I+1, A(I,I), LDA, T(1,I), LDT, IINFO )
  181. IF( I+IB.LE.M ) THEN
  182. *
  183. * Update by applying H**T to A(I:M,I+IB:N) from the right
  184. *
  185. CALL CLARFB( 'R', 'N', 'F', 'R', M-I-IB+1, N-I+1, IB,
  186. $ A( I, I ), LDA, T( 1, I ), LDT,
  187. $ A( I+IB, I ), LDA, WORK , M-I-IB+1 )
  188. END IF
  189. END DO
  190. RETURN
  191. *
  192. * End of CGELQT
  193. *
  194. END