1 | *> \brief \b SLAQPS |
---|
2 | * |
---|
3 | * =========== DOCUMENTATION =========== |
---|
4 | * |
---|
5 | * Online html documentation available at |
---|
6 | * http://www.netlib.org/lapack/explore-html/ |
---|
7 | * |
---|
8 | *> \htmlonly |
---|
9 | *> Download SLAQPS + dependencies |
---|
10 | *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaqps.f"> |
---|
11 | *> [TGZ]</a> |
---|
12 | *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaqps.f"> |
---|
13 | *> [ZIP]</a> |
---|
14 | *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaqps.f"> |
---|
15 | *> [TXT]</a> |
---|
16 | *> \endhtmlonly |
---|
17 | * |
---|
18 | * Definition: |
---|
19 | * =========== |
---|
20 | * |
---|
21 | * SUBROUTINE SLAQPS( M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU, VN1, |
---|
22 | * VN2, AUXV, F, LDF ) |
---|
23 | * |
---|
24 | * .. Scalar Arguments .. |
---|
25 | * INTEGER KB, LDA, LDF, M, N, NB, OFFSET |
---|
26 | * .. |
---|
27 | * .. Array Arguments .. |
---|
28 | * INTEGER JPVT( * ) |
---|
29 | * REAL A( LDA, * ), AUXV( * ), F( LDF, * ), TAU( * ), |
---|
30 | * $ VN1( * ), VN2( * ) |
---|
31 | * .. |
---|
32 | * |
---|
33 | * |
---|
34 | *> \par Purpose: |
---|
35 | * ============= |
---|
36 | *> |
---|
37 | *> \verbatim |
---|
38 | *> |
---|
39 | *> SLAQPS computes a step of QR factorization with column pivoting |
---|
40 | *> of a real M-by-N matrix A by using Blas-3. It tries to factorize |
---|
41 | *> NB columns from A starting from the row OFFSET+1, and updates all |
---|
42 | *> of the matrix with Blas-3 xGEMM. |
---|
43 | *> |
---|
44 | *> In some cases, due to catastrophic cancellations, it cannot |
---|
45 | *> factorize NB columns. Hence, the actual number of factorized |
---|
46 | *> columns is returned in KB. |
---|
47 | *> |
---|
48 | *> Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized. |
---|
49 | *> \endverbatim |
---|
50 | * |
---|
51 | * Arguments: |
---|
52 | * ========== |
---|
53 | * |
---|
54 | *> \param[in] M |
---|
55 | *> \verbatim |
---|
56 | *> M is INTEGER |
---|
57 | *> The number of rows of the matrix A. M >= 0. |
---|
58 | *> \endverbatim |
---|
59 | *> |
---|
60 | *> \param[in] N |
---|
61 | *> \verbatim |
---|
62 | *> N is INTEGER |
---|
63 | *> The number of columns of the matrix A. N >= 0 |
---|
64 | *> \endverbatim |
---|
65 | *> |
---|
66 | *> \param[in] OFFSET |
---|
67 | *> \verbatim |
---|
68 | *> OFFSET is INTEGER |
---|
69 | *> The number of rows of A that have been factorized in |
---|
70 | *> previous steps. |
---|
71 | *> \endverbatim |
---|
72 | *> |
---|
73 | *> \param[in] NB |
---|
74 | *> \verbatim |
---|
75 | *> NB is INTEGER |
---|
76 | *> The number of columns to factorize. |
---|
77 | *> \endverbatim |
---|
78 | *> |
---|
79 | *> \param[out] KB |
---|
80 | *> \verbatim |
---|
81 | *> KB is INTEGER |
---|
82 | *> The number of columns actually factorized. |
---|
83 | *> \endverbatim |
---|
84 | *> |
---|
85 | *> \param[in,out] A |
---|
86 | *> \verbatim |
---|
87 | *> A is REAL array, dimension (LDA,N) |
---|
88 | *> On entry, the M-by-N matrix A. |
---|
89 | *> On exit, block A(OFFSET+1:M,1:KB) is the triangular |
---|
90 | *> factor obtained and block A(1:OFFSET,1:N) has been |
---|
91 | *> accordingly pivoted, but no factorized. |
---|
92 | *> The rest of the matrix, block A(OFFSET+1:M,KB+1:N) has |
---|
93 | *> been updated. |
---|
94 | *> \endverbatim |
---|
95 | *> |
---|
96 | *> \param[in] LDA |
---|
97 | *> \verbatim |
---|
98 | *> LDA is INTEGER |
---|
99 | *> The leading dimension of the array A. LDA >= max(1,M). |
---|
100 | *> \endverbatim |
---|
101 | *> |
---|
102 | *> \param[in,out] JPVT |
---|
103 | *> \verbatim |
---|
104 | *> JPVT is INTEGER array, dimension (N) |
---|
105 | *> JPVT(I) = K <==> Column K of the full matrix A has been |
---|
106 | *> permuted into position I in AP. |
---|
107 | *> \endverbatim |
---|
108 | *> |
---|
109 | *> \param[out] TAU |
---|
110 | *> \verbatim |
---|
111 | *> TAU is REAL array, dimension (KB) |
---|
112 | *> The scalar factors of the elementary reflectors. |
---|
113 | *> \endverbatim |
---|
114 | *> |
---|
115 | *> \param[in,out] VN1 |
---|
116 | *> \verbatim |
---|
117 | *> VN1 is REAL array, dimension (N) |
---|
118 | *> The vector with the partial column norms. |
---|
119 | *> \endverbatim |
---|
120 | *> |
---|
121 | *> \param[in,out] VN2 |
---|
122 | *> \verbatim |
---|
123 | *> VN2 is REAL array, dimension (N) |
---|
124 | *> The vector with the exact column norms. |
---|
125 | *> \endverbatim |
---|
126 | *> |
---|
127 | *> \param[in,out] AUXV |
---|
128 | *> \verbatim |
---|
129 | *> AUXV is REAL array, dimension (NB) |
---|
130 | *> Auxiliar vector. |
---|
131 | *> \endverbatim |
---|
132 | *> |
---|
133 | *> \param[in,out] F |
---|
134 | *> \verbatim |
---|
135 | *> F is REAL array, dimension (LDF,NB) |
---|
136 | *> Matrix F**T = L*Y**T*A. |
---|
137 | *> \endverbatim |
---|
138 | *> |
---|
139 | *> \param[in] LDF |
---|
140 | *> \verbatim |
---|
141 | *> LDF is INTEGER |
---|
142 | *> The leading dimension of the array F. LDF >= max(1,N). |
---|
143 | *> \endverbatim |
---|
144 | * |
---|
145 | * Authors: |
---|
146 | * ======== |
---|
147 | * |
---|
148 | *> \author Univ. of Tennessee |
---|
149 | *> \author Univ. of California Berkeley |
---|
150 | *> \author Univ. of Colorado Denver |
---|
151 | *> \author NAG Ltd. |
---|
152 | * |
---|
153 | *> \date November 2011 |
---|
154 | * |
---|
155 | *> \ingroup realOTHERauxiliary |
---|
156 | * |
---|
157 | *> \par Contributors: |
---|
158 | * ================== |
---|
159 | *> |
---|
160 | *> G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain |
---|
161 | *> X. Sun, Computer Science Dept., Duke University, USA |
---|
162 | *> |
---|
163 | *> \n |
---|
164 | *> Partial column norm updating strategy modified on April 2011 |
---|
165 | *> Z. Drmac and Z. Bujanovic, Dept. of Mathematics, |
---|
166 | *> University of Zagreb, Croatia. |
---|
167 | * |
---|
168 | *> \par References: |
---|
169 | * ================ |
---|
170 | *> |
---|
171 | *> LAPACK Working Note 176 |
---|
172 | * |
---|
173 | *> \htmlonly |
---|
174 | *> <a href="http://www.netlib.org/lapack/lawnspdf/lawn176.pdf">[PDF]</a> |
---|
175 | *> \endhtmlonly |
---|
176 | * |
---|
177 | * ===================================================================== |
---|
178 | SUBROUTINE SLAQPS( M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU, VN1, |
---|
179 | $ VN2, AUXV, F, LDF ) |
---|
180 | * |
---|
181 | * -- LAPACK auxiliary routine (version 3.4.0) -- |
---|
182 | * -- LAPACK is a software package provided by Univ. of Tennessee, -- |
---|
183 | * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
---|
184 | * November 2011 |
---|
185 | * |
---|
186 | * .. Scalar Arguments .. |
---|
187 | INTEGER KB, LDA, LDF, M, N, NB, OFFSET |
---|
188 | * .. |
---|
189 | * .. Array Arguments .. |
---|
190 | INTEGER JPVT( * ) |
---|
191 | REAL A( LDA, * ), AUXV( * ), F( LDF, * ), TAU( * ), |
---|
192 | $ VN1( * ), VN2( * ) |
---|
193 | * .. |
---|
194 | * |
---|
195 | * ===================================================================== |
---|
196 | * |
---|
197 | * .. Parameters .. |
---|
198 | REAL ZERO, ONE |
---|
199 | PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) |
---|
200 | * .. |
---|
201 | * .. Local Scalars .. |
---|
202 | INTEGER ITEMP, J, K, LASTRK, LSTICC, PVT, RK |
---|
203 | REAL AKK, TEMP, TEMP2, TOL3Z |
---|
204 | * .. |
---|
205 | * .. External Subroutines .. |
---|
206 | EXTERNAL SGEMM, SGEMV, SLARFG, SSWAP |
---|
207 | * .. |
---|
208 | * .. Intrinsic Functions .. |
---|
209 | INTRINSIC ABS, MAX, MIN, NINT, REAL, SQRT |
---|
210 | * .. |
---|
211 | * .. External Functions .. |
---|
212 | INTEGER ISAMAX |
---|
213 | REAL SLAMCH, SNRM2 |
---|
214 | EXTERNAL ISAMAX, SLAMCH, SNRM2 |
---|
215 | * .. |
---|
216 | * .. Executable Statements .. |
---|
217 | * |
---|
218 | LASTRK = MIN( M, N+OFFSET ) |
---|
219 | LSTICC = 0 |
---|
220 | K = 0 |
---|
221 | TOL3Z = SQRT(SLAMCH('Epsilon')) |
---|
222 | * |
---|
223 | * Beginning of while loop. |
---|
224 | * |
---|
225 | 10 CONTINUE |
---|
226 | IF( ( K.LT.NB ) .AND. ( LSTICC.EQ.0 ) ) THEN |
---|
227 | K = K + 1 |
---|
228 | RK = OFFSET + K |
---|
229 | * |
---|
230 | * Determine ith pivot column and swap if necessary |
---|
231 | * |
---|
232 | PVT = ( K-1 ) + ISAMAX( N-K+1, VN1( K ), 1 ) |
---|
233 | IF( PVT.NE.K ) THEN |
---|
234 | CALL SSWAP( M, A( 1, PVT ), 1, A( 1, K ), 1 ) |
---|
235 | CALL SSWAP( K-1, F( PVT, 1 ), LDF, F( K, 1 ), LDF ) |
---|
236 | ITEMP = JPVT( PVT ) |
---|
237 | JPVT( PVT ) = JPVT( K ) |
---|
238 | JPVT( K ) = ITEMP |
---|
239 | VN1( PVT ) = VN1( K ) |
---|
240 | VN2( PVT ) = VN2( K ) |
---|
241 | END IF |
---|
242 | * |
---|
243 | * Apply previous Householder reflectors to column K: |
---|
244 | * A(RK:M,K) := A(RK:M,K) - A(RK:M,1:K-1)*F(K,1:K-1)**T. |
---|
245 | * |
---|
246 | IF( K.GT.1 ) THEN |
---|
247 | CALL SGEMV( 'No transpose', M-RK+1, K-1, -ONE, A( RK, 1 ), |
---|
248 | $ LDA, F( K, 1 ), LDF, ONE, A( RK, K ), 1 ) |
---|
249 | END IF |
---|
250 | * |
---|
251 | * Generate elementary reflector H(k). |
---|
252 | * |
---|
253 | IF( RK.LT.M ) THEN |
---|
254 | CALL SLARFG( M-RK+1, A( RK, K ), A( RK+1, K ), 1, TAU( K ) ) |
---|
255 | ELSE |
---|
256 | CALL SLARFG( 1, A( RK, K ), A( RK, K ), 1, TAU( K ) ) |
---|
257 | END IF |
---|
258 | * |
---|
259 | AKK = A( RK, K ) |
---|
260 | A( RK, K ) = ONE |
---|
261 | * |
---|
262 | * Compute Kth column of F: |
---|
263 | * |
---|
264 | * Compute F(K+1:N,K) := tau(K)*A(RK:M,K+1:N)**T*A(RK:M,K). |
---|
265 | * |
---|
266 | IF( K.LT.N ) THEN |
---|
267 | CALL SGEMV( 'Transpose', M-RK+1, N-K, TAU( K ), |
---|
268 | $ A( RK, K+1 ), LDA, A( RK, K ), 1, ZERO, |
---|
269 | $ F( K+1, K ), 1 ) |
---|
270 | END IF |
---|
271 | * |
---|
272 | * Padding F(1:K,K) with zeros. |
---|
273 | * |
---|
274 | DO 20 J = 1, K |
---|
275 | F( J, K ) = ZERO |
---|
276 | 20 CONTINUE |
---|
277 | * |
---|
278 | * Incremental updating of F: |
---|
279 | * F(1:N,K) := F(1:N,K) - tau(K)*F(1:N,1:K-1)*A(RK:M,1:K-1)**T |
---|
280 | * *A(RK:M,K). |
---|
281 | * |
---|
282 | IF( K.GT.1 ) THEN |
---|
283 | CALL SGEMV( 'Transpose', M-RK+1, K-1, -TAU( K ), A( RK, 1 ), |
---|
284 | $ LDA, A( RK, K ), 1, ZERO, AUXV( 1 ), 1 ) |
---|
285 | * |
---|
286 | CALL SGEMV( 'No transpose', N, K-1, ONE, F( 1, 1 ), LDF, |
---|
287 | $ AUXV( 1 ), 1, ONE, F( 1, K ), 1 ) |
---|
288 | END IF |
---|
289 | * |
---|
290 | * Update the current row of A: |
---|
291 | * A(RK,K+1:N) := A(RK,K+1:N) - A(RK,1:K)*F(K+1:N,1:K)**T. |
---|
292 | * |
---|
293 | IF( K.LT.N ) THEN |
---|
294 | CALL SGEMV( 'No transpose', N-K, K, -ONE, F( K+1, 1 ), LDF, |
---|
295 | $ A( RK, 1 ), LDA, ONE, A( RK, K+1 ), LDA ) |
---|
296 | END IF |
---|
297 | * |
---|
298 | * Update partial column norms. |
---|
299 | * |
---|
300 | IF( RK.LT.LASTRK ) THEN |
---|
301 | DO 30 J = K + 1, N |
---|
302 | IF( VN1( J ).NE.ZERO ) THEN |
---|
303 | * |
---|
304 | * NOTE: The following 4 lines follow from the analysis in |
---|
305 | * Lapack Working Note 176. |
---|
306 | * |
---|
307 | TEMP = ABS( A( RK, J ) ) / VN1( J ) |
---|
308 | TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) ) |
---|
309 | TEMP2 = TEMP*( VN1( J ) / VN2( J ) )**2 |
---|
310 | IF( TEMP2 .LE. TOL3Z ) THEN |
---|
311 | VN2( J ) = REAL( LSTICC ) |
---|
312 | LSTICC = J |
---|
313 | ELSE |
---|
314 | VN1( J ) = VN1( J )*SQRT( TEMP ) |
---|
315 | END IF |
---|
316 | END IF |
---|
317 | 30 CONTINUE |
---|
318 | END IF |
---|
319 | * |
---|
320 | A( RK, K ) = AKK |
---|
321 | * |
---|
322 | * End of while loop. |
---|
323 | * |
---|
324 | GO TO 10 |
---|
325 | END IF |
---|
326 | KB = K |
---|
327 | RK = OFFSET + KB |
---|
328 | * |
---|
329 | * Apply the block reflector to the rest of the matrix: |
---|
330 | * A(OFFSET+KB+1:M,KB+1:N) := A(OFFSET+KB+1:M,KB+1:N) - |
---|
331 | * A(OFFSET+KB+1:M,1:KB)*F(KB+1:N,1:KB)**T. |
---|
332 | * |
---|
333 | IF( KB.LT.MIN( N, M-OFFSET ) ) THEN |
---|
334 | CALL SGEMM( 'No transpose', 'Transpose', M-RK, N-KB, KB, -ONE, |
---|
335 | $ A( RK+1, 1 ), LDA, F( KB+1, 1 ), LDF, ONE, |
---|
336 | $ A( RK+1, KB+1 ), LDA ) |
---|
337 | END IF |
---|
338 | * |
---|
339 | * Recomputation of difficult columns. |
---|
340 | * |
---|
341 | 40 CONTINUE |
---|
342 | IF( LSTICC.GT.0 ) THEN |
---|
343 | ITEMP = NINT( VN2( LSTICC ) ) |
---|
344 | VN1( LSTICC ) = SNRM2( M-RK, A( RK+1, LSTICC ), 1 ) |
---|
345 | * |
---|
346 | * NOTE: The computation of VN1( LSTICC ) relies on the fact that |
---|
347 | * SNRM2 does not fail on vectors with norm below the value of |
---|
348 | * SQRT(DLAMCH('S')) |
---|
349 | * |
---|
350 | VN2( LSTICC ) = VN1( LSTICC ) |
---|
351 | LSTICC = ITEMP |
---|
352 | GO TO 40 |
---|
353 | END IF |
---|
354 | * |
---|
355 | RETURN |
---|
356 | * |
---|
357 | * End of SLAQPS |
---|
358 | * |
---|
359 | END |
---|