001:       SUBROUTINE DGEQPF( M, N, A, LDA, JPVT, TAU, WORK, 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:       DOUBLE PRECISION   A( LDA, * ), TAU( * ), WORK( * )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  This routine is deprecated and has been replaced by routine DGEQP3.
020: *
021: *  DGEQPF computes a QR factorization with column pivoting of a
022: *  real M-by-N matrix A: A*P = Q*R.
023: *
024: *  Arguments
025: *  =========
026: *
027: *  M       (input) INTEGER
028: *          The number of rows of the matrix A. M >= 0.
029: *
030: *  N       (input) INTEGER
031: *          The number of columns of the matrix A. N >= 0
032: *
033: *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
034: *          On entry, the M-by-N matrix A.
035: *          On exit, the upper triangle of the array contains the
036: *          min(M,N)-by-N upper triangular matrix R; the elements
037: *          below the diagonal, together with the array TAU,
038: *          represent the orthogonal matrix Q as a product of
039: *          min(m,n) elementary reflectors.
040: *
041: *  LDA     (input) INTEGER
042: *          The leading dimension of the array A. LDA >= max(1,M).
043: *
044: *  JPVT    (input/output) INTEGER array, dimension (N)
045: *          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
046: *          to the front of A*P (a leading column); if JPVT(i) = 0,
047: *          the i-th column of A is a free column.
048: *          On exit, if JPVT(i) = k, then the i-th column of A*P
049: *          was the k-th column of A.
050: *
051: *  TAU     (output) DOUBLE PRECISION array, dimension (min(M,N))
052: *          The scalar factors of the elementary reflectors.
053: *
054: *  WORK    (workspace) DOUBLE PRECISION array, dimension (3*N)
055: *
056: *  INFO    (output) INTEGER
057: *          = 0:  successful exit
058: *          < 0:  if INFO = -i, the i-th argument had an illegal value
059: *
060: *  Further Details
061: *  ===============
062: *
063: *  The matrix Q is represented as a product of elementary reflectors
064: *
065: *     Q = H(1) H(2) . . . H(n)
066: *
067: *  Each H(i) has the form
068: *
069: *     H = I - tau * v * v'
070: *
071: *  where tau is a real scalar, and v is a real vector with
072: *  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i).
073: *
074: *  The matrix P is represented in jpvt as follows: If
075: *     jpvt(j) = i
076: *  then the jth column of P is the ith canonical unit vector.
077: *
078: *  Partial column norm updating strategy modified by
079: *    Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
080: *    University of Zagreb, Croatia.
081: *    June 2006.
082: *  For more details see LAPACK Working Note 176.
083: *
084: *  =====================================================================
085: *
086: *     .. Parameters ..
087:       DOUBLE PRECISION   ZERO, ONE
088:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
089: *     ..
090: *     .. Local Scalars ..
091:       INTEGER            I, ITEMP, J, MA, MN, PVT
092:       DOUBLE PRECISION   AII, TEMP, TEMP2, TOL3Z
093: *     ..
094: *     .. External Subroutines ..
095:       EXTERNAL           DGEQR2, DLARF, DLARFP, DORM2R, DSWAP, XERBLA
096: *     ..
097: *     .. Intrinsic Functions ..
098:       INTRINSIC          ABS, MAX, MIN, SQRT
099: *     ..
100: *     .. External Functions ..
101:       INTEGER            IDAMAX
102:       DOUBLE PRECISION   DLAMCH, DNRM2
103:       EXTERNAL           IDAMAX, DLAMCH, DNRM2
104: *     ..
105: *     .. Executable Statements ..
106: *
107: *     Test the input arguments
108: *
109:       INFO = 0
110:       IF( M.LT.0 ) THEN
111:          INFO = -1
112:       ELSE IF( N.LT.0 ) THEN
113:          INFO = -2
114:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
115:          INFO = -4
116:       END IF
117:       IF( INFO.NE.0 ) THEN
118:          CALL XERBLA( 'DGEQPF', -INFO )
119:          RETURN
120:       END IF
121: *
122:       MN = MIN( M, N )
123:       TOL3Z = SQRT(DLAMCH('Epsilon'))
124: *
125: *     Move initial columns up front
126: *
127:       ITEMP = 1
128:       DO 10 I = 1, N
129:          IF( JPVT( I ).NE.0 ) THEN
130:             IF( I.NE.ITEMP ) THEN
131:                CALL DSWAP( M, A( 1, I ), 1, A( 1, ITEMP ), 1 )
132:                JPVT( I ) = JPVT( ITEMP )
133:                JPVT( ITEMP ) = I
134:             ELSE
135:                JPVT( I ) = I
136:             END IF
137:             ITEMP = ITEMP + 1
138:          ELSE
139:             JPVT( I ) = I
140:          END IF
141:    10 CONTINUE
142:       ITEMP = ITEMP - 1
143: *
144: *     Compute the QR factorization and update remaining columns
145: *
146:       IF( ITEMP.GT.0 ) THEN
147:          MA = MIN( ITEMP, M )
148:          CALL DGEQR2( M, MA, A, LDA, TAU, WORK, INFO )
149:          IF( MA.LT.N ) THEN
150:             CALL DORM2R( 'Left', 'Transpose', M, N-MA, MA, A, LDA, TAU,
151:      $                   A( 1, MA+1 ), LDA, WORK, INFO )
152:          END IF
153:       END IF
154: *
155:       IF( ITEMP.LT.MN ) THEN
156: *
157: *        Initialize partial column norms. The first n elements of
158: *        work store the exact column norms.
159: *
160:          DO 20 I = ITEMP + 1, N
161:             WORK( I ) = DNRM2( M-ITEMP, A( ITEMP+1, I ), 1 )
162:             WORK( N+I ) = WORK( I )
163:    20    CONTINUE
164: *
165: *        Compute factorization
166: *
167:          DO 40 I = ITEMP + 1, MN
168: *
169: *           Determine ith pivot column and swap if necessary
170: *
171:             PVT = ( I-1 ) + IDAMAX( N-I+1, WORK( I ), 1 )
172: *
173:             IF( PVT.NE.I ) THEN
174:                CALL DSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 )
175:                ITEMP = JPVT( PVT )
176:                JPVT( PVT ) = JPVT( I )
177:                JPVT( I ) = ITEMP
178:                WORK( PVT ) = WORK( I )
179:                WORK( N+PVT ) = WORK( N+I )
180:             END IF
181: *
182: *           Generate elementary reflector H(i)
183: *
184:             IF( I.LT.M ) THEN
185:                CALL DLARFP( M-I+1, A( I, I ), A( I+1, I ), 1, TAU( I ) )
186:             ELSE
187:                CALL DLARFP( 1, A( M, M ), A( M, M ), 1, TAU( M ) )
188:             END IF
189: *
190:             IF( I.LT.N ) THEN
191: *
192: *              Apply H(i) to A(i:m,i+1:n) from the left
193: *
194:                AII = A( I, I )
195:                A( I, I ) = ONE
196:                CALL DLARF( 'LEFT', M-I+1, N-I, A( I, I ), 1, TAU( I ),
197:      $                     A( I, I+1 ), LDA, WORK( 2*N+1 ) )
198:                A( I, I ) = AII
199:             END IF
200: *
201: *           Update partial column norms
202: *
203:             DO 30 J = I + 1, N
204:                IF( WORK( J ).NE.ZERO ) THEN
205: *
206: *                 NOTE: The following 4 lines follow from the analysis in
207: *                 Lapack Working Note 176.
208: *                 
209:                   TEMP = ABS( A( I, J ) ) / WORK( J )
210:                   TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) )
211:                   TEMP2 = TEMP*( WORK( J ) / WORK( N+J ) )**2
212:                   IF( TEMP2 .LE. TOL3Z ) THEN 
213:                      IF( M-I.GT.0 ) THEN
214:                         WORK( J ) = DNRM2( M-I, A( I+1, J ), 1 )
215:                         WORK( N+J ) = WORK( J )
216:                      ELSE
217:                         WORK( J ) = ZERO
218:                         WORK( N+J ) = ZERO
219:                      END IF
220:                   ELSE
221:                      WORK( J ) = WORK( J )*SQRT( TEMP )
222:                   END IF
223:                END IF
224:    30       CONTINUE
225: *
226:    40    CONTINUE
227:       END IF
228:       RETURN
229: *
230: *     End of DGEQPF
231: *
232:       END
233: