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.

sgelqt.f 5.2 kB

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