LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
cgeqpf.f
Go to the documentation of this file.
1 *> \brief \b CGEQPF
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CGEQPF + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgeqpf.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgeqpf.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgeqpf.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CGEQPF( M, N, A, LDA, JPVT, TAU, WORK, RWORK, INFO )
22 *
23 * .. Scalar Arguments ..
24 * INTEGER INFO, LDA, M, N
25 * ..
26 * .. Array Arguments ..
27 * INTEGER JPVT( * )
28 * REAL RWORK( * )
29 * COMPLEX A( LDA, * ), TAU( * ), WORK( * )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> This routine is deprecated and has been replaced by routine CGEQP3.
39 *>
40 *> CGEQPF 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 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 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 array, dimension (N)
95 *> \endverbatim
96 *>
97 *> \param[out] RWORK
98 *> \verbatim
99 *> RWORK is REAL 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 *
112 *> \author Univ. of Tennessee
113 *> \author Univ. of California Berkeley
114 *> \author Univ. of Colorado Denver
115 *> \author NAG Ltd.
116 *
117 *> \date November 2011
118 *
119 *> \ingroup complexGEcomputational
120 *
121 *> \par Further Details:
122 * =====================
123 *>
124 *> \verbatim
125 *>
126 *> The matrix Q is represented as a product of elementary reflectors
127 *>
128 *> Q = H(1) H(2) . . . H(n)
129 *>
130 *> Each H(i) has the form
131 *>
132 *> H = I - tau * v * v**H
133 *>
134 *> where tau is a complex scalar, and v is a complex vector with
135 *> v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i).
136 *>
137 *> The matrix P is represented in jpvt as follows: If
138 *> jpvt(j) = i
139 *> then the jth column of P is the ith canonical unit vector.
140 *>
141 *> Partial column norm updating strategy modified by
142 *> Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
143 *> University of Zagreb, Croatia.
144 *> -- April 2011 --
145 *> For more details see LAPACK Working Note 176.
146 *> \endverbatim
147 *>
148 * =====================================================================
149  SUBROUTINE cgeqpf( M, N, A, LDA, JPVT, TAU, WORK, RWORK, INFO )
150 *
151 * -- LAPACK computational routine (version 3.4.0) --
152 * -- LAPACK is a software package provided by Univ. of Tennessee, --
153 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
154 * November 2011
155 *
156 * .. Scalar Arguments ..
157  INTEGER info, lda, m, n
158 * ..
159 * .. Array Arguments ..
160  INTEGER jpvt( * )
161  REAL rwork( * )
162  COMPLEX a( lda, * ), tau( * ), work( * )
163 * ..
164 *
165 * =====================================================================
166 *
167 * .. Parameters ..
168  REAL zero, one
169  parameter( zero = 0.0e+0, one = 1.0e+0 )
170 * ..
171 * .. Local Scalars ..
172  INTEGER i, itemp, j, ma, mn, pvt
173  REAL temp, temp2, tol3z
174  COMPLEX aii
175 * ..
176 * .. External Subroutines ..
177  EXTERNAL cgeqr2, clarf, clarfg, cswap, cunm2r, xerbla
178 * ..
179 * .. Intrinsic Functions ..
180  INTRINSIC abs, cmplx, conjg, max, min, sqrt
181 * ..
182 * .. External Functions ..
183  INTEGER isamax
184  REAL scnrm2, slamch
185  EXTERNAL isamax, scnrm2, slamch
186 * ..
187 * .. Executable Statements ..
188 *
189 * Test the input arguments
190 *
191  info = 0
192  IF( m.LT.0 ) THEN
193  info = -1
194  ELSE IF( n.LT.0 ) THEN
195  info = -2
196  ELSE IF( lda.LT.max( 1, m ) ) THEN
197  info = -4
198  END IF
199  IF( info.NE.0 ) THEN
200  CALL xerbla( 'CGEQPF', -info )
201  RETURN
202  END IF
203 *
204  mn = min( m, n )
205  tol3z = sqrt(slamch('Epsilon'))
206 *
207 * Move initial columns up front
208 *
209  itemp = 1
210  DO 10 i = 1, n
211  IF( jpvt( i ).NE.0 ) THEN
212  IF( i.NE.itemp ) THEN
213  CALL cswap( m, a( 1, i ), 1, a( 1, itemp ), 1 )
214  jpvt( i ) = jpvt( itemp )
215  jpvt( itemp ) = i
216  ELSE
217  jpvt( i ) = i
218  END IF
219  itemp = itemp + 1
220  ELSE
221  jpvt( i ) = i
222  END IF
223  10 CONTINUE
224  itemp = itemp - 1
225 *
226 * Compute the QR factorization and update remaining columns
227 *
228  IF( itemp.GT.0 ) THEN
229  ma = min( itemp, m )
230  CALL cgeqr2( m, ma, a, lda, tau, work, info )
231  IF( ma.LT.n ) THEN
232  CALL cunm2r( 'Left', 'Conjugate transpose', m, n-ma, ma, a,
233  $ lda, tau, a( 1, ma+1 ), lda, work, info )
234  END IF
235  END IF
236 *
237  IF( itemp.LT.mn ) THEN
238 *
239 * Initialize partial column norms. The first n elements of
240 * work store the exact column norms.
241 *
242  DO 20 i = itemp + 1, n
243  rwork( i ) = scnrm2( m-itemp, a( itemp+1, i ), 1 )
244  rwork( n+i ) = rwork( i )
245  20 CONTINUE
246 *
247 * Compute factorization
248 *
249  DO 40 i = itemp + 1, mn
250 *
251 * Determine ith pivot column and swap if necessary
252 *
253  pvt = ( i-1 ) + isamax( n-i+1, rwork( i ), 1 )
254 *
255  IF( pvt.NE.i ) THEN
256  CALL cswap( m, a( 1, pvt ), 1, a( 1, i ), 1 )
257  itemp = jpvt( pvt )
258  jpvt( pvt ) = jpvt( i )
259  jpvt( i ) = itemp
260  rwork( pvt ) = rwork( i )
261  rwork( n+pvt ) = rwork( n+i )
262  END IF
263 *
264 * Generate elementary reflector H(i)
265 *
266  aii = a( i, i )
267  CALL clarfg( m-i+1, aii, a( min( i+1, m ), i ), 1,
268  $ tau( i ) )
269  a( i, i ) = aii
270 *
271  IF( i.LT.n ) THEN
272 *
273 * Apply H(i) to A(i:m,i+1:n) from the left
274 *
275  aii = a( i, i )
276  a( i, i ) = cmplx( one )
277  CALL clarf( 'Left', m-i+1, n-i, a( i, i ), 1,
278  $ conjg( tau( i ) ), a( i, i+1 ), lda, work )
279  a( i, i ) = aii
280  END IF
281 *
282 * Update partial column norms
283 *
284  DO 30 j = i + 1, n
285  IF( rwork( j ).NE.zero ) THEN
286 *
287 * NOTE: The following 4 lines follow from the analysis in
288 * Lapack Working Note 176.
289 *
290  temp = abs( a( i, j ) ) / rwork( j )
291  temp = max( zero, ( one+temp )*( one-temp ) )
292  temp2 = temp*( rwork( j ) / rwork( n+j ) )**2
293  IF( temp2 .LE. tol3z ) THEN
294  IF( m-i.GT.0 ) THEN
295  rwork( j ) = scnrm2( m-i, a( i+1, j ), 1 )
296  rwork( n+j ) = rwork( j )
297  ELSE
298  rwork( j ) = zero
299  rwork( n+j ) = zero
300  END IF
301  ELSE
302  rwork( j ) = rwork( j )*sqrt( temp )
303  END IF
304  END IF
305  30 CONTINUE
306 *
307  40 CONTINUE
308  END IF
309  RETURN
310 *
311 * End of CGEQPF
312 *
313  END