1 SUBROUTINE DGEQPF( M, N, A, LDA, JPVT, TAU, WORK, INFO )
2 *
3 * -- LAPACK deprecated computational routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
7 *
8 * .. Scalar Arguments ..
9 INTEGER INFO, LDA, M, N
10 * ..
11 * .. Array Arguments ..
12 INTEGER JPVT( * )
13 DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * This routine is deprecated and has been replaced by routine DGEQP3.
20 *
21 * DGEQPF computes a QR factorization with column pivoting of a
22 * real M-by-N matrix A: A*P = Q*R.
23 *
24 * Arguments
25 * =========
26 *
27 * M (input) INTEGER
28 * The number of rows of the matrix A. M >= 0.
29 *
30 * N (input) INTEGER
31 * The number of columns of the matrix A. N >= 0
32 *
33 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
34 * On entry, the M-by-N matrix A.
35 * On exit, the upper triangle of the array contains the
36 * min(M,N)-by-N upper triangular matrix R; the elements
37 * below the diagonal, together with the array TAU,
38 * represent the orthogonal matrix Q as a product of
39 * min(m,n) elementary reflectors.
40 *
41 * LDA (input) INTEGER
42 * The leading dimension of the array A. LDA >= max(1,M).
43 *
44 * JPVT (input/output) INTEGER array, dimension (N)
45 * On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
46 * to the front of A*P (a leading column); if JPVT(i) = 0,
47 * the i-th column of A is a free column.
48 * On exit, if JPVT(i) = k, then the i-th column of A*P
49 * was the k-th column of A.
50 *
51 * TAU (output) DOUBLE PRECISION array, dimension (min(M,N))
52 * The scalar factors of the elementary reflectors.
53 *
54 * WORK (workspace) DOUBLE PRECISION array, dimension (3*N)
55 *
56 * INFO (output) INTEGER
57 * = 0: successful exit
58 * < 0: if INFO = -i, the i-th argument had an illegal value
59 *
60 * Further Details
61 * ===============
62 *
63 * The matrix Q is represented as a product of elementary reflectors
64 *
65 * Q = H(1) H(2) . . . H(n)
66 *
67 * Each H(i) has the form
68 *
69 * H = I - tau * v * v**T
70 *
71 * where tau is a real scalar, and v is a real vector with
72 * v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i).
73 *
74 * The matrix P is represented in jpvt as follows: If
75 * jpvt(j) = i
76 * then the jth column of P is the ith canonical unit vector.
77 *
78 * Partial column norm updating strategy modified by
79 * Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
80 * University of Zagreb, Croatia.
81 * -- April 2011 --
82 * For more details see LAPACK Working Note 176.
83 *
84 * =====================================================================
85 *
86 * .. Parameters ..
87 DOUBLE PRECISION ZERO, ONE
88 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
89 * ..
90 * .. Local Scalars ..
91 INTEGER I, ITEMP, J, MA, MN, PVT
92 DOUBLE PRECISION AII, TEMP, TEMP2, TOL3Z
93 * ..
94 * .. External Subroutines ..
95 EXTERNAL DGEQR2, DLARF, DLARFG, DORM2R, DSWAP, XERBLA
96 * ..
97 * .. Intrinsic Functions ..
98 INTRINSIC ABS, MAX, MIN, SQRT
99 * ..
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 DLARFG( M-I+1, A( I, I ), A( I+1, I ), 1, TAU( I ) )
186 ELSE
187 CALL DLARFG( 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
2 *
3 * -- LAPACK deprecated computational routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
7 *
8 * .. Scalar Arguments ..
9 INTEGER INFO, LDA, M, N
10 * ..
11 * .. Array Arguments ..
12 INTEGER JPVT( * )
13 DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * This routine is deprecated and has been replaced by routine DGEQP3.
20 *
21 * DGEQPF computes a QR factorization with column pivoting of a
22 * real M-by-N matrix A: A*P = Q*R.
23 *
24 * Arguments
25 * =========
26 *
27 * M (input) INTEGER
28 * The number of rows of the matrix A. M >= 0.
29 *
30 * N (input) INTEGER
31 * The number of columns of the matrix A. N >= 0
32 *
33 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
34 * On entry, the M-by-N matrix A.
35 * On exit, the upper triangle of the array contains the
36 * min(M,N)-by-N upper triangular matrix R; the elements
37 * below the diagonal, together with the array TAU,
38 * represent the orthogonal matrix Q as a product of
39 * min(m,n) elementary reflectors.
40 *
41 * LDA (input) INTEGER
42 * The leading dimension of the array A. LDA >= max(1,M).
43 *
44 * JPVT (input/output) INTEGER array, dimension (N)
45 * On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
46 * to the front of A*P (a leading column); if JPVT(i) = 0,
47 * the i-th column of A is a free column.
48 * On exit, if JPVT(i) = k, then the i-th column of A*P
49 * was the k-th column of A.
50 *
51 * TAU (output) DOUBLE PRECISION array, dimension (min(M,N))
52 * The scalar factors of the elementary reflectors.
53 *
54 * WORK (workspace) DOUBLE PRECISION array, dimension (3*N)
55 *
56 * INFO (output) INTEGER
57 * = 0: successful exit
58 * < 0: if INFO = -i, the i-th argument had an illegal value
59 *
60 * Further Details
61 * ===============
62 *
63 * The matrix Q is represented as a product of elementary reflectors
64 *
65 * Q = H(1) H(2) . . . H(n)
66 *
67 * Each H(i) has the form
68 *
69 * H = I - tau * v * v**T
70 *
71 * where tau is a real scalar, and v is a real vector with
72 * v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i).
73 *
74 * The matrix P is represented in jpvt as follows: If
75 * jpvt(j) = i
76 * then the jth column of P is the ith canonical unit vector.
77 *
78 * Partial column norm updating strategy modified by
79 * Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
80 * University of Zagreb, Croatia.
81 * -- April 2011 --
82 * For more details see LAPACK Working Note 176.
83 *
84 * =====================================================================
85 *
86 * .. Parameters ..
87 DOUBLE PRECISION ZERO, ONE
88 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
89 * ..
90 * .. Local Scalars ..
91 INTEGER I, ITEMP, J, MA, MN, PVT
92 DOUBLE PRECISION AII, TEMP, TEMP2, TOL3Z
93 * ..
94 * .. External Subroutines ..
95 EXTERNAL DGEQR2, DLARF, DLARFG, DORM2R, DSWAP, XERBLA
96 * ..
97 * .. Intrinsic Functions ..
98 INTRINSIC ABS, MAX, MIN, SQRT
99 * ..
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 DLARFG( M-I+1, A( I, I ), A( I+1, I ), 1, TAU( I ) )
186 ELSE
187 CALL DLARFG( 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