1 SUBROUTINE DLAQPS( M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU, VN1,
2 $ VN2, AUXV, F, LDF )
3 *
4 * -- LAPACK auxiliary routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * .. Scalar Arguments ..
10 INTEGER KB, LDA, LDF, M, N, NB, OFFSET
11 * ..
12 * .. Array Arguments ..
13 INTEGER JPVT( * )
14 DOUBLE PRECISION A( LDA, * ), AUXV( * ), F( LDF, * ), TAU( * ),
15 $ VN1( * ), VN2( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DLAQPS computes a step of QR factorization with column pivoting
22 * of a real M-by-N matrix A by using Blas-3. It tries to factorize
23 * NB columns from A starting from the row OFFSET+1, and updates all
24 * of the matrix with Blas-3 xGEMM.
25 *
26 * In some cases, due to catastrophic cancellations, it cannot
27 * factorize NB columns. Hence, the actual number of factorized
28 * columns is returned in KB.
29 *
30 * Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.
31 *
32 * Arguments
33 * =========
34 *
35 * M (input) INTEGER
36 * The number of rows of the matrix A. M >= 0.
37 *
38 * N (input) INTEGER
39 * The number of columns of the matrix A. N >= 0
40 *
41 * OFFSET (input) INTEGER
42 * The number of rows of A that have been factorized in
43 * previous steps.
44 *
45 * NB (input) INTEGER
46 * The number of columns to factorize.
47 *
48 * KB (output) INTEGER
49 * The number of columns actually factorized.
50 *
51 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
52 * On entry, the M-by-N matrix A.
53 * On exit, block A(OFFSET+1:M,1:KB) is the triangular
54 * factor obtained and block A(1:OFFSET,1:N) has been
55 * accordingly pivoted, but no factorized.
56 * The rest of the matrix, block A(OFFSET+1:M,KB+1:N) has
57 * been updated.
58 *
59 * LDA (input) INTEGER
60 * The leading dimension of the array A. LDA >= max(1,M).
61 *
62 * JPVT (input/output) INTEGER array, dimension (N)
63 * JPVT(I) = K <==> Column K of the full matrix A has been
64 * permuted into position I in AP.
65 *
66 * TAU (output) DOUBLE PRECISION array, dimension (KB)
67 * The scalar factors of the elementary reflectors.
68 *
69 * VN1 (input/output) DOUBLE PRECISION array, dimension (N)
70 * The vector with the partial column norms.
71 *
72 * VN2 (input/output) DOUBLE PRECISION array, dimension (N)
73 * The vector with the exact column norms.
74 *
75 * AUXV (input/output) DOUBLE PRECISION array, dimension (NB)
76 * Auxiliar vector.
77 *
78 * F (input/output) DOUBLE PRECISION array, dimension (LDF,NB)
79 * Matrix F**T = L*Y**T*A.
80 *
81 * LDF (input) INTEGER
82 * The leading dimension of the array F. LDF >= max(1,N).
83 *
84 * Further Details
85 * ===============
86 *
87 * Based on contributions by
88 * G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
89 * X. Sun, Computer Science Dept., Duke University, USA
90 *
91 * Partial column norm updating strategy modified by
92 * Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
93 * University of Zagreb, Croatia.
94 * -- April 2011 --
95 * For more details see LAPACK Working Note 176.
96 * =====================================================================
97 *
98 * .. Parameters ..
99 DOUBLE PRECISION ZERO, ONE
100 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
101 * ..
102 * .. Local Scalars ..
103 INTEGER ITEMP, J, K, LASTRK, LSTICC, PVT, RK
104 DOUBLE PRECISION AKK, TEMP, TEMP2, TOL3Z
105 * ..
106 * .. External Subroutines ..
107 EXTERNAL DGEMM, DGEMV, DLARFG, DSWAP
108 * ..
109 * .. Intrinsic Functions ..
110 INTRINSIC ABS, DBLE, MAX, MIN, NINT, SQRT
111 * ..
112 * .. External Functions ..
113 INTEGER IDAMAX
114 DOUBLE PRECISION DLAMCH, DNRM2
115 EXTERNAL IDAMAX, DLAMCH, DNRM2
116 * ..
117 * .. Executable Statements ..
118 *
119 LASTRK = MIN( M, N+OFFSET )
120 LSTICC = 0
121 K = 0
122 TOL3Z = SQRT(DLAMCH('Epsilon'))
123 *
124 * Beginning of while loop.
125 *
126 10 CONTINUE
127 IF( ( K.LT.NB ) .AND. ( LSTICC.EQ.0 ) ) THEN
128 K = K + 1
129 RK = OFFSET + K
130 *
131 * Determine ith pivot column and swap if necessary
132 *
133 PVT = ( K-1 ) + IDAMAX( N-K+1, VN1( K ), 1 )
134 IF( PVT.NE.K ) THEN
135 CALL DSWAP( M, A( 1, PVT ), 1, A( 1, K ), 1 )
136 CALL DSWAP( K-1, F( PVT, 1 ), LDF, F( K, 1 ), LDF )
137 ITEMP = JPVT( PVT )
138 JPVT( PVT ) = JPVT( K )
139 JPVT( K ) = ITEMP
140 VN1( PVT ) = VN1( K )
141 VN2( PVT ) = VN2( K )
142 END IF
143 *
144 * Apply previous Householder reflectors to column K:
145 * A(RK:M,K) := A(RK:M,K) - A(RK:M,1:K-1)*F(K,1:K-1)**T.
146 *
147 IF( K.GT.1 ) THEN
148 CALL DGEMV( 'No transpose', M-RK+1, K-1, -ONE, A( RK, 1 ),
149 $ LDA, F( K, 1 ), LDF, ONE, A( RK, K ), 1 )
150 END IF
151 *
152 * Generate elementary reflector H(k).
153 *
154 IF( RK.LT.M ) THEN
155 CALL DLARFG( M-RK+1, A( RK, K ), A( RK+1, K ), 1, TAU( K ) )
156 ELSE
157 CALL DLARFG( 1, A( RK, K ), A( RK, K ), 1, TAU( K ) )
158 END IF
159 *
160 AKK = A( RK, K )
161 A( RK, K ) = ONE
162 *
163 * Compute Kth column of F:
164 *
165 * Compute F(K+1:N,K) := tau(K)*A(RK:M,K+1:N)**T*A(RK:M,K).
166 *
167 IF( K.LT.N ) THEN
168 CALL DGEMV( 'Transpose', M-RK+1, N-K, TAU( K ),
169 $ A( RK, K+1 ), LDA, A( RK, K ), 1, ZERO,
170 $ F( K+1, K ), 1 )
171 END IF
172 *
173 * Padding F(1:K,K) with zeros.
174 *
175 DO 20 J = 1, K
176 F( J, K ) = ZERO
177 20 CONTINUE
178 *
179 * Incremental updating of F:
180 * F(1:N,K) := F(1:N,K) - tau(K)*F(1:N,1:K-1)*A(RK:M,1:K-1)**T
181 * *A(RK:M,K).
182 *
183 IF( K.GT.1 ) THEN
184 CALL DGEMV( 'Transpose', M-RK+1, K-1, -TAU( K ), A( RK, 1 ),
185 $ LDA, A( RK, K ), 1, ZERO, AUXV( 1 ), 1 )
186 *
187 CALL DGEMV( 'No transpose', N, K-1, ONE, F( 1, 1 ), LDF,
188 $ AUXV( 1 ), 1, ONE, F( 1, K ), 1 )
189 END IF
190 *
191 * Update the current row of A:
192 * A(RK,K+1:N) := A(RK,K+1:N) - A(RK,1:K)*F(K+1:N,1:K)**T.
193 *
194 IF( K.LT.N ) THEN
195 CALL DGEMV( 'No transpose', N-K, K, -ONE, F( K+1, 1 ), LDF,
196 $ A( RK, 1 ), LDA, ONE, A( RK, K+1 ), LDA )
197 END IF
198 *
199 * Update partial column norms.
200 *
201 IF( RK.LT.LASTRK ) THEN
202 DO 30 J = K + 1, N
203 IF( VN1( 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( RK, J ) ) / VN1( J )
209 TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) )
210 TEMP2 = TEMP*( VN1( J ) / VN2( J ) )**2
211 IF( TEMP2 .LE. TOL3Z ) THEN
212 VN2( J ) = DBLE( LSTICC )
213 LSTICC = J
214 ELSE
215 VN1( J ) = VN1( J )*SQRT( TEMP )
216 END IF
217 END IF
218 30 CONTINUE
219 END IF
220 *
221 A( RK, K ) = AKK
222 *
223 * End of while loop.
224 *
225 GO TO 10
226 END IF
227 KB = K
228 RK = OFFSET + KB
229 *
230 * Apply the block reflector to the rest of the matrix:
231 * A(OFFSET+KB+1:M,KB+1:N) := A(OFFSET+KB+1:M,KB+1:N) -
232 * A(OFFSET+KB+1:M,1:KB)*F(KB+1:N,1:KB)**T.
233 *
234 IF( KB.LT.MIN( N, M-OFFSET ) ) THEN
235 CALL DGEMM( 'No transpose', 'Transpose', M-RK, N-KB, KB, -ONE,
236 $ A( RK+1, 1 ), LDA, F( KB+1, 1 ), LDF, ONE,
237 $ A( RK+1, KB+1 ), LDA )
238 END IF
239 *
240 * Recomputation of difficult columns.
241 *
242 40 CONTINUE
243 IF( LSTICC.GT.0 ) THEN
244 ITEMP = NINT( VN2( LSTICC ) )
245 VN1( LSTICC ) = DNRM2( M-RK, A( RK+1, LSTICC ), 1 )
246 *
247 * NOTE: The computation of VN1( LSTICC ) relies on the fact that
248 * SNRM2 does not fail on vectors with norm below the value of
249 * SQRT(DLAMCH('S'))
250 *
251 VN2( LSTICC ) = VN1( LSTICC )
252 LSTICC = ITEMP
253 GO TO 40
254 END IF
255 *
256 RETURN
257 *
258 * End of DLAQPS
259 *
260 END
2 $ VN2, AUXV, F, LDF )
3 *
4 * -- LAPACK auxiliary routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * .. Scalar Arguments ..
10 INTEGER KB, LDA, LDF, M, N, NB, OFFSET
11 * ..
12 * .. Array Arguments ..
13 INTEGER JPVT( * )
14 DOUBLE PRECISION A( LDA, * ), AUXV( * ), F( LDF, * ), TAU( * ),
15 $ VN1( * ), VN2( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DLAQPS computes a step of QR factorization with column pivoting
22 * of a real M-by-N matrix A by using Blas-3. It tries to factorize
23 * NB columns from A starting from the row OFFSET+1, and updates all
24 * of the matrix with Blas-3 xGEMM.
25 *
26 * In some cases, due to catastrophic cancellations, it cannot
27 * factorize NB columns. Hence, the actual number of factorized
28 * columns is returned in KB.
29 *
30 * Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.
31 *
32 * Arguments
33 * =========
34 *
35 * M (input) INTEGER
36 * The number of rows of the matrix A. M >= 0.
37 *
38 * N (input) INTEGER
39 * The number of columns of the matrix A. N >= 0
40 *
41 * OFFSET (input) INTEGER
42 * The number of rows of A that have been factorized in
43 * previous steps.
44 *
45 * NB (input) INTEGER
46 * The number of columns to factorize.
47 *
48 * KB (output) INTEGER
49 * The number of columns actually factorized.
50 *
51 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
52 * On entry, the M-by-N matrix A.
53 * On exit, block A(OFFSET+1:M,1:KB) is the triangular
54 * factor obtained and block A(1:OFFSET,1:N) has been
55 * accordingly pivoted, but no factorized.
56 * The rest of the matrix, block A(OFFSET+1:M,KB+1:N) has
57 * been updated.
58 *
59 * LDA (input) INTEGER
60 * The leading dimension of the array A. LDA >= max(1,M).
61 *
62 * JPVT (input/output) INTEGER array, dimension (N)
63 * JPVT(I) = K <==> Column K of the full matrix A has been
64 * permuted into position I in AP.
65 *
66 * TAU (output) DOUBLE PRECISION array, dimension (KB)
67 * The scalar factors of the elementary reflectors.
68 *
69 * VN1 (input/output) DOUBLE PRECISION array, dimension (N)
70 * The vector with the partial column norms.
71 *
72 * VN2 (input/output) DOUBLE PRECISION array, dimension (N)
73 * The vector with the exact column norms.
74 *
75 * AUXV (input/output) DOUBLE PRECISION array, dimension (NB)
76 * Auxiliar vector.
77 *
78 * F (input/output) DOUBLE PRECISION array, dimension (LDF,NB)
79 * Matrix F**T = L*Y**T*A.
80 *
81 * LDF (input) INTEGER
82 * The leading dimension of the array F. LDF >= max(1,N).
83 *
84 * Further Details
85 * ===============
86 *
87 * Based on contributions by
88 * G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
89 * X. Sun, Computer Science Dept., Duke University, USA
90 *
91 * Partial column norm updating strategy modified by
92 * Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
93 * University of Zagreb, Croatia.
94 * -- April 2011 --
95 * For more details see LAPACK Working Note 176.
96 * =====================================================================
97 *
98 * .. Parameters ..
99 DOUBLE PRECISION ZERO, ONE
100 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
101 * ..
102 * .. Local Scalars ..
103 INTEGER ITEMP, J, K, LASTRK, LSTICC, PVT, RK
104 DOUBLE PRECISION AKK, TEMP, TEMP2, TOL3Z
105 * ..
106 * .. External Subroutines ..
107 EXTERNAL DGEMM, DGEMV, DLARFG, DSWAP
108 * ..
109 * .. Intrinsic Functions ..
110 INTRINSIC ABS, DBLE, MAX, MIN, NINT, SQRT
111 * ..
112 * .. External Functions ..
113 INTEGER IDAMAX
114 DOUBLE PRECISION DLAMCH, DNRM2
115 EXTERNAL IDAMAX, DLAMCH, DNRM2
116 * ..
117 * .. Executable Statements ..
118 *
119 LASTRK = MIN( M, N+OFFSET )
120 LSTICC = 0
121 K = 0
122 TOL3Z = SQRT(DLAMCH('Epsilon'))
123 *
124 * Beginning of while loop.
125 *
126 10 CONTINUE
127 IF( ( K.LT.NB ) .AND. ( LSTICC.EQ.0 ) ) THEN
128 K = K + 1
129 RK = OFFSET + K
130 *
131 * Determine ith pivot column and swap if necessary
132 *
133 PVT = ( K-1 ) + IDAMAX( N-K+1, VN1( K ), 1 )
134 IF( PVT.NE.K ) THEN
135 CALL DSWAP( M, A( 1, PVT ), 1, A( 1, K ), 1 )
136 CALL DSWAP( K-1, F( PVT, 1 ), LDF, F( K, 1 ), LDF )
137 ITEMP = JPVT( PVT )
138 JPVT( PVT ) = JPVT( K )
139 JPVT( K ) = ITEMP
140 VN1( PVT ) = VN1( K )
141 VN2( PVT ) = VN2( K )
142 END IF
143 *
144 * Apply previous Householder reflectors to column K:
145 * A(RK:M,K) := A(RK:M,K) - A(RK:M,1:K-1)*F(K,1:K-1)**T.
146 *
147 IF( K.GT.1 ) THEN
148 CALL DGEMV( 'No transpose', M-RK+1, K-1, -ONE, A( RK, 1 ),
149 $ LDA, F( K, 1 ), LDF, ONE, A( RK, K ), 1 )
150 END IF
151 *
152 * Generate elementary reflector H(k).
153 *
154 IF( RK.LT.M ) THEN
155 CALL DLARFG( M-RK+1, A( RK, K ), A( RK+1, K ), 1, TAU( K ) )
156 ELSE
157 CALL DLARFG( 1, A( RK, K ), A( RK, K ), 1, TAU( K ) )
158 END IF
159 *
160 AKK = A( RK, K )
161 A( RK, K ) = ONE
162 *
163 * Compute Kth column of F:
164 *
165 * Compute F(K+1:N,K) := tau(K)*A(RK:M,K+1:N)**T*A(RK:M,K).
166 *
167 IF( K.LT.N ) THEN
168 CALL DGEMV( 'Transpose', M-RK+1, N-K, TAU( K ),
169 $ A( RK, K+1 ), LDA, A( RK, K ), 1, ZERO,
170 $ F( K+1, K ), 1 )
171 END IF
172 *
173 * Padding F(1:K,K) with zeros.
174 *
175 DO 20 J = 1, K
176 F( J, K ) = ZERO
177 20 CONTINUE
178 *
179 * Incremental updating of F:
180 * F(1:N,K) := F(1:N,K) - tau(K)*F(1:N,1:K-1)*A(RK:M,1:K-1)**T
181 * *A(RK:M,K).
182 *
183 IF( K.GT.1 ) THEN
184 CALL DGEMV( 'Transpose', M-RK+1, K-1, -TAU( K ), A( RK, 1 ),
185 $ LDA, A( RK, K ), 1, ZERO, AUXV( 1 ), 1 )
186 *
187 CALL DGEMV( 'No transpose', N, K-1, ONE, F( 1, 1 ), LDF,
188 $ AUXV( 1 ), 1, ONE, F( 1, K ), 1 )
189 END IF
190 *
191 * Update the current row of A:
192 * A(RK,K+1:N) := A(RK,K+1:N) - A(RK,1:K)*F(K+1:N,1:K)**T.
193 *
194 IF( K.LT.N ) THEN
195 CALL DGEMV( 'No transpose', N-K, K, -ONE, F( K+1, 1 ), LDF,
196 $ A( RK, 1 ), LDA, ONE, A( RK, K+1 ), LDA )
197 END IF
198 *
199 * Update partial column norms.
200 *
201 IF( RK.LT.LASTRK ) THEN
202 DO 30 J = K + 1, N
203 IF( VN1( 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( RK, J ) ) / VN1( J )
209 TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) )
210 TEMP2 = TEMP*( VN1( J ) / VN2( J ) )**2
211 IF( TEMP2 .LE. TOL3Z ) THEN
212 VN2( J ) = DBLE( LSTICC )
213 LSTICC = J
214 ELSE
215 VN1( J ) = VN1( J )*SQRT( TEMP )
216 END IF
217 END IF
218 30 CONTINUE
219 END IF
220 *
221 A( RK, K ) = AKK
222 *
223 * End of while loop.
224 *
225 GO TO 10
226 END IF
227 KB = K
228 RK = OFFSET + KB
229 *
230 * Apply the block reflector to the rest of the matrix:
231 * A(OFFSET+KB+1:M,KB+1:N) := A(OFFSET+KB+1:M,KB+1:N) -
232 * A(OFFSET+KB+1:M,1:KB)*F(KB+1:N,1:KB)**T.
233 *
234 IF( KB.LT.MIN( N, M-OFFSET ) ) THEN
235 CALL DGEMM( 'No transpose', 'Transpose', M-RK, N-KB, KB, -ONE,
236 $ A( RK+1, 1 ), LDA, F( KB+1, 1 ), LDF, ONE,
237 $ A( RK+1, KB+1 ), LDA )
238 END IF
239 *
240 * Recomputation of difficult columns.
241 *
242 40 CONTINUE
243 IF( LSTICC.GT.0 ) THEN
244 ITEMP = NINT( VN2( LSTICC ) )
245 VN1( LSTICC ) = DNRM2( M-RK, A( RK+1, LSTICC ), 1 )
246 *
247 * NOTE: The computation of VN1( LSTICC ) relies on the fact that
248 * SNRM2 does not fail on vectors with norm below the value of
249 * SQRT(DLAMCH('S'))
250 *
251 VN2( LSTICC ) = VN1( LSTICC )
252 LSTICC = ITEMP
253 GO TO 40
254 END IF
255 *
256 RETURN
257 *
258 * End of DLAQPS
259 *
260 END