|
2329
|
1 SUBROUTINE ZGEQPF( M, N, A, LDA, JPVT, TAU, WORK, RWORK, INFO ) |
|
|
2 * |
|
3333
|
3 * -- LAPACK auxiliary routine (version 3.0) -- |
|
2329
|
4 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., |
|
|
5 * Courant Institute, Argonne National Lab, and Rice University |
|
3333
|
6 * June 30, 1999 |
|
2329
|
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 * |
|
3333
|
20 * This routine is deprecated and has been replaced by routine ZGEQP3. |
|
|
21 * |
|
2329
|
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, |
|
3333
|
39 * represent the unitary matrix Q as a product of |
|
2329
|
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 * ===================================================================== |
|
|
82 * |
|
|
83 * .. Parameters .. |
|
|
84 DOUBLE PRECISION ZERO, ONE |
|
|
85 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) |
|
|
86 * .. |
|
|
87 * .. Local Scalars .. |
|
|
88 INTEGER I, ITEMP, J, MA, MN, PVT |
|
|
89 DOUBLE PRECISION TEMP, TEMP2 |
|
|
90 COMPLEX*16 AII |
|
|
91 * .. |
|
|
92 * .. External Subroutines .. |
|
|
93 EXTERNAL XERBLA, ZGEQR2, ZLARF, ZLARFG, ZSWAP, ZUNM2R |
|
|
94 * .. |
|
|
95 * .. Intrinsic Functions .. |
|
|
96 INTRINSIC ABS, DCMPLX, DCONJG, MAX, MIN, SQRT |
|
|
97 * .. |
|
|
98 * .. External Functions .. |
|
|
99 INTEGER IDAMAX |
|
|
100 DOUBLE PRECISION DZNRM2 |
|
|
101 EXTERNAL IDAMAX, DZNRM2 |
|
|
102 * .. |
|
|
103 * .. Executable Statements .. |
|
|
104 * |
|
|
105 * Test the input arguments |
|
|
106 * |
|
|
107 INFO = 0 |
|
|
108 IF( M.LT.0 ) THEN |
|
|
109 INFO = -1 |
|
|
110 ELSE IF( N.LT.0 ) THEN |
|
|
111 INFO = -2 |
|
|
112 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN |
|
|
113 INFO = -4 |
|
|
114 END IF |
|
|
115 IF( INFO.NE.0 ) THEN |
|
|
116 CALL XERBLA( 'ZGEQPF', -INFO ) |
|
|
117 RETURN |
|
|
118 END IF |
|
|
119 * |
|
|
120 MN = MIN( M, N ) |
|
|
121 * |
|
|
122 * Move initial columns up front |
|
|
123 * |
|
|
124 ITEMP = 1 |
|
|
125 DO 10 I = 1, N |
|
|
126 IF( JPVT( I ).NE.0 ) THEN |
|
|
127 IF( I.NE.ITEMP ) THEN |
|
|
128 CALL ZSWAP( M, A( 1, I ), 1, A( 1, ITEMP ), 1 ) |
|
|
129 JPVT( I ) = JPVT( ITEMP ) |
|
|
130 JPVT( ITEMP ) = I |
|
|
131 ELSE |
|
|
132 JPVT( I ) = I |
|
|
133 END IF |
|
|
134 ITEMP = ITEMP + 1 |
|
|
135 ELSE |
|
|
136 JPVT( I ) = I |
|
|
137 END IF |
|
|
138 10 CONTINUE |
|
|
139 ITEMP = ITEMP - 1 |
|
|
140 * |
|
|
141 * Compute the QR factorization and update remaining columns |
|
|
142 * |
|
|
143 IF( ITEMP.GT.0 ) THEN |
|
|
144 MA = MIN( ITEMP, M ) |
|
|
145 CALL ZGEQR2( M, MA, A, LDA, TAU, WORK, INFO ) |
|
|
146 IF( MA.LT.N ) THEN |
|
|
147 CALL ZUNM2R( 'Left', 'Conjugate transpose', M, N-MA, MA, A, |
|
|
148 $ LDA, TAU, A( 1, MA+1 ), LDA, WORK, INFO ) |
|
|
149 END IF |
|
|
150 END IF |
|
|
151 * |
|
|
152 IF( ITEMP.LT.MN ) THEN |
|
|
153 * |
|
|
154 * Initialize partial column norms. The first n elements of |
|
|
155 * work store the exact column norms. |
|
|
156 * |
|
|
157 DO 20 I = ITEMP + 1, N |
|
|
158 RWORK( I ) = DZNRM2( M-ITEMP, A( ITEMP+1, I ), 1 ) |
|
|
159 RWORK( N+I ) = RWORK( I ) |
|
|
160 20 CONTINUE |
|
|
161 * |
|
|
162 * Compute factorization |
|
|
163 * |
|
|
164 DO 40 I = ITEMP + 1, MN |
|
|
165 * |
|
|
166 * Determine ith pivot column and swap if necessary |
|
|
167 * |
|
|
168 PVT = ( I-1 ) + IDAMAX( N-I+1, RWORK( I ), 1 ) |
|
|
169 * |
|
|
170 IF( PVT.NE.I ) THEN |
|
|
171 CALL ZSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 ) |
|
|
172 ITEMP = JPVT( PVT ) |
|
|
173 JPVT( PVT ) = JPVT( I ) |
|
|
174 JPVT( I ) = ITEMP |
|
|
175 RWORK( PVT ) = RWORK( I ) |
|
|
176 RWORK( N+PVT ) = RWORK( N+I ) |
|
|
177 END IF |
|
|
178 * |
|
|
179 * Generate elementary reflector H(i) |
|
|
180 * |
|
|
181 AII = A( I, I ) |
|
|
182 CALL ZLARFG( M-I+1, AII, A( MIN( I+1, M ), I ), 1, |
|
|
183 $ TAU( I ) ) |
|
|
184 A( I, I ) = AII |
|
|
185 * |
|
|
186 IF( I.LT.N ) THEN |
|
|
187 * |
|
|
188 * Apply H(i) to A(i:m,i+1:n) from the left |
|
|
189 * |
|
|
190 AII = A( I, I ) |
|
|
191 A( I, I ) = DCMPLX( ONE ) |
|
|
192 CALL ZLARF( 'Left', M-I+1, N-I, A( I, I ), 1, |
|
|
193 $ DCONJG( TAU( I ) ), A( I, I+1 ), LDA, WORK ) |
|
|
194 A( I, I ) = AII |
|
|
195 END IF |
|
|
196 * |
|
|
197 * Update partial column norms |
|
|
198 * |
|
|
199 DO 30 J = I + 1, N |
|
|
200 IF( RWORK( J ).NE.ZERO ) THEN |
|
|
201 TEMP = ONE - ( ABS( A( I, J ) ) / RWORK( J ) )**2 |
|
|
202 TEMP = MAX( TEMP, ZERO ) |
|
|
203 TEMP2 = ONE + 0.05D0*TEMP* |
|
|
204 $ ( RWORK( J ) / RWORK( N+J ) )**2 |
|
|
205 IF( TEMP2.EQ.ONE ) THEN |
|
|
206 IF( M-I.GT.0 ) THEN |
|
|
207 RWORK( J ) = DZNRM2( M-I, A( I+1, J ), 1 ) |
|
|
208 RWORK( N+J ) = RWORK( J ) |
|
|
209 ELSE |
|
|
210 RWORK( J ) = ZERO |
|
|
211 RWORK( N+J ) = ZERO |
|
|
212 END IF |
|
|
213 ELSE |
|
|
214 RWORK( J ) = RWORK( J )*SQRT( TEMP ) |
|
|
215 END IF |
|
|
216 END IF |
|
|
217 30 CONTINUE |
|
|
218 * |
|
|
219 40 CONTINUE |
|
|
220 END IF |
|
|
221 RETURN |
|
|
222 * |
|
|
223 * End of ZGEQPF |
|
|
224 * |
|
|
225 END |
