Annotation of rpl/lapack/lapack/zgeqp3.f, revision 1.5
1.1 bertrand 1: SUBROUTINE ZGEQP3( M, N, A, LDA, JPVT, TAU, WORK, LWORK, RWORK,
2: $ INFO )
3: *
4: * -- LAPACK routine (version 3.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
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'
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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgeqp3.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack