LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
cgeqp3.f
Go to the documentation of this file.
1*> \brief \b CGEQP3
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CGEQP3 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgeqp3.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgeqp3.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgeqp3.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CGEQP3( M, N, A, LDA, JPVT, TAU, WORK, LWORK, RWORK,
22* INFO )
23*
24* .. Scalar Arguments ..
25* INTEGER INFO, LDA, LWORK, M, N
26* ..
27* .. Array Arguments ..
28* INTEGER JPVT( * )
29* REAL RWORK( * )
30* COMPLEX A( LDA, * ), TAU( * ), WORK( * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> CGEQP3 computes a QR factorization with column pivoting of a
40*> matrix A: A*P = Q*R using Level 3 BLAS.
41*> \endverbatim
42*
43* Arguments:
44* ==========
45*
46*> \param[in] M
47*> \verbatim
48*> M is INTEGER
49*> The number of rows of the matrix A. M >= 0.
50*> \endverbatim
51*>
52*> \param[in] N
53*> \verbatim
54*> N is INTEGER
55*> The number of columns of the matrix A. N >= 0.
56*> \endverbatim
57*>
58*> \param[in,out] A
59*> \verbatim
60*> A is COMPLEX array, dimension (LDA,N)
61*> On entry, the M-by-N matrix A.
62*> On exit, the upper triangle of the array contains the
63*> min(M,N)-by-N upper trapezoidal matrix R; the elements below
64*> the diagonal, together with the array TAU, represent the
65*> unitary matrix Q as a product of min(M,N) elementary
66*> reflectors.
67*> \endverbatim
68*>
69*> \param[in] LDA
70*> \verbatim
71*> LDA is INTEGER
72*> The leading dimension of the array A. LDA >= max(1,M).
73*> \endverbatim
74*>
75*> \param[in,out] JPVT
76*> \verbatim
77*> JPVT is INTEGER array, dimension (N)
78*> On entry, if JPVT(J).ne.0, the J-th column of A is permuted
79*> to the front of A*P (a leading column); if JPVT(J)=0,
80*> the J-th column of A is a free column.
81*> On exit, if JPVT(J)=K, then the J-th column of A*P was the
82*> the K-th column of A.
83*> \endverbatim
84*>
85*> \param[out] TAU
86*> \verbatim
87*> TAU is COMPLEX array, dimension (min(M,N))
88*> The scalar factors of the elementary reflectors.
89*> \endverbatim
90*>
91*> \param[out] WORK
92*> \verbatim
93*> WORK is COMPLEX array, dimension (MAX(1,LWORK))
94*> On exit, if INFO=0, WORK(1) returns the optimal LWORK.
95*> \endverbatim
96*>
97*> \param[in] LWORK
98*> \verbatim
99*> LWORK is INTEGER
100*> The dimension of the array WORK. LWORK >= N+1.
101*> For optimal performance LWORK >= ( N+1 )*NB, where NB
102*> is the optimal blocksize.
103*>
104*> If LWORK = -1, then a workspace query is assumed; the routine
105*> only calculates the optimal size of the WORK array, returns
106*> this value as the first entry of the WORK array, and no error
107*> message related to LWORK is issued by XERBLA.
108*> \endverbatim
109*>
110*> \param[out] RWORK
111*> \verbatim
112*> RWORK is REAL array, dimension (2*N)
113*> \endverbatim
114*>
115*> \param[out] INFO
116*> \verbatim
117*> INFO is INTEGER
118*> = 0: successful exit.
119*> < 0: if INFO = -i, the i-th argument had an illegal value.
120*> \endverbatim
121*
122* Authors:
123* ========
124*
125*> \author Univ. of Tennessee
126*> \author Univ. of California Berkeley
127*> \author Univ. of Colorado Denver
128*> \author NAG Ltd.
129*
130*> \ingroup complexGEcomputational
131*
132*> \par Further Details:
133* =====================
134*>
135*> \verbatim
136*>
137*> The matrix Q is represented as a product of elementary reflectors
138*>
139*> Q = H(1) H(2) . . . H(k), where k = min(m,n).
140*>
141*> Each H(i) has the form
142*>
143*> H(i) = I - tau * v * v**H
144*>
145*> where tau is a complex scalar, and v is a real/complex vector
146*> with v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in
147*> A(i+1:m,i), and tau in TAU(i).
148*> \endverbatim
149*
150*> \par Contributors:
151* ==================
152*>
153*> G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
154*> X. Sun, Computer Science Dept., Duke University, USA
155*>
156* =====================================================================
157 SUBROUTINE cgeqp3( M, N, A, LDA, JPVT, TAU, WORK, LWORK, RWORK,
158 $ INFO )
159*
160* -- LAPACK computational routine --
161* -- LAPACK is a software package provided by Univ. of Tennessee, --
162* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
163*
164* .. Scalar Arguments ..
165 INTEGER INFO, LDA, LWORK, M, N
166* ..
167* .. Array Arguments ..
168 INTEGER JPVT( * )
169 REAL RWORK( * )
170 COMPLEX A( LDA, * ), TAU( * ), WORK( * )
171* ..
172*
173* =====================================================================
174*
175* .. Parameters ..
176 INTEGER INB, INBMIN, IXOVER
177 parameter( inb = 1, inbmin = 2, ixover = 3 )
178* ..
179* .. Local Scalars ..
180 LOGICAL LQUERY
181 INTEGER FJB, IWS, J, JB, LWKOPT, MINMN, MINWS, NA, NB,
182 $ nbmin, nfxd, nx, sm, sminmn, sn, topbmn
183* ..
184* .. External Subroutines ..
185 EXTERNAL cgeqrf, claqp2, claqps, cswap, cunmqr, xerbla
186* ..
187* .. External Functions ..
188 INTEGER ILAENV
189 REAL SCNRM2
190 EXTERNAL ilaenv, scnrm2
191* ..
192* .. Intrinsic Functions ..
193 INTRINSIC int, max, min
194* ..
195* .. Executable Statements ..
196*
197* Test input arguments
198* ====================
199*
200 info = 0
201 lquery = ( lwork.EQ.-1 )
202 IF( m.LT.0 ) THEN
203 info = -1
204 ELSE IF( n.LT.0 ) THEN
205 info = -2
206 ELSE IF( lda.LT.max( 1, m ) ) THEN
207 info = -4
208 END IF
209*
210 IF( info.EQ.0 ) THEN
211 minmn = min( m, n )
212 IF( minmn.EQ.0 ) THEN
213 iws = 1
214 lwkopt = 1
215 ELSE
216 iws = n + 1
217 nb = ilaenv( inb, 'CGEQRF', ' ', m, n, -1, -1 )
218 lwkopt = ( n + 1 )*nb
219 END IF
220 work( 1 ) = cmplx( lwkopt )
221*
222 IF( ( lwork.LT.iws ) .AND. .NOT.lquery ) THEN
223 info = -8
224 END IF
225 END IF
226*
227 IF( info.NE.0 ) THEN
228 CALL xerbla( 'CGEQP3', -info )
229 RETURN
230 ELSE IF( lquery ) THEN
231 RETURN
232 END IF
233*
234* Move initial columns up front.
235*
236 nfxd = 1
237 DO 10 j = 1, n
238 IF( jpvt( j ).NE.0 ) THEN
239 IF( j.NE.nfxd ) THEN
240 CALL cswap( m, a( 1, j ), 1, a( 1, nfxd ), 1 )
241 jpvt( j ) = jpvt( nfxd )
242 jpvt( nfxd ) = j
243 ELSE
244 jpvt( j ) = j
245 END IF
246 nfxd = nfxd + 1
247 ELSE
248 jpvt( j ) = j
249 END IF
250 10 CONTINUE
251 nfxd = nfxd - 1
252*
253* Factorize fixed columns
254* =======================
255*
256* Compute the QR factorization of fixed columns and update
257* remaining columns.
258*
259 IF( nfxd.GT.0 ) THEN
260 na = min( m, nfxd )
261*CC CALL CGEQR2( M, NA, A, LDA, TAU, WORK, INFO )
262 CALL cgeqrf( m, na, a, lda, tau, work, lwork, info )
263 iws = max( iws, int( work( 1 ) ) )
264 IF( na.LT.n ) THEN
265*CC CALL CUNM2R( 'Left', 'Conjugate Transpose', M, N-NA,
266*CC $ NA, A, LDA, TAU, A( 1, NA+1 ), LDA, WORK,
267*CC $ INFO )
268 CALL cunmqr( 'Left', 'Conjugate Transpose', m, n-na, na, a,
269 $ lda, tau, a( 1, na+1 ), lda, work, lwork,
270 $ info )
271 iws = max( iws, int( work( 1 ) ) )
272 END IF
273 END IF
274*
275* Factorize free columns
276* ======================
277*
278 IF( nfxd.LT.minmn ) THEN
279*
280 sm = m - nfxd
281 sn = n - nfxd
282 sminmn = minmn - nfxd
283*
284* Determine the block size.
285*
286 nb = ilaenv( inb, 'CGEQRF', ' ', sm, sn, -1, -1 )
287 nbmin = 2
288 nx = 0
289*
290 IF( ( nb.GT.1 ) .AND. ( nb.LT.sminmn ) ) THEN
291*
292* Determine when to cross over from blocked to unblocked code.
293*
294 nx = max( 0, ilaenv( ixover, 'CGEQRF', ' ', sm, sn, -1,
295 $ -1 ) )
296*
297*
298 IF( nx.LT.sminmn ) THEN
299*
300* Determine if workspace is large enough for blocked code.
301*
302 minws = ( sn+1 )*nb
303 iws = max( iws, minws )
304 IF( lwork.LT.minws ) THEN
305*
306* Not enough workspace to use optimal NB: Reduce NB and
307* determine the minimum value of NB.
308*
309 nb = lwork / ( sn+1 )
310 nbmin = max( 2, ilaenv( inbmin, 'CGEQRF', ' ', sm, sn,
311 $ -1, -1 ) )
312*
313*
314 END IF
315 END IF
316 END IF
317*
318* Initialize partial column norms. The first N elements of work
319* store the exact column norms.
320*
321 DO 20 j = nfxd + 1, n
322 rwork( j ) = scnrm2( sm, a( nfxd+1, j ), 1 )
323 rwork( n+j ) = rwork( j )
324 20 CONTINUE
325*
326 IF( ( nb.GE.nbmin ) .AND. ( nb.LT.sminmn ) .AND.
327 $ ( nx.LT.sminmn ) ) THEN
328*
329* Use blocked code initially.
330*
331 j = nfxd + 1
332*
333* Compute factorization: while loop.
334*
335*
336 topbmn = minmn - nx
337 30 CONTINUE
338 IF( j.LE.topbmn ) THEN
339 jb = min( nb, topbmn-j+1 )
340*
341* Factorize JB columns among columns J:N.
342*
343 CALL claqps( m, n-j+1, j-1, jb, fjb, a( 1, j ), lda,
344 $ jpvt( j ), tau( j ), rwork( j ),
345 $ rwork( n+j ), work( 1 ), work( jb+1 ),
346 $ n-j+1 )
347*
348 j = j + fjb
349 GO TO 30
350 END IF
351 ELSE
352 j = nfxd + 1
353 END IF
354*
355* Use unblocked code to factor the last or only block.
356*
357*
358 IF( j.LE.minmn )
359 $ CALL claqp2( m, n-j+1, j-1, a( 1, j ), lda, jpvt( j ),
360 $ tau( j ), rwork( j ), rwork( n+j ), work( 1 ) )
361*
362 END IF
363*
364 work( 1 ) = cmplx( lwkopt )
365 RETURN
366*
367* End of CGEQP3
368*
369 END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:81
subroutine cgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGEQRF
Definition: cgeqrf.f:146
subroutine cgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, RWORK, INFO)
CGEQP3
Definition: cgeqp3.f:159
subroutine claqp2(M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2, WORK)
CLAQP2 computes a QR factorization with column pivoting of the matrix block.
Definition: claqp2.f:149
subroutine claqps(M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU, VN1, VN2, AUXV, F, LDF)
CLAQPS computes a step of QR factorization with column pivoting of a real m-by-n matrix A by using BL...
Definition: claqps.f:178
subroutine cunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMQR
Definition: cunmqr.f:168