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