001:       SUBROUTINE ZGEQP3( M, N, A, LDA, JPVT, TAU, WORK, LWORK, RWORK,
002:      $                   INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       INTEGER            INFO, LDA, LWORK, M, N
011: *     ..
012: *     .. Array Arguments ..
013:       INTEGER            JPVT( * )
014:       DOUBLE PRECISION   RWORK( * )
015:       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  ZGEQP3 computes a QR factorization with column pivoting of a
022: *  matrix A:  A*P = Q*R  using Level 3 BLAS.
023: *
024: *  Arguments
025: *  =========
026: *
027: *  M       (input) INTEGER
028: *          The number of rows of the matrix A. M >= 0.
029: *
030: *  N       (input) INTEGER
031: *          The number of columns of the matrix A.  N >= 0.
032: *
033: *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
034: *          On entry, the M-by-N matrix A.
035: *          On exit, the upper triangle of the array contains the
036: *          min(M,N)-by-N upper trapezoidal matrix R; the elements below
037: *          the diagonal, together with the array TAU, represent the
038: *          unitary matrix Q as a product of min(M,N) elementary
039: *          reflectors.
040: *
041: *  LDA     (input) INTEGER
042: *          The leading dimension of the array A. LDA >= max(1,M).
043: *
044: *  JPVT    (input/output) INTEGER array, dimension (N)
045: *          On entry, if JPVT(J).ne.0, the J-th column of A is permuted
046: *          to the front of A*P (a leading column); if JPVT(J)=0,
047: *          the J-th column of A is a free column.
048: *          On exit, if JPVT(J)=K, then the J-th column of A*P was the
049: *          the K-th column of A.
050: *
051: *  TAU     (output) COMPLEX*16 array, dimension (min(M,N))
052: *          The scalar factors of the elementary reflectors.
053: *
054: *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
055: *          On exit, if INFO=0, WORK(1) returns the optimal LWORK.
056: *
057: *  LWORK   (input) INTEGER
058: *          The dimension of the array WORK. LWORK >= N+1.
059: *          For optimal performance LWORK >= ( N+1 )*NB, where NB
060: *          is the optimal blocksize.
061: *
062: *          If LWORK = -1, then a workspace query is assumed; the routine
063: *          only calculates the optimal size of the WORK array, returns
064: *          this value as the first entry of the WORK array, and no error
065: *          message related to LWORK is issued by XERBLA.
066: *
067: *  RWORK   (workspace) DOUBLE PRECISION array, dimension (2*N)
068: *
069: *  INFO    (output) INTEGER
070: *          = 0: successful exit.
071: *          < 0: if INFO = -i, the i-th argument had an illegal value.
072: *
073: *  Further Details
074: *  ===============
075: *
076: *  The matrix Q is represented as a product of elementary reflectors
077: *
078: *     Q = H(1) H(2) . . . H(k), where k = min(m,n).
079: *
080: *  Each H(i) has the form
081: *
082: *     H(i) = I - tau * v * v'
083: *
084: *  where tau is a real/complex scalar, and v is a real/complex vector
085: *  with v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in
086: *  A(i+1:m,i), and tau in TAU(i).
087: *
088: *  Based on contributions by
089: *    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
090: *    X. Sun, Computer Science Dept., Duke University, USA
091: *
092: *  =====================================================================
093: *
094: *     .. Parameters ..
095:       INTEGER            INB, INBMIN, IXOVER
096:       PARAMETER          ( INB = 1, INBMIN = 2, IXOVER = 3 )
097: *     ..
098: *     .. Local Scalars ..
099:       LOGICAL            LQUERY
100:       INTEGER            FJB, IWS, J, JB, LWKOPT, MINMN, MINWS, NA, NB,
101:      $                   NBMIN, NFXD, NX, SM, SMINMN, SN, TOPBMN
102: *     ..
103: *     .. External Subroutines ..
104:       EXTERNAL           XERBLA, ZGEQRF, ZLAQP2, ZLAQPS, ZSWAP, ZUNMQR
105: *     ..
106: *     .. External Functions ..
107:       INTEGER            ILAENV
108:       DOUBLE PRECISION   DZNRM2
109:       EXTERNAL           ILAENV, DZNRM2
110: *     ..
111: *     .. Intrinsic Functions ..
112:       INTRINSIC          INT, MAX, MIN
113: *     ..
114: *     .. Executable Statements ..
115: *
116: *     Test input arguments
117: *     ====================
118: *
119:       INFO = 0
120:       LQUERY = ( LWORK.EQ.-1 )
121:       IF( M.LT.0 ) THEN
122:          INFO = -1
123:       ELSE IF( N.LT.0 ) THEN
124:          INFO = -2
125:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
126:          INFO = -4
127:       END IF
128: *
129:       IF( INFO.EQ.0 ) THEN
130:          MINMN = MIN( M, N )
131:          IF( MINMN.EQ.0 ) THEN
132:             IWS = 1
133:             LWKOPT = 1
134:          ELSE
135:             IWS = N + 1
136:             NB = ILAENV( INB, 'ZGEQRF', ' ', M, N, -1, -1 )
137:             LWKOPT = ( N + 1 )*NB
138:          END IF
139:          WORK( 1 ) = LWKOPT
140: *
141:          IF( ( LWORK.LT.IWS ) .AND. .NOT.LQUERY ) THEN
142:             INFO = -8
143:          END IF
144:       END IF
145: *
146:       IF( INFO.NE.0 ) THEN
147:          CALL XERBLA( 'ZGEQP3', -INFO )
148:          RETURN
149:       ELSE IF( LQUERY ) THEN
150:          RETURN
151:       END IF
152: *
153: *     Quick return if possible.
154: *
155:       IF( MINMN.EQ.0 ) THEN
156:          RETURN
157:       END IF
158: *
159: *     Move initial columns up front.
160: *
161:       NFXD = 1
162:       DO 10 J = 1, N
163:          IF( JPVT( J ).NE.0 ) THEN
164:             IF( J.NE.NFXD ) THEN
165:                CALL ZSWAP( M, A( 1, J ), 1, A( 1, NFXD ), 1 )
166:                JPVT( J ) = JPVT( NFXD )
167:                JPVT( NFXD ) = J
168:             ELSE
169:                JPVT( J ) = J
170:             END IF
171:             NFXD = NFXD + 1
172:          ELSE
173:             JPVT( J ) = J
174:          END IF
175:    10 CONTINUE
176:       NFXD = NFXD - 1
177: *
178: *     Factorize fixed columns
179: *     =======================
180: *
181: *     Compute the QR factorization of fixed columns and update
182: *     remaining columns.
183: *
184:       IF( NFXD.GT.0 ) THEN
185:          NA = MIN( M, NFXD )
186: *CC      CALL ZGEQR2( M, NA, A, LDA, TAU, WORK, INFO )
187:          CALL ZGEQRF( M, NA, A, LDA, TAU, WORK, LWORK, INFO )
188:          IWS = MAX( IWS, INT( WORK( 1 ) ) )
189:          IF( NA.LT.N ) THEN
190: *CC         CALL ZUNM2R( 'Left', 'Conjugate Transpose', M, N-NA,
191: *CC  $                   NA, A, LDA, TAU, A( 1, NA+1 ), LDA, WORK,
192: *CC  $                   INFO )
193:             CALL ZUNMQR( 'Left', 'Conjugate Transpose', M, N-NA, NA, A,
194:      $                   LDA, TAU, A( 1, NA+1 ), LDA, WORK, LWORK,
195:      $                   INFO )
196:             IWS = MAX( IWS, INT( WORK( 1 ) ) )
197:          END IF
198:       END IF
199: *
200: *     Factorize free columns
201: *     ======================
202: *
203:       IF( NFXD.LT.MINMN ) THEN
204: *
205:          SM = M - NFXD
206:          SN = N - NFXD
207:          SMINMN = MINMN - NFXD
208: *
209: *        Determine the block size.
210: *
211:          NB = ILAENV( INB, 'ZGEQRF', ' ', SM, SN, -1, -1 )
212:          NBMIN = 2
213:          NX = 0
214: *
215:          IF( ( NB.GT.1 ) .AND. ( NB.LT.SMINMN ) ) THEN
216: *
217: *           Determine when to cross over from blocked to unblocked code.
218: *
219:             NX = MAX( 0, ILAENV( IXOVER, 'ZGEQRF', ' ', SM, SN, -1,
220:      $           -1 ) )
221: *
222: *
223:             IF( NX.LT.SMINMN ) THEN
224: *
225: *              Determine if workspace is large enough for blocked code.
226: *
227:                MINWS = ( SN+1 )*NB
228:                IWS = MAX( IWS, MINWS )
229:                IF( LWORK.LT.MINWS ) THEN
230: *
231: *                 Not enough workspace to use optimal NB: Reduce NB and
232: *                 determine the minimum value of NB.
233: *
234:                   NB = LWORK / ( SN+1 )
235:                   NBMIN = MAX( 2, ILAENV( INBMIN, 'ZGEQRF', ' ', SM, SN,
236:      $                    -1, -1 ) )
237: *
238: *
239:                END IF
240:             END IF
241:          END IF
242: *
243: *        Initialize partial column norms. The first N elements of work
244: *        store the exact column norms.
245: *
246:          DO 20 J = NFXD + 1, N
247:             RWORK( J ) = DZNRM2( SM, A( NFXD+1, J ), 1 )
248:             RWORK( N+J ) = RWORK( J )
249:    20    CONTINUE
250: *
251:          IF( ( NB.GE.NBMIN ) .AND. ( NB.LT.SMINMN ) .AND.
252:      $       ( NX.LT.SMINMN ) ) THEN
253: *
254: *           Use blocked code initially.
255: *
256:             J = NFXD + 1
257: *
258: *           Compute factorization: while loop.
259: *
260: *
261:             TOPBMN = MINMN - NX
262:    30       CONTINUE
263:             IF( J.LE.TOPBMN ) THEN
264:                JB = MIN( NB, TOPBMN-J+1 )
265: *
266: *              Factorize JB columns among columns J:N.
267: *
268:                CALL ZLAQPS( M, N-J+1, J-1, JB, FJB, A( 1, J ), LDA,
269:      $                      JPVT( J ), TAU( J ), RWORK( J ),
270:      $                      RWORK( N+J ), WORK( 1 ), WORK( JB+1 ),
271:      $                      N-J+1 )
272: *
273:                J = J + FJB
274:                GO TO 30
275:             END IF
276:          ELSE
277:             J = NFXD + 1
278:          END IF
279: *
280: *        Use unblocked code to factor the last or only block.
281: *
282: *
283:          IF( J.LE.MINMN )
284:      $      CALL ZLAQP2( M, N-J+1, J-1, A( 1, J ), LDA, JPVT( J ),
285:      $                   TAU( J ), RWORK( J ), RWORK( N+J ), WORK( 1 ) )
286: *
287:       END IF
288: *
289:       WORK( 1 ) = IWS
290:       RETURN
291: *
292: *     End of ZGEQP3
293: *
294:       END
295: