Annotation of rpl/lapack/lapack/zgeqpf.f, revision 1.19
1.10 bertrand 1: *> \brief \b ZGEQPF
2: *
3: * =========== DOCUMENTATION ===========
4: *
1.16 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.10 bertrand 7: *
8: *> \htmlonly
1.16 bertrand 9: *> Download ZGEQPF + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgeqpf.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgeqpf.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgeqpf.f">
1.10 bertrand 15: *> [TXT]</a>
1.16 bertrand 16: *> \endhtmlonly
1.10 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZGEQPF( M, N, A, LDA, JPVT, TAU, WORK, RWORK, INFO )
1.16 bertrand 22: *
1.10 bertrand 23: * .. Scalar Arguments ..
24: * INTEGER INFO, LDA, M, N
25: * ..
26: * .. Array Arguments ..
27: * INTEGER JPVT( * )
28: * DOUBLE PRECISION RWORK( * )
29: * COMPLEX*16 A( LDA, * ), TAU( * ), WORK( * )
30: * ..
1.16 bertrand 31: *
1.10 bertrand 32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> This routine is deprecated and has been replaced by routine ZGEQP3.
39: *>
40: *> ZGEQPF computes a QR factorization with column pivoting of a
41: *> complex M-by-N matrix A: A*P = Q*R.
42: *> \endverbatim
43: *
44: * Arguments:
45: * ==========
46: *
47: *> \param[in] M
48: *> \verbatim
49: *> M is INTEGER
50: *> The number of rows of the matrix A. M >= 0.
51: *> \endverbatim
52: *>
53: *> \param[in] N
54: *> \verbatim
55: *> N is INTEGER
56: *> The number of columns of the matrix A. N >= 0
57: *> \endverbatim
58: *>
59: *> \param[in,out] A
60: *> \verbatim
61: *> A is COMPLEX*16 array, dimension (LDA,N)
62: *> On entry, the M-by-N matrix A.
63: *> On exit, the upper triangle of the array contains the
64: *> min(M,N)-by-N upper triangular matrix R; the elements
65: *> below the diagonal, together with the array TAU,
66: *> represent the unitary matrix Q as a product of
67: *> min(m,n) elementary reflectors.
68: *> \endverbatim
69: *>
70: *> \param[in] LDA
71: *> \verbatim
72: *> LDA is INTEGER
73: *> The leading dimension of the array A. LDA >= max(1,M).
74: *> \endverbatim
75: *>
76: *> \param[in,out] JPVT
77: *> \verbatim
78: *> JPVT is INTEGER array, dimension (N)
79: *> On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
80: *> to the front of A*P (a leading column); if JPVT(i) = 0,
81: *> the i-th column of A is a free column.
82: *> On exit, if JPVT(i) = k, then the i-th column of A*P
83: *> was the k-th column of A.
84: *> \endverbatim
85: *>
86: *> \param[out] TAU
87: *> \verbatim
88: *> TAU is COMPLEX*16 array, dimension (min(M,N))
89: *> The scalar factors of the elementary reflectors.
90: *> \endverbatim
91: *>
92: *> \param[out] WORK
93: *> \verbatim
94: *> WORK is COMPLEX*16 array, dimension (N)
95: *> \endverbatim
96: *>
97: *> \param[out] RWORK
98: *> \verbatim
99: *> RWORK is DOUBLE PRECISION array, dimension (2*N)
100: *> \endverbatim
101: *>
102: *> \param[out] INFO
103: *> \verbatim
104: *> INFO is INTEGER
105: *> = 0: successful exit
106: *> < 0: if INFO = -i, the i-th argument had an illegal value
107: *> \endverbatim
108: *
109: * Authors:
110: * ========
111: *
1.16 bertrand 112: *> \author Univ. of Tennessee
113: *> \author Univ. of California Berkeley
114: *> \author Univ. of Colorado Denver
115: *> \author NAG Ltd.
1.10 bertrand 116: *
117: *> \ingroup complex16GEcomputational
118: *
119: *> \par Further Details:
120: * =====================
121: *>
122: *> \verbatim
123: *>
124: *> The matrix Q is represented as a product of elementary reflectors
125: *>
126: *> Q = H(1) H(2) . . . H(n)
127: *>
128: *> Each H(i) has the form
129: *>
130: *> H = I - tau * v * v**H
131: *>
132: *> where tau is a complex scalar, and v is a complex vector with
133: *> v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i).
134: *>
135: *> The matrix P is represented in jpvt as follows: If
136: *> jpvt(j) = i
137: *> then the jth column of P is the ith canonical unit vector.
138: *>
139: *> Partial column norm updating strategy modified by
140: *> Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
141: *> University of Zagreb, Croatia.
142: *> -- April 2011 --
143: *> For more details see LAPACK Working Note 176.
144: *> \endverbatim
145: *>
146: * =====================================================================
1.1 bertrand 147: SUBROUTINE ZGEQPF( M, N, A, LDA, JPVT, TAU, WORK, RWORK, INFO )
148: *
1.19 ! bertrand 149: * -- LAPACK computational routine --
1.1 bertrand 150: * -- LAPACK is a software package provided by Univ. of Tennessee, --
151: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
152: *
153: * .. Scalar Arguments ..
154: INTEGER INFO, LDA, M, N
155: * ..
156: * .. Array Arguments ..
157: INTEGER JPVT( * )
158: DOUBLE PRECISION RWORK( * )
159: COMPLEX*16 A( LDA, * ), TAU( * ), WORK( * )
160: * ..
161: *
162: * =====================================================================
163: *
164: * .. Parameters ..
165: DOUBLE PRECISION ZERO, ONE
166: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
167: * ..
168: * .. Local Scalars ..
169: INTEGER I, ITEMP, J, MA, MN, PVT
170: DOUBLE PRECISION TEMP, TEMP2, TOL3Z
171: COMPLEX*16 AII
172: * ..
173: * .. External Subroutines ..
1.5 bertrand 174: EXTERNAL XERBLA, ZGEQR2, ZLARF, ZLARFG, ZSWAP, ZUNM2R
1.1 bertrand 175: * ..
176: * .. Intrinsic Functions ..
177: INTRINSIC ABS, DCMPLX, DCONJG, MAX, MIN, SQRT
178: * ..
179: * .. External Functions ..
180: INTEGER IDAMAX
181: DOUBLE PRECISION DLAMCH, DZNRM2
182: EXTERNAL IDAMAX, DLAMCH, DZNRM2
183: * ..
184: * .. Executable Statements ..
185: *
186: * Test the input arguments
187: *
188: INFO = 0
189: IF( M.LT.0 ) THEN
190: INFO = -1
191: ELSE IF( N.LT.0 ) THEN
192: INFO = -2
193: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
194: INFO = -4
195: END IF
196: IF( INFO.NE.0 ) THEN
197: CALL XERBLA( 'ZGEQPF', -INFO )
198: RETURN
199: END IF
200: *
201: MN = MIN( M, N )
202: TOL3Z = SQRT(DLAMCH('Epsilon'))
203: *
204: * Move initial columns up front
205: *
206: ITEMP = 1
207: DO 10 I = 1, N
208: IF( JPVT( I ).NE.0 ) THEN
209: IF( I.NE.ITEMP ) THEN
210: CALL ZSWAP( M, A( 1, I ), 1, A( 1, ITEMP ), 1 )
211: JPVT( I ) = JPVT( ITEMP )
212: JPVT( ITEMP ) = I
213: ELSE
214: JPVT( I ) = I
215: END IF
216: ITEMP = ITEMP + 1
217: ELSE
218: JPVT( I ) = I
219: END IF
220: 10 CONTINUE
221: ITEMP = ITEMP - 1
222: *
223: * Compute the QR factorization and update remaining columns
224: *
225: IF( ITEMP.GT.0 ) THEN
226: MA = MIN( ITEMP, M )
227: CALL ZGEQR2( M, MA, A, LDA, TAU, WORK, INFO )
228: IF( MA.LT.N ) THEN
229: CALL ZUNM2R( 'Left', 'Conjugate transpose', M, N-MA, MA, A,
230: $ LDA, TAU, A( 1, MA+1 ), LDA, WORK, INFO )
231: END IF
232: END IF
233: *
234: IF( ITEMP.LT.MN ) THEN
235: *
236: * Initialize partial column norms. The first n elements of
237: * work store the exact column norms.
238: *
239: DO 20 I = ITEMP + 1, N
240: RWORK( I ) = DZNRM2( M-ITEMP, A( ITEMP+1, I ), 1 )
241: RWORK( N+I ) = RWORK( I )
242: 20 CONTINUE
243: *
244: * Compute factorization
245: *
246: DO 40 I = ITEMP + 1, MN
247: *
248: * Determine ith pivot column and swap if necessary
249: *
250: PVT = ( I-1 ) + IDAMAX( N-I+1, RWORK( I ), 1 )
251: *
252: IF( PVT.NE.I ) THEN
253: CALL ZSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 )
254: ITEMP = JPVT( PVT )
255: JPVT( PVT ) = JPVT( I )
256: JPVT( I ) = ITEMP
257: RWORK( PVT ) = RWORK( I )
258: RWORK( N+PVT ) = RWORK( N+I )
259: END IF
260: *
261: * Generate elementary reflector H(i)
262: *
263: AII = A( I, I )
1.5 bertrand 264: CALL ZLARFG( M-I+1, AII, A( MIN( I+1, M ), I ), 1,
1.1 bertrand 265: $ TAU( I ) )
266: A( I, I ) = AII
267: *
268: IF( I.LT.N ) THEN
269: *
270: * Apply H(i) to A(i:m,i+1:n) from the left
271: *
272: AII = A( I, I )
273: A( I, I ) = DCMPLX( ONE )
274: CALL ZLARF( 'Left', M-I+1, N-I, A( I, I ), 1,
275: $ DCONJG( TAU( I ) ), A( I, I+1 ), LDA, WORK )
276: A( I, I ) = AII
277: END IF
278: *
279: * Update partial column norms
280: *
281: DO 30 J = I + 1, N
282: IF( RWORK( J ).NE.ZERO ) THEN
283: *
284: * NOTE: The following 4 lines follow from the analysis in
285: * Lapack Working Note 176.
1.16 bertrand 286: *
1.1 bertrand 287: TEMP = ABS( A( I, J ) ) / RWORK( J )
288: TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) )
289: TEMP2 = TEMP*( RWORK( J ) / RWORK( N+J ) )**2
1.16 bertrand 290: IF( TEMP2 .LE. TOL3Z ) THEN
1.1 bertrand 291: IF( M-I.GT.0 ) THEN
292: RWORK( J ) = DZNRM2( M-I, A( I+1, J ), 1 )
293: RWORK( N+J ) = RWORK( J )
294: ELSE
295: RWORK( J ) = ZERO
296: RWORK( N+J ) = ZERO
297: END IF
298: ELSE
299: RWORK( J ) = RWORK( J )*SQRT( TEMP )
300: END IF
301: END IF
302: 30 CONTINUE
303: *
304: 40 CONTINUE
305: END IF
306: RETURN
307: *
308: * End of ZGEQPF
309: *
310: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgeqpf.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack