001:       SUBROUTINE CGEQPF( M, N, A, LDA, JPVT, TAU, WORK, RWORK, INFO )
002: *
003: *  -- LAPACK deprecated driver routine (version 3.2) --
004: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
005: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       INTEGER            INFO, LDA, M, N
010: *     ..
011: *     .. Array Arguments ..
012:       INTEGER            JPVT( * )
013:       REAL               RWORK( * )
014:       COMPLEX            A( LDA, * ), TAU( * ), WORK( * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  This routine is deprecated and has been replaced by routine CGEQP3.
021: *
022: *  CGEQPF computes a QR factorization with column pivoting of a
023: *  complex M-by-N matrix A: A*P = Q*R.
024: *
025: *  Arguments
026: *  =========
027: *
028: *  M       (input) INTEGER
029: *          The number of rows of the matrix A. M >= 0.
030: *
031: *  N       (input) INTEGER
032: *          The number of columns of the matrix A. N >= 0
033: *
034: *  A       (input/output) COMPLEX array, dimension (LDA,N)
035: *          On entry, the M-by-N matrix A.
036: *          On exit, the upper triangle of the array contains the
037: *          min(M,N)-by-N upper triangular matrix R; the elements
038: *          below the diagonal, together with the array TAU,
039: *          represent the unitary matrix Q as a product of
040: *          min(m,n) elementary reflectors.
041: *
042: *  LDA     (input) INTEGER
043: *          The leading dimension of the array A. LDA >= max(1,M).
044: *
045: *  JPVT    (input/output) INTEGER array, dimension (N)
046: *          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
047: *          to the front of A*P (a leading column); if JPVT(i) = 0,
048: *          the i-th column of A is a free column.
049: *          On exit, if JPVT(i) = k, then the i-th column of A*P
050: *          was the k-th column of A.
051: *
052: *  TAU     (output) COMPLEX array, dimension (min(M,N))
053: *          The scalar factors of the elementary reflectors.
054: *
055: *  WORK    (workspace) COMPLEX array, dimension (N)
056: *
057: *  RWORK   (workspace) REAL array, dimension (2*N)
058: *
059: *  INFO    (output) INTEGER
060: *          = 0:  successful exit
061: *          < 0:  if INFO = -i, the i-th argument had an illegal value
062: *
063: *  Further Details
064: *  ===============
065: *
066: *  The matrix Q is represented as a product of elementary reflectors
067: *
068: *     Q = H(1) H(2) . . . H(n)
069: *
070: *  Each H(i) has the form
071: *
072: *     H = I - tau * v * v'
073: *
074: *  where tau is a complex scalar, and v is a complex vector with
075: *  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i).
076: *
077: *  The matrix P is represented in jpvt as follows: If
078: *     jpvt(j) = i
079: *  then the jth column of P is the ith canonical unit vector.
080: *
081: *  Partial column norm updating strategy modified by
082: *    Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
083: *    University of Zagreb, Croatia.
084: *    June 2006.
085: *  For more details see LAPACK Working Note 176.
086: *
087: *  =====================================================================
088: *
089: *     .. Parameters ..
090:       REAL               ZERO, ONE
091:       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
092: *     ..
093: *     .. Local Scalars ..
094:       INTEGER            I, ITEMP, J, MA, MN, PVT
095:       REAL               TEMP, TEMP2, TOL3Z
096:       COMPLEX            AII
097: *     ..
098: *     .. External Subroutines ..
099:       EXTERNAL           CGEQR2, CLARF, CLARFP, CSWAP, CUNM2R, XERBLA
100: *     ..
101: *     .. Intrinsic Functions ..
102:       INTRINSIC          ABS, CMPLX, CONJG, MAX, MIN, SQRT
103: *     ..
104: *     .. External Functions ..
105:       INTEGER            ISAMAX
106:       REAL               SCNRM2, SLAMCH
107:       EXTERNAL           ISAMAX, SCNRM2, SLAMCH
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( 'CGEQPF', -INFO )
123:          RETURN
124:       END IF
125: *
126:       MN = MIN( M, N )
127:       TOL3Z = SQRT(SLAMCH('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 CSWAP( 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 CGEQR2( M, MA, A, LDA, TAU, WORK, INFO )
153:          IF( MA.LT.N ) THEN
154:             CALL CUNM2R( '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 ) = SCNRM2( 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 ) + ISAMAX( N-I+1, RWORK( I ), 1 )
176: *
177:             IF( PVT.NE.I ) THEN
178:                CALL CSWAP( 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 CLARFP( 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 ) = CMPLX( ONE )
199:                CALL CLARF( 'Left', M-I+1, N-I, A( I, I ), 1,
200:      $                     CONJG( 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 ) = SCNRM2( 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 CGEQPF
234: *
235:       END
236: