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.

sorgbr.f 9.8 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337
  1. *> \brief \b SORGBR
  2. *
  3. * =========== DOCUMENTATION ===========
  4. *
  5. * Online html documentation available at
  6. * http://www.netlib.org/lapack/explore-html/
  7. *
  8. *> \htmlonly
  9. *> Download SORGBR + dependencies
  10. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sorgbr.f">
  11. *> [TGZ]</a>
  12. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sorgbr.f">
  13. *> [ZIP]</a>
  14. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sorgbr.f">
  15. *> [TXT]</a>
  16. *> \endhtmlonly
  17. *
  18. * Definition:
  19. * ===========
  20. *
  21. * SUBROUTINE SORGBR( VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO )
  22. *
  23. * .. Scalar Arguments ..
  24. * CHARACTER VECT
  25. * INTEGER INFO, K, LDA, LWORK, M, N
  26. * ..
  27. * .. Array Arguments ..
  28. * REAL A( LDA, * ), TAU( * ), WORK( * )
  29. * ..
  30. *
  31. *
  32. *> \par Purpose:
  33. * =============
  34. *>
  35. *> \verbatim
  36. *>
  37. *> SORGBR generates one of the real orthogonal matrices Q or P**T
  38. *> determined by SGEBRD when reducing a real matrix A to bidiagonal
  39. *> form: A = Q * B * P**T. Q and P**T are defined as products of
  40. *> elementary reflectors H(i) or G(i) respectively.
  41. *>
  42. *> If VECT = 'Q', A is assumed to have been an M-by-K matrix, and Q
  43. *> is of order M:
  44. *> if m >= k, Q = H(1) H(2) . . . H(k) and SORGBR returns the first n
  45. *> columns of Q, where m >= n >= k;
  46. *> if m < k, Q = H(1) H(2) . . . H(m-1) and SORGBR returns Q as an
  47. *> M-by-M matrix.
  48. *>
  49. *> If VECT = 'P', A is assumed to have been a K-by-N matrix, and P**T
  50. *> is of order N:
  51. *> if k < n, P**T = G(k) . . . G(2) G(1) and SORGBR returns the first m
  52. *> rows of P**T, where n >= m >= k;
  53. *> if k >= n, P**T = G(n-1) . . . G(2) G(1) and SORGBR returns P**T as
  54. *> an N-by-N matrix.
  55. *> \endverbatim
  56. *
  57. * Arguments:
  58. * ==========
  59. *
  60. *> \param[in] VECT
  61. *> \verbatim
  62. *> VECT is CHARACTER*1
  63. *> Specifies whether the matrix Q or the matrix P**T is
  64. *> required, as defined in the transformation applied by SGEBRD:
  65. *> = 'Q': generate Q;
  66. *> = 'P': generate P**T.
  67. *> \endverbatim
  68. *>
  69. *> \param[in] M
  70. *> \verbatim
  71. *> M is INTEGER
  72. *> The number of rows of the matrix Q or P**T to be returned.
  73. *> M >= 0.
  74. *> \endverbatim
  75. *>
  76. *> \param[in] N
  77. *> \verbatim
  78. *> N is INTEGER
  79. *> The number of columns of the matrix Q or P**T to be returned.
  80. *> N >= 0.
  81. *> If VECT = 'Q', M >= N >= min(M,K);
  82. *> if VECT = 'P', N >= M >= min(N,K).
  83. *> \endverbatim
  84. *>
  85. *> \param[in] K
  86. *> \verbatim
  87. *> K is INTEGER
  88. *> If VECT = 'Q', the number of columns in the original M-by-K
  89. *> matrix reduced by SGEBRD.
  90. *> If VECT = 'P', the number of rows in the original K-by-N
  91. *> matrix reduced by SGEBRD.
  92. *> K >= 0.
  93. *> \endverbatim
  94. *>
  95. *> \param[in,out] A
  96. *> \verbatim
  97. *> A is REAL array, dimension (LDA,N)
  98. *> On entry, the vectors which define the elementary reflectors,
  99. *> as returned by SGEBRD.
  100. *> On exit, the M-by-N matrix Q or P**T.
  101. *> \endverbatim
  102. *>
  103. *> \param[in] LDA
  104. *> \verbatim
  105. *> LDA is INTEGER
  106. *> The leading dimension of the array A. LDA >= max(1,M).
  107. *> \endverbatim
  108. *>
  109. *> \param[in] TAU
  110. *> \verbatim
  111. *> TAU is REAL array, dimension
  112. *> (min(M,K)) if VECT = 'Q'
  113. *> (min(N,K)) if VECT = 'P'
  114. *> TAU(i) must contain the scalar factor of the elementary
  115. *> reflector H(i) or G(i), which determines Q or P**T, as
  116. *> returned by SGEBRD in its array argument TAUQ or TAUP.
  117. *> \endverbatim
  118. *>
  119. *> \param[out] WORK
  120. *> \verbatim
  121. *> WORK is REAL array, dimension (MAX(1,LWORK))
  122. *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  123. *> \endverbatim
  124. *>
  125. *> \param[in] LWORK
  126. *> \verbatim
  127. *> LWORK is INTEGER
  128. *> The dimension of the array WORK. LWORK >= max(1,min(M,N)).
  129. *> For optimum performance LWORK >= min(M,N)*NB, where NB
  130. *> is the optimal blocksize.
  131. *>
  132. *> If LWORK = -1, then a workspace query is assumed; the routine
  133. *> only calculates the optimal size of the WORK array, returns
  134. *> this value as the first entry of the WORK array, and no error
  135. *> message related to LWORK is issued by XERBLA.
  136. *> \endverbatim
  137. *>
  138. *> \param[out] INFO
  139. *> \verbatim
  140. *> INFO is INTEGER
  141. *> = 0: successful exit
  142. *> < 0: if INFO = -i, the i-th argument had an illegal value
  143. *> \endverbatim
  144. *
  145. * Authors:
  146. * ========
  147. *
  148. *> \author Univ. of Tennessee
  149. *> \author Univ. of California Berkeley
  150. *> \author Univ. of Colorado Denver
  151. *> \author NAG Ltd.
  152. *
  153. *> \date April 2012
  154. *
  155. *> \ingroup realGBcomputational
  156. *
  157. * =====================================================================
  158. SUBROUTINE SORGBR( VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO )
  159. *
  160. * -- LAPACK computational routine (version 3.7.0) --
  161. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  162. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  163. * April 2012
  164. *
  165. * .. Scalar Arguments ..
  166. CHARACTER VECT
  167. INTEGER INFO, K, LDA, LWORK, M, N
  168. * ..
  169. * .. Array Arguments ..
  170. REAL A( LDA, * ), TAU( * ), WORK( * )
  171. * ..
  172. *
  173. * =====================================================================
  174. *
  175. * .. Parameters ..
  176. REAL ZERO, ONE
  177. PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
  178. * ..
  179. * .. Local Scalars ..
  180. LOGICAL LQUERY, WANTQ
  181. INTEGER I, IINFO, J, LWKOPT, MN
  182. * ..
  183. * .. External Functions ..
  184. LOGICAL LSAME
  185. EXTERNAL LSAME
  186. * ..
  187. * .. External Subroutines ..
  188. EXTERNAL SORGLQ, SORGQR, XERBLA
  189. * ..
  190. * .. Intrinsic Functions ..
  191. INTRINSIC MAX, MIN
  192. * ..
  193. * .. Executable Statements ..
  194. *
  195. * Test the input arguments
  196. *
  197. INFO = 0
  198. WANTQ = LSAME( VECT, 'Q' )
  199. MN = MIN( M, N )
  200. LQUERY = ( LWORK.EQ.-1 )
  201. IF( .NOT.WANTQ .AND. .NOT.LSAME( VECT, 'P' ) ) THEN
  202. INFO = -1
  203. ELSE IF( M.LT.0 ) THEN
  204. INFO = -2
  205. ELSE IF( N.LT.0 .OR. ( WANTQ .AND. ( N.GT.M .OR. N.LT.MIN( M,
  206. $ K ) ) ) .OR. ( .NOT.WANTQ .AND. ( M.GT.N .OR. M.LT.
  207. $ MIN( N, K ) ) ) ) THEN
  208. INFO = -3
  209. ELSE IF( K.LT.0 ) THEN
  210. INFO = -4
  211. ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  212. INFO = -6
  213. ELSE IF( LWORK.LT.MAX( 1, MN ) .AND. .NOT.LQUERY ) THEN
  214. INFO = -9
  215. END IF
  216. *
  217. IF( INFO.EQ.0 ) THEN
  218. WORK( 1 ) = 1
  219. IF( WANTQ ) THEN
  220. IF( M.GE.K ) THEN
  221. CALL SORGQR( M, N, K, A, LDA, TAU, WORK, -1, IINFO )
  222. ELSE
  223. IF( M.GT.1 ) THEN
  224. CALL SORGQR( M-1, M-1, M-1, A( 2, 2 ), LDA, TAU, WORK,
  225. $ -1, IINFO )
  226. END IF
  227. END IF
  228. ELSE
  229. IF( K.LT.N ) THEN
  230. CALL SORGLQ( M, N, K, A, LDA, TAU, WORK, -1, IINFO )
  231. ELSE
  232. IF( N.GT.1 ) THEN
  233. CALL SORGLQ( N-1, N-1, N-1, A( 2, 2 ), LDA, TAU, WORK,
  234. $ -1, IINFO )
  235. END IF
  236. END IF
  237. END IF
  238. LWKOPT = WORK( 1 )
  239. LWKOPT = MAX (LWKOPT, MN)
  240. END IF
  241. *
  242. IF( INFO.NE.0 ) THEN
  243. CALL XERBLA( 'SORGBR', -INFO )
  244. RETURN
  245. ELSE IF( LQUERY ) THEN
  246. WORK( 1 ) = LWKOPT
  247. RETURN
  248. END IF
  249. *
  250. * Quick return if possible
  251. *
  252. IF( M.EQ.0 .OR. N.EQ.0 ) THEN
  253. WORK( 1 ) = 1
  254. RETURN
  255. END IF
  256. *
  257. IF( WANTQ ) THEN
  258. *
  259. * Form Q, determined by a call to SGEBRD to reduce an m-by-k
  260. * matrix
  261. *
  262. IF( M.GE.K ) THEN
  263. *
  264. * If m >= k, assume m >= n >= k
  265. *
  266. CALL SORGQR( M, N, K, A, LDA, TAU, WORK, LWORK, IINFO )
  267. *
  268. ELSE
  269. *
  270. * If m < k, assume m = n
  271. *
  272. * Shift the vectors which define the elementary reflectors one
  273. * column to the right, and set the first row and column of Q
  274. * to those of the unit matrix
  275. *
  276. DO 20 J = M, 2, -1
  277. A( 1, J ) = ZERO
  278. DO 10 I = J + 1, M
  279. A( I, J ) = A( I, J-1 )
  280. 10 CONTINUE
  281. 20 CONTINUE
  282. A( 1, 1 ) = ONE
  283. DO 30 I = 2, M
  284. A( I, 1 ) = ZERO
  285. 30 CONTINUE
  286. IF( M.GT.1 ) THEN
  287. *
  288. * Form Q(2:m,2:m)
  289. *
  290. CALL SORGQR( M-1, M-1, M-1, A( 2, 2 ), LDA, TAU, WORK,
  291. $ LWORK, IINFO )
  292. END IF
  293. END IF
  294. ELSE
  295. *
  296. * Form P**T, determined by a call to SGEBRD to reduce a k-by-n
  297. * matrix
  298. *
  299. IF( K.LT.N ) THEN
  300. *
  301. * If k < n, assume k <= m <= n
  302. *
  303. CALL SORGLQ( M, N, K, A, LDA, TAU, WORK, LWORK, IINFO )
  304. *
  305. ELSE
  306. *
  307. * If k >= n, assume m = n
  308. *
  309. * Shift the vectors which define the elementary reflectors one
  310. * row downward, and set the first row and column of P**T to
  311. * those of the unit matrix
  312. *
  313. A( 1, 1 ) = ONE
  314. DO 40 I = 2, N
  315. A( I, 1 ) = ZERO
  316. 40 CONTINUE
  317. DO 60 J = 2, N
  318. DO 50 I = J - 1, 2, -1
  319. A( I, J ) = A( I-1, J )
  320. 50 CONTINUE
  321. A( 1, J ) = ZERO
  322. 60 CONTINUE
  323. IF( N.GT.1 ) THEN
  324. *
  325. * Form P**T(2:n,2:n)
  326. *
  327. CALL SORGLQ( N-1, N-1, N-1, A( 2, 2 ), LDA, TAU, WORK,
  328. $ LWORK, IINFO )
  329. END IF
  330. END IF
  331. END IF
  332. WORK( 1 ) = LWKOPT
  333. RETURN
  334. *
  335. * End of SORGBR
  336. *
  337. END