1 SUBROUTINE ZGEQPF( M, N, A, LDA, JPVT, TAU, WORK, RWORK, INFO )
2 *
3 * -- LAPACK deprecated computational routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
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**H
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.
84 * -- April 2011 --
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 ..
99 EXTERNAL XERBLA, ZGEQR2, ZLARF, ZLARFG, ZSWAP, ZUNM2R
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 )
189 CALL ZLARFG( M-I+1, AII, A( MIN( I+1, M ), I ), 1,
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
2 *
3 * -- LAPACK deprecated computational routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
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**H
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.
84 * -- April 2011 --
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 ..
99 EXTERNAL XERBLA, ZGEQR2, ZLARF, ZLARFG, ZSWAP, ZUNM2R
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 )
189 CALL ZLARFG( M-I+1, AII, A( MIN( I+1, M ), I ), 1,
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