Annotation of rpl/lapack/lapack/zgeqpf.f, revision 1.6
1.1 bertrand 1: SUBROUTINE ZGEQPF( M, N, A, LDA, JPVT, TAU, WORK, RWORK, INFO )
2: *
1.5 bertrand 3: * -- LAPACK deprecated computational routine (version 3.2.2) --
1.1 bertrand 4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.5 bertrand 6: * June 2010
1.1 bertrand 7: *
8: * .. Scalar Arguments ..
9: INTEGER INFO, LDA, M, N
10: * ..
11: * .. Array Arguments ..
12: INTEGER JPVT( * )
13: DOUBLE PRECISION RWORK( * )
14: COMPLEX*16 A( LDA, * ), TAU( * ), WORK( * )
15: * ..
16: *
17: * Purpose
18: * =======
19: *
20: * This routine is deprecated and has been replaced by routine ZGEQP3.
21: *
22: * ZGEQPF computes a QR factorization with column pivoting of a
23: * complex M-by-N matrix A: A*P = Q*R.
24: *
25: * Arguments
26: * =========
27: *
28: * M (input) INTEGER
29: * The number of rows of the matrix A. M >= 0.
30: *
31: * N (input) INTEGER
32: * The number of columns of the matrix A. N >= 0
33: *
34: * A (input/output) COMPLEX*16 array, dimension (LDA,N)
35: * On entry, the M-by-N matrix A.
36: * On exit, the upper triangle of the array contains the
37: * min(M,N)-by-N upper triangular matrix R; the elements
38: * below the diagonal, together with the array TAU,
39: * represent the unitary matrix Q as a product of
40: * min(m,n) elementary reflectors.
41: *
42: * LDA (input) INTEGER
43: * The leading dimension of the array A. LDA >= max(1,M).
44: *
45: * JPVT (input/output) INTEGER array, dimension (N)
46: * On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
47: * to the front of A*P (a leading column); if JPVT(i) = 0,
48: * the i-th column of A is a free column.
49: * On exit, if JPVT(i) = k, then the i-th column of A*P
50: * was the k-th column of A.
51: *
52: * TAU (output) COMPLEX*16 array, dimension (min(M,N))
53: * The scalar factors of the elementary reflectors.
54: *
55: * WORK (workspace) COMPLEX*16 array, dimension (N)
56: *
57: * RWORK (workspace) DOUBLE PRECISION array, dimension (2*N)
58: *
59: * INFO (output) INTEGER
60: * = 0: successful exit
61: * < 0: if INFO = -i, the i-th argument had an illegal value
62: *
63: * Further Details
64: * ===============
65: *
66: * The matrix Q is represented as a product of elementary reflectors
67: *
68: * Q = H(1) H(2) . . . H(n)
69: *
70: * Each H(i) has the form
71: *
72: * H = I - tau * v * v'
73: *
74: * where tau is a complex scalar, and v is a complex vector with
75: * v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i).
76: *
77: * The matrix P is represented in jpvt as follows: If
78: * jpvt(j) = i
79: * then the jth column of P is the ith canonical unit vector.
80: *
81: * Partial column norm updating strategy modified by
82: * Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
83: * University of Zagreb, Croatia.
1.5 bertrand 84: * June 2010
1.1 bertrand 85: * For more details see LAPACK Working Note 176.
86: *
87: * =====================================================================
88: *
89: * .. Parameters ..
90: DOUBLE PRECISION ZERO, ONE
91: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
92: * ..
93: * .. Local Scalars ..
94: INTEGER I, ITEMP, J, MA, MN, PVT
95: DOUBLE PRECISION TEMP, TEMP2, TOL3Z
96: COMPLEX*16 AII
97: * ..
98: * .. External Subroutines ..
1.5 bertrand 99: EXTERNAL XERBLA, ZGEQR2, ZLARF, ZLARFG, ZSWAP, ZUNM2R
1.1 bertrand 100: * ..
101: * .. Intrinsic Functions ..
102: INTRINSIC ABS, DCMPLX, DCONJG, MAX, MIN, SQRT
103: * ..
104: * .. External Functions ..
105: INTEGER IDAMAX
106: DOUBLE PRECISION DLAMCH, DZNRM2
107: EXTERNAL IDAMAX, DLAMCH, DZNRM2
108: * ..
109: * .. Executable Statements ..
110: *
111: * Test the input arguments
112: *
113: INFO = 0
114: IF( M.LT.0 ) THEN
115: INFO = -1
116: ELSE IF( N.LT.0 ) THEN
117: INFO = -2
118: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
119: INFO = -4
120: END IF
121: IF( INFO.NE.0 ) THEN
122: CALL XERBLA( 'ZGEQPF', -INFO )
123: RETURN
124: END IF
125: *
126: MN = MIN( M, N )
127: TOL3Z = SQRT(DLAMCH('Epsilon'))
128: *
129: * Move initial columns up front
130: *
131: ITEMP = 1
132: DO 10 I = 1, N
133: IF( JPVT( I ).NE.0 ) THEN
134: IF( I.NE.ITEMP ) THEN
135: CALL ZSWAP( M, A( 1, I ), 1, A( 1, ITEMP ), 1 )
136: JPVT( I ) = JPVT( ITEMP )
137: JPVT( ITEMP ) = I
138: ELSE
139: JPVT( I ) = I
140: END IF
141: ITEMP = ITEMP + 1
142: ELSE
143: JPVT( I ) = I
144: END IF
145: 10 CONTINUE
146: ITEMP = ITEMP - 1
147: *
148: * Compute the QR factorization and update remaining columns
149: *
150: IF( ITEMP.GT.0 ) THEN
151: MA = MIN( ITEMP, M )
152: CALL ZGEQR2( M, MA, A, LDA, TAU, WORK, INFO )
153: IF( MA.LT.N ) THEN
154: CALL ZUNM2R( 'Left', 'Conjugate transpose', M, N-MA, MA, A,
155: $ LDA, TAU, A( 1, MA+1 ), LDA, WORK, INFO )
156: END IF
157: END IF
158: *
159: IF( ITEMP.LT.MN ) THEN
160: *
161: * Initialize partial column norms. The first n elements of
162: * work store the exact column norms.
163: *
164: DO 20 I = ITEMP + 1, N
165: RWORK( I ) = DZNRM2( M-ITEMP, A( ITEMP+1, I ), 1 )
166: RWORK( N+I ) = RWORK( I )
167: 20 CONTINUE
168: *
169: * Compute factorization
170: *
171: DO 40 I = ITEMP + 1, MN
172: *
173: * Determine ith pivot column and swap if necessary
174: *
175: PVT = ( I-1 ) + IDAMAX( N-I+1, RWORK( I ), 1 )
176: *
177: IF( PVT.NE.I ) THEN
178: CALL ZSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 )
179: ITEMP = JPVT( PVT )
180: JPVT( PVT ) = JPVT( I )
181: JPVT( I ) = ITEMP
182: RWORK( PVT ) = RWORK( I )
183: RWORK( N+PVT ) = RWORK( N+I )
184: END IF
185: *
186: * Generate elementary reflector H(i)
187: *
188: AII = A( I, I )
1.5 bertrand 189: CALL ZLARFG( M-I+1, AII, A( MIN( I+1, M ), I ), 1,
1.1 bertrand 190: $ TAU( I ) )
191: A( I, I ) = AII
192: *
193: IF( I.LT.N ) THEN
194: *
195: * Apply H(i) to A(i:m,i+1:n) from the left
196: *
197: AII = A( I, I )
198: A( I, I ) = DCMPLX( ONE )
199: CALL ZLARF( 'Left', M-I+1, N-I, A( I, I ), 1,
200: $ DCONJG( TAU( I ) ), A( I, I+1 ), LDA, WORK )
201: A( I, I ) = AII
202: END IF
203: *
204: * Update partial column norms
205: *
206: DO 30 J = I + 1, N
207: IF( RWORK( J ).NE.ZERO ) THEN
208: *
209: * NOTE: The following 4 lines follow from the analysis in
210: * Lapack Working Note 176.
211: *
212: TEMP = ABS( A( I, J ) ) / RWORK( J )
213: TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) )
214: TEMP2 = TEMP*( RWORK( J ) / RWORK( N+J ) )**2
215: IF( TEMP2 .LE. TOL3Z ) THEN
216: IF( M-I.GT.0 ) THEN
217: RWORK( J ) = DZNRM2( M-I, A( I+1, J ), 1 )
218: RWORK( N+J ) = RWORK( J )
219: ELSE
220: RWORK( J ) = ZERO
221: RWORK( N+J ) = ZERO
222: END IF
223: ELSE
224: RWORK( J ) = RWORK( J )*SQRT( TEMP )
225: END IF
226: END IF
227: 30 CONTINUE
228: *
229: 40 CONTINUE
230: END IF
231: RETURN
232: *
233: * End of ZGEQPF
234: *
235: 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