1 SUBROUTINE ZGEQP3( M, N, A, LDA, JPVT, TAU, WORK, LWORK, RWORK,
2 $ INFO )
3 *
4 * -- LAPACK routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * .. Scalar Arguments ..
10 INTEGER INFO, LDA, LWORK, M, N
11 * ..
12 * .. Array Arguments ..
13 INTEGER JPVT( * )
14 DOUBLE PRECISION RWORK( * )
15 COMPLEX*16 A( LDA, * ), TAU( * ), WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZGEQP3 computes a QR factorization with column pivoting of a
22 * matrix A: A*P = Q*R using Level 3 BLAS.
23 *
24 * Arguments
25 * =========
26 *
27 * M (input) INTEGER
28 * The number of rows of the matrix A. M >= 0.
29 *
30 * N (input) INTEGER
31 * The number of columns of the matrix A. N >= 0.
32 *
33 * A (input/output) COMPLEX*16 array, dimension (LDA,N)
34 * On entry, the M-by-N matrix A.
35 * On exit, the upper triangle of the array contains the
36 * min(M,N)-by-N upper trapezoidal matrix R; the elements below
37 * the diagonal, together with the array TAU, represent the
38 * unitary matrix Q as a product of min(M,N) elementary
39 * reflectors.
40 *
41 * LDA (input) INTEGER
42 * The leading dimension of the array A. LDA >= max(1,M).
43 *
44 * JPVT (input/output) INTEGER array, dimension (N)
45 * On entry, if JPVT(J).ne.0, the J-th column of A is permuted
46 * to the front of A*P (a leading column); if JPVT(J)=0,
47 * the J-th column of A is a free column.
48 * On exit, if JPVT(J)=K, then the J-th column of A*P was the
49 * the K-th column of A.
50 *
51 * TAU (output) COMPLEX*16 array, dimension (min(M,N))
52 * The scalar factors of the elementary reflectors.
53 *
54 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
55 * On exit, if INFO=0, WORK(1) returns the optimal LWORK.
56 *
57 * LWORK (input) INTEGER
58 * The dimension of the array WORK. LWORK >= N+1.
59 * For optimal performance LWORK >= ( N+1 )*NB, where NB
60 * is the optimal blocksize.
61 *
62 * If LWORK = -1, then a workspace query is assumed; the routine
63 * only calculates the optimal size of the WORK array, returns
64 * this value as the first entry of the WORK array, and no error
65 * message related to LWORK is issued by XERBLA.
66 *
67 * RWORK (workspace) DOUBLE PRECISION array, dimension (2*N)
68 *
69 * INFO (output) INTEGER
70 * = 0: successful exit.
71 * < 0: if INFO = -i, the i-th argument had an illegal value.
72 *
73 * Further Details
74 * ===============
75 *
76 * The matrix Q is represented as a product of elementary reflectors
77 *
78 * Q = H(1) H(2) . . . H(k), where k = min(m,n).
79 *
80 * Each H(i) has the form
81 *
82 * H(i) = I - tau * v * v**H
83 *
84 * where tau is a real/complex scalar, and v is a real/complex vector
85 * with v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in
86 * A(i+1:m,i), and tau in TAU(i).
87 *
88 * Based on contributions by
89 * G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
90 * X. Sun, Computer Science Dept., Duke University, USA
91 *
92 * =====================================================================
93 *
94 * .. Parameters ..
95 INTEGER INB, INBMIN, IXOVER
96 PARAMETER ( INB = 1, INBMIN = 2, IXOVER = 3 )
97 * ..
98 * .. Local Scalars ..
99 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
2 $ INFO )
3 *
4 * -- LAPACK routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * .. Scalar Arguments ..
10 INTEGER INFO, LDA, LWORK, M, N
11 * ..
12 * .. Array Arguments ..
13 INTEGER JPVT( * )
14 DOUBLE PRECISION RWORK( * )
15 COMPLEX*16 A( LDA, * ), TAU( * ), WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZGEQP3 computes a QR factorization with column pivoting of a
22 * matrix A: A*P = Q*R using Level 3 BLAS.
23 *
24 * Arguments
25 * =========
26 *
27 * M (input) INTEGER
28 * The number of rows of the matrix A. M >= 0.
29 *
30 * N (input) INTEGER
31 * The number of columns of the matrix A. N >= 0.
32 *
33 * A (input/output) COMPLEX*16 array, dimension (LDA,N)
34 * On entry, the M-by-N matrix A.
35 * On exit, the upper triangle of the array contains the
36 * min(M,N)-by-N upper trapezoidal matrix R; the elements below
37 * the diagonal, together with the array TAU, represent the
38 * unitary matrix Q as a product of min(M,N) elementary
39 * reflectors.
40 *
41 * LDA (input) INTEGER
42 * The leading dimension of the array A. LDA >= max(1,M).
43 *
44 * JPVT (input/output) INTEGER array, dimension (N)
45 * On entry, if JPVT(J).ne.0, the J-th column of A is permuted
46 * to the front of A*P (a leading column); if JPVT(J)=0,
47 * the J-th column of A is a free column.
48 * On exit, if JPVT(J)=K, then the J-th column of A*P was the
49 * the K-th column of A.
50 *
51 * TAU (output) COMPLEX*16 array, dimension (min(M,N))
52 * The scalar factors of the elementary reflectors.
53 *
54 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
55 * On exit, if INFO=0, WORK(1) returns the optimal LWORK.
56 *
57 * LWORK (input) INTEGER
58 * The dimension of the array WORK. LWORK >= N+1.
59 * For optimal performance LWORK >= ( N+1 )*NB, where NB
60 * is the optimal blocksize.
61 *
62 * If LWORK = -1, then a workspace query is assumed; the routine
63 * only calculates the optimal size of the WORK array, returns
64 * this value as the first entry of the WORK array, and no error
65 * message related to LWORK is issued by XERBLA.
66 *
67 * RWORK (workspace) DOUBLE PRECISION array, dimension (2*N)
68 *
69 * INFO (output) INTEGER
70 * = 0: successful exit.
71 * < 0: if INFO = -i, the i-th argument had an illegal value.
72 *
73 * Further Details
74 * ===============
75 *
76 * The matrix Q is represented as a product of elementary reflectors
77 *
78 * Q = H(1) H(2) . . . H(k), where k = min(m,n).
79 *
80 * Each H(i) has the form
81 *
82 * H(i) = I - tau * v * v**H
83 *
84 * where tau is a real/complex scalar, and v is a real/complex vector
85 * with v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in
86 * A(i+1:m,i), and tau in TAU(i).
87 *
88 * Based on contributions by
89 * G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
90 * X. Sun, Computer Science Dept., Duke University, USA
91 *
92 * =====================================================================
93 *
94 * .. Parameters ..
95 INTEGER INB, INBMIN, IXOVER
96 PARAMETER ( INB = 1, INBMIN = 2, IXOVER = 3 )
97 * ..
98 * .. Local Scalars ..
99 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