LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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 getsqrhrt
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 ..
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 xerbla(srname, info)
Definition cblat2.f:3285
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
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:169
subroutine cungtsqr_row(m, n, mb, nb, a, lda, t, ldt, work, lwork, info)
CUNGTSQR_ROW
subroutine cunhr_col(m, n, nb, a, lda, t, ldt, d, info)
CUNHR_COL
Definition cunhr_col.f:259