LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
sgetsqrhrt.f
Go to the documentation of this file.
1*> \brief \b SGETSQRHRT
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download SGETSQRHRT + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgetsqrhrt.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgetsqrhrt.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgetsqrhrt.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SGETSQRHRT( 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* REAL A( LDA, * ), T( LDT, * ), WORK( * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> SGETSQRHRT 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 SGEQRT
52*> (Q is in blocked compact WY-representation). See the documentation
53*> of SGEQRT 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 REAL 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 REAL 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) REAL 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 sgetsqrhrt( 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 REAL A( LDA, * ), T( LDT, * ), WORK( * )
190* ..
191*
192* =====================================================================
193*
194* .. Parameters ..
195 REAL ONE
196 parameter( one = 1.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 Functions ..
204 REAL SROUNDUP_LWORK
205 EXTERNAL sroundup_lwork
206* ..
207* .. External Subroutines ..
209 $ xerbla
210* ..
211* .. Intrinsic Functions ..
212 INTRINSIC ceiling, max, min
213* ..
214* .. Executable Statements ..
215*
216* Test the input arguments
217*
218 info = 0
219 lquery = lwork.EQ.-1
220 IF( m.LT.0 ) THEN
221 info = -1
222 ELSE IF( n.LT.0 .OR. m.LT.n ) THEN
223 info = -2
224 ELSE IF( mb1.LE.n ) THEN
225 info = -3
226 ELSE IF( nb1.LT.1 ) THEN
227 info = -4
228 ELSE IF( nb2.LT.1 ) THEN
229 info = -5
230 ELSE IF( lda.LT.max( 1, m ) ) THEN
231 info = -7
232 ELSE IF( ldt.LT.max( 1, min( nb2, n ) ) ) THEN
233 info = -9
234 ELSE
235*
236* Test the input LWORK for the dimension of the array WORK.
237* This workspace is used to store array:
238* a) Matrix T and WORK for SLATSQR;
239* b) N-by-N upper-triangular factor R_tsqr;
240* c) Matrix T and array WORK for SORGTSQR_ROW;
241* d) Diagonal D for SORHR_COL.
242*
243 IF( lwork.LT.n*n+1 .AND. .NOT.lquery ) THEN
244 info = -11
245 ELSE
246*
247* Set block size for column blocks
248*
249 nb1local = min( nb1, n )
250*
251 num_all_row_blocks = max( 1,
252 $ ceiling( real( m - n ) / real( mb1 - n ) ) )
253*
254* Length and leading dimension of WORK array to place
255* T array in TSQR.
256*
257 lwt = num_all_row_blocks * n * nb1local
258
259 ldwt = nb1local
260*
261* Length of TSQR work array
262*
263 lw1 = nb1local * n
264*
265* Length of SORGTSQR_ROW work array.
266*
267 lw2 = nb1local * max( nb1local, ( n - nb1local ) )
268*
269 lworkopt = max( lwt + lw1, max( lwt+n*n+lw2, lwt+n*n+n ) )
270*
271 IF( ( lwork.LT.max( 1, 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( 'SGETSQRHRT', -info )
282 RETURN
283 ELSE IF ( lquery ) THEN
284 work( 1 ) = sroundup_lwork( lworkopt )
285 RETURN
286 END IF
287*
288* Quick return if possible
289*
290 IF( min( m, n ).EQ.0 ) THEN
291 work( 1 ) = sroundup_lwork( 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 slatsqr( 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 scopy( 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 sorgtsqr_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 sorhr_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 SORHR_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 scopy( n-i+1, work(lwt+n*(i-1)+i), n, a( i, i ), lda )
345 END IF
346 END DO
347*
348 work( 1 ) = sroundup_lwork( lworkopt )
349 RETURN
350*
351* End of SGETSQRHRT
352*
353 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
subroutine sgetsqrhrt(m, n, mb1, nb1, nb2, a, lda, t, ldt, work, lwork, info)
SGETSQRHRT
Definition sgetsqrhrt.f:179
subroutine slatsqr(m, n, mb, nb, a, lda, t, ldt, work, lwork, info)
SLATSQR
Definition slatsqr.f:169
subroutine sorgtsqr_row(m, n, mb, nb, a, lda, t, ldt, work, lwork, info)
SORGTSQR_ROW
subroutine sorhr_col(m, n, nb, a, lda, t, ldt, d, info)
SORHR_COL
Definition sorhr_col.f:259