1*> \brief \b ZGEQP3
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download ZGEQP3 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgeqp3.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgeqp3.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgeqp3.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE ZGEQP3( 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*       DOUBLE PRECISION   RWORK( * )
30*       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
31*       ..
32*
33*
34*> \par Purpose:
35*  =============
36*>
37*> \verbatim
38*>
39*> ZGEQP3 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*16 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*16 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*16 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 DOUBLE PRECISION 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 complex16GEcomputational
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 ZGEQP3( 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      DOUBLE PRECISION   RWORK( * )
170      COMPLEX*16         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           XERBLA, ZGEQRF, ZLAQP2, ZLAQPS, ZSWAP, ZUNMQR
186*     ..
187*     .. External Functions ..
188      INTEGER            ILAENV
189      DOUBLE PRECISION   DZNRM2
190      EXTERNAL           ILAENV, DZNRM2
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, 'ZGEQRF', ' ', M, N, -1, -1 )
218            LWKOPT = ( N + 1 )*NB
219         END IF
220         WORK( 1 ) = DCMPLX( 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( 'ZGEQP3', -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 ZSWAP( 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 ZGEQR2( M, NA, A, LDA, TAU, WORK, INFO )
262         CALL ZGEQRF( M, NA, A, LDA, TAU, WORK, LWORK, INFO )
263         IWS = MAX( IWS, INT( WORK( 1 ) ) )
264         IF( NA.LT.N ) THEN
265*CC         CALL ZUNM2R( 'Left', 'Conjugate Transpose', M, N-NA,
266*CC  $                   NA, A, LDA, TAU, A( 1, NA+1 ), LDA, WORK,
267*CC  $                   INFO )
268            CALL ZUNMQR( '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, 'ZGEQRF', ' ', 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, 'ZGEQRF', ' ', 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, 'ZGEQRF', ' ', 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 ) = DZNRM2( 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 ZLAQPS( 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 ZLAQP2( 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 ) = DCMPLX( LWKOPT )
365      RETURN
366*
367*     End of ZGEQP3
368*
369      END
370