LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
cgetsqrhrt.f
Go to the documentation of this file.
1 *> \brief \b CGETSQRHRT
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CGETSQRHRT + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgetsqrhrt.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgetsqrhrt.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgetsqrhrt.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CGETSQRHRT( M, N, MB1, NB1, NB2, A, LDA, T, LDT, WORK,
22 * $ LWORK, INFO )
23 * IMPLICIT NONE
24 *
25 * .. Scalar Arguments ..
26 * INTEGER INFO, LDA, LDT, LWORK, M, N, NB1, NB2, MB1
27 * ..
28 * .. Array Arguments ..
29 * COMPLEX*16 A( LDA, * ), T( LDT, * ), WORK( * )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> CGETSQRHRT computes a NB2-sized column blocked QR-factorization
39 *> of a complex M-by-N matrix A with M >= N,
40 *>
41 *> A = Q * R.
42 *>
43 *> The routine uses internally a NB1-sized column blocked and MB1-sized
44 *> row blocked TSQR-factorization and perfors the reconstruction
45 *> of the Householder vectors from the TSQR output. The routine also
46 *> converts the R_tsqr factor from the TSQR-factorization output into
47 *> the R factor that corresponds to the Householder QR-factorization,
48 *>
49 *> A = Q_tsqr * R_tsqr = Q * R.
50 *>
51 *> The output Q and R factors are stored in the same format as in CGEQRT
52 *> (Q is in blocked compact WY-representation). See the documentation
53 *> of CGEQRT for more details on the format.
54 *> \endverbatim
55 *
56 * Arguments:
57 * ==========
58 *
59 *> \param[in] M
60 *> \verbatim
61 *> M is INTEGER
62 *> The number of rows of the matrix A. M >= 0.
63 *> \endverbatim
64 *>
65 *> \param[in] N
66 *> \verbatim
67 *> N is INTEGER
68 *> The number of columns of the matrix A. M >= N >= 0.
69 *> \endverbatim
70 *>
71 *> \param[in] MB1
72 *> \verbatim
73 *> MB1 is INTEGER
74 *> The row block size to be used in the blocked TSQR.
75 *> MB1 > N.
76 *> \endverbatim
77 *>
78 *> \param[in] NB1
79 *> \verbatim
80 *> NB1 is INTEGER
81 *> The column block size to be used in the blocked TSQR.
82 *> N >= NB1 >= 1.
83 *> \endverbatim
84 *>
85 *> \param[in] NB2
86 *> \verbatim
87 *> NB2 is INTEGER
88 *> The block size to be used in the blocked QR that is
89 *> output. NB2 >= 1.
90 *> \endverbatim
91 *>
92 *> \param[in,out] A
93 *> \verbatim
94 *> A is COMPLEX*16 array, dimension (LDA,N)
95 *>
96 *> On entry: an M-by-N matrix A.
97 *>
98 *> On exit:
99 *> a) the elements on and above the diagonal
100 *> of the array contain the N-by-N upper-triangular
101 *> matrix R corresponding to the Householder QR;
102 *> b) the elements below the diagonal represent Q by
103 *> the columns of blocked V (compact WY-representation).
104 *> \endverbatim
105 *>
106 *> \param[in] LDA
107 *> \verbatim
108 *> LDA is INTEGER
109 *> The leading dimension of the array A. LDA >= max(1,M).
110 *> \endverbatim
111 *>
112 *> \param[out] T
113 *> \verbatim
114 *> T is COMPLEX array, dimension (LDT,N))
115 *> The upper triangular block reflectors stored in compact form
116 *> as a sequence of upper triangular blocks.
117 *> \endverbatim
118 *>
119 *> \param[in] LDT
120 *> \verbatim
121 *> LDT is INTEGER
122 *> The leading dimension of the array T. LDT >= NB2.
123 *> \endverbatim
124 *>
125 *> \param[out] WORK
126 *> \verbatim
127 *> (workspace) COMPLEX array, dimension (MAX(1,LWORK))
128 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
129 *> \endverbatim
130 *>
131 *> \param[in] LWORK
132 *> \verbatim
133 *> The dimension of the array WORK.
134 *> LWORK >= MAX( LWT + LW1, MAX( LWT+N*N+LW2, LWT+N*N+N ) ),
135 *> where
136 *> NUM_ALL_ROW_BLOCKS = CEIL((M-N)/(MB1-N)),
137 *> NB1LOCAL = MIN(NB1,N).
138 *> LWT = NUM_ALL_ROW_BLOCKS * N * NB1LOCAL,
139 *> LW1 = NB1LOCAL * N,
140 *> LW2 = NB1LOCAL * MAX( NB1LOCAL, ( N - NB1LOCAL ) ),
141 *> If LWORK = -1, then a workspace query is assumed.
142 *> The routine only calculates the optimal size of the WORK
143 *> array, returns this value as the first entry of the WORK
144 *> array, and no error message related to LWORK is issued
145 *> by XERBLA.
146 *> \endverbatim
147 *>
148 *> \param[out] INFO
149 *> \verbatim
150 *> INFO is INTEGER
151 *> = 0: successful exit
152 *> < 0: if INFO = -i, the i-th argument had an illegal value
153 *> \endverbatim
154 *
155 * Authors:
156 * ========
157 *
158 *> \author Univ. of Tennessee
159 *> \author Univ. of California Berkeley
160 *> \author Univ. of Colorado Denver
161 *> \author NAG Ltd.
162 *
163 *> \ingroup comlpexOTHERcomputational
164 *
165 *> \par Contributors:
166 * ==================
167 *>
168 *> \verbatim
169 *>
170 *> November 2020, Igor Kozachenko,
171 *> Computer Science Division,
172 *> University of California, Berkeley
173 *>
174 *> \endverbatim
175 *>
176 * =====================================================================
177  SUBROUTINE cgetsqrhrt( M, N, MB1, NB1, NB2, A, LDA, T, LDT, WORK,
178  $ LWORK, INFO )
179  IMPLICIT NONE
180 *
181 * -- LAPACK computational routine --
182 * -- LAPACK is a software package provided by Univ. of Tennessee, --
183 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
184 *
185 * .. Scalar Arguments ..
186  INTEGER INFO, LDA, LDT, LWORK, M, N, NB1, NB2, MB1
187 * ..
188 * .. Array Arguments ..
189  COMPLEX A( LDA, * ), T( LDT, * ), WORK( * )
190 * ..
191 *
192 * =====================================================================
193 *
194 * .. Parameters ..
195  COMPLEX CONE
196  parameter( cone = ( 1.0e+0, 0.0e+0 ) )
197 * ..
198 * .. Local Scalars ..
199  LOGICAL LQUERY
200  INTEGER I, IINFO, J, LW1, LW2, LWT, LDWT, LWORKOPT,
201  $ nb1local, nb2local, num_all_row_blocks
202 * ..
203 * .. External Subroutines ..
204  EXTERNAL ccopy, clatsqr, cungtsqr_row, cunhr_col,
205  $ xerbla
206 * ..
207 * .. Intrinsic Functions ..
208  INTRINSIC ceiling, real, cmplx, max, min
209 * ..
210 * .. Executable Statements ..
211 *
212 * Test the input arguments
213 *
214  info = 0
215  lquery = lwork.EQ.-1
216  IF( m.LT.0 ) THEN
217  info = -1
218  ELSE IF( n.LT.0 .OR. m.LT.n ) THEN
219  info = -2
220  ELSE IF( mb1.LE.n ) THEN
221  info = -3
222  ELSE IF( nb1.LT.1 ) THEN
223  info = -4
224  ELSE IF( nb2.LT.1 ) THEN
225  info = -5
226  ELSE IF( lda.LT.max( 1, m ) ) THEN
227  info = -7
228  ELSE IF( ldt.LT.max( 1, min( nb2, n ) ) ) THEN
229  info = -9
230  ELSE
231 *
232 * Test the input LWORK for the dimension of the array WORK.
233 * This workspace is used to store array:
234 * a) Matrix T and WORK for CLATSQR;
235 * b) N-by-N upper-triangular factor R_tsqr;
236 * c) Matrix T and array WORK for CUNGTSQR_ROW;
237 * d) Diagonal D for CUNHR_COL.
238 *
239  IF( lwork.LT.n*n+1 .AND. .NOT.lquery ) THEN
240  info = -11
241  ELSE
242 *
243 * Set block size for column blocks
244 *
245  nb1local = min( nb1, n )
246 *
247  num_all_row_blocks = max( 1,
248  $ ceiling( real( m - n ) / real( mb1 - n ) ) )
249 *
250 * Length and leading dimension of WORK array to place
251 * T array in TSQR.
252 *
253  lwt = num_all_row_blocks * n * nb1local
254 
255  ldwt = nb1local
256 *
257 * Length of TSQR work array
258 *
259  lw1 = nb1local * n
260 *
261 * Length of CUNGTSQR_ROW work array.
262 *
263  lw2 = nb1local * max( nb1local, ( n - nb1local ) )
264 *
265  lworkopt = max( lwt + lw1, max( lwt+n*n+lw2, lwt+n*n+n ) )
266 *
267  IF( ( lwork.LT.max( 1, lworkopt ) ).AND.(.NOT.lquery) ) THEN
268  info = -11
269  END IF
270 *
271  END IF
272  END IF
273 *
274 * Handle error in the input parameters and return workspace query.
275 *
276  IF( info.NE.0 ) THEN
277  CALL xerbla( 'CGETSQRHRT', -info )
278  RETURN
279  ELSE IF ( lquery ) THEN
280  work( 1 ) = cmplx( lworkopt )
281  RETURN
282  END IF
283 *
284 * Quick return if possible
285 *
286  IF( min( m, n ).EQ.0 ) THEN
287  work( 1 ) = cmplx( lworkopt )
288  RETURN
289  END IF
290 *
291  nb2local = min( nb2, n )
292 *
293 *
294 * (1) Perform TSQR-factorization of the M-by-N matrix A.
295 *
296  CALL clatsqr( m, n, mb1, nb1local, a, lda, work, ldwt,
297  $ work(lwt+1), lw1, iinfo )
298 *
299 * (2) Copy the factor R_tsqr stored in the upper-triangular part
300 * of A into the square matrix in the work array
301 * WORK(LWT+1:LWT+N*N) column-by-column.
302 *
303  DO j = 1, n
304  CALL ccopy( j, a( 1, j ), 1, work( lwt + n*(j-1)+1 ), 1 )
305  END DO
306 *
307 * (3) Generate a M-by-N matrix Q with orthonormal columns from
308 * the result stored below the diagonal in the array A in place.
309 *
310 
311  CALL cungtsqr_row( m, n, mb1, nb1local, a, lda, work, ldwt,
312  $ work( lwt+n*n+1 ), lw2, iinfo )
313 *
314 * (4) Perform the reconstruction of Householder vectors from
315 * the matrix Q (stored in A) in place.
316 *
317  CALL cunhr_col( m, n, nb2local, a, lda, t, ldt,
318  $ work( lwt+n*n+1 ), iinfo )
319 *
320 * (5) Copy the factor R_tsqr stored in the square matrix in the
321 * work array WORK(LWT+1:LWT+N*N) into the upper-triangular
322 * part of A.
323 *
324 * (6) Compute from R_tsqr the factor R_hr corresponding to
325 * the reconstructed Householder vectors, i.e. R_hr = S * R_tsqr.
326 * This multiplication by the sign matrix S on the left means
327 * changing the sign of I-th row of the matrix R_tsqr according
328 * to sign of the I-th diagonal element DIAG(I) of the matrix S.
329 * DIAG is stored in WORK( LWT+N*N+1 ) from the CUNHR_COL output.
330 *
331 * (5) and (6) can be combined in a single loop, so the rows in A
332 * are accessed only once.
333 *
334  DO i = 1, n
335  IF( work( lwt+n*n+i ).EQ.-cone ) THEN
336  DO j = i, n
337  a( i, j ) = -cone * work( lwt+n*(j-1)+i )
338  END DO
339  ELSE
340  CALL ccopy( n-i+1, work(lwt+n*(i-1)+i), n, a( i, i ), lda )
341  END IF
342  END DO
343 *
344  work( 1 ) = cmplx( lworkopt )
345  RETURN
346 *
347 * End of CGETSQRHRT
348 *
349  END
subroutine cgetsqrhrt(M, N, MB1, NB1, NB2, A, LDA, T, LDT, WORK, LWORK, INFO)
CGETSQRHRT
Definition: cgetsqrhrt.f:179
subroutine clatsqr(M, N, MB, NB, A, LDA, T, LDT, WORK, LWORK, INFO)
CLATSQR
Definition: clatsqr.f:166
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine cunhr_col(M, N, NB, A, LDA, T, LDT, D, INFO)
CUNHR_COL
Definition: cunhr_col.f:259
subroutine cungtsqr_row(M, N, MB, NB, A, LDA, T, LDT, WORK, LWORK, INFO)
CUNGTSQR_ROW
Definition: cungtsqr_row.f:188