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